From 918b389c392f4d920c9ab3214e8161e72f9eb46a Mon Sep 17 00:00:00 2001 From: Henri Lefebvre Date: Wed, 18 Sep 2024 12:08:14 +0200 Subject: [PATCH] write KKT --- examples/bilevel-optimization/KKT.cpp | 280 ++++++++++++------ examples/bilevel-optimization/KKT.h | 52 ++-- examples/bilevel-optimization/kkt.example.cpp | 31 +- 3 files changed, 223 insertions(+), 140 deletions(-) diff --git a/examples/bilevel-optimization/KKT.cpp b/examples/bilevel-optimization/KKT.cpp index cd10bdc7..1318d3fd 100644 --- a/examples/bilevel-optimization/KKT.cpp +++ b/examples/bilevel-optimization/KKT.cpp @@ -6,204 +6,288 @@ #include "idol/modeling/expressions/operations/operators.h" idol::Reformulators::KKT::KKT(const idol::Model &t_src_model, - const std::function &t_indicator_for_vars, - const std::function &t_indicator_for_ctrs) - : m_src_model(t_src_model), - m_indicator_for_vars(t_indicator_for_vars), - m_indicator_for_ctrs(t_indicator_for_ctrs) -{ - - create_primal_variables(); - create_primal_constraints(); + const idol::Bilevel::LowerLevelDescription &t_lower_level_description) + : m_src_model(t_src_model), m_description(t_lower_level_description) { create_dual_variables_for_constraints(); - // create_dual_variables_for_upper_bounds(); - // create_dual_variables_for_lower_bounds(); - // create_dual_constraints(); + create_dual_variables_for_bounds(); + create_dual_constraints(); + create_complementarity_constraints(); + create_dual_objective(); + +} - // create_complementary_for_constraints(); - // create_complementary_for_upper_bounds(); - // create_complementary_for_lower_bounds(); +void idol::Reformulators::KKT::add_primal_variables(idol::Model &t_destination) const { + + for (const auto& var : m_src_model.vars()) { + + if (m_description.is_follower(var)) { + t_destination.add(var, TempVar(-Inf, Inf, Continuous, Column())); + continue; + } + + const auto lb = m_src_model.get_var_lb(var); + const auto ub = m_src_model.get_var_ub(var); + const auto type = m_src_model.get_var_type(var); + t_destination.add(var, TempVar(lb, ub, type, Column())); + } } -void idol::Reformulators::KKT::create_primal_variables() { +void idol::Reformulators::KKT::add_primal_constraints(Model &t_destination) const { - const unsigned int n_variables = m_src_model.vars().size(); - auto& env = m_src_model.env(); - m_primal_vars.resize(n_variables); + for (const auto& ctr : m_src_model.ctrs()) { + const auto& row = m_src_model.get_ctr_row(ctr); + const auto type = m_src_model.get_ctr_type(ctr); + t_destination.add(ctr, TempCtr(Row(row), type)); + } for (const auto& var : m_src_model.vars()) { - if (!m_indicator_for_vars(var)) { + if (m_description.is_leader(var)) { continue; } - const unsigned int index = m_src_model.get_var_index(var); + const auto lb = m_src_model.get_var_lb(var); + const auto ub = m_src_model.get_var_ub(var); - if (m_keep_primal_objects) { - m_primal_vars[index] = var; - continue; + if (!is_pos_inf(ub)) { + t_destination.add_ctr(var <= ub, "ub_" + var.name()); } - m_primal_vars[index] = Var(env, 0, Inf, Continuous, "primal_" + var.name()); + if (!is_neg_inf(lb)) { + t_destination.add_ctr(var >= lb, "lb_" + var.name()); + } } } -void idol::Reformulators::KKT::create_primal_constraints() { +void idol::Reformulators::KKT::create_dual_variables_for_constraints() { - const unsigned int n_constraints = m_src_model.ctrs().size(); auto& env = m_src_model.env(); - m_primal_ctrs.resize(n_constraints); + m_dual_variables_for_constraints.resize(m_src_model.ctrs().size()); for (const auto& ctr : m_src_model.ctrs()) { - if (!m_indicator_for_ctrs(ctr)) { + if (m_description.is_leader(ctr)) { continue; } - if (m_keep_primal_objects) { - m_primal_ctrs[m_src_model.get_ctr_index(ctr)] = ctr; - continue; + const auto type = m_src_model.get_ctr_type(ctr); + const auto index = m_src_model.get_ctr_index(ctr); + + double lb, ub; + switch (type) { + case LessOrEqual: lb = -Inf; ub = 0; break; + case Equal: lb = -Inf; ub = Inf; break; + case GreaterOrEqual: lb = 0; ub = Inf; break; } - const unsigned int index = m_src_model.get_ctr_index(ctr); - m_primal_ctrs[index] = Ctr(env, TempCtr(), "primal_" + ctr.name()); + m_dual_variables_for_constraints[index] = Var(env, lb, ub, Continuous, "dual_" + ctr.name()); } } -void idol::Reformulators::KKT::create_dual_variables_for_constraints() { +void idol::Reformulators::KKT::create_dual_variables_for_bounds() { - const unsigned int n_constraints = m_src_model.ctrs().size(); auto& env = m_src_model.env(); + m_dual_variables_for_lower_bounds.resize(m_src_model.vars().size()); + m_dual_variables_for_upper_bounds.resize(m_src_model.vars().size()); - m_dual_vars_for_constraints.resize(n_constraints); - - for (const auto& ctr : m_src_model.ctrs()) { + for (const auto& var : m_src_model.vars()) { - if (!m_indicator_for_ctrs(ctr)) { + if (m_description.is_leader(var)) { continue; } - const unsigned int index = m_src_model.get_ctr_index(ctr); - m_dual_vars_for_constraints[index] = Var(env, 0, idol::Inf, idol::Continuous, "dual_" + ctr.name()); + const auto lb = m_src_model.get_var_lb(var); + const auto ub = m_src_model.get_var_ub(var); + const auto index = m_src_model.get_var_index(var); + + if (!is_pos_inf(ub)) { + m_dual_variables_for_upper_bounds[index] = Var(env, -Inf, 0, Continuous, "dual_ub_" + var.name()); + } + + if (!is_neg_inf(lb)) { + m_dual_variables_for_lower_bounds[index] = Var(env, 0, Inf, Continuous, "dual_lb_" + var.name()); + } } } -const idol::Var &idol::Reformulators::KKT::to_destination_space(const idol::Var &t_src_var) const { - return m_primal_vars[m_src_model.get_var_index(t_src_var)].value(); -} +void idol::Reformulators::KKT::add_dual_variables(idol::Model &t_destination) const { -const idol::Var &idol::Reformulators::KKT::dual(const idol::Ctr &t_ctr) const { - return m_dual_vars_for_constraints[m_src_model.get_ctr_index(t_ctr)].value(); -} + for (const auto& opt_var : m_dual_variables_for_constraints) { + if (opt_var.has_value()) { + t_destination.add(opt_var.value()); + } + } -void idol::Reformulators::KKT::add_dual_variables(Model& t_destination) { + for (const auto& opt_var : m_dual_variables_for_lower_bounds) { + if (opt_var.has_value()) { + t_destination.add(opt_var.value()); + } + } - add_variables(t_destination, m_dual_vars_for_constraints); - add_variables(t_destination, m_dual_vars_for_upper_bounds); - add_variables(t_destination, m_dual_vars_for_lower_bounds); + for (const auto& opt_var : m_dual_variables_for_upper_bounds) { + if (opt_var.has_value()) { + t_destination.add(opt_var.value()); + } + } } -void idol::Reformulators::KKT::add_primal_variables(idol::Model &t_destination) { +void idol::Reformulators::KKT::create_dual_constraints() { + + m_dual_constraints.resize(m_src_model.vars().size()); - for (const auto& opt_var : m_primal_vars) { + for (const auto& var : m_src_model.vars()) { - if (!opt_var.has_value()) { + if (m_description.is_leader(var)) { continue; } - const auto& var = opt_var.value(); - const double lb = m_src_model.get_var_lb(var); - const double ub = m_src_model.get_var_ub(var); - const auto type = m_src_model.get_var_type(var); - - t_destination.add(var, TempVar(lb, ub, type, Column())); + create_dual_constraint(var); } } -void -idol::Reformulators::KKT::add_variables(idol::Model &t_destination, const std::vector> &t_vars) { +void idol::Reformulators::KKT::create_dual_constraint(const idol::Var &t_var) { - for (const auto& opt_var : t_vars) { - if (opt_var.has_value()) { - t_destination.add(opt_var.value()); + auto& env = m_src_model.env(); + const auto index = m_src_model.get_var_index(t_var); + const auto& col = m_src_model.get_var_column(t_var); + + Expr expr; + + for (const auto& [ctr, coeff] : col.linear()) { + + if (m_description.is_leader(ctr)) { + continue; } + + const auto index_ctr = m_src_model.get_ctr_index(ctr); + const auto& dual_var = *m_dual_variables_for_constraints[index_ctr]; + expr += coeff * dual_var; + } -} + for (const auto& [ctr, var, coeff] : col.quadratic()) { -void -idol::Reformulators::KKT::add_constraints(idol::Model &t_destination, const std::vector> &t_ctrs) { + const auto index_ctr = m_src_model.get_ctr_index(ctr); + const auto& dual_var = *m_dual_variables_for_constraints[index_ctr]; + expr += coeff * dual_var * var; - for (const auto& opt_ctr : t_ctrs) { - if (opt_ctr.has_value()) { - t_destination.add(opt_ctr.value()); - } } + if (const auto dual_var = m_dual_variables_for_lower_bounds[index]; dual_var.has_value()) { + expr += t_var * dual_var.value(); + } + + if (const auto dual_var = m_dual_variables_for_upper_bounds[index]; dual_var.has_value()) { + expr += t_var * dual_var.value(); + } + + Expr obj = m_description.follower_obj().linear().get(t_var); + + m_dual_constraints[index] = Ctr(env, expr == obj, "dual_" + t_var.name()); + } -void idol::Reformulators::KKT::add_primal_constraints(idol::Model &t_destination) { +void idol::Reformulators::KKT::create_complementarity_constraints() { - for (const auto& opt_ctr : m_primal_ctrs) { + m_complementarity_constraints.resize(m_src_model.ctrs().size()); + auto& env = m_src_model.env(); + + for (const auto& ctr : m_src_model.ctrs()) { - if (!opt_ctr.has_value()) { + if (m_description.is_leader(ctr)) { continue; } - const auto& ctr = opt_ctr.value(); const auto type = m_src_model.get_ctr_type(ctr); - const auto& src_row = m_src_model.get_ctr_row(ctr); - auto row = to_destination_space(src_row); - t_destination.add(ctr, TempCtr(std::move(row), type)); + if (type == Equal) { + continue; + } + + const auto index = m_src_model.get_ctr_index(ctr); + const auto& dual_var = *m_dual_variables_for_constraints[index]; + const auto& row = m_src_model.get_ctr_row(ctr); + + if (!row.quadratic().empty()) { + throw Exception("Cannot write complementarity constraints for quadratic constraints."); + } + + Expr expr; + + for (const auto& [var, constant] : row.linear()) { + expr += constant.as_numerical() * var * dual_var; + } + + expr -= row.rhs() * dual_var; + + m_complementarity_constraints[index] = Ctr(env, expr == 0, "complementarity_" + ctr.name()); } } -void idol::Reformulators::KKT::add_dual_constraints(idol::Model &t_destination) { +void idol::Reformulators::KKT::add_dual_constraints(idol::Model &t_destination) const { - add_constraints(t_destination, m_dual_ctrs); + for (const auto& opt_ctr : m_dual_constraints) { + if (opt_ctr.has_value()) { + t_destination.add(opt_ctr.value()); + } + } } -idol::Row idol::Reformulators::KKT::to_destination_space(const idol::Row &t_src_row) const { +void idol::Reformulators::KKT::add_complementarity_constraints(idol::Model &t_destination) const { - Row result; - result.rhs() = t_src_row.rhs(); + for (const auto& opt_ctr : m_complementarity_constraints) { + if (opt_ctr.has_value()) { + t_destination.add(opt_ctr.value()); + } + } - for (const auto& [var, coeff] : t_src_row.linear()) { +} - if (!m_indicator_for_vars(var)) { +void idol::Reformulators::KKT::add_leader_objective(idol::Model &t_destination) const { + t_destination.set_obj_expr(m_src_model.get_obj_expr()); +} - if (m_keep_original_variables) { - result.linear() += coeff * var; - } else { - result.rhs() -= coeff.as_numerical() * !var; - } +void idol::Reformulators::KKT::add_strong_duality_constraint(idol::Model &t_destination) const { + + t_destination.add_ctr(m_description.follower_obj() == m_dual_objective, "strong_duality"); + +} + +void idol::Reformulators::KKT::create_dual_objective() { - } else { + for (const auto& ctr : m_src_model.ctrs()) { + + if (m_description.is_leader(ctr)) { + continue; + } - result.linear() += coeff * to_destination_space(var); + const auto index = m_src_model.get_ctr_index(ctr); + const auto& dual_var = *m_dual_variables_for_constraints[index]; + const auto& row = m_src_model.get_ctr_row(ctr); + for (const auto& [var, constant] : row.linear()) { + if (m_description.is_follower(var)) { + continue; + } + m_dual_objective -= constant.as_numerical() * var * dual_var; } - } + m_dual_objective += row.rhs() * dual_var; - if (!t_src_row.quadratic().empty()) { - throw Exception("Not implemented."); } - return result; } diff --git a/examples/bilevel-optimization/KKT.h b/examples/bilevel-optimization/KKT.h index 91ca6e7f..6b6b3954 100644 --- a/examples/bilevel-optimization/KKT.h +++ b/examples/bilevel-optimization/KKT.h @@ -6,6 +6,7 @@ #define IDOL_KKT_H #include "idol/modeling/models/Model.h" +#include "idol/modeling/bilevel-optimization/LowerLevelDescription.h" namespace idol::Reformulators { class KKT; @@ -14,41 +15,32 @@ namespace idol::Reformulators { class idol::Reformulators::KKT { const Model& m_src_model; - std::function m_indicator_for_vars; - std::function m_indicator_for_ctrs; - const bool m_keep_original_variables = true; // If false, original variables are seen as constants in the KKT model. - const bool m_keep_primal_objects = true; // If false, new objects (variables and constraints) are created for the KKT model. + const Bilevel::LowerLevelDescription& m_description; - std::vector> m_primal_vars; - std::vector> m_primal_ctrs; - - std::vector> m_dual_vars_for_constraints; - std::vector> m_dual_vars_for_upper_bounds; - std::vector> m_dual_vars_for_lower_bounds; - std::vector> m_dual_ctrs; - - const Var& to_destination_space(const Var& t_src_var) const; // Returns the primal variable used to represent the given primal variable in the source model. - Row to_destination_space(const Row& t_src_row) const; - - const Var& dual(const Ctr& t_ctr) const; // Returns the dual variable associated with the given constraint. - - void create_primal_variables(); - void create_primal_constraints(); + std::vector> m_dual_variables_for_constraints; + std::vector> m_dual_variables_for_lower_bounds; + std::vector> m_dual_variables_for_upper_bounds; + std::vector> m_dual_constraints; + std::vector> m_complementarity_constraints; + Expr m_dual_objective; void create_dual_variables_for_constraints(); - - static void add_variables(Model& t_destination, const std::vector>& t_vars); - static void add_constraints(Model& t_destination, const std::vector>& t_ctrs); + void create_dual_variables_for_bounds(); + void create_dual_constraints(); + void create_dual_constraint(const Var& t_var); + void create_complementarity_constraints(); + void create_dual_objective(); public: KKT(const Model& t_src_model, - const std::function& t_indicator_for_vars, - const std::function& t_indicator_for_ctrs); - - void add_primal_variables(Model& t_destination); - void add_primal_constraints(Model& t_destination); - - void add_dual_variables(Model& t_destination); - void add_dual_constraints(Model& t_destination); + const Bilevel::LowerLevelDescription& t_lower_level_description); + + void add_primal_variables(idol::Model &t_destination) const; + void add_primal_constraints(Model &t_destination) const; + void add_dual_variables(Model &t_destination) const; + void add_dual_constraints(Model &t_destination) const; + void add_complementarity_constraints(Model& t_destination) const; + void add_strong_duality_constraint(Model& t_destination) const; + void add_leader_objective(Model& t_destination) const; }; diff --git a/examples/bilevel-optimization/kkt.example.cpp b/examples/bilevel-optimization/kkt.example.cpp index 7f876c0d..3c7b0ae9 100644 --- a/examples/bilevel-optimization/kkt.example.cpp +++ b/examples/bilevel-optimization/kkt.example.cpp @@ -22,7 +22,7 @@ int main(int t_argc, const char** t_argv) { y in argmin { -y : 2 x - y >= 0, - -x - y >= 0, + -x - y >= -6, -x + 6 y >= -3, x + 3 y >= 3. } @@ -41,7 +41,7 @@ int main(int t_argc, const char** t_argv) { high_point_relaxation.set_obj_expr(x + 6 * y); auto follower_c1 = high_point_relaxation.add_ctr(2 * x - y >= 0, "f1"); - auto follower_c2 = high_point_relaxation.add_ctr(-x - y >= 0, "f2"); + auto follower_c2 = high_point_relaxation.add_ctr(-x - y >= -6, "f2"); auto follower_c3 = high_point_relaxation.add_ctr(-x + 6 * y >= -3, "f3"); auto follower_c4 = high_point_relaxation.add_ctr(x + 3 * y >= 3, "f4"); high_point_relaxation.add_ctr(-x + 5 * y <= 12.5); @@ -55,18 +55,25 @@ int main(int t_argc, const char** t_argv) { description.make_follower_ctr(follower_c3); description.make_follower_ctr(follower_c4); - Reformulators::KKT kkt(high_point_relaxation, - [&](const Var& t_var) { return t_var.get(description.follower_vars()) != MasterId; }, - [&](const Ctr& t_ctr) { return t_ctr.get(description.follower_ctrs()) != MasterId; }); + Reformulators::KKT kkt(high_point_relaxation, description); - Model dual(env); - dual.add(x); - kkt.add_primal_variables(dual); - kkt.add_primal_constraints(dual); - kkt.add_dual_variables(dual); - kkt.add_dual_constraints(dual); + Model reformulation(env); + kkt.add_primal_variables(reformulation); + kkt.add_primal_constraints(reformulation); + kkt.add_dual_variables(reformulation); + kkt.add_dual_constraints(reformulation); + kkt.add_complementarity_constraints(reformulation); + kkt.add_strong_duality_constraint(reformulation); + kkt.add_leader_objective(reformulation); - std::cout << dual << std::endl; + std::cout << high_point_relaxation << std::endl; + std::cout << reformulation << std::endl; + + reformulation.use(Gurobi().with_logs(true).with_presolve(false)); + + reformulation.optimize(); + + std::cout << save_primal(reformulation) << std::endl; return 0; }