implement restart and adaptive penalty updates in padm

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

@@ -5,6 +5,8 @@
 #ifndef IDOL_ADM_FORMULATION_H
 #define IDOL_ADM_FORMULATION_H
 
+#include <utility>
+
 #include "idol/modeling/models/Model.h"
 
 namespace idol {
@@ -38,23 +40,23 @@ 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_in_sub_problem[t_sub_problem_id]); }
+    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(); }
 
     void fix_sub_problem(unsigned int t_sub_problem_id, const std::vector<Solution::Primal>& t_primals);
 
     void initialize_penalty_parameters(double t_value);
 
-    void update_penalty_parameters(const std::vector<Solution::Primal>& t_primals, PenaltyUpdate& t_penalty_update);
+    bool update_penalty_parameters(const std::vector<Solution::Primal>& t_primals, PenaltyUpdate& t_penalty_update); // Returns true if penalty parameters have been resacled
 
     struct CurrentPenalty {
         const Ctr constraint;
         const Var variable;
         const double max_violation;
         double penalty;
         CurrentPenalty(Ctr t_constraint, Var t_variable, double t_max_violation, double t_penalty)
-            : constraint(std::move(t_constraint)), variable(t_variable), max_violation(t_max_violation), penalty(t_penalty) {}
+            : constraint(std::move(t_constraint)), variable(std::move(t_variable)), max_violation(t_max_violation), penalty(t_penalty) {}
     };
 
 private:
@@ -82,7 +84,7 @@ class idol::ADM::Formulation {
     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(std::list<CurrentPenalty>& t_penalties);
+    bool rescale_penalty_parameters(std::list<CurrentPenalty>& t_penalties);
 
     double fix(const Constant& t_constant, const std::vector<Solution::Primal>& t_primals);
 };

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

@@ -96,6 +96,8 @@ class idol::Optimizers::PADM : public Algorithm {
     double infeasibility_linf(unsigned int t_sub_problem_id, const Solution::Primal& t_solution) const;
     double infeasibility_l1(unsigned int t_sub_problem_id, const Solution::Primal& t_solution) const;
     void detect_stagnation(bool t_feasibility_has_changed);
+    void detect_stagnation_due_to_rescaling();
+    void restart();
 
     void check_feasibility();
     void check_time_limit();
@@ -108,15 +110,19 @@ class idol::Optimizers::PADM : public Algorithm {
     const std::vector<idol::ADM::SubProblem> m_sub_problem_specs;
     const std::unique_ptr<PenaltyUpdate> m_penalty_update;
     const double m_initial_penalty_parameter;
-    const unsigned int m_max_inner_loop_iterations = 1000;
+    const unsigned int m_max_inner_loop_iterations = std::numeric_limits<unsigned int>::max();
     const SolutionStatus m_feasible_solution_status;
-    const unsigned int m_max_iterations_without_feasibility_change = 500;
+    const unsigned int m_max_iterations_without_feasibility_change = 1000;
 
     std::optional<unsigned int> m_last_iteration_with_no_feasibility_change;
+    std::optional<double> m_last_objective_value_when_rescaled;
 
+    bool m_first_run = true;
+    unsigned int m_n_restart = 0;
     unsigned int m_outer_loop_iteration = 0;
     unsigned int m_inner_loop_iterations = 0;
     std::vector<Solution::Primal> m_last_solutions;
+    double m_current_initial_penalty_parameter = m_initial_penalty_parameter;
 
     struct IterationLog {
         unsigned int outer_iteration;

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

@@ -6,6 +6,8 @@
 #define IDOL_PENALTYUPDATES_H
 
 #include <list>
+#include <vector>
+#include <ostream>
 #include "Formulation.h"
 
 namespace idol {
@@ -14,7 +16,10 @@ namespace idol {
     namespace PenaltyUpdates {
         class Additive;
         class Multiplicative;
+        class Adaptive;
     }
+
+    std::ostream &operator<<(std::ostream &t_os, const PenaltyUpdate &t_penalty_update);
 }
 
 class idol::PenaltyUpdate {
@@ -25,7 +30,11 @@ class idol::PenaltyUpdate {
 
     virtual void operator()(std::list<ADM::Formulation::CurrentPenalty>& t_current_penalties);
 
-    virtual PenaltyUpdate* clone() const = 0;
+    virtual bool diversify() { return false; }
+
+    virtual std::ostream &describe(std::ostream &t_os) const = 0;
+
+    [[nodiscard]] virtual PenaltyUpdate* clone() const = 0;
 };
 
 class idol::PenaltyUpdates::Additive : public PenaltyUpdate {
@@ -37,27 +46,63 @@ class idol::PenaltyUpdates::Additive : public PenaltyUpdate {
         return t_current_penalty + m_increment;
     }
 
-    PenaltyUpdate* clone() const override {
+    std::ostream &describe(std::ostream &t_os) const override {
+        return t_os << "Additive(" << m_increment << ")";
+    }
+
+    [[nodiscard]] PenaltyUpdate* clone() const override {
         return new Additive(*this);
     }
 };
 
 class idol::PenaltyUpdates::Multiplicative : public PenaltyUpdate {
-    double m_factor;
+    std::vector<double> m_factor;
+    unsigned int m_current_factor_index = 0;
     bool m_normalized;
 public:
-    explicit Multiplicative(double t_factor, bool t_normalized = false) : m_factor(t_factor), m_normalized(t_normalized) {}
+    explicit Multiplicative(double t_factor, bool t_normalized = false) : m_factor({ t_factor }), m_normalized(t_normalized) {}
+
+    explicit Multiplicative(std::vector<double> t_factor, bool t_normalized = false) : m_factor(std::move(t_factor)), m_normalized(t_normalized) {}
 
     double operator()(double t_current_penalty) override {
-        return t_current_penalty * m_factor;
+        return t_current_penalty * m_factor[m_current_factor_index];
     }
 
     void operator()(std::list<ADM::Formulation::CurrentPenalty> &t_current_penalties) override;
 
-    PenaltyUpdate* clone() const override {
+    bool diversify() override;
+
+    std::ostream &describe(std::ostream &t_os) const override;
+
+    [[nodiscard]] PenaltyUpdate* clone() const override {
         return new Multiplicative(*this);
     }
 };
 
+class idol::PenaltyUpdates::Adaptive : public PenaltyUpdate {
+    std::vector<std::unique_ptr<PenaltyUpdate>> m_penalty_updates;
+    unsigned int m_current_penalty_update_index = 0;
+public:
+    explicit Adaptive(const std::vector<PenaltyUpdate*>& t_penalty_updates);
+
+    Adaptive(const Adaptive& t_src);
+
+    double operator()(double t_current_penalty) override {
+        return m_penalty_updates[m_current_penalty_update_index]->operator()(t_current_penalty);
+    }
+
+    void operator()(std::list<ADM::Formulation::CurrentPenalty> &t_current_penalties) override {
+        m_penalty_updates[m_current_penalty_update_index]->operator()(t_current_penalties);
+    }
+
+    bool diversify() override;
+
+    std::ostream &describe(std::ostream &t_os) const override;
+
+    [[nodiscard]] PenaltyUpdate* clone() const override {
+        return new Adaptive(*this);
+    }
+};
+
 
 #endif //IDOL_PENALTYUPDATES_H

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

@@ -301,15 +301,15 @@ const idol::Model &idol::ADM::Formulation::sub_problem(const idol::Var &t_var) c
     return m_sub_problems[t_var.get(m_decomposition)];
 }
 
-void
+bool
 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;
+        return false;
     }
 
     std::list<CurrentPenalty> current_penalties;
@@ -340,14 +340,16 @@ idol::ADM::Formulation::update_penalty_parameters(const std::vector<Solution::Pr
 
     t_penalty_update(current_penalties);
 
+    bool has_rescaled = false;
     if (m_rescaling.first) {
-        rescale_penalty_parameters(current_penalties);
+        has_rescaled = rescale_penalty_parameters(current_penalties);
     }
 
     for (const auto& [ctr, var, violation, penalty] : current_penalties) {
         set_penalty_in_all_sub_problems(var, penalty);
     }
 
+    return has_rescaled;
 }
 
 idol::Var idol::ADM::Formulation::get_or_create_l1_var(const idol::Ctr &t_ctr) {
@@ -412,15 +414,15 @@ void idol::ADM::Formulation::update_penalty_parameters_independently(const std::
 
 }
 
-void idol::ADM::Formulation::rescale_penalty_parameters(std::list<CurrentPenalty>& t_penalties) {
+bool idol::ADM::Formulation::rescale_penalty_parameters(std::list<CurrentPenalty>& t_penalties) {
 
     double max = 0;
     for (const auto& penalty : t_penalties) {
         max = std::max(max, penalty.penalty);
     }
 
     if (max < m_rescaling.second) {
-        return;
+        return false;
     }
 
     const auto sigmoid = [&](double t_val) {
@@ -431,10 +433,13 @@ void idol::ADM::Formulation::rescale_penalty_parameters(std::list<CurrentPenalty
         return beta * (t_val - sigma) / (alpha + std::abs(t_val - sigma)) + omega;
     };
 
+    max = 0;
     for (auto& penalty: t_penalties) {
         penalty.penalty = sigmoid(penalty.penalty);
+        max = std::max(max, penalty.penalty);
     }
 
+    return true;
 }
 
 unsigned int idol::ADM::Formulation::sub_problem_id(const idol::Var &t_var) const {

lib/src/optimizers/mixed-integer-optimization/padm/Optimizers_PADM.cpp

@@ -170,8 +170,12 @@ void idol::Optimizers::PADM::hook_before_optimize() {
 
     m_outer_loop_iteration = 0;
     m_inner_loop_iterations = 0;
+    m_first_run = true;
+    m_n_restart = 0;
     m_last_iteration_with_no_feasibility_change.reset();
+    m_last_objective_value_when_rescaled.reset();
     m_last_solutions = std::vector<Solution::Primal>(n_sub_problems);
+    m_current_initial_penalty_parameter = m_initial_penalty_parameter;
 
     for (unsigned int i = 0 ; i < n_sub_problems ; ++i) {
         auto& model = m_formulation.sub_problem(i);
@@ -304,6 +308,7 @@ void idol::Optimizers::PADM::run_inner_loop() {
             feasibility_has_changed |= feas;
         }
 
+        m_first_run = false;
         ++m_inner_loop_iterations;
 
         check_time_limit();
@@ -336,7 +341,11 @@ void idol::Optimizers::PADM::update_penalty_parameters() {
         return;
     }
 
-    m_formulation.update_penalty_parameters(m_last_solutions, *m_penalty_update);
+    bool has_rescaled = m_formulation.update_penalty_parameters(m_last_solutions, *m_penalty_update);
+
+    if (has_rescaled) {
+        detect_stagnation_due_to_rescaling();
+    }
 
 }
 
@@ -371,8 +380,8 @@ idol::Optimizers::PADM::solve_sub_problem(unsigned int t_sub_problem_id) {
 
     auto current_solution = save_primal(model);
     // bool obj_has_changed = m_inner_loop_iterations == 0 || (m_last_solutions[t_sub_problem_id] + -1. * current_solution).norm(2) > 1e-4;
-    const bool obj_has_changed = m_inner_loop_iterations == 0 || std::abs(m_last_solutions[t_sub_problem_id].objective_value() - current_solution.objective_value()) > 1e-4;
-    const bool feas_has_changed = m_inner_loop_iterations == 0 || std::abs(infeasibility_l1(t_sub_problem_id, m_last_solutions[t_sub_problem_id]) - infeasibility_l1(t_sub_problem_id, current_solution)) > 1e-5;
+    const bool obj_has_changed = m_first_run || std::abs(m_last_solutions[t_sub_problem_id].objective_value() - current_solution.objective_value()) > 1e-4;
+    const bool feas_has_changed = m_first_run || std::abs(infeasibility_l1(t_sub_problem_id, m_last_solutions[t_sub_problem_id]) - infeasibility_l1(t_sub_problem_id, current_solution)) > 1e-5;
     m_last_solutions[t_sub_problem_id] = std::move(current_solution);
 
     return { obj_has_changed, feas_has_changed };
@@ -411,7 +420,7 @@ void idol::Optimizers::PADM::log_inner_loop(unsigned int t_inner_loop_iteration)
 
     for (unsigned int i = 0 ; i < n_sub_problems ; ++i) {
         std::cout << std::setw(12) << m_last_solutions[i].status() << '\t';
-        std::cout << std::setw(12) << std::fixed << std::setprecision(3) << infeasibility_linf(i, m_last_solutions[i]) << '\t';
+        std::cout << std::setw(12) << std::fixed << std::setprecision(3) << infeasibility_l1(i, m_last_solutions[i]) << '\t';
     }
 
     std::cout << std::endl;
@@ -538,6 +547,52 @@ void idol::Optimizers::PADM::detect_stagnation(bool t_feasibility_has_changed) {
         return;
     }
 
-    std::cout << "Stagnation detected." << std::endl;
+    restart();
+
+}
+
+void idol::Optimizers::PADM::detect_stagnation_due_to_rescaling() {
+
+    const auto sum_sub_problems = std::accumulate(m_formulation.sub_problems().begin(), m_formulation.sub_problems().end(), 0., [](double t_acc, const Model& t_model) {
+        return t_acc + t_model.get_best_obj();
+    });
+
+    if (!m_last_objective_value_when_rescaled.has_value()) {
+        m_last_objective_value_when_rescaled = sum_sub_problems;
+        return;
+    }
+
+    if (std::abs(sum_sub_problems - m_last_objective_value_when_rescaled.value()) > 1e-4) {
+        m_last_objective_value_when_rescaled = sum_sub_problems;
+        return;
+    }
+
+    const bool has_diversified = m_penalty_update->diversify();
+
+    if(!has_diversified) {
+        restart();
+    }
+
+}
+
+void idol::Optimizers::PADM::restart() {
+
+    if (m_n_restart >= 1) {
+        set_status(Fail);
+        set_reason(IterLimit);
+        terminate();
+        return;
+    }
+
+    std::cout << "Restarting with inverse initial penalty parameter..." << std::endl;
+
+    m_last_objective_value_when_rescaled.reset();
+    m_last_iteration_with_no_feasibility_change.reset();
+    while(m_penalty_update->diversify());
+    m_current_initial_penalty_parameter = 1. / m_initial_penalty_parameter;
+    m_first_run = true;
+
+    m_formulation.initialize_penalty_parameters(m_current_initial_penalty_parameter);
 
+    ++m_n_restart;
 }