From 447a6d95ec7561907ef2762eeb2fc5a6eaa9ae1b Mon Sep 17 00:00:00 2001 From: Henri Lefebvre Date: Wed, 18 Sep 2024 18:08:37 +0200 Subject: [PATCH] first success with PADM --- .../bilevel-optimization/padm.example.cpp | 23 +- lib/CMakeLists.txt | 10 +- .../padm/AlternatingDirectionMethod.cpp | 78 +++++ .../padm/AlternatingDirectionMethod.h | 47 +++ .../padm/Formulation.cpp | 310 ++++++++++++++++++ .../padm/Formulation.h | 65 ++++ .../Optimizers_AlternatingDirectionMethod.cpp | 236 +++++++++++++ .../Optimizers_AlternatingDirectionMethod.h | 95 ++++++ .../mixed-integer-optimization/padm/PADM.cpp | 53 --- .../mixed-integer-optimization/padm/PADM.h | 57 ---- .../padm/SubProblem.cpp | 29 ++ .../padm/SubProblem.h | 32 ++ 12 files changed, 916 insertions(+), 119 deletions(-) create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.cpp create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.h create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.cpp create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.h create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.cpp create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.h delete mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.cpp delete mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.h create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.cpp create mode 100644 lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.h diff --git a/examples/bilevel-optimization/padm.example.cpp b/examples/bilevel-optimization/padm.example.cpp index 3452133f..4580e8a3 100644 --- a/examples/bilevel-optimization/padm.example.cpp +++ b/examples/bilevel-optimization/padm.example.cpp @@ -6,7 +6,8 @@ #include #include #include -#include "idol/optimizers/mixed-integer-optimization/padm/PADM.h" +#include "idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.h" +#include "idol/optimizers/mixed-integer-optimization/padm/SubProblem.h" int main(int t_argc, const char** t_argv) { @@ -40,17 +41,25 @@ int main(int t_argc, const char** t_argv) { var.set(decomposition, var.id() == delta.id()); } + Annotation penalize(env, "penalize"); + for (const auto& ctr : single_level.ctrs()) { + ctr.set(penalize, true); + } + single_level.use( - PADM(decomposition) - .with_default_sub_problem_spec( - PADM::SubProblem() - .with_optimizer(Gurobi()) - ) + AlternatingDirectionMethod(decomposition) + .with_penalization(penalize) + .with_default_sub_problem_spec( + AlternatingDirection::SubProblem() + .with_optimizer(Gurobi()) + ) ); single_level.optimize(); - //std::cout << save_primal(single_level) << std::endl; + std::cout << single_level.get_status() << std::endl; + + std::cout << save_primal(single_level) << std::endl; return 0; } \ No newline at end of file diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 351d7393..8369a5a7 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -250,8 +250,14 @@ add_library(idol STATIC src/optimizers/robust-optimization/convexification/Convexification.cpp include/idol/modeling/models/KKT.h src/modeling/models/KKT.cpp - include/idol/optimizers/mixed-integer-optimization/padm/PADM.cpp - include/idol/optimizers/mixed-integer-optimization/padm/PADM.h + include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.cpp + include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.h + include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.cpp + include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.h + include/idol/optimizers/mixed-integer-optimization/padm/Formulation.cpp + include/idol/optimizers/mixed-integer-optimization/padm/Formulation.h + include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.cpp + include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.h ) find_package(OpenMP REQUIRED) diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.cpp b/lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.cpp new file mode 100644 index 00000000..6b3f24c9 --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.cpp @@ -0,0 +1,78 @@ +// +// Created by henri on 18.09.24. +// + +#include "AlternatingDirectionMethod.h" +#include "Optimizers_AlternatingDirectionMethod.h" + +idol::AlternatingDirectionMethod::AlternatingDirectionMethod(idol::Annotation t_decomposition) + : m_decomposition(t_decomposition) { + +} + +idol::AlternatingDirectionMethod &idol::AlternatingDirectionMethod::with_default_sub_problem_spec(idol::AlternatingDirection::SubProblem t_sub_problem) { + + if (m_default_sub_problem_spec) { + throw Exception("The default sub-problem has already been set."); + } + + m_default_sub_problem_spec = std::move(t_sub_problem); + + return *this; + +} + +idol::Optimizer *idol::AlternatingDirectionMethod::operator()(const idol::Model &t_model) const { + + if (!m_decomposition) { + throw Exception("The decomposition has not been set."); + } + + AlternatingDirection::Formulation formulation(t_model, + *m_decomposition, + m_penalized_constraints); + + // create sub-problems specs + auto sub_problem_specs = create_sub_problem_specs(t_model, formulation); + + auto* result = new Optimizers::AlternatingDirectionMethod( + t_model, + std::move(formulation), + std::move(sub_problem_specs) + ); + + handle_default_parameters(result); + + return result; +} + +idol::OptimizerFactory *idol::AlternatingDirectionMethod::clone() const { + return new AlternatingDirectionMethod(*this); +} + +idol::AlternatingDirectionMethod & +idol::AlternatingDirectionMethod::with_penalization(const idol::Annotation& t_penalized_constraints) { + + if (m_penalized_constraints) { + throw Exception("The penalized constraints have already been set."); + } + + m_penalized_constraints = t_penalized_constraints; + + return *this; +} + +std::vector +idol::AlternatingDirectionMethod::create_sub_problem_specs(const idol::Model &t_model, + const idol::AlternatingDirection::Formulation &t_formulation) const { + + const unsigned int n_sub_problem = t_formulation.n_sub_problems(); + + auto result = std::vector(n_sub_problem, *m_default_sub_problem_spec); + + for (const auto& [sub_problem_id, sub_problem_spec] : m_sub_problem_specs) { + result[sub_problem_id] = AlternatingDirection::SubProblem(sub_problem_spec); + } + + return result; +} diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.h b/lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.h new file mode 100644 index 00000000..4abc7d59 --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/AlternatingDirectionMethod.h @@ -0,0 +1,47 @@ +// +// Created by henri on 18.09.24. +// + +#ifndef IDOL_ALTERNATINGDIRECTIONMETHOD_H +#define IDOL_ALTERNATINGDIRECTIONMETHOD_H + +#include "idol/optimizers/OptimizerFactory.h" +#include "idol/modeling/annotations/Annotation.h" +#include "idol/containers/Map.h" +#include "SubProblem.h" +#include "Formulation.h" +#include + +namespace idol { + class AlternatingDirectionMethod; +} + +class idol::AlternatingDirectionMethod : public OptimizerFactoryWithDefaultParameters { +public: + explicit AlternatingDirectionMethod(Annotation t_decomposition); + + AlternatingDirectionMethod(const AlternatingDirectionMethod& t_src) = default; + AlternatingDirectionMethod(AlternatingDirectionMethod&&) = default; + + AlternatingDirectionMethod& operator=(const AlternatingDirectionMethod&) = default; + AlternatingDirectionMethod& operator=(AlternatingDirectionMethod&&) = default; + + AlternatingDirectionMethod& with_default_sub_problem_spec(AlternatingDirection::SubProblem t_sub_problem); + + AlternatingDirectionMethod& with_penalization(const Annotation& t_penalized_constraints); + + Optimizer *operator()(const Model &t_model) const override; + + OptimizerFactory *clone() const override; + +private: + std::optional> m_decomposition; + std::optional> m_penalized_constraints; + std::optional m_default_sub_problem_spec; + Map m_sub_problem_specs; + + std::vector create_sub_problem_specs(const Model& t_model, const AlternatingDirection::Formulation& t_formulation) const; +}; + + +#endif //IDOL_ALTERNATINGDIRECTIONMETHOD_H diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.cpp b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.cpp new file mode 100644 index 00000000..355f8727 --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.cpp @@ -0,0 +1,310 @@ +// +// Created by henri on 18.09.24. +// + +#include "Formulation.h" +#include "idol/modeling/objects/Versions.h" +#include "idol/modeling/expressions/operations/operators.h" + +#include + +idol::AlternatingDirection::Formulation::Formulation(const Model& t_src_model, + Annotation t_decomposition, + std::optional> t_penalized_constraints) + : m_decomposition(std::move(t_decomposition)), + m_penalized_constraints(std::move(t_penalized_constraints)) { + + const auto n_sub_problems = compute_n_sub_problems(t_src_model); + + initialize_sub_problems(t_src_model, n_sub_problems); + initialize_patterns(t_src_model, n_sub_problems); + initialize_slacks(t_src_model, n_sub_problems); + + dispatch_vars(t_src_model); + dispatch_ctrs(t_src_model); + dispatch_obj(t_src_model); + +} + +unsigned int idol::AlternatingDirection::Formulation::compute_n_sub_problems(const idol::Model &t_src_model) const { + unsigned int result = 0; + for (const auto& var : t_src_model.vars()) { + result = std::max(result, var.get(m_decomposition)); + } + return result + 1; +} + +void idol::AlternatingDirection::Formulation::initialize_sub_problems(const idol::Model &t_src_model, + unsigned int n_sub_problems) { + + auto& env = t_src_model.env(); + m_sub_problems.reserve(n_sub_problems); + + for (unsigned int i = 0 ; i < n_sub_problems ; ++i) { + m_sub_problems.emplace_back(env); + } + +} + +void idol::AlternatingDirection::Formulation::initialize_patterns(const idol::Model &t_src_model, + unsigned int n_sub_problems) { + + m_objective_patterns.resize(n_sub_problems); + m_constraint_patterns.resize(n_sub_problems); + +} + +void idol::AlternatingDirection::Formulation::initialize_slacks(const idol::Model &t_src_model, + unsigned int n_sub_problems) { + m_l1_vars.resize(n_sub_problems); +} + +void idol::AlternatingDirection::Formulation::dispatch_vars(const idol::Model &t_src_model) { + + for (const auto& var : t_src_model.vars()) { + + const unsigned int sub_problem_id = var.get(m_decomposition); + const auto lb = t_src_model.get_var_lb(var); + const auto ub = t_src_model.get_var_ub(var); + const auto type = t_src_model.get_var_type(var); + + m_sub_problems[sub_problem_id].add(var, TempVar(lb, ub, type, Column())); + + } + +} + +void idol::AlternatingDirection::Formulation::dispatch_ctrs(const idol::Model &t_src_model) { + + const auto n_sub_problems = m_sub_problems.size(); + + for (const auto& ctr : t_src_model.ctrs()) { + for (unsigned int i = 0 ; i < n_sub_problems ; ++i) { + dispatch_ctr(t_src_model, ctr, i); + } + } + +} + +void idol::AlternatingDirection::Formulation::dispatch_ctr(const idol::Model &t_src_model, const idol::Ctr &t_ctr, unsigned int t_sub_problem_id) { + + const auto& row = t_src_model.get_ctr_row(t_ctr); + const auto type = t_src_model.get_ctr_type(t_ctr); + + auto [pattern, is_pure] = dispatch(t_src_model, row.linear(), row.quadratic(), t_sub_problem_id); + pattern.constant() -= row.rhs(); + + if (pattern.linear().empty() && pattern.quadratic().empty()) { + return; + } + + if (m_penalized_constraints && t_ctr.get(*m_penalized_constraints)) { + 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()); + pattern.linear() += t_coefficient * var; + m_l1_vars[t_sub_problem_id].emplace_back(var); + }; + + switch (type) { + case Equal: + add_l1_var(+1); + add_l1_var(-1); + break; + case LessOrEqual: + add_l1_var(-1); + break; + case GreaterOrEqual: + add_l1_var(+1); + break; + } + } + + if (is_pure) { + m_sub_problems[t_sub_problem_id].add(t_ctr, TempCtr(Row(std::move(pattern), 0), type)); + } else { + m_sub_problems[t_sub_problem_id].add(t_ctr, TempCtr(Row(), type)); + m_constraint_patterns[t_sub_problem_id].emplace_back(t_ctr, std::move(pattern)); + } + +} + +void +idol::AlternatingDirection::Formulation::dispatch_obj(const Model &t_src_model) { + + const unsigned int n_sub_problems = m_sub_problems.size(); + + for (unsigned int i = 0 ; i < n_sub_problems ; ++i) { + dispatch_obj(t_src_model, i); + } + +} + +std::pair, bool> idol::AlternatingDirection::Formulation::dispatch(const idol::Model &t_src_model, + const idol::LinExpr &t_lin_expr, + const idol::QuadExpr &t_quad_expr, + unsigned int t_sub_problem_id) { + + bool is_pure = true; // true if the row only has variables from the same sub-problem + + Expr pattern; + + for (const auto& [var, coefficient] : t_lin_expr) { + + const unsigned int var_sub_problem_id = var.get(m_decomposition); + + if (var_sub_problem_id != t_sub_problem_id) { + is_pure = false; + pattern.constant() += coefficient.as_numerical() * !var; + continue; + } + + pattern.linear() += coefficient * var; + } + + for (const auto& [var1, var2, constant] : t_quad_expr) { + + const unsigned int var1_sub_problem_id = var1.get(m_decomposition); + const unsigned int var2_sub_problem_id = var2.get(m_decomposition); + + if (var1_sub_problem_id != t_sub_problem_id && var2_sub_problem_id != t_sub_problem_id) { + is_pure = false; + pattern.constant() += constant.as_numerical() * (!var1 * !var2); + continue; + } + + if (var1_sub_problem_id != t_sub_problem_id) { + is_pure = false; + pattern.linear() += constant.as_numerical() * !var1 * var2; + continue; + } + + if (var2_sub_problem_id != t_sub_problem_id) { + is_pure = false; + pattern.linear() += constant.as_numerical() * !var2 * var1; + continue; + } + + pattern.quadratic() += constant * var1 * var2; + + } + + return { + std::move(pattern), + is_pure + }; +} + +void +idol::AlternatingDirection::Formulation::dispatch_obj(const Model &t_src_model, unsigned int t_sub_problem_id) { + + const auto& obj = t_src_model.get_obj_expr(); + auto [pattern, is_pure] = dispatch(t_src_model, obj.linear(), obj.quadratic(), t_sub_problem_id); + pattern += obj.constant(); + + if (pattern.linear().empty() && pattern.quadratic().empty()) { + return; + } + + if (is_pure) { + m_sub_problems[t_sub_problem_id].set_obj_expr(pattern); + return; + } + + m_objective_patterns[t_sub_problem_id] = std::move(pattern); + +} + +void idol::AlternatingDirection::Formulation::fix_sub_problem(unsigned int t_sub_problem_id, + const std::vector &t_primals) { + + // Constraints + for (const auto& [ctr, pattern] : m_constraint_patterns[t_sub_problem_id]) { + + Expr lhs = fix(pattern.constant(), t_primals); + + for (const auto& [var, coefficient] : pattern.linear()) { + lhs += fix(coefficient, t_primals) * var; + } + + for (const auto& [var1, var2, coefficient] : pattern.quadratic()) { + lhs += fix(coefficient, t_primals) * var1 * var2; + } + + m_sub_problems[t_sub_problem_id].set_ctr_row(ctr, Row(std::move(lhs), 0)); + } + + // Objective + if (m_objective_patterns[t_sub_problem_id]) { + const auto& obj_pattern = *m_objective_patterns[t_sub_problem_id]; + Expr obj = fix(obj_pattern.constant(), t_primals); + + for (const auto& [var, coefficient] : obj_pattern.linear()) { + obj += fix(coefficient, t_primals) * var; + } + + for (const auto& [var1, var2, coefficient] : obj_pattern.quadratic()) { + obj += fix(coefficient, t_primals) * var1 * var2; + } + + m_sub_problems[t_sub_problem_id].set_obj_expr(std::move(obj)); + + } + +} + +double idol::AlternatingDirection::Formulation::fix(const idol::Constant &t_constant, + const std::vector &t_primals) { + double result = t_constant.numerical(); + + for (const auto& [param, coefficient] : t_constant.linear()) { + const auto& var = param.as(); + const auto var_sub_problem_id = var.get(m_decomposition); + const auto& solution = t_primals[var_sub_problem_id]; + result += coefficient * solution.get(var); + } + + for (const auto& [params, coefficient] : t_constant.quadratic()) { + const auto& var1 = params.first.as(); + const auto& var2 = params.second.as(); + const auto var1_sub_problem_id = var1.get(m_decomposition); + const auto var2_sub_problem_id = var2.get(m_decomposition); + const auto& solution1 = t_primals[var1_sub_problem_id]; + const auto& solution2 = t_primals[var2_sub_problem_id]; + result += coefficient * solution1.get(var1) * solution2.get(var2); + } + + return result; +} + +idol::Model &idol::AlternatingDirection::Formulation::sub_problem(const idol::Var &t_var) { + return m_sub_problems[t_var.get(m_decomposition)]; +} + +const idol::Model &idol::AlternatingDirection::Formulation::sub_problem(const idol::Var &t_var) const { + return m_sub_problems[t_var.get(m_decomposition)]; +} + +void +idol::AlternatingDirection::Formulation::update_penalty_parameters(const std::vector &t_primals) { + + 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]) { + + double current_penalty = model.get_var_column(var).obj().as_numerical(); + + if (current_penalty < 1e-3) { + model.set_var_obj(var, 1); + continue; + } + + if (t_primals[i].get(var) > 1e-4) { + model.set_var_obj(var, current_penalty + 1); + } + + } + } + +} diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.h b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.h new file mode 100644 index 00000000..d42f38a3 --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Formulation.h @@ -0,0 +1,65 @@ +// +// Created by henri on 18.09.24. +// + +#ifndef IDOL_ADM_FORMULATION_H +#define IDOL_ADM_FORMULATION_H + +#include "idol/modeling/models/Model.h" + +namespace idol::AlternatingDirection { + class Formulation; +} + +class idol::AlternatingDirection::Formulation { +public: + Formulation(const Model& t_src_model, + Annotation t_decomposition, + std::optional> t_penalized_constraints); + + Model& sub_problem(const Var& t_var); + + const Model& sub_problem(const Var& t_var) const; + + Model& sub_problem(unsigned int t_sub_problem_id) { return m_sub_problems[t_sub_problem_id]; } + + const Model& sub_problem(unsigned int t_sub_problem_id) const { return m_sub_problems[t_sub_problem_id]; } + + unsigned int n_sub_problems() const { return m_sub_problems.size(); } + + auto sub_problems() { return IteratorForward(m_sub_problems); } + + 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]); } + + bool has_penalized_constraints() const { return m_penalized_constraints.has_value(); } + + void fix_sub_problem(unsigned int t_sub_problem_id, const std::vector& t_primals); + + void update_penalty_parameters(const std::vector& t_primals); +private: + Annotation m_decomposition; + std::optional> m_penalized_constraints; + + std::vector m_sub_problems; + std::vector>> m_objective_patterns; + std::vector>>> m_constraint_patterns; // as ctr: lhs <= 0 + std::vector> 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); + void initialize_patterns(const Model& t_src_model, unsigned int n_sub_problems); + void initialize_slacks(const Model& t_src_model, unsigned int n_sub_problems); + void dispatch_vars(const Model& t_src_model); + void dispatch_ctrs(const Model& t_src_model); + void dispatch_ctr(const Model& t_src_model, const Ctr& t_ctr, unsigned int t_sub_problem_id); + void dispatch_obj(const Model& t_src_model); + void dispatch_obj(const Model& t_src_model, unsigned int t_sub_problem_id); + std::pair, bool> dispatch(const Model& t_src_model, const LinExpr& t_lin_expr, const QuadExpr& t_quad_expr, unsigned int t_sub_problem_id); + + double fix(const Constant& t_constant, const std::vector& t_primals); +}; + + +#endif //IDOL_ADM_FORMULATION_H diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.cpp b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.cpp new file mode 100644 index 00000000..67f4cc11 --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.cpp @@ -0,0 +1,236 @@ +// +// Created by henri on 18.09.24. +// + +#include "Optimizers_AlternatingDirectionMethod.h" + +idol::Optimizers::AlternatingDirectionMethod::AlternatingDirectionMethod(const Model& t_model, + AlternatingDirection::Formulation t_formulation, + std::vector&& t_sub_problem_specs) + : Algorithm(t_model), + m_formulation(std::move(t_formulation)), + m_sub_problem_specs(std::move(t_sub_problem_specs)) { + +} + +double idol::Optimizers::AlternatingDirectionMethod::get_var_primal(const idol::Var &t_var) const { + return m_formulation.sub_problem(t_var).get_var_primal(t_var); +} + +double idol::Optimizers::AlternatingDirectionMethod::get_var_reduced_cost(const idol::Var &t_var) const { + return m_formulation.sub_problem(t_var).get_var_reduced_cost(t_var); +} + +double idol::Optimizers::AlternatingDirectionMethod::get_var_ray(const idol::Var &t_var) const { + return m_formulation.sub_problem(t_var).get_var_ray(t_var); +} + +double idol::Optimizers::AlternatingDirectionMethod::get_ctr_dual(const idol::Ctr &t_ctr) const { + throw Exception("Not implemented"); +} + +double idol::Optimizers::AlternatingDirectionMethod::get_ctr_farkas(const idol::Ctr &t_ctr) const { + throw Exception("Not implemented"); +} + +unsigned int idol::Optimizers::AlternatingDirectionMethod::get_n_solutions() const { + return Algorithm::get_status() == Feasible ? 1 : 0; +} + +unsigned int idol::Optimizers::AlternatingDirectionMethod::get_solution_index() const { + return 0; +} + +void idol::Optimizers::AlternatingDirectionMethod::add(const idol::Var &t_var) { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::add(const idol::Ctr &t_ctr) { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::remove(const idol::Var &t_var) { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::remove(const idol::Ctr &t_ctr) { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update() { + +} + +void idol::Optimizers::AlternatingDirectionMethod::write(const std::string &t_name) { + throw Exception("Not implemented"); +} + + +void idol::Optimizers::AlternatingDirectionMethod::hook_before_optimize() { + Optimizer::hook_before_optimize(); + + set_status(Loaded); + set_reason(NotSpecified); + set_best_obj(Inf); + set_best_bound(-Inf); + + const unsigned int n_sub_problems = m_formulation.n_sub_problems(); + + m_outer_loop_iteration = 0; + m_inner_loop_iterations = 0; + m_last_solutions = std::vector(n_sub_problems); + + for (unsigned int i = 0 ; i < n_sub_problems ; ++i) { + auto& model = m_formulation.sub_problem(i); + if (!model.has_optimizer()) { + model.use(m_sub_problem_specs[i].optimizer_factory()); + } + } +} + + +void idol::Optimizers::AlternatingDirectionMethod::hook_optimize() { + + for (unsigned int max_outer_loop = get_param_iteration_limit() ; m_outer_loop_iteration < max_outer_loop ; ++m_outer_loop_iteration) { + + update_penalty_parameters(); + + run_inner_loop(); + + if (is_feasible()) { + set_status(Feasible); + set_reason(Proved); + std::cout << "DONE mega DONE!" << std::endl; + break; + } + + } + +} + +void idol::Optimizers::AlternatingDirectionMethod::set_solution_index(unsigned int t_index) { + if (t_index != 0) { + throw Exception("Invalid solution index"); + } +} + +void idol::Optimizers::AlternatingDirectionMethod::update_obj_sense() { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_obj() { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_rhs() { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_obj_constant() { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_mat_coeff(const idol::Ctr &t_ctr, const idol::Var &t_var) { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_ctr_type(const idol::Ctr &t_ctr) { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_ctr_rhs(const idol::Ctr &t_ctr) { + throw Exception("Not implemented"); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_var_type(const idol::Var &t_var) { + const auto type = parent().get_var_type(t_var); + m_formulation.sub_problem(t_var).set_var_type(t_var, type); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_var_lb(const idol::Var &t_var) { + const auto lb = parent().get_var_lb(t_var); + m_formulation.sub_problem(t_var).set_var_lb(t_var, lb); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_var_ub(const idol::Var &t_var) { + const auto ub = parent().get_var_ub(t_var); + m_formulation.sub_problem(t_var).set_var_ub(t_var, ub); +} + +void idol::Optimizers::AlternatingDirectionMethod::update_var_obj(const idol::Var &t_var) { + throw Exception("Not implemented"); +} + +bool idol::Optimizers::AlternatingDirectionMethod::is_feasible() const { + for (unsigned int i = 0, n = m_formulation.n_sub_problems() ; i < n ; ++i) { + if (!is_feasible(i)) { + return false; + } + } + return true; +} + +bool idol::Optimizers::AlternatingDirectionMethod::is_feasible(unsigned int t_sub_problem_id) const { + + std::cout << m_last_solutions[t_sub_problem_id] << std::endl; + + for (const auto& var : m_formulation.l1_vars(t_sub_problem_id)) { + if (m_last_solutions[t_sub_problem_id].get(var) > 1e-4) { + return false; + } + } + + return true; +} + +void idol::Optimizers::AlternatingDirectionMethod::run_inner_loop() { + + for (unsigned int inner_loop_iteration = 0 ; inner_loop_iteration < m_max_inner_loop_iterations ; ++inner_loop_iteration) { + + bool has_changed = false; + + for (unsigned int i = 0, n = m_formulation.n_sub_problems() ; i < n ; ++i) { + has_changed |= solve_sub_problem(i); + } + + if (!has_changed || is_terminated()) { + break; + } + + ++m_inner_loop_iterations; + } + +} + +void idol::Optimizers::AlternatingDirectionMethod::update_penalty_parameters() { + + m_formulation.update_penalty_parameters(m_last_solutions); + +} + +bool idol::Optimizers::AlternatingDirectionMethod::solve_sub_problem(unsigned int t_sub_problem_id) { + + std::cout << "Solving sub-problem " << t_sub_problem_id << std::endl; + + m_formulation.fix_sub_problem(t_sub_problem_id, m_last_solutions); + + auto& model = m_formulation.sub_problem(t_sub_problem_id); + + model.optimize(); + + const auto status = model.get_status(); + + std::cout << "Sub-problem status: " << status << std::endl; + + if (status != Optimal && status != Feasible) { + set_status(status); + set_reason(NotSpecified); + terminate(); + return true; + } + + bool has_changed = m_inner_loop_iterations == 0 || std::abs(m_last_solutions[t_sub_problem_id].objective_value() - model.get_best_obj()) > 1e-4; + m_last_solutions[t_sub_problem_id] = save_primal(model); + + return has_changed; +} diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.h b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.h new file mode 100644 index 00000000..05491e09 --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/Optimizers_AlternatingDirectionMethod.h @@ -0,0 +1,95 @@ +// +// Created by henri on 18.09.24. +// + +#ifndef IDOL_OPTIMIZERS_ALTERNATINGDIRECTIONMETHOD_H +#define IDOL_OPTIMIZERS_ALTERNATINGDIRECTIONMETHOD_H + +#include "idol/optimizers/Algorithm.h" +#include "idol/optimizers/mixed-integer-optimization/padm/SubProblem.h" +#include "Formulation.h" + +namespace idol::Optimizers { + class AlternatingDirectionMethod; +} + +class idol::Optimizers::AlternatingDirectionMethod : public Algorithm { +public: + AlternatingDirectionMethod(const Model& t_model, + AlternatingDirection::Formulation t_formulation, + std::vector&& t_sub_problem_specs); + + std::string name() const override { return "AlternatingDirectionMethod"; } + + double get_var_primal(const Var &t_var) const override; + + double get_var_reduced_cost(const Var &t_var) const override; + + double get_var_ray(const Var &t_var) const override; + + double get_ctr_dual(const Ctr &t_ctr) const override; + + double get_ctr_farkas(const Ctr &t_ctr) const override; + + unsigned int get_n_solutions() const override; + + unsigned int get_solution_index() const override; + +protected: + void add(const Var &t_var) override; + + void add(const Ctr &t_ctr) override; + + void remove(const Var &t_var) override; + + void remove(const Ctr &t_ctr) override; + + void update() override; + + void write(const std::string &t_name) override; + + void hook_before_optimize() override; + + void hook_optimize() override; + + void set_solution_index(unsigned int t_index) override; + + void update_obj_sense() override; + + void update_obj() override; + + void update_rhs() override; + + void update_obj_constant() override; + + void update_mat_coeff(const Ctr &t_ctr, const Var &t_var) override; + + void update_ctr_type(const Ctr &t_ctr) override; + + void update_ctr_rhs(const Ctr &t_ctr) override; + + void update_var_type(const Var &t_var) override; + + void update_var_lb(const Var &t_var) override; + + void update_var_ub(const Var &t_var) override; + + void update_var_obj(const Var &t_var) override; + + void update_penalty_parameters(); + void run_inner_loop(); + bool is_feasible() const; + bool is_feasible(unsigned int t_sub_problem_id) const; + bool solve_sub_problem(unsigned int t_sub_problem_id); +private: + AlternatingDirection::Formulation m_formulation; + std::vector m_sub_problem_specs; + unsigned int m_max_inner_loop_iterations = 1000; + + unsigned int m_outer_loop_iteration = 0; + unsigned int m_inner_loop_iterations = 0; + std::vector m_last_solutions; +}; + + +#endif //IDOL_OPTIMIZERS_ALTERNATINGDIRECTIONMETHOD_H diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.cpp b/lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.cpp deleted file mode 100644 index cb95cb70..00000000 --- a/lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.cpp +++ /dev/null @@ -1,53 +0,0 @@ -// -// Created by henri on 18.09.24. -// - -#include "PADM.h" - -idol::PADM::PADM(idol::Annotation t_decomposition) - : m_decomposition(t_decomposition) { - -} - -idol::PADM::PADM(idol::Annotation t_decomposition, - idol::Annotation t_penalized_constraints) - : m_decomposition(t_decomposition), - m_penalized_constraints(t_penalized_constraints) { - -} - -idol::PADM::SubProblem &idol::PADM::SubProblem::with_optimizer(const idol::OptimizerFactory &t_optimizer_factory) { - - if (m_optimizer_factory) { - throw Exception("The optimizer has already been set."); - } - - m_optimizer_factory.reset(t_optimizer_factory.clone()); - - return *this; -} - -idol::PADM::SubProblem::SubProblem(const idol::PADM::SubProblem & t_src) - : m_optimizer_factory(t_src.m_optimizer_factory ? t_src.m_optimizer_factory->clone() : nullptr) { - -} - -idol::PADM &idol::PADM::with_default_sub_problem_spec(idol::PADM::SubProblem t_sub_problem) { - - if (m_default_sub_problem_spec) { - throw Exception("The default sub-problem has already been set."); - } - - m_default_sub_problem_spec = std::move(t_sub_problem); - - return *this; - -} - -idol::Optimizer *idol::PADM::operator()(const idol::Model &t_model) const { - throw Exception("Not implemented."); -} - -idol::OptimizerFactory *idol::PADM::clone() const { - return new PADM(*this); -} diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.h b/lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.h deleted file mode 100644 index 6ba4c438..00000000 --- a/lib/include/idol/optimizers/mixed-integer-optimization/padm/PADM.h +++ /dev/null @@ -1,57 +0,0 @@ -// -// Created by henri on 18.09.24. -// - -#ifndef IDOL_PADM_H -#define IDOL_PADM_H - -#include "idol/optimizers/OptimizerFactory.h" -#include "idol/modeling/annotations/Annotation.h" -#include "idol/containers/Map.h" -#include - -namespace idol { - class PADM; -} - -class idol::PADM : public OptimizerFactoryWithDefaultParameters { -public: - explicit PADM(Annotation t_decomposition); - - PADM(Annotation t_decomposition, Annotation t_penalized_constraints); - - PADM(const PADM& t_src) = default; - PADM(PADM&&) = default; - - PADM& operator=(const PADM&) = default; - PADM& operator=(PADM&&) = default; - - class SubProblem { - std::unique_ptr m_optimizer_factory; - public: - SubProblem() = default; - - SubProblem(const SubProblem&); - SubProblem(SubProblem&&) = default; - - SubProblem& operator=(const SubProblem&) = delete; - SubProblem& operator=(SubProblem&&) = default; - - SubProblem& with_optimizer(const OptimizerFactory& t_optimizer_factory); - }; - - PADM& with_default_sub_problem_spec(SubProblem t_sub_problem); - - Optimizer *operator()(const Model &t_model) const override; - - OptimizerFactory *clone() const override; - -private: - std::optional> m_decomposition; - std::optional> m_penalized_constraints; - std::optional m_default_sub_problem_spec; - Map m_sub_problem_specs; -}; - - -#endif //IDOL_PADM_H diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.cpp b/lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.cpp new file mode 100644 index 00000000..166ba51d --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.cpp @@ -0,0 +1,29 @@ +// +// Created by henri on 18.09.24. +// + +#include "SubProblem.h" + +idol::AlternatingDirection::SubProblem &idol::AlternatingDirection::SubProblem::with_optimizer(const idol::OptimizerFactory &t_optimizer_factory) { + + if (m_optimizer_factory) { + throw Exception("The optimizer has already been set."); + } + + m_optimizer_factory.reset(t_optimizer_factory.clone()); + + return *this; +} + +idol::AlternatingDirection::SubProblem::SubProblem(const idol::AlternatingDirection::SubProblem & t_src) + : m_optimizer_factory(t_src.m_optimizer_factory ? t_src.m_optimizer_factory->clone() : nullptr) { + +} + +const idol::OptimizerFactory &idol::AlternatingDirection::SubProblem::optimizer_factory() const { + if (!m_optimizer_factory) { + throw Exception("The optimizer has not been set."); + } + + return *m_optimizer_factory; +} diff --git a/lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.h b/lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.h new file mode 100644 index 00000000..4edad825 --- /dev/null +++ b/lib/include/idol/optimizers/mixed-integer-optimization/padm/SubProblem.h @@ -0,0 +1,32 @@ +// +// Created by henri on 18.09.24. +// + +#ifndef IDOL_ADM_SUBPROBLEM_H +#define IDOL_ADM_SUBPROBLEM_H + + +#include "idol/optimizers/OptimizerFactory.h" + +namespace idol::AlternatingDirection { + class SubProblem; +} + +class idol::AlternatingDirection::SubProblem { + std::unique_ptr m_optimizer_factory; +public: + SubProblem() = default; + + SubProblem(const SubProblem&); + SubProblem(SubProblem&&) = default; + + SubProblem& operator=(const SubProblem&) = delete; + SubProblem& operator=(SubProblem&&) = default; + + SubProblem& with_optimizer(const OptimizerFactory& t_optimizer_factory); + + const OptimizerFactory& optimizer_factory() const; +}; + + +#endif //IDOL_ADM_SUBPROBLEM_H