use Ctr if strong duality constraint is linear

dev/main.cpp

@@ -20,45 +20,52 @@
 #include "idol/mixed-integer/modeling/expressions/QuadExpr.h"
 #include "idol/general/utils/GenerationPattern.h"
 #include "idol/mixed-integer/modeling/models/KKT.h"
+#include "idol/mixed-integer/problems/facility-location-problem/FLP_Instance.h"
 
 using namespace idol;
 
 int main(int t_argc, const char** t_argv) {
 
     Env env;
-    Model primal(env);
-
-    const auto x = primal.add_vars(Dim<1>(10), 0, 1, Continuous, -1, "x");
-    primal.add_ctr(idol_Sum(i, Range(10), i * x[i]) <= 5);
-
-    Model dual(env);
-    Model strong_duality(env);
-    Model kkt(env);
-
-    KKT dualizer(primal, primal.get_obj_expr());
-    dualizer.add_dual(dual);
-    dualizer.add_strong_duality_reformulation(strong_duality);
-    dualizer.add_kkt_reformulation(kkt);
-
-    std::cout << primal << std::endl;
-    std::cout << dual << std::endl;
-    std::cout << strong_duality << std::endl;
-    std::cout << kkt << std::endl;
-
-    primal.use(Gurobi());
-    dual.use(Gurobi());
-    strong_duality.use(Gurobi());
-    kkt.use(Gurobi());
-
-    primal.optimize();
-    dual.optimize();
-    strong_duality.optimize();
-    kkt.optimize();
-
-    std::cout << "Primal solution: " << primal.get_best_obj() << std::endl;
-    std::cout << "Dual solution: " << dual.get_best_obj() << std::endl;
-    std::cout << "Strong duality solution: " << evaluate(dualizer.get_dual_obj_expr(), save_primal(strong_duality)) << std::endl;
-    std::cout << "KKT solution: " << evaluate(dualizer.get_dual_obj_expr(), save_primal(kkt)) << std::endl;
+
+    // Generate parameters
+    const auto filename = "/home/henri/Research/idol/tests/data/facility-location-problem/instance_F10_C20__0.txt";
+
+    // Read instance
+    const auto instance = idol::Problems::FLP::read_instance_1991_Cornuejols_et_al(filename);
+    const unsigned int n_customers = instance.n_customers();
+    const unsigned int n_facilities = instance.n_facilities();
+
+    // Make model
+    auto x = Var::make_vector(env, Dim<1>(n_facilities), 0., 1., Continuous, 0., "x");
+    auto y = Var::make_vector(env, Dim<2>(n_facilities, n_customers), 0., 1., Continuous, 0., "y");
+
+    Model model(env);
+
+    model.add_vector<Var, 1>(x);
+    model.add_vector<Var, 2>(y);
+
+    for (unsigned int i = 0 ; i < n_facilities ; ++i) {
+        model.add(Ctr(env, idol_Sum(j, Range(n_customers), instance.demand(j) * y[i][j]) <= instance.capacity(i) * x[i]));
+    }
+
+    for (unsigned int j = 0 ; j < n_customers ; ++j) {
+        model.add(Ctr(env, idol_Sum(i, Range(n_facilities), y[i][j]) == 1));
+    }
+
+    model.set_obj_expr(idol_Sum(i, Range(n_facilities), instance.fixed_cost(i) * x[i] + idol_Sum(j, Range(n_customers),
+                                                                                                 instance.per_unit_transportation_cost(
+                                                                                                         i, j) *
+                                                                                                 instance.demand(j) *
+                                                                                                 y[i][j])));
+
+    Model reformulated(env);
+
+    Reformulators::KKT reformulator(model, model.get_obj_expr());
+
+    reformulator.add_strong_duality_reformulation(reformulated);
+
+    std::cout << reformulated << std::endl;
 
     return 0;
 }

lib/include/idol/general/utils/SparseVector.h

@@ -82,7 +82,11 @@ class idol::SparseVector {
 public:
     SparseVector() = default;
 
-    SparseVector(const IndexT& t_index, const ValueT& t_value) : m_map({ { t_index, t_value} }) {}
+    SparseVector(const IndexT& t_index, const ValueT& t_value) : m_map({ { t_index, t_value} }) {
+        if (::idol::is_zero(t_value, Tolerance::Sparsity)) {
+            m_map.clear();
+        }
+    }
 
     SparseVector(const SparseVector&) = default;
     SparseVector(SparseVector&&) = default;

lib/src/mixed-integer/modeling/models/KKT.cpp

@@ -530,7 +530,13 @@ void idol::Reformulators::KKT::add_strong_duality_reformulation(idol::Model &t_d
     add_primal_constraints(t_destination);
     add_dual_constraints(t_destination);
 
-    t_destination.add_qctr(m_primal_objective - m_dual_objective, LessOrEqual);
+    auto duality_gap = m_primal_objective - m_dual_objective;
+
+    if (duality_gap.has_quadratic()) {
+        t_destination.add_qctr(std::move(duality_gap), LessOrEqual);
+    } else {
+        t_destination.add_ctr(std::move(duality_gap.affine()) <= 0);
+    }
 }
 
 void idol::Reformulators::KKT::add_kkt_reformulation(idol::Model &t_destination) {

tests/mixed-integer/dualizer/location.cpp

@@ -76,17 +76,17 @@ TEST_CASE("Reformulate Facility Location Problem instances using optimality cond
 
     Model reformulated(env);
 
-    KKT dualizer(model, model.get_obj_expr());
+    Reformulators::KKT reformulator(model, model.get_obj_expr());
 
     switch (reformulation) {
         case Dual:
-            dualizer.add_dual(reformulated);
+            reformulator.add_dual(reformulated);
             break;
         case StrongDuality:
-            dualizer.add_strong_duality_reformulation(reformulated);
+            reformulator.add_strong_duality_reformulation(reformulated);
             break;
         case KKT:
-            dualizer.add_kkt_reformulation(reformulated);
+            reformulator.add_kkt_reformulation(reformulated);
             break;
         default: throw std::runtime_error("Invalid reformulation.");
     }
@@ -106,7 +106,7 @@ TEST_CASE("Reformulate Facility Location Problem instances using optimality cond
         AND_THEN("The objective value is " + std::to_string(objective_value)) {
 
             const auto primal_solution = save_primal(reformulated);
-            double value = evaluate(dualizer.get_dual_obj_expr(), primal_solution);
+            double value = evaluate(reformulator.get_dual_obj_expr(), primal_solution);
 
             CHECK(value == Catch::Approx(objective_value));