allow for independent updates

examples/bilevel-optimization/padm.example.cpp

@@ -50,7 +50,7 @@ int main(int t_argc, const char** t_argv) {
     single_level.use(
             PADM(decomposition, penalized_constraints)
                 .with_default_sub_problem_spec(ADM::SubProblem().with_optimizer(Gurobi()))
-                .with_penalty_update(PenaltyUpdates::Additive(1))
+                .with_penalty_update(PenaltyUpdates::Multiplicative(2))
                 .with_rescaling(true, 1e3)
     );
 

lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.cpp

@@ -10,9 +10,13 @@
 
 idol::ADM::Formulation::Formulation(const Model& t_src_model,
                                     Annotation<Var, unsigned int> t_decomposition,
-                                    std::optional<Annotation<Ctr, bool>> t_penalized_constraints)
-                                                     : m_decomposition(std::move(t_decomposition)),
-                                                       m_penalized_constraints(std::move(t_penalized_constraints)) {
+                                    std::optional<Annotation<Ctr, bool>> t_penalized_constraints,
+                                    bool t_independent_penalty_update,
+                                    std::pair<bool, double> t_rescaling)
+                                 : m_decomposition(std::move(t_decomposition)),
+                                   m_penalized_constraints(std::move(t_penalized_constraints)),
+                                   m_independent_penalty_update(t_independent_penalty_update),
+                                   m_rescaling(std::move(t_rescaling)) {
 
     const auto n_sub_problems = compute_n_sub_problems(t_src_model);
 
@@ -56,7 +60,7 @@ void idol::ADM::Formulation::initialize_patterns(const idol::Model &t_src_model,
 
 void idol::ADM::Formulation::initialize_slacks(const idol::Model &t_src_model,
                                                unsigned int n_sub_problems) {
-    m_l1_vars.resize(n_sub_problems);
+    m_l1_vars_in_sub_problem.resize(n_sub_problems);
 }
 
 void idol::ADM::Formulation::dispatch_vars(const idol::Model &t_src_model) {
@@ -102,16 +106,22 @@ void idol::ADM::Formulation::dispatch_ctr(const idol::Model &t_src_model, const
         auto& model = m_sub_problems[t_sub_problem_id];
 
         const auto add_l1_var = [&](double t_coefficient) {
-            auto var = model.add_var(0, Inf, Continuous, Column(), "l1_norm_" + std::to_string(t_coefficient < 0) + "_" + t_ctr.name());
+            auto var = get_or_create_l1_var(t_ctr);
+            model.add(var);
             pattern.linear() += t_coefficient * var;
-            m_l1_vars[t_sub_problem_id].emplace_back(var);
+            m_l1_vars_in_sub_problem[t_sub_problem_id].emplace_back(var);
+            return var;
         };
 
         switch (type) {
-            case Equal:
-                add_l1_var(+1);
-                add_l1_var(-1);
+            case Equal: {
+                auto var = add_l1_var(0);
+                auto var_minus = model.add_var(0, Inf, Continuous, var.name() + "_minus");
+                auto var_plus = model.add_var(0, Inf, Continuous, var.name() + "_plus");
+                model.add_ctr(var == var_plus + var_minus);
+                pattern.linear() += var_plus - var_minus;
                 break;
+            }
             case LessOrEqual:
                 add_l1_var(-1);
                 break;
@@ -290,9 +300,83 @@ void
 idol::ADM::Formulation::update_penalty_parameters(const std::vector<Solution::Primal> &t_primals,
                                                   PenaltyUpdate& t_penalty_update) {
 
+    const unsigned int n_sub_problems = m_sub_problems.size();
+
+    if (m_independent_penalty_update) {
+        update_penalty_parameters_independently(t_primals, t_penalty_update);
+        return;
+    }
+
+    for (auto &[ctr, penalty]: m_l1_vars) {
+
+        auto &[var, value] = penalty;
+
+        if (value < 1e-3) {
+            value = 1;
+            set_penalty_in_all_sub_problems(var, 1);
+            continue;
+        }
+
+        unsigned int argmax = -1;
+        double max = 0.;
+
+        for (unsigned int i = 0; i < n_sub_problems; ++i) {
+
+            if (const double val = t_primals[i].get(var); val > max) {
+                max = val;
+                argmax = i;
+            }
+
+        }
+
+        if (argmax == -1 || max <= 1e-4) {
+            continue;
+        }
+
+        value = t_penalty_update(value);
+        set_penalty_in_all_sub_problems(var, value);
+
+    }
+
+    if (m_rescaling.first) {
+        rescale_penalty_parameters();
+    }
+
+}
+
+idol::Var idol::ADM::Formulation::get_or_create_l1_var(const idol::Ctr &t_ctr) {
+
+    auto it = std::lower_bound(m_l1_vars.begin(), m_l1_vars.end(), t_ctr, [](const auto& t_lhs, const auto& t_rhs) {
+        return t_lhs.first.id() < t_rhs.id();
+    });
+
+    if (it != m_l1_vars.end() && it->first.id() == t_ctr.id()) {
+        return it->second.first;
+    }
+
+    auto& env = m_sub_problems.front().env();
+    Var var (env, 0, Inf, Continuous, Column(), "l1_norm_" + t_ctr.name());
+    m_l1_vars.emplace_hint(it, t_ctr, std::make_pair( var, 0. ));
+
+    return var;
+}
+
+void idol::ADM::Formulation::set_penalty_in_all_sub_problems(const Var &t_var, double t_value) {
+
+    for (auto& model : m_sub_problems) {
+        if (model.has(t_var)) {
+            model.set_var_obj(t_var, t_value);
+        }
+    }
+
+}
+
+void idol::ADM::Formulation::update_penalty_parameters_independently(const std::vector<Solution::Primal> &t_primals,
+                                                                     idol::PenaltyUpdate &t_penalty_update) {
+
     for (unsigned int i = 0, n_sub_problems = m_sub_problems.size() ; i < n_sub_problems ; ++i) {
         auto& model = m_sub_problems[i];
-        for (const auto& var : m_l1_vars[i]) {
+        for (const auto& var : m_l1_vars_in_sub_problem[i]) {
 
             double current_penalty = model.get_var_column(var).obj().as_numerical();
 
@@ -309,3 +393,35 @@ idol::ADM::Formulation::update_penalty_parameters(const std::vector<Solution::Pr
     }
 
 }
+
+void idol::ADM::Formulation::rescale_penalty_parameters() {
+
+    double max = 0;
+    for (const auto& [ctr, penalty] : m_l1_vars) {
+        max = std::max(max, penalty.second);
+    }
+
+    if (max < m_rescaling.second) {
+        return;
+    }
+
+    const auto sigmoid = [&](double t_val) {
+        const double alpha = max;
+        constexpr double beta = 5;
+        const double sigma = max;
+        constexpr double omega = 5;
+        return beta * (t_val - sigma) / (alpha + std::abs(t_val - sigma)) + omega;
+    };
+
+    const unsigned int n_sub_problems = m_sub_problems.size();
+
+    for (unsigned int i = 0 ; i < n_sub_problems ; ++i) {
+
+        for (const auto& var : m_l1_vars_in_sub_problem[i]) {
+            const double current_penalty = m_sub_problems[i].get_var_column(var).obj().as_numerical();
+            m_sub_problems[i].set_var_obj(var, sigmoid(current_penalty));
+        }
+
+    }
+
+}

lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.h

@@ -16,7 +16,9 @@ class idol::ADM::Formulation {
 public:
     Formulation(const Model& t_src_model,
                 Annotation<Var, unsigned int> t_decomposition,
-                std::optional<Annotation<Ctr, bool>> t_penalized_constraints);
+                std::optional<Annotation<Ctr, bool>> t_penalized_constraints,
+                bool t_independent_penalty_update,
+                std::pair<bool, double> t_rescaling);
 
     Model& sub_problem(const Var& t_var);
 
@@ -32,7 +34,7 @@ class idol::ADM::Formulation {
 
     auto sub_problems() const { return ConstIteratorForward(m_sub_problems); }
 
-    const auto l1_vars(unsigned int t_sub_problem_id) const { return ConstIteratorForward(m_l1_vars[t_sub_problem_id]); }
+    const auto l1_vars(unsigned int t_sub_problem_id) const { return ConstIteratorForward(m_l1_vars_in_sub_problem[t_sub_problem_id]); }
 
     bool has_penalized_constraints() const { return m_penalized_constraints.has_value(); }
 
@@ -42,11 +44,14 @@ class idol::ADM::Formulation {
 private:
     Annotation<Var, unsigned int> m_decomposition;
     std::optional<Annotation<Ctr, bool>> m_penalized_constraints;
+    bool m_independent_penalty_update;
+    std::pair<bool, double> m_rescaling;
 
     std::vector<Model> m_sub_problems;
     std::vector<std::optional<Expr<Var, Var>>> m_objective_patterns;
     std::vector<std::list<std::pair<Ctr, Expr<Var, Var>>>> m_constraint_patterns; // as ctr: lhs <= 0
-    std::vector<std::list<Var>> m_l1_vars;
+    std::vector<std::list<Var>> m_l1_vars_in_sub_problem;
+    Map<Ctr, std::pair<Var, double>> m_l1_vars;
 
     unsigned int compute_n_sub_problems(const Model& t_src_model) const;
     void initialize_sub_problems(const Model& t_src_model, unsigned int n_sub_problems);
@@ -58,6 +63,10 @@ class idol::ADM::Formulation {
     void dispatch_obj(const Model& t_src_model);
     void dispatch_obj(const Model& t_src_model, unsigned int t_sub_problem_id);
     std::pair<Expr<Var, Var>, bool> dispatch(const Model& t_src_model, const LinExpr<Var>& t_lin_expr, const QuadExpr<Var, Var>& t_quad_expr, unsigned int t_sub_problem_id);
+    Var get_or_create_l1_var(const Ctr& t_ctr);
+    void set_penalty_in_all_sub_problems(const Var& t_var, double t_value);
+    void update_penalty_parameters_independently(const std::vector<Solution::Primal>& t_primals, PenaltyUpdate& t_penalty_update);
+    void rescale_penalty_parameters();
 
     double fix(const Constant& t_constant, const std::vector<Solution::Primal>& t_primals);
 };

lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_PADM.cpp

@@ -9,12 +9,10 @@
 idol::Optimizers::PADM::PADM(const Model& t_model,
                              ADM::Formulation t_formulation,
                              std::vector<idol::ADM::SubProblem>&& t_sub_problem_specs,
-                             std::pair<bool, double> t_rescaling,
                              PenaltyUpdate* t_penalty_update)
                             : Algorithm(t_model),
                               m_formulation(std::move(t_formulation)),
                               m_sub_problem_specs(std::move(t_sub_problem_specs)),
-                              m_rescaling(std::move(t_rescaling)),
                               m_penalty_update(t_penalty_update) {
 
 }

lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_PADM.h

@@ -19,7 +19,6 @@ class idol::Optimizers::PADM : public Algorithm {
     PADM(const Model& t_model,
          ADM::Formulation t_formulation,
          std::vector<idol::ADM::SubProblem>&& t_sub_problem_specs,
-         std::pair<bool, double> t_rescaling,
          PenaltyUpdate* t_penalty_update);
 
     std::string name() const override { return "PADM"; }
@@ -87,7 +86,6 @@ class idol::Optimizers::PADM : public Algorithm {
 private:
     ADM::Formulation m_formulation;
     std::vector<idol::ADM::SubProblem> m_sub_problem_specs;
-    std::pair<bool, double> m_rescaling;
     std::unique_ptr<PenaltyUpdate> m_penalty_update;
     unsigned int m_max_inner_loop_iterations = 1000;
 

lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.cpp

@@ -36,13 +36,19 @@ idol::Optimizer *idol::PADM::operator()(const idol::Model &t_model) const {
         throw Exception("The decomposition has not been set.");
     }
 
-    if (!m_penalized_constraints && (m_rescaling || m_penalty_update)) {
+    if (!m_penalized_constraints && (m_rescaling || m_penalty_update || m_independent_penalty_update)) {
         std::cout << "Warning: The penalized constraints have not been set. The rescaling and penalty update will be ignored." << std::endl;
     }
 
+    if (m_rescaling && m_rescaling->first && m_independent_penalty_update && m_independent_penalty_update.value()) {
+        throw Exception("The rescaling and independent penalty update cannot be used together.");
+    }
+
     ADM::Formulation formulation(t_model,
                                  *m_decomposition,
-                                 m_penalized_constraints);
+                                 m_penalized_constraints,
+                                 m_independent_penalty_update && *m_independent_penalty_update,
+                                 m_rescaling ? *m_rescaling : std::make_pair(false, 0.));
 
     auto sub_problem_specs = create_sub_problem_specs(t_model, formulation);
 
@@ -55,7 +61,6 @@ idol::Optimizer *idol::PADM::operator()(const idol::Model &t_model) const {
                 t_model,
                 std::move(formulation),
                 std::move(sub_problem_specs),
-                m_rescaling ? *m_rescaling : std::make_pair(false, 0.),
                 penalty_update
             );
 
@@ -115,3 +120,14 @@ idol::PADM::PADM(const idol::PADM &t_src)
       m_penalty_update(t_src.m_penalty_update ? t_src.m_penalty_update->clone() : nullptr) {
 
 }
+
+idol::PADM &idol::PADM::with_independent_penalty_update(bool t_value) {
+
+    if (m_independent_penalty_update) {
+        throw Exception("The independent penalty update has already been set.");
+    }
+
+    m_independent_penalty_update = t_value;
+
+    return *this;
+}

lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.h

@@ -35,6 +35,8 @@ class idol::PADM : public OptimizerFactoryWithDefaultParameters<PADM> {
 
     PADM& with_penalty_update(const PenaltyUpdate& t_penalty_update);
 
+    PADM& with_independent_penalty_update(bool t_value);
+
     Optimizer *operator()(const Model &t_model) const override;
 
     OptimizerFactory *clone() const override;
@@ -46,6 +48,7 @@ class idol::PADM : public OptimizerFactoryWithDefaultParameters<PADM> {
     Map<unsigned int, ADM::SubProblem> m_sub_problem_specs;
     std::optional<std::pair<bool, double>> m_rescaling;
     std::unique_ptr<PenaltyUpdate> m_penalty_update;
+    std::optional<bool> m_independent_penalty_update;
 
     std::vector<ADM::SubProblem> create_sub_problem_specs(const Model& t_model, const ADM::Formulation& t_formulation) const;
 };