use more than one separator in CCG

examples/robust/ccg-discrete-uncertainty.example.cpp

@@ -23,6 +23,10 @@
 #include "idol/robust/optimizers/column-and-constraint-generation/ColumnAndConstraintGeneration.h"
 #include "idol/bilevel/optimizers/wrappers/MibS/MibS.h"
 #include "idol/bilevel/optimizers/KKT/KKT.h"
+#include "idol/mixed-integer/optimizers/padm/SubProblem.h"
+#include "idol/mixed-integer/optimizers/padm/PenaltyUpdates.h"
+#include "idol/bilevel/optimizers/StrongDuality/StrongDuality.h"
+#include "idol/mixed-integer/optimizers/padm/PADM.h"
 
 using namespace idol;
 
@@ -52,7 +56,7 @@ int main(int t_argc, const char** t_argv) {
     Bilevel::Description bilevel_description(env);
 
     auto x = model.add_vars(Dim<1>(n_facilities), 0., 1., Binary, 0., "x");
-    auto y = model.add_vars(Dim<2>(n_facilities, n_customers), 0., 1., Binary, 0., "y");
+    auto y = model.add_vars(Dim<2>(n_facilities, n_customers), 0., 1., Continuous, 0., "y");
 
     for (unsigned int i = 0 ; i < n_facilities ; ++i) {
         const auto c = model.add_ctr(
@@ -97,6 +101,27 @@ int main(int t_argc, const char** t_argv) {
         }
     }
 
+    Annotation<unsigned int> decomposition(env, "sub_problem", 0);
+    for (unsigned int i = 0 ; i < n_facilities ; ++i) {
+        for (unsigned int j = 0; j < n_customers; ++j) {
+            xi[i][j].set(decomposition, 1);
+        }
+    }
+
+    Annotation<double> initial_penalties(env, "initial_penalties", 1e1);
+
+    const auto padm = Bilevel::StrongDuality()
+                .with_single_level_optimizer(
+                    PADM(decomposition, initial_penalties)
+                        .with_default_sub_problem_spec(
+                            ADM::SubProblem().with_optimizer(Gurobi())
+                        )
+                        .with_penalty_update(PenaltyUpdates::Multiplicative(2))
+                        .with_rescaling_threshold(1e4)
+                        .with_logs(false)
+                )
+            ;
+
     const auto mibs = Bilevel::MibS()
             .with_cplex_for_feasibility(true)
             .with_file_interface(false)
@@ -107,8 +132,13 @@ int main(int t_argc, const char** t_argv) {
             Robust::ColumnAndConstraintGeneration(robust_description, bilevel_description)
                     //.with_initial_scenario_by_maximization(Gurobi())
                     //.with_initial_scenario_by_minimization(Gurobi())
+
                     .with_master_optimizer(Gurobi())
+
+                    .add_feasibility_separation_optimizer(padm)
                     .add_feasibility_separation_optimizer(mibs)
+
+                    .add_optimality_separation_optimizer(padm)
                     .add_optimality_separation_optimizer(mibs)
                     .with_logs(true)
     );

lib/include/idol/robust/optimizers/column-and-constraint-generation/ColumnAndConstraintGeneration.h

@@ -24,8 +24,8 @@ class idol::Robust::ColumnAndConstraintGeneration : public OptimizerFactoryWithD
     std::vector<Point<Var>> m_initial_scenarios;
 
     // Separation
-    std::unique_ptr<OptimizerFactory> m_optimizer_feasibility_separation;
-    std::unique_ptr<OptimizerFactory> m_optimizer_optimality_separation;
+    std::list<std::unique_ptr<OptimizerFactory>> m_optimizer_feasibility_separation;
+    std::list<std::unique_ptr<OptimizerFactory>> m_optimizer_optimality_separation;
 public:
     ColumnAndConstraintGeneration(const Robust::Description& t_robust_description,
                                   const Bilevel::Description& t_bilevel_description);

lib/include/idol/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.h

@@ -29,10 +29,12 @@ class idol::Optimizers::Robust::ColumnAndConstraintGeneration : public Algorithm
     std::unique_ptr<OptimizerFactory> m_initial_scenario_by_maximization;
 
     // Feasibility Separation
-    std::unique_ptr<OptimizerFactory> m_optimizer_feasibility_separation;
+    std::vector<std::unique_ptr<OptimizerFactory>> m_optimizer_feasibility_separation;
+    unsigned int m_index_feasibility_separation = 0;
 
     // Optimality Separation
-    std::unique_ptr<OptimizerFactory> m_optimizer_optimality_separation;
+    std::vector<std::unique_ptr<OptimizerFactory>> m_optimizer_optimality_separation;
+    unsigned int m_index_optimality_separation = 0;
 public:
     ColumnAndConstraintGeneration(const Model& t_parent,
                                   const ::idol::Robust::Description &t_robust_description,
@@ -41,8 +43,8 @@ class idol::Optimizers::Robust::ColumnAndConstraintGeneration : public Algorithm
                                   std::vector<Point<Var>> t_initial_scenarios,
                                   OptimizerFactory* t_initial_scenario_by_minimization,
                                   OptimizerFactory* t_initial_scenario_by_maximization,
-                                  OptimizerFactory* t_optimizer_feasibility_separation,
-                                  OptimizerFactory* t_optimizer_optimality_separation);
+                                  const std::list<std::unique_ptr<OptimizerFactory>>& t_optimizer_feasibility_separation,
+                                  const std::list<std::unique_ptr<OptimizerFactory>>& t_optimizer_optimality_separation);
 
     [[nodiscard]] std::string name() const override;
 

lib/src/bilevel/optimizers/StrongDuality/Optimizers_StrongDuality.cpp

@@ -15,7 +15,7 @@ idol::Optimizers::Bilevel::StrongDuality::StrongDuality(const Model& t_parent,
 }
 
 std::string idol::Optimizers::Bilevel::StrongDuality::name() const {
-    return "deterministic";
+    return "strong-duality reformulation";
 }
 
 double idol::Optimizers::Bilevel::StrongDuality::get_var_primal(const idol::Var &t_var) const {

lib/src/mixed-integer/optimizers/padm/Formulation.cpp

@@ -435,6 +435,9 @@ void idol::ADM::Formulation::dispatch_qctr(const idol::Model &t_src_model,
 
         if (full_row->has_quadratic()) {
             // we have to add a quadratic constraint here
+            for (const auto& [var, expr] : *full_row) {
+                std::cout << var << " <-> " << expr << std::endl;
+            }
             throw Exception("Quadratic constraints in fixed problems are not implemented");
         }
 

lib/src/robust/optimizers/column-and-constraint-generation/ColumnAndConstraintGeneration.cpp

@@ -19,9 +19,16 @@ idol::Robust::ColumnAndConstraintGeneration::ColumnAndConstraintGeneration(
           m_bilevel_description(t_src.m_bilevel_description),
           m_initial_scenarios(t_src.m_initial_scenarios),
           m_initial_scenario_by_minimization(t_src.m_initial_scenario_by_minimization ? t_src.m_initial_scenario_by_minimization->clone() : nullptr),
-          m_initial_scenario_by_maximization(t_src.m_initial_scenario_by_maximization ? t_src.m_initial_scenario_by_maximization->clone() : nullptr),
-          m_optimizer_feasibility_separation(t_src.m_optimizer_feasibility_separation ? t_src.m_optimizer_feasibility_separation->clone() : nullptr),
-          m_optimizer_optimality_separation(t_src.m_optimizer_optimality_separation ? t_src.m_optimizer_optimality_separation->clone() : nullptr) {
+          m_initial_scenario_by_maximization(t_src.m_initial_scenario_by_maximization ? t_src.m_initial_scenario_by_maximization->clone() : nullptr)
+          {
+
+    for (const auto& optimizer : t_src.m_optimizer_feasibility_separation) {
+        m_optimizer_feasibility_separation.emplace_back(optimizer->clone());
+    }
+
+    for (const auto& optimizer : t_src.m_optimizer_optimality_separation) {
+        m_optimizer_optimality_separation.emplace_back(optimizer->clone());
+    }
 
 }
 
@@ -47,7 +54,7 @@ idol::Optimizer *idol::Robust::ColumnAndConstraintGeneration::operator()(const i
         throw Exception("Master optimizer not set");
     }
 
-    if (!m_optimizer_feasibility_separation && !m_optimizer_optimality_separation) {
+    if (m_optimizer_feasibility_separation.empty() && m_optimizer_optimality_separation.empty()) {
         throw Exception("At least one of feasibility or optimality separation optimizers must be set");
     }
 
@@ -58,8 +65,8 @@ idol::Optimizer *idol::Robust::ColumnAndConstraintGeneration::operator()(const i
                                                                          m_initial_scenarios,
                                                                          m_initial_scenario_by_minimization ? m_initial_scenario_by_minimization->clone() : nullptr,
                                                                          m_initial_scenario_by_maximization ? m_initial_scenario_by_maximization->clone() : nullptr,
-                                                                         m_optimizer_feasibility_separation ? m_optimizer_feasibility_separation->clone() : nullptr,
-                                                                         m_optimizer_optimality_separation ? m_optimizer_optimality_separation->clone() : nullptr
+                                                                         m_optimizer_feasibility_separation,
+                                                                         m_optimizer_optimality_separation
                                                                          );
 
     handle_default_parameters(result);
@@ -107,15 +114,11 @@ idol::Robust::ColumnAndConstraintGeneration &
 idol::Robust::ColumnAndConstraintGeneration::add_feasibility_separation_optimizer(
         const OptimizerFactory &t_optimizer) {
 
-    if (m_optimizer_feasibility_separation) {
-        throw Exception("Feasibility separation optimizer already set");
-    }
-
     if (!t_optimizer.is<Bilevel::OptimizerInterface>()) {
         throw Exception("Feasibility separation optimizer must be a bilevel optimizer");
     }
 
-    m_optimizer_feasibility_separation.reset(t_optimizer.clone());
+    m_optimizer_feasibility_separation.emplace_back(t_optimizer.clone());
 
     return *this;
 }
@@ -124,15 +127,11 @@ idol::Robust::ColumnAndConstraintGeneration &
 idol::Robust::ColumnAndConstraintGeneration::add_optimality_separation_optimizer(
         const OptimizerFactory &t_optimizer) {
 
-    if (m_optimizer_optimality_separation) {
-        throw Exception("Optimality separation optimizer already set");
-    }
-
     if (!t_optimizer.is<Bilevel::OptimizerInterface>()) {
         throw Exception("Optimality separation optimizer must be a bilevel optimizer");
     }
 
-    m_optimizer_optimality_separation.reset(t_optimizer.clone());
+    m_optimizer_optimality_separation.emplace_back(t_optimizer.clone());
 
     return *this;
 }

lib/src/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.cpp

@@ -14,17 +14,25 @@ idol::Optimizers::Robust::ColumnAndConstraintGeneration::ColumnAndConstraintGene
                                                                                        std::vector<Point<Var>> t_initial_scenarios,
                                                                                        OptimizerFactory* t_initial_scenario_by_minimization,
                                                                                        OptimizerFactory* t_initial_scenario_by_maximization,
-                                                                                       OptimizerFactory* t_optimizer_feasibility_separation,
-                                                                                       OptimizerFactory* t_optimizer_optimality_separation)
+                                                                                       const std::list<std::unique_ptr<OptimizerFactory>>& t_optimizer_feasibility_separation,
+                                                                                       const std::list<std::unique_ptr<OptimizerFactory>>& t_optimizer_optimality_separation)
                                                                                        : Algorithm(t_parent),
                                                                                          m_robust_description(t_robust_description),
                                                                                          m_bilevel_description(t_bilevel_description),
                                                                                          m_master_optimizer(t_master_optimizer.clone()),
                                                                                          m_initial_scenario_by_minimization(t_initial_scenario_by_minimization),
                                                                                          m_initial_scenario_by_maximization(t_initial_scenario_by_maximization),
-                                                                                         m_initial_scenarios(std::move(t_initial_scenarios)),
-                                                                                         m_optimizer_feasibility_separation(t_optimizer_feasibility_separation),
-                                                                                         m_optimizer_optimality_separation(t_optimizer_optimality_separation) {
+                                                                                         m_initial_scenarios(std::move(t_initial_scenarios)) {
+
+    m_optimizer_feasibility_separation.reserve(t_optimizer_feasibility_separation.size());
+    for (const auto& optimizer : t_optimizer_feasibility_separation) {
+        m_optimizer_feasibility_separation.emplace_back(optimizer->clone());
+    }
+
+    m_optimizer_optimality_separation.reserve(t_optimizer_optimality_separation.size());
+    for (const auto& optimizer : t_optimizer_optimality_separation) {
+        m_optimizer_optimality_separation.emplace_back(optimizer->clone());
+    }
 
 }
 
@@ -100,6 +108,9 @@ void idol::Optimizers::Robust::ColumnAndConstraintGeneration::hook_before_optimi
                                                        m_bilevel_description);
     m_formulation->master().use(*m_master_optimizer);
 
+    m_index_feasibility_separation = 0;
+    m_index_optimality_separation = 0;
+
 }
 
 void idol::Optimizers::Robust::ColumnAndConstraintGeneration::hook_optimize() {
@@ -284,14 +295,24 @@ void idol::Optimizers::Robust::ColumnAndConstraintGeneration::log_iteration() {
 
 void idol::Optimizers::Robust::ColumnAndConstraintGeneration::solve_adversarial_problem() {
 
-    if (m_optimizer_feasibility_separation) {
-        const bool n_added_scenario = solve_feasibility_adversarial_problem();
+    const unsigned int n_feasibility_separation_optimizers = m_optimizer_feasibility_separation.size();
+
+    do {
+        const unsigned int n_added_scenario = solve_feasibility_adversarial_problem();
         if (n_added_scenario > 0) {
             return;
         }
-    }
+    } while (m_index_feasibility_separation + 1 != n_feasibility_separation_optimizers);
 
-    solve_optimality_adversarial_problem();
+
+    const unsigned int n_optimality_separation_optimizers = m_optimizer_optimality_separation.size();
+
+    do {
+        const unsigned int n_added_scenario = solve_optimality_adversarial_problem();
+        if (n_added_scenario > 0) {
+            return;
+        }
+    } while (m_index_optimality_separation + 1 != n_optimality_separation_optimizers);
 
 }
 
@@ -304,18 +325,33 @@ unsigned int idol::Optimizers::Robust::ColumnAndConstraintGeneration::solve_feas
 
     // Set bilevel description
     const auto& separation_bilevel_description = m_formulation->separation_bilevel_description();
-    m_optimizer_feasibility_separation->as<Bilevel::OptimizerInterface>().set_bilevel_description(separation_bilevel_description);
+    m_optimizer_feasibility_separation[m_index_feasibility_separation]->as<Bilevel::OptimizerInterface>().set_bilevel_description(separation_bilevel_description);
 
     // Set optimizer
-    high_point_relaxation.use(*m_optimizer_feasibility_separation);
+    high_point_relaxation.use(*m_optimizer_feasibility_separation[m_index_feasibility_separation]);
     high_point_relaxation.optimizer().set_param_time_limit(get_remaining_time());
+
+    // Solve adversarial problem
+    std::cout << "Solving feasibility adversarial problem with " << high_point_relaxation.optimizer().name() << std::endl;
     high_point_relaxation.optimize();
-    Bilevel::write_to_file(high_point_relaxation, separation_bilevel_description, "feasibility_separation");
 
     // Analyze results
     const auto status = high_point_relaxation.get_status();
 
+    const bool is_last_optimizer = m_index_feasibility_separation == m_optimizer_feasibility_separation.size() - 1;
+
     if (status != Optimal && status != Feasible) {
+        if (!is_last_optimizer) { // If we can, skip this optimizer
+            ++m_index_feasibility_separation;
+        } else { // otherwise, it's a fail
+            set_status(Fail);
+            set_reason(high_point_relaxation.get_reason());
+            terminate();
+            return 1;
+        }
+    }
+
+    if (is_last_optimizer && status != Optimal) { // if the last optimizer is not reporting optimal, it's a fail
         set_status(Fail);
         set_reason(high_point_relaxation.get_reason());
         terminate();
@@ -330,6 +366,10 @@ unsigned int idol::Optimizers::Robust::ColumnAndConstraintGeneration::solve_feas
         return 1;
     }
 
+    if (status == Feasible) {
+        ++m_index_feasibility_separation;
+    }
+
     return 0;
 }
 
@@ -371,17 +411,33 @@ idol::Optimizers::Robust::ColumnAndConstraintGeneration::solve_optimality_advers
 
     // Set bilevel description
     const auto& separation_bilevel_description = m_formulation->separation_bilevel_description();
-    m_optimizer_optimality_separation->as<Bilevel::OptimizerInterface>().set_bilevel_description(separation_bilevel_description);
+    m_optimizer_optimality_separation[m_index_optimality_separation]->as<Bilevel::OptimizerInterface>().set_bilevel_description(separation_bilevel_description);
 
     // Set optimizer
-    high_point_relaxation.use(*m_optimizer_optimality_separation);
+    high_point_relaxation.use(*m_optimizer_optimality_separation[m_index_optimality_separation]);
     high_point_relaxation.optimizer().set_param_time_limit(get_remaining_time());
+
+    // Solve adversarial problem
+    std::cout << "Solving optimality adversarial problem with " << high_point_relaxation.optimizer().name() << std::endl;
     high_point_relaxation.optimize();
 
     // Analyze results
     const auto status = high_point_relaxation.get_status();
 
+    const bool is_last_optimizer = m_index_optimality_separation == m_optimizer_optimality_separation.size() - 1;
+
     if (status != Optimal && status != Feasible) {
+        if (!is_last_optimizer) { // If we can, skip this optimizer
+            ++m_index_feasibility_separation;
+        } else { // otherwise, it's a fail
+            set_status(Fail);
+            set_reason(high_point_relaxation.get_reason());
+            terminate();
+            return 1;
+        }
+    }
+
+    if (is_last_optimizer && status != Optimal) { // if the last optimizer is not reporting optimal, it's a fail
         set_status(Fail);
         set_reason(high_point_relaxation.get_reason());
         terminate();
@@ -390,7 +446,9 @@ idol::Optimizers::Robust::ColumnAndConstraintGeneration::solve_optimality_advers
 
     bool add_scenario;
     if (t_coupling_constraint_index == 0) {
-        add_scenario = true;
+        const double LB = m_formulation->master().get_best_obj();
+        const double UB = -high_point_relaxation.get_best_obj();
+        add_scenario = absolute_gap(LB, UB) > Tolerance::MIPAbsoluteGap || relative_gap(LB, UB) > Tolerance::MIPRelativeGap;
         if (status == Optimal) {
             set_best_obj(std::min(get_best_obj(), -high_point_relaxation.get_best_obj()));
         }
@@ -405,5 +463,9 @@ idol::Optimizers::Robust::ColumnAndConstraintGeneration::solve_optimality_advers
         return 1;
     }
 
+    if (status == Feasible) {
+        ++m_index_optimality_separation;
+    }
+
     return 0;
 }