first version of Dualizer

hlefebvr · Nov 24, 2024 · 382cc89 · 382cc89
1 parent 7a9260e
commit 382cc89
Show file tree

Hide file tree

Showing 5 changed files with 599 additions and 11 deletions.
diff --git a/dev/main.cpp b/dev/main.cpp
@@ -26,11 +26,39 @@ using namespace idol;
 int main(int t_argc, const char** t_argv) {
 
     Env env;
-    Model model(env);
+    Model primal(env);
 
-    const auto x = model.add_vars(Dim<1>(10), 0, 1, Binary, 1, "x");
+    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);
 
-    Dualizer dualizer(model);
+    Model dual(env);
+    Model strong_duality(env);
+    Model kkt(env);
+
+    Dualizer 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;
 
     return 0;
 }
diff --git a/lib/include/idol/mixed-integer/modeling/expressions/operations/operators.h b/lib/include/idol/mixed-integer/modeling/expressions/operations/operators.h
@@ -127,7 +127,17 @@ DEFINE_COMMUTATIVE_OPERATOR(*, *=, double, QuadExpr<T>, QuadExpr<T>) \
 static QuadExpr<T> operator*(T t_a, T t_b) { return { t_a, t_b }; }   \
 \
 static LinExpr<T, double> operator+(const T& t_a) { return t_a; } \
-static LinExpr<T, double> operator-(const T& t_a) { return { -1., t_a }; }
+static LinExpr<T, double> operator-(const T& t_a) { return { -1., t_a }; }    \
+static QuadExpr<T> operator*(const T& t_a, const LinExpr<T>& t_b) {           \
+    QuadExpr<T> result;      \
+    for (const auto& [var, constant] : t_b) {                      \
+        result += constant * (var * t_a);\
+    }                        \
+    return result;           \
+}                            \
+static QuadExpr<T> operator*(const LinExpr<T>& t_a, const T& t_b) { return t_b * t_a; }       \
+static QuadExpr<T> operator*(const T& t_a, const AffExpr<T>& t_b) { return t_a * t_b.constant() + t_a * t_b.linear(); } \
+static QuadExpr<T> operator*(const AffExpr<T>& t_a, const T& t_b) { return t_b * t_a; }
 
 namespace idol {
     DEFINE_OPERATIONS(Var)

diff --git a/lib/include/idol/mixed-integer/modeling/expressions/operations/operators_utils.h b/lib/include/idol/mixed-integer/modeling/expressions/operations/operators_utils.h
@@ -24,6 +24,7 @@ class idol::Range {
         explicit iterator(NumT t_index) : m_index(t_index) {}
     public:
         bool operator!=(const iterator& t_it) const { return t_it.m_index != m_index; }
+        bool operator==(const iterator& t_it) const { return t_it.m_index == m_index; }
         iterator& operator++() { ++m_index; return *this; }
         iterator operator++(int) & { auto result = *this; ++m_index; return *this; }
         NumT operator*() const { return m_index; }
@@ -34,7 +35,6 @@ class idol::Range {
     [[nodiscard]] iterator end() const { return iterator(m_end); }
 };
 
-
 #define idol_Sum(INDEX, ITERABLE, EXPR) \
 [&]() {                              \
     const auto compute_sum = [&]() { \

diff --git a/lib/include/idol/mixed-integer/modeling/models/Dualizer.h b/lib/include/idol/mixed-integer/modeling/models/Dualizer.h
@@ -10,19 +10,58 @@
 #include "idol/mixed-integer/modeling/variables/Var.h"
 #include "idol/mixed-integer/modeling/constraints/Ctr.h"
 #include "idol/mixed-integer/modeling/constraints/QCtr.h"
+#include "idol/mixed-integer/modeling/expressions/QuadExpr.h"
 
 namespace idol {
     class Dualizer;
     class Model;
 }
 
 class idol::Dualizer {
-    const Model& m_parent;
+    const Model& m_primal;
+    const QuadExpr<Var, double>& m_primal_objective;
 
-    std::vector<Var> m_dual_variables; // constraint - qconstraint index <-> dual variable
-    std::vector<std::variant<Ctr, QCtr>> m_dual_constraints; // variable index <-> dual constraint
+    std::function<bool(const Var&)> m_primal_variable_indicator;
+    std::function<bool(const Ctr&)> m_primal_constraint_indicator;
+    std::function<bool(const QCtr&)> m_primal_qconstraint_indicator;
+
+    QuadExpr<Var, double> m_dual_objective;
+    std::vector<std::optional<Var>> m_dual_variables_for_constraints; // constraint index <-> dual variable
+    std::vector<std::optional<Var>> m_dual_variables_for_qconstraints; // qconstraint index <-> dual variable
+    std::vector<std::optional<Var>> m_dual_variables_for_lower_bounds; // variable index <-> dual variable
+    std::vector<std::optional<Var>> m_dual_variables_for_upper_bounds; // variable index <-> dual variable
+    std::vector<std::optional<std::variant<Ctr, QCtr>>> m_dual_constraints; // variable index <-> dual constraint
+
+    [[nodiscard]] bool primal_is_a_linear_problem() const;
+    [[nodiscard]] bool is_linear(const QuadExpr<Var, double>& t_expr) const;
+    void create_dual_variables();
+    [[nodiscard]] std::pair<double, double> bounds_for_dual_variable(idol::CtrType t_type) const;
+    void create_dual_constraints();
 public:
-    explicit Dualizer(const Model& t_parent);
+    explicit Dualizer(const Model& t_parent,
+                      const QuadExpr<Var, double>& t_primal_objective,
+                      const std::function<bool(const Var&)>& t_primal_variable_indicator = [](const Var&) { return true; },
+                      const std::function<bool(const Ctr&)>& t_primal_constraint_indicator = [](const Ctr&) { return true; },
+                      const std::function<bool(const QCtr&)>& t_primal_qconstraint_indicator = [](const QCtr&) { return true; });
+
+    const QuadExpr<Var, double>& get_dual_obj_expr() const { return m_dual_objective; }
+
+    void add_coupling(Model& t_destination);
+    void add_coupling_variables(Model& t_destination);
+    void add_coupling_constraints(Model& t_destination);
+
+    void add_dual(Model& t_destination);
+    void add_dual_variables(Model& t_destination);
+    void add_dual_constraints(Model& t_destination);
+    void add_dual_objective(Model& t_destination);
+
+    void add_primal(Model& t_destination);
+    void add_primal_variables(idol::Model &t_destination);
+    void add_primal_constraints(idol::Model &t_destination);
+    void add_primal_objective(idol::Model &t_destination);
+
+    void add_strong_duality_reformulation(Model& t_destination);
+    void add_kkt_reformulation(idol::Model &t_destination);
 };
 
 #endif //IDOL_DUALIZER_H