From 71492ae6e907c81dfa155919d388dba7a995a43d Mon Sep 17 00:00:00 2001 From: Hassan Date: Thu, 15 Feb 2018 17:58:21 -0700 Subject: [PATCH] implementing lazy constraint generation, SDPOPF on 2383wp_k takes 407 secs on my mac. --- .../NeuralNet/NeuralNet_main.cpp | 4 +- examples/MINLP/Power/ACOPF/ACOPF_main.cpp | 8 +- examples/MINLP/Power/PowerNet.cpp | 3 + examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 15 +- examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp | 114 ++-- include/gravity/constraint.h | 14 + include/gravity/func.h | 2 - include/gravity/model.h | 9 +- include/gravity/solver.h | 2 +- src/GurobiProgram.cpp | 120 ++-- src/IpoptProgram.cpp | 6 +- src/constraint.cpp | 19 + src/func.cpp | 24 +- src/model.cpp | 594 ++++++++++++------ src/solver.cpp | 16 +- 15 files changed, 596 insertions(+), 354 deletions(-) diff --git a/examples/Classification/NeuralNet/NeuralNet_main.cpp b/examples/Classification/NeuralNet/NeuralNet_main.cpp index 3fc5d2a5d..51a0ee0ee 100644 --- a/examples/Classification/NeuralNet/NeuralNet_main.cpp +++ b/examples/Classification/NeuralNet/NeuralNet_main.cpp @@ -78,11 +78,11 @@ int main (int argc, const char * argv[]) SolverType stype = cplex; double wall0 = get_wall_time(); double cpu0 = get_cpu_time(); - cout << "Running the NONCONVEX model\n"; + //s.run(); double wall1 = get_wall_time(); double cpu1 = get_cpu_time(); - cout << "Done running the NONCONVEX model\n"; + cout << "Wall clock computing time = " << wall1 - wall0 << "\n"; cout << "CPU computing time = " << cpu1 - cpu0 << "\n"; return 0; diff --git a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp index 4ff2c32f9..f6bc0116d 100644 --- a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp +++ b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp @@ -242,20 +242,20 @@ int main (int argc, char * argv[]) DebugOff(grid.th_min.to_str(true) << endl); DebugOff(grid.th_max.to_str(true) << endl); } - ACOPF.add_constraint(PAD_UB.in(bus_pairs) <= 0); - ACOPF.add_constraint(PAD_LB.in(bus_pairs) >= 0); +// ACOPF.add_constraint(PAD_UB.in(bus_pairs) <= 0); +// ACOPF.add_constraint(PAD_LB.in(bus_pairs) >= 0); /* Thermal Limit Constraints */ Constraint Thermal_Limit_from("Thermal_Limit_from"); Thermal_Limit_from += power(Pf_from, 2) + power(Qf_from, 2); Thermal_Limit_from -= power(grid.S_max, 2); - ACOPF.add_constraint(Thermal_Limit_from.in(grid.arcs) <= 0); +// ACOPF.add_constraint(Thermal_Limit_from.in(grid.arcs) <= 0); Constraint Thermal_Limit_to("Thermal_Limit_to"); Thermal_Limit_to += power(Pf_to, 2) + power(Qf_to, 2); Thermal_Limit_to -= power(grid.S_max,2); - ACOPF.add_constraint(Thermal_Limit_to.in(grid.arcs) <= 0); +// ACOPF.add_constraint(Thermal_Limit_to.in(grid.arcs) <= 0); DebugOff(grid.S_max.in(grid.arcs).to_str(true) << endl); solver OPF(ACOPF,ipopt); diff --git a/examples/MINLP/Power/PowerNet.cpp b/examples/MINLP/Power/PowerNet.cpp index 7afe99431..d3fcbe199 100755 --- a/examples/MINLP/Power/PowerNet.cpp +++ b/examples/MINLP/Power/PowerNet.cpp @@ -554,6 +554,9 @@ int PowerNet::readgrid(const char* fname) { DebugOff(tr.to_str(true) << endl); file.close(); + if (nodes.size()>2000) { + add_3d_nlin = false; + } return 0; } diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index 47b6c2144..55ac8acf1 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -189,6 +189,7 @@ int main (int argc, char * argv[]) { /** Constraints */ if(grid.add_3d_nlin) { + DebugOn("Adding 3d determinant polynomial cuts\n"); auto R_Wij_ = R_Wij.pairs_in_directed(grid, grid._bags, 3); auto Im_Wij_ = Im_Wij.pairs_in_directed(grid, grid._bags, 3); auto Wii_ = Wii.in(grid._bags, 3); @@ -216,7 +217,7 @@ int main (int argc, char * argv[]) { SDP3 -= (power(R_Wij_[2], 2) + power(Im_Wij_[2], 2)) * Wii_[1]; SDP3 += Wii_[0] * Wii_[1] * Wii_[2]; DebugOff("\nsdp nb inst = " << SDP3.get_nb_instances() << endl); - SDP.add_constraint(SDP3 >= 0); + SDP.add(SDP3 >= 0); // SDP3.print_expanded(); } @@ -257,24 +258,24 @@ int main (int argc, char * argv[]) { Constraint PAD_UB("PAD_UB"); PAD_UB = Im_Wij; PAD_UB <= grid.tan_th_max*R_Wij; - SDP.add_constraint(PAD_UB.in(bus_pairs)); + SDP.add_lazy(PAD_UB.in(bus_pairs)); Constraint PAD_LB("PAD_LB"); PAD_LB = Im_Wij; PAD_LB >= grid.tan_th_min*R_Wij; - SDP.add_constraint(PAD_LB.in(bus_pairs)); + SDP.add_lazy(PAD_LB.in(bus_pairs)); /* Thermal Limit Constraints */ Constraint Thermal_Limit_from("Thermal_Limit_from"); Thermal_Limit_from = power(Pf_from, 2) + power(Qf_from, 2); Thermal_Limit_from <= power(grid.S_max,2); - SDP.add_constraint(Thermal_Limit_from.in(grid.arcs)); + SDP.add_lazy(Thermal_Limit_from.in(grid.arcs)); Constraint Thermal_Limit_to("Thermal_Limit_to"); Thermal_Limit_to = power(Pf_to, 2) + power(Qf_to, 2); Thermal_Limit_to <= power(grid.S_max,2); - SDP.add_constraint(Thermal_Limit_to.in(grid.arcs)); + SDP.add_lazy(Thermal_Limit_to.in(grid.arcs)); /* Lifted Nonlinear Cuts */ Constraint LNC1("LNC1"); @@ -282,14 +283,14 @@ int main (int argc, char * argv[]) { LNC1 -= grid.v_max.to()*cos(0.5*(grid.th_max-grid.th_min))*(grid.v_min.to()+grid.v_max.to())*Wii.from(); LNC1 -= grid.v_max.from()*cos(0.5*(grid.th_max-grid.th_min))*(grid.v_min.from()+grid.v_max.from())*Wii.to(); LNC1 -= grid.v_max.from()*grid.v_max.to()*cos(0.5*(grid.th_max-grid.th_min))*(grid.v_min.from()*grid.v_min.to() - grid.v_max.from()*grid.v_max.to()); - SDP.add_constraint(LNC1.in(bus_pairs) >= 0); + SDP.add_lazy(LNC1.in(bus_pairs) >= 0); Constraint LNC2("LNC2"); LNC2 = (grid.v_min.from()+grid.v_max.from())*(grid.v_min.to()+grid.v_max.to())*(sin(0.5*(grid.th_max+grid.th_min))*Im_Wij + cos(0.5*(grid.th_max+grid.th_min))*R_Wij); LNC2 -= grid.v_min.to()*cos(0.5*(grid.th_max-grid.th_min))*(grid.v_min.to()+grid.v_max.to())*Wii.from(); LNC2 -= grid.v_min.from()*cos(0.5*(grid.th_max-grid.th_min))*(grid.v_min.from()+grid.v_max.from())*Wii.to(); LNC2 += grid.v_min.from()*grid.v_min.to()*cos(0.5*(grid.th_max-grid.th_min))*(grid.v_min.from()*grid.v_min.to() - grid.v_max.from()*grid.v_max.to()); - SDP.add_constraint(LNC2.in(bus_pairs) >= 0); + SDP.add_lazy(LNC2.in(bus_pairs) >= 0); vector bags; int n3 = 0; diff --git a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp index 3b56240e9..b4979c2a7 100644 --- a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp +++ b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp @@ -62,15 +62,15 @@ int main (int argc, char * argv[]) use_cplex = true; } double total_time_start = get_wall_time(); - PowerNet* grid = new PowerNet(); - grid->readgrid(fname.c_str()); + PowerNet grid; + grid.readgrid(fname.c_str()); /* Grid Parameters */ - auto bus_pairs = grid->get_bus_pairs(); - auto nb_bus_pairs = grid->get_nb_active_bus_pairs(); - auto nb_gen = grid->get_nb_active_gens(); - auto nb_lines = grid->get_nb_active_arcs(); - auto nb_buses = grid->get_nb_active_nodes(); + auto bus_pairs = grid.get_bus_pairs(); + auto nb_bus_pairs = grid.get_nb_active_bus_pairs(); + auto nb_gen = grid.get_nb_active_gens(); + auto nb_lines = grid.get_nb_active_arcs(); + auto nb_buses = grid.get_nb_active_nodes(); DebugOn("nb gens = " << nb_gen << endl); DebugOn("nb lines = " << nb_lines << endl); DebugOn("nb buses = " << nb_buses << endl); @@ -81,29 +81,29 @@ int main (int argc, char * argv[]) /** Variables */ /* power generation variables */ - var Pg("Pg", grid->pg_min, grid->pg_max); - var Qg ("Qg", grid->qg_min, grid->qg_max); - SOCP.add_var(Pg.in(grid->gens)); - SOCP.add_var(Qg.in(grid->gens)); + var Pg("Pg", grid.pg_min, grid.pg_max); + var Qg ("Qg", grid.qg_min, grid.qg_max); + SOCP.add_var(Pg.in(grid.gens)); + SOCP.add_var(Qg.in(grid.gens)); /* power flow variables */ - var Pf_from("Pf_from", grid->S_max); - var Qf_from("Qf_from", grid->S_max); - var Pf_to("Pf_to", grid->S_max); - var Qf_to("Qf_to", grid->S_max); - SOCP.add_var(Pf_from.in(grid->arcs)); - SOCP.add_var(Qf_from.in(grid->arcs)); - SOCP.add_var(Pf_to.in(grid->arcs)); - SOCP.add_var(Qf_to.in(grid->arcs)); + var Pf_from("Pf_from", grid.S_max); + var Qf_from("Qf_from", grid.S_max); + var Pf_to("Pf_to", grid.S_max); + var Qf_to("Qf_to", grid.S_max); + SOCP.add_var(Pf_from.in(grid.arcs)); + SOCP.add_var(Qf_from.in(grid.arcs)); + SOCP.add_var(Pf_to.in(grid.arcs)); + SOCP.add_var(Qf_to.in(grid.arcs)); /* Real part of Wij = ViVj */ - var R_Wij("R_Wij", grid->wr_min, grid->wr_max); + var R_Wij("R_Wij", grid.wr_min, grid.wr_max); /* Imaginary part of Wij = ViVj */ - var Im_Wij("Im_Wij", grid->wi_min, grid->wi_max); + var Im_Wij("Im_Wij", grid.wi_min, grid.wi_max); /* Magnitude of Wii = Vi^2 */ - var Wii("Wii", grid->w_min, grid->w_max); - SOCP.add_var(Wii.in(grid->nodes)); + var Wii("Wii", grid.w_min, grid.w_max); + SOCP.add_var(Wii.in(grid.nodes)); SOCP.add_var(R_Wij.in(bus_pairs)); SOCP.add_var(Im_Wij.in(bus_pairs)); @@ -112,8 +112,8 @@ int main (int argc, char * argv[]) Wii.initialize_all(1.001); /** Objective */ - auto obj = product(grid->c1, Pg) + product(grid->c2, power(Pg,2)) + sum(grid->c0); - SOCP.min(obj.in(grid->gens)); + auto obj = product(grid.c1, Pg) + product(grid.c2, power(Pg,2)) + sum(grid.c0); + SOCP.min(obj.in(grid.gens)); /** Constraints */ /* Second-order cone constraints */ @@ -123,67 +123,69 @@ int main (int argc, char * argv[]) /* Flow conservation */ Constraint KCL_P("KCL_P"); - KCL_P = sum(Pf_from.out_arcs()) + sum(Pf_to.in_arcs()) + grid->pl - sum(Pg.in_gens()) + grid->gs*Wii; - SOCP.add_constraint(KCL_P.in(grid->nodes) == 0); + KCL_P = sum(Pf_from.out_arcs()) + sum(Pf_to.in_arcs()) + grid.pl - sum(Pg.in_gens()) + grid.gs*Wii; + SOCP.add_constraint(KCL_P.in(grid.nodes) == 0); Constraint KCL_Q("KCL_Q"); - KCL_Q = sum(Qf_from.out_arcs()) + sum(Qf_to.in_arcs()) + grid->ql - sum(Qg.in_gens()) - grid->bs*Wii; - SOCP.add_constraint(KCL_Q.in(grid->nodes) == 0); + KCL_Q = sum(Qf_from.out_arcs()) + sum(Qf_to.in_arcs()) + grid.ql - sum(Qg.in_gens()) - grid.bs*Wii; + SOCP.add_constraint(KCL_Q.in(grid.nodes) == 0); /* AC Power Flow */ Constraint Flow_P_From("Flow_P_From"); - Flow_P_From = Pf_from - (grid->g_ff*Wii.from() + grid->g_ft*R_Wij.in_pairs() + grid->b_ft*Im_Wij.in_pairs()); - SOCP.add_constraint(Flow_P_From.in(grid->arcs) == 0); + Flow_P_From = Pf_from - (grid.g_ff*Wii.from() + grid.g_ft*R_Wij.in_pairs() + grid.b_ft*Im_Wij.in_pairs()); + SOCP.add_constraint(Flow_P_From.in(grid.arcs) == 0); Constraint Flow_P_To("Flow_P_To"); - Flow_P_To = Pf_to - (grid->g_tt*Wii.to() + grid->g_tf*R_Wij.in_pairs() - grid->b_tf*Im_Wij.in_pairs()); - SOCP.add_constraint(Flow_P_To.in(grid->arcs) == 0); + Flow_P_To = Pf_to - (grid.g_tt*Wii.to() + grid.g_tf*R_Wij.in_pairs() - grid.b_tf*Im_Wij.in_pairs()); + SOCP.add_constraint(Flow_P_To.in(grid.arcs) == 0); Constraint Flow_Q_From("Flow_Q_From"); - Flow_Q_From = Qf_from - (grid->g_ft*Im_Wij.in_pairs() - grid->b_ff*Wii.from() - grid->b_ft*R_Wij.in_pairs()); - SOCP.add_constraint(Flow_Q_From.in(grid->arcs) == 0); + Flow_Q_From = Qf_from - (grid.g_ft*Im_Wij.in_pairs() - grid.b_ff*Wii.from() - grid.b_ft*R_Wij.in_pairs()); + SOCP.add_constraint(Flow_Q_From.in(grid.arcs) == 0); Constraint Flow_Q_To("Flow_Q_To"); - Flow_Q_To = Qf_to + (grid->b_tt*Wii.to() + grid->b_tf*R_Wij.in_pairs() + grid->g_tf*Im_Wij.in_pairs()); - SOCP.add_constraint(Flow_Q_To.in(grid->arcs) == 0); + Flow_Q_To = Qf_to + (grid.b_tt*Wii.to() + grid.b_tf*R_Wij.in_pairs() + grid.g_tf*Im_Wij.in_pairs()); + SOCP.add_constraint(Flow_Q_To.in(grid.arcs) == 0); /* Phase Angle Bounds constraints */ Constraint PAD_UB("PAD_UB"); PAD_UB = Im_Wij; - PAD_UB <= grid->tan_th_max*R_Wij; - SOCP.add_constraint(PAD_UB.in(bus_pairs)); + PAD_UB <= grid.tan_th_max*R_Wij; + SOCP.add_lazy(PAD_UB.in(bus_pairs)); Constraint PAD_LB("PAD_LB"); PAD_LB = Im_Wij; - PAD_LB >= grid->tan_th_min*R_Wij; - SOCP.add_constraint(PAD_LB.in(bus_pairs)); + PAD_LB >= grid.tan_th_min*R_Wij; + SOCP.add_lazy(PAD_LB.in(bus_pairs)); /* Thermal Limit Constraints */ Constraint Thermal_Limit_from("Thermal_Limit_from"); Thermal_Limit_from = power(Pf_from, 2) + power(Qf_from, 2); - Thermal_Limit_from <= power(grid->S_max,2); - SOCP.add_constraint(Thermal_Limit_from.in(grid->arcs)); + Thermal_Limit_from <= power(grid.S_max,2); + SOCP.add_lazy(Thermal_Limit_from.in(grid.arcs)); +// SOCP.add(Thermal_Limit_from.in(grid.arcs)); Constraint Thermal_Limit_to("Thermal_Limit_to"); Thermal_Limit_to = power(Pf_to, 2) + power(Qf_to, 2); - Thermal_Limit_to <= power(grid->S_max,2); - SOCP.add_constraint(Thermal_Limit_to.in(grid->arcs)); + Thermal_Limit_to <= power(grid.S_max,2); + SOCP.add_lazy(Thermal_Limit_to.in(grid.arcs)); +// SOCP.add(Thermal_Limit_to.in(grid.arcs)); /* Lifted Nonlinear Cuts */ Constraint LNC1("LNC1"); - LNC1 += (grid->v_min.from()+grid->v_max.from())*(grid->v_min.to()+grid->v_max.to())*(grid->sphi*Im_Wij + grid->cphi*R_Wij); - LNC1 -= grid->v_max.to()*grid->cos_d*(grid->v_min.to()+grid->v_max.to())*Wii.from(); - LNC1 -= grid->v_max.from()*grid->cos_d*(grid->v_min.from()+grid->v_max.from())*Wii.to(); - LNC1 -= grid->v_max.from()*grid->v_max.to()*grid->cos_d*(grid->v_min.from()*grid->v_min.to() - grid->v_max.from()*grid->v_max.to()); - SOCP.add_constraint(LNC1.in(bus_pairs) >= 0); + LNC1 += (grid.v_min.from()+grid.v_max.from())*(grid.v_min.to()+grid.v_max.to())*(grid.sphi*Im_Wij + grid.cphi*R_Wij); + LNC1 -= grid.v_max.to()*grid.cos_d*(grid.v_min.to()+grid.v_max.to())*Wii.from(); + LNC1 -= grid.v_max.from()*grid.cos_d*(grid.v_min.from()+grid.v_max.from())*Wii.to(); + LNC1 -= grid.v_max.from()*grid.v_max.to()*grid.cos_d*(grid.v_min.from()*grid.v_min.to() - grid.v_max.from()*grid.v_max.to()); + SOCP.add_lazy(LNC1.in(bus_pairs) >= 0); Constraint LNC2("LNC2"); - LNC2 += (grid->v_min.from()+grid->v_max.from())*(grid->v_min.to()+grid->v_max.to())*(grid->sphi*Im_Wij + grid->cphi*R_Wij); - LNC2 -= grid->v_min.to()*grid->cos_d*(grid->v_min.to()+grid->v_max.to())*Wii.from(); - LNC2 -= grid->v_min.from()*grid->cos_d*(grid->v_min.from()+grid->v_max.from())*Wii.to(); - LNC2 += grid->v_min.from()*grid->v_min.to()*grid->cos_d*(grid->v_min.from()*grid->v_min.to() - grid->v_max.from()*grid->v_max.to()); - SOCP.add_constraint(LNC2.in(bus_pairs) >= 0); + LNC2 += (grid.v_min.from()+grid.v_max.from())*(grid.v_min.to()+grid.v_max.to())*(grid.sphi*Im_Wij + grid.cphi*R_Wij); + LNC2 -= grid.v_min.to()*grid.cos_d*(grid.v_min.to()+grid.v_max.to())*Wii.from(); + LNC2 -= grid.v_min.from()*grid.cos_d*(grid.v_min.from()+grid.v_max.from())*Wii.to(); + LNC2 += grid.v_min.from()*grid.v_min.to()*grid.cos_d*(grid.v_min.from()*grid.v_min.to() - grid.v_max.from()*grid.v_max.to()); + SOCP.add_lazy(LNC2.in(bus_pairs) >= 0); /* Solver selection */ @@ -217,7 +219,7 @@ int main (int argc, char * argv[]) } /** Uncomment next line to print expanded model */ /* SOCP.print_expanded(); */ - string out = "DATA_OPF, " + grid->_name + ", " + to_string(nb_buses) + ", " + to_string(nb_lines) +", " + to_string(SOCP._obj_val) + ", " + to_string(-numeric_limits::infinity()) + ", " + to_string(solve_time) + ", LocalOptimal, " + to_string(total_time); + string out = "DATA_OPF, " + grid._name + ", " + to_string(nb_buses) + ", " + to_string(nb_lines) +", " + to_string(SOCP._obj_val) + ", " + to_string(-numeric_limits::infinity()) + ", " + to_string(solve_time) + ", LocalOptimal, " + to_string(total_time); DebugOn(out < _dual ; /**< Lagrange multipliers at a KKT point */ + bool _all_active = true; + vector _active; + shared_ptr _all_lazy; + vector _lazy; + bool _all_satisfied = true; + vector _violated; + /** Constructor */ //@{ @@ -53,18 +60,25 @@ namespace gravity { Constraint& operator =(const func_& rhs); /* Accessors */ + size_t get_nb_instances() const; string get_name() const; int get_type() const; double get_rhs() const; bool is_active(unsigned inst = 0) const; bool is_convex() const; bool is_concave() const; + bool is_ineq() const; size_t get_id_inst(size_t ind) const; /* Modifiers */ + void make_lazy() { + *_all_lazy = true; + _lazy.resize(_nb_instances,true); + } + Constraint& in(const node_pairs& np){ this->func_::in(np); return *this; diff --git a/include/gravity/func.h b/include/gravity/func.h index 9b56a4afb..f1f1e40c8 100755 --- a/include/gravity/func.h +++ b/include/gravity/func.h @@ -475,8 +475,6 @@ namespace gravity { size_t _nb_instances = 1; /**< Number of different instances this constraint has (different indices, constant coefficients and bounds, but same structure).>>**/ - bool _all_active = true; - vector _active; bool _new = true; /**< Will become false once this function is added to a program. Can be useful for iterative model solving. */ bool _is_constraint = false; bool _is_hessian = false; diff --git a/include/gravity/model.h b/include/gravity/model.h index 003d202a0..33b6c268d 100755 --- a/include/gravity/model.h +++ b/include/gravity/model.h @@ -44,6 +44,7 @@ namespace gravity { void add_param(param_* v); //Add variables without reallocating memory public: + bool _has_lazy = false; /*<< Has lazy constraints. */ bool _built = false; /* Indicates if this model has been already built */ bool _first_run = true; /* Indicates if this model has been already */ @@ -100,6 +101,8 @@ namespace gravity { size_t get_nb_cons() const; + size_t get_nb_ineq() const; + size_t get_nb_nnz_g(); size_t get_nb_nnz_h(); @@ -204,8 +207,9 @@ namespace gravity { void del_param(const param_& v); - - void add_constraint(const Constraint& c); + void add_lazy(const Constraint& c); + void add(const Constraint& c); + shared_ptr add_constraint(const Constraint& c); shared_ptr embed(shared_ptr f);/**< Transfer all variables and parameters to the model, useful for a centralized memory management. */ void embed(expr& e);/**< Transfer all variables and parameters to the model, useful for a centralized memory management. */ @@ -216,6 +220,7 @@ namespace gravity { void set_objective(pair p); void set_objective_type(ObjectiveType); void init_indices();// Initialize the indices of all variables involved in the model + void reindex(); /*<< Reindexes the constraints after violated ones have been detected and aded to the formulation */ bool has_violated_constraints(double tol); /*<< Returns true if some constraints are violated by the current solution with tolerance tol */ void check_feasible(const double* x); void reset_funcs(); diff --git a/include/gravity/solver.h b/include/gravity/solver.h index fc27e5ec8..e4e0f116a 100755 --- a/include/gravity/solver.h +++ b/include/gravity/solver.h @@ -36,7 +36,7 @@ namespace gravity { public: Model* _model = nullptr; Program* _prog = nullptr; - SolverType _stype; + SolverType _stype; double _tol = 1e-6; /*<< Solver tolerance. */ /** Constructor */ diff --git a/src/GurobiProgram.cpp b/src/GurobiProgram.cpp index 909e34505..fe5f8fb11 100755 --- a/src/GurobiProgram.cpp +++ b/src/GurobiProgram.cpp @@ -201,7 +201,7 @@ void GurobiProgram::create_grb_constraints(){ double coeff; for(auto& p: _model->_cons){ auto c = p.second; - if (!c->_new) { + if (!c->_new && c->_all_satisfied) { continue; } c->_new = false; @@ -224,63 +224,76 @@ void GurobiProgram::create_grb_constraints(){ nb_inst = c->_nb_instances; if (c->is_linear()) { for (int i = 0; i< nb_inst; i++){ - linlhs = 0; - for (auto& it1: c->get_lterms()) { - lterm = 0; - if (it1.second._coef->_is_transposed) { - auto dim =it1.second._p->get_dim(i); - for (int j = 0; jget_id() + it1.second._p->get_id_inst(i,j)]; + if (c->_violated[i]) { + linlhs = 0; + for (auto& it1: c->get_lterms()) { + lterm = 0; + if (it1.second._coef->_is_transposed) { + auto dim =it1.second._p->get_dim(i); + for (int j = 0; jget_id() + it1.second._p->get_id_inst(i,j)]; + lterm += coeff*gvar1; + } + } + else { + coeff = poly_eval(it1.second._coef,i); + gvar1 = _grb_vars[it1.second._p->get_id() + it1.second._p->get_id_inst(i)]; lterm += coeff*gvar1; } + if (!it1.second._sign) { + lterm *= -1; + } + linlhs += lterm; } - else { - coeff = poly_eval(it1.second._coef,i); - gvar1 = _grb_vars[it1.second._p->get_id() + it1.second._p->get_id_inst(i)]; - lterm += coeff*gvar1; - } - if (!it1.second._sign) { - lterm *= -1; - } - linlhs += lterm; + linlhs += poly_eval(c->get_cst(), i); + grb_mod->addConstr(linlhs,sense,c->get_rhs(),c->get_name()+"_"+to_string(i)); } - linlhs += poly_eval(c->get_cst(), i); - grb_mod->addConstr(linlhs,sense,c->get_rhs(),c->get_name()+"_"+to_string(i)); } } else { for (int i = 0; i< nb_inst; i++){ - quadlhs = 0; - for (auto& it1: c->get_lterms()) { - lterm = 0; - if (it1.second._coef->_is_transposed) { - auto dim =it1.second._p->get_dim(i); - for (int j = 0; jget_id() + it1.second._p->get_id_inst(i,j)]; + if (c->_violated[i]) { + quadlhs = 0; + for (auto& it1: c->get_lterms()) { + lterm = 0; + if (it1.second._coef->_is_transposed) { + auto dim =it1.second._p->get_dim(i); + for (int j = 0; jget_id() + it1.second._p->get_id_inst(i,j)]; + lterm += coeff*gvar1; + } + } + else { + coeff = poly_eval(it1.second._coef,i); + gvar1 = _grb_vars[it1.second._p->get_id() + it1.second._p->get_id_inst(i)]; lterm += coeff*gvar1; } + if (!it1.second._sign) { + lterm *= -1; + } + quadlhs += lterm; } - else { - coeff = poly_eval(it1.second._coef,i); - gvar1 = _grb_vars[it1.second._p->get_id() + it1.second._p->get_id_inst(i)]; - lterm += coeff*gvar1; - } - if (!it1.second._sign) { - lterm *= -1; - } - quadlhs += lterm; - } - for (auto& it1: c->get_qterms()) { - gvar1 = _grb_vars[it1.second._p->first->get_id() + it1.second._p->first->get_id_inst(i)]; - gvar2 = _grb_vars[it1.second._p->second->get_id() + it1.second._p->second->get_id_inst(i)]; - if (it1.second._coef->_is_transposed) { - auto dim =it1.second._p->first->get_dim(i); - for (int j = 0; jfirst->get_id() + it1.second._p->first->get_id_inst(i,j)]; - gvar2 = _grb_vars[it1.second._p->second->get_id() + it1.second._p->second->get_id_inst(i,j)]; + for (auto& it1: c->get_qterms()) { + gvar1 = _grb_vars[it1.second._p->first->get_id() + it1.second._p->first->get_id_inst(i)]; + gvar2 = _grb_vars[it1.second._p->second->get_id() + it1.second._p->second->get_id_inst(i)]; + if (it1.second._coef->_is_transposed) { + auto dim =it1.second._p->first->get_dim(i); + for (int j = 0; jfirst->get_id() + it1.second._p->first->get_id_inst(i,j)]; + gvar2 = _grb_vars[it1.second._p->second->get_id() + it1.second._p->second->get_id_inst(i,j)]; + if (!it1.second._sign) { + quadlhs += -1*coeff*gvar1*gvar2; + } + else { + quadlhs += coeff*gvar1*gvar2; + } + } + } + else { + coeff = poly_eval(it1.second._coef,i); if (!it1.second._sign) { quadlhs += -1*coeff*gvar1*gvar2; } @@ -289,18 +302,9 @@ void GurobiProgram::create_grb_constraints(){ } } } - else { - coeff = poly_eval(it1.second._coef,i); - if (!it1.second._sign) { - quadlhs += -1*coeff*gvar1*gvar2; - } - else { - quadlhs += coeff*gvar1*gvar2; - } - } + quadlhs += poly_eval(c->get_cst(), i); + grb_mod->addQConstr(quadlhs,sense,c->get_rhs(),c->get_name()+"_"+to_string(i)); } - quadlhs += poly_eval(c->get_cst(), i); - grb_mod->addQConstr(quadlhs,sense,c->get_rhs(),c->get_name()+"_"+to_string(i)); } } } diff --git a/src/IpoptProgram.cpp b/src/IpoptProgram.cpp index adf8b95f2..f8dada4dc 100755 --- a/src/IpoptProgram.cpp +++ b/src/IpoptProgram.cpp @@ -43,6 +43,7 @@ void IpoptProgram::finalize_solution(Ipopt::SolverReturn status , ) { _model->set_x(x); + _model->compute_funcs(); // _model->check_feasible(x); if(_model->_objt==maximize){ _model->_obj *= -1.; @@ -50,8 +51,11 @@ void IpoptProgram::finalize_solution(Ipopt::SolverReturn status , _model->_obj_val = _model->_obj.eval(); for (auto &cp: _model->_cons) { cp.second->_dual.resize(cp.second->_nb_instances); + auto idx = 0; for (unsigned inst = 0; inst < cp.second->_nb_instances; inst++) { - cp.second->_dual[inst] = lambda[cp.second->_id + inst]; + if (!*cp.second->_all_lazy || !cp.second->_lazy[inst]) { + cp.second->_dual[inst] = lambda[cp.second->_id + idx++]; + } } } for (auto &vp: _model->_vars) { diff --git a/src/constraint.cpp b/src/constraint.cpp index 61ef5d303..fdc3a0057 100755 --- a/src/constraint.cpp +++ b/src/constraint.cpp @@ -19,6 +19,7 @@ Constraint::Constraint(string name, ConstraintType ctype){ _ctype = ctype; _rhs = 0; _is_constraint = true; + _all_lazy = make_shared(false); }; Constraint::Constraint(const Constraint& c):func_(c){ @@ -27,6 +28,7 @@ Constraint::Constraint(const Constraint& c):func_(c){ _rhs = c._rhs; _id = c._id; _is_constraint = true; + _all_lazy = make_shared(false); } //@} @@ -40,6 +42,19 @@ Constraint::~Constraint(){}; /* Accessors */ +size_t Constraint::get_nb_instances() const{ + size_t nb = 0; + if (!*_all_lazy) { + return _nb_instances; + } + for (unsigned i = 0; i<_nb_instances; i++) { + if (!_lazy[i]) { + nb++; + } + } + return nb; +} + double Constraint::get_rhs() const{ return _rhs; }; @@ -180,6 +195,10 @@ bool Constraint::is_concave() const{ return ((_all_convexity==convex_ &&_ctype==geq) || (_all_convexity==concave_ &&_ctype==leq)); } +bool Constraint::is_ineq() const{ + return (_ctype==leq || _ctype==geq); +} + void Constraint::print(){ cout << _name; if (is_convex()) { diff --git a/src/func.cpp b/src/func.cpp index 9f0ac51d7..ced753537 100755 --- a/src/func.cpp +++ b/src/func.cpp @@ -2059,10 +2059,10 @@ namespace gravity{ if (is_number() && all_zeros(poly_to_str(this))){ return true; } -// if (is_param() || is_var()) { -// auto p_c = (param_*)this; -// return p_c->get_all_sign()==zero_; -// } + if (is_param() || is_var()) { + auto p_c = (param_*)this; + return p_c->get_all_sign()==zero_; + } // if (is_uexpr() || is_bexpr()) { // auto e_p = (expr*)this; // return e_p->get_all_sign()==zero_; @@ -7372,14 +7372,14 @@ namespace gravity{ if(is_number() && !_is_vector && !_is_transposed && poly_eval(this)==0){ return true; } -// if (_ftype==const_ && _cst->is_zero()){ -// for (auto& it:*_params) { -// if (!it.second.first->is_zero()) { -// return false; -// } -// } -// return true; -// } + if (_ftype==const_ && _cst->is_zero()){ + for (auto& it:*_params) { + if (!it.second.first->is_zero()) { + return false; + } + } + return true; + } return false; } diff --git a/src/model.cpp b/src/model.cpp index 7db6143c0..1dd33b26e 100755 --- a/src/model.cpp +++ b/src/model.cpp @@ -42,7 +42,17 @@ size_t Model::get_nb_vars() const{ size_t Model::get_nb_cons() const{ size_t n = 0; for (auto &cp:_cons) { - n += cp.second->_nb_instances; + n += cp.second->get_nb_instances(); + } + return n; +}; + +size_t Model::get_nb_ineq() const{ + size_t n = 0; + for (auto &cp:_cons) { + if (cp.second->is_ineq()) { + n += cp.second->get_nb_instances(); + } } return n; }; @@ -54,7 +64,9 @@ size_t Model::get_nb_nnz_g(){ auto c = cp.second; auto nb_inst = c->_nb_instances; for (unsigned inst = 0; instget_nb_vars(inst); + if (!*c->_all_lazy || !c->_lazy[inst]) { + _nnz_g += c->get_nb_vars(inst); + } } } return _nnz_g; @@ -65,29 +77,35 @@ size_t Model::get_nb_nnz_g(){ size_t Model::get_nb_nnz_h(){ size_t idx = 0; bool idx_inc = false; + Constraint* c = nullptr; for (auto &pairs: _hess_link) { auto f_pair = *pairs.second.begin(); // for (auto &f_pair:pairs.second) { idx_inc = false; auto f = f_pair.first; + if (f->_is_constraint) { + c = (Constraint*)f; + } auto d2f = f_pair.second; size_t nb_inst = f->_nb_instances; for (unsigned inst = 0; inst_is_matrix) { - for (unsigned i = 0; i < d2f->_dim[0]; i++) { - for (unsigned j = i; j < d2f->_dim[1]; j++) { + if (!(f->_is_constraint && *c->_all_lazy && c->_lazy[inst])) { + if (d2f->_is_matrix) { + for (unsigned i = 0; i < d2f->_dim[0]; i++) { + for (unsigned j = i; j < d2f->_dim[1]; j++) { + idx++; + } + } + } + else if(d2f->_is_vector){ + for (unsigned j = 0; j < d2f->_dim[0]; j++) { idx++; } } - } - else if(d2f->_is_vector){ - for (unsigned j = 0; j < d2f->_dim[0]; j++) { + else { idx++; } } - else { - idx++; - } } } // } @@ -101,7 +119,7 @@ void Model::add_indices(const string& constr_name, const node_pairs& indices){ map> new_cons; c->add_indices_in(indices); for (auto &cc : _cons) { - if (cc.second->_id > c->_id) { + if (cc.second->_id > c->_id) {//Ids of all constraints coming after c will be shifted by the size of indices cc.second->_id += indices._keys.size(); } new_cons[cc.second->_id] = cc.second; @@ -122,6 +140,25 @@ param_* Model::get_var(const string& vname) const{ /* Modifiers */ +void Model::reindex(){ + unsigned cid = 0, new_cid = 0, nb_inst = 0; + shared_ptr c = nullptr; + map> new_cons; + for(auto& c_p: _cons_name) + { + c = c_p.second; + cid = c->_id; + if (cid!=new_cid) { + c->_id = new_cid; + } + new_cons[c->_id] = c; + nb_inst = c->get_nb_instances(); + new_cid = c->_id+nb_inst; + } + _cons = new_cons; + _nb_cons = get_nb_cons(); +} + void Model::init_indices(){// Initialize the indices of all variables involved in the model param_* v= nullptr; size_t idx = 0; @@ -195,7 +232,16 @@ void Model::del_param(const param_& v){ } }; -void Model::add_constraint(const Constraint& c){ +void Model::add(const Constraint& c){ + add_constraint(c); +} +void Model::add_lazy(const Constraint& c){ + auto newc = add_constraint(c); + newc->make_lazy(); + _has_lazy = true; +} + +shared_ptr Model::add_constraint(const Constraint& c){ if (_cons_name.count(c.get_name())==0) { auto newc = make_shared(c); if (newc->is_constant()) { @@ -219,7 +265,7 @@ void Model::add_constraint(const Constraint& c){ break; } DebugOn("WARNING: Adding redundant constant constraint, Gravity will be ignoring it.\n"); - return; + return newc; } newc->update_to_str(); @@ -393,30 +439,36 @@ void Model::add_constraint(const Constraint& c){ } // embed(); // } + if(*newc->_all_lazy){ + newc->_lazy.resize(newc->_nb_instances,true); + } + newc->_violated.resize(newc->_nb_instances,true); newc->_id = get_nb_cons(); // embed(newc); _cons_name[c.get_name()] = newc; _cons[newc->_id] = newc; _built = false; // newc->print_expanded(); + int nb_inst = c._nb_instances; + _nb_cons += nb_inst; + for (unsigned inst = 0; inst c = nullptr; - for(auto& c_p: _cons) + bool violated = false; + for(auto& c_p: _cons_name) { c = c_p.second; - cid = c->_id; +// cid = c->_id; + nb_inst = c->_nb_instances; + nb_viol = 0; + nb_active = 0; + c->_all_satisfied = true; + c->_violated.resize(nb_inst); + c->_active.resize(nb_inst); switch (c->get_type()) { case eq: - nb_inst = c->_nb_instances; for (unsigned inst=0; insteval(inst) - c->_rhs); if(diff > tol) { DebugOff("Violated equation: "); DebugOff(c->to_str(inst)); DebugOff(", violation = "<< diff << endl); - return true; + nb_viol++; + if (*c->_all_lazy) { + c->_all_satisfied = false; + c->_violated[inst] = true; + violated = true; + c->_lazy[inst] = false; + } + else { +// throw runtime_error("Non-lazy constraint is violated, solution declared optimal by solver!\n" + c->to_str(inst)); + } + c->_violated[inst] = true; } + else { + c->_violated[inst] = false; + } +// nb_active++; } break; case leq: - nb_inst = c->_nb_instances; for (unsigned inst=0; insteval(inst) - c->_rhs); + c->_violated[inst] = false; + diff = c->eval(inst) - c->_rhs; if(diff > tol) { DebugOff("Violated inequality: "); DebugOff(c->to_str(inst)); DebugOff(", violation = "<< diff << endl); - return true; + nb_viol++; + if (*c->_all_lazy) { + c->_all_satisfied = false; + c->_violated[inst] = true; + violated = true; + c->_lazy[inst] = false; + } + else { +// throw runtime_error("Non-lazy constraint is violated, solution declared optimal by solver!\n" + c->to_str(inst)); + } + } + else if (fabs(diff)>tol) { + c->_active[inst] = false; +// if (*c->_all_lazy) { +// c->_lazy[inst] = true; +// } + } + else { + nb_active++; } } break; case geq: - nb_inst = c->_nb_instances; for (unsigned inst=0; insteval(inst) - c->_rhs); + c->_violated[inst] = false; + diff = c->eval(inst) - c->_rhs; if(diff < -tol) { DebugOff("Violated inequality: "); DebugOff(c->to_str(inst)); DebugOff(", violation = "<< diff << endl); - return true; + nb_viol++; + if (*c->_all_lazy) { + c->_all_satisfied = false; + c->_violated[inst] = true; + violated = true; + c->_lazy[inst] = false; + } + else { +// throw runtime_error("Non-lazy constraint is violated, solution declared optimal by solver!\n" + c->to_str(inst)); + } + } + else if (fabs(diff)> tol) { + c->_active[inst] = false; +// if (*c->_all_lazy) { +// c->_lazy[inst] = true; +// } + } + else { + nb_active++; } } break; @@ -516,8 +628,19 @@ bool Model::has_violated_constraints(double tol){ default: break; } + nb_viol_all += nb_viol; + nb_active_all += nb_active; + if (nb_viol>0) { + DebugOn("Percentage of violated constraints for " << c->get_name() << " = " << to_string(100.*nb_viol/nb_inst) << "%\n"); + } + if (c->get_type()!=eq) { + DebugOff("Percentage of active constraints for " << c->get_name() << " = " << to_string(100.*nb_active/nb_inst) << "%\n"); + } } - return false; + DebugOn("Total percentage of violated constraints = " << to_string(100.*nb_viol_all/_nb_cons) << "%\n"); + auto nb_ineq = get_nb_ineq(); + DebugOn("Total percentage of active constraints = " << to_string(100.*nb_active_all/nb_ineq) << "%\n"); + return violated; } void Model::check_feasible(const double* x){ @@ -913,11 +1036,14 @@ void compute_constrs(vector>& v, double* res, int i, int for (unsigned idx = i; idx < j; idx++) { auto c = v[idx]; auto nb_ins = c->_nb_instances; + auto id = 0; for (int inst = 0; inst< nb_ins; inst++){ - res[c->_id+inst] = c->eval(inst); - DebugOff("Accessing res at position " << c->_id+inst << endl); + if (!*c->_all_lazy || !c->_lazy[inst]) { + res[c->_id+id++] = c->eval(inst); + DebugOff("Accessing res at position " << c->_id+inst << endl); // _cons_vals[index++] = res[c->_id+inst]; - DebugOff("g[" << to_string(c->_id+inst) << "] = " << to_string(res[c->_id+inst]) << endl); + DebugOff("g[" << to_string(c->_id+inst) << "] = " << to_string(res[c->_id+inst]) << endl); + } } } } @@ -995,7 +1121,7 @@ void Model::fill_in_cstr(const double* x , double* res, bool new_x){ void Model::fill_in_jac_nnz(int* iRow , int* jCol){ - size_t idx=0; + size_t idx=0, id = 0; size_t cid = 0; size_t vid = 0; Constraint* c = NULL; @@ -1009,21 +1135,24 @@ void Model::fill_in_jac_nnz(int* iRow , int* jCol){ for (auto &v_p: c->get_vars()){ v = v_p.second.first.get(); vid = v->get_id(); + id = 0; for (int inst = 0; inst< nb_ins; inst++){ - cid = c->_id+inst; - if (v->_is_vector) { - auto dim = v->get_dim(inst); - for (int j = 0; j_all_lazy || !c->_lazy[inst]) { + cid = c->_id+id++; + if (v->_is_vector) { + auto dim = v->get_dim(inst); + for (int j = 0; jget_id_inst(inst, j); + idx++; + } + } + else { iRow[idx] = cid; - jCol[idx] = vid + v->get_id_inst(inst, j); + jCol[idx] = vid + v->get_id_inst(inst); idx++; } } - else { - iRow[idx] = cid; - jCol[idx] = vid + v->get_id_inst(inst); - idx++; - } } } } @@ -1034,7 +1163,7 @@ void Model::fill_in_jac_nnz(int* iRow , int* jCol){ void compute_jac(vector>& vec, double* res, int i, int j, bool first_call, vector& jac_vals){ - size_t cid = 0; + size_t cid = 0, id_inst = 0; unique_id vid; shared_ptr c = NULL; param_* v = NULL; @@ -1043,13 +1172,16 @@ void compute_jac(vector>& vec, double* res, int i, int j, for (unsigned id = i; id < j; id++) { c = vec[id]; auto nb_ins = c->_nb_instances; + id_inst = 0; if (c->is_linear() && !first_call) { // if (false) { DebugOff("Linear constraint, using stored jacobian!\n"); for (unsigned i = 0; iget_nb_vars(i); j++) { - res[idx] = jac_vals[idx]; - idx++; + if (!*c->_all_lazy || !c->_lazy[i]) { + for (unsigned j = 0; jget_nb_vars(i); j++) { + res[idx] = jac_vals[idx]; + idx++; + } } } } @@ -1058,21 +1190,24 @@ void compute_jac(vector>& vec, double* res, int i, int j, v = v_p.second.first.get(); vid = v->_unique_id; dfdx = c->get_stored_derivative(vid); + id_inst = 0; for (int inst = 0; inst< nb_ins; inst++){ - cid = c->_id+inst; - if (v->_is_vector) { - auto dim = v->get_dim(inst); - for (int j = 0; jeval(inst,j); + if (!*c->_all_lazy || !c->_lazy[inst]) { + cid = c->_id+id_inst++; + if (v->_is_vector) { + auto dim = v->get_dim(inst); + for (int j = 0; jeval(inst,j); + jac_vals[idx] = res[idx]; + idx++; + } + } + else { + res[idx] = dfdx->eval(inst); jac_vals[idx] = res[idx]; idx++; } } - else { - res[idx] = dfdx->eval(inst); - jac_vals[idx] = res[idx]; - idx++; - } } } } @@ -1129,13 +1264,16 @@ void Model::fill_in_jac(const double* x , double* res, bool new_x){ { c = c_p.second.get(); auto nb_ins = c->_nb_instances; + unsigned id = 0; if (c->is_linear() && !_first_call_jac) { // if (false) { DebugOff("Linear constraint, using stored jacobian!\n"); for (unsigned i = 0; iget_nb_vars(i); j++) { - res[idx] = _jac_vals[idx]; - idx++; + if (!*c->_all_lazy || !c->_lazy[i]) { + for (unsigned j = 0; jget_nb_vars(i); j++) { + res[idx] = _jac_vals[idx]; + idx++; + } } } } @@ -1147,47 +1285,51 @@ void Model::fill_in_jac(const double* x , double* res, bool new_x){ dfdx = c->get_stored_derivative(vid); if (dfdx->is_number()) { for (int inst = 0; inst< nb_ins; inst++){ - cid = c->_id+inst; - if (v->_is_vector) { - auto dim = v->get_dim(inst); - for (int j = 0; j_all_lazy || !c->_lazy[inst]) { + cid = c->_id+id++; + if (v->_is_vector) { + auto dim = v->get_dim(inst); + for (int j = 0; j_val->at(0); + _jac_vals[idx] = res[idx]; + idx++; + } + } + else { res[idx] = dfdx->_val->at(0); _jac_vals[idx] = res[idx]; idx++; - } - } - else { - res[idx] = dfdx->_val->at(0); - _jac_vals[idx] = res[idx]; - idx++; + } } } } else { for (int inst = 0; inst< nb_ins; inst++){ - cid = c->_id+inst; - if (v->_is_vector) { - auto dim = v->get_dim(inst); - if (dfdx->_is_matrix) { - for (int j = 0; jget_val(j,inst); - _jac_vals[idx] = res[idx]; - idx++; + if (!*c->_all_lazy || !c->_lazy[inst]) { + cid = c->_id+id++; + if (v->_is_vector) { + auto dim = v->get_dim(inst); + if (dfdx->_is_matrix) { + for (int j = 0; jget_val(j,inst); + _jac_vals[idx] = res[idx]; + idx++; + } + } + else { + for (int j = 0; jget_val(j);//TODO check double indexed funcs + _jac_vals[idx] = res[idx]; + idx++; + } } } else { - for (int j = 0; jget_val(j);//TODO check double indexed funcs - _jac_vals[idx] = res[idx]; - idx++; - } + res[idx] = dfdx->get_val(inst); + _jac_vals[idx] = res[idx]; + idx++; } } - else { - res[idx] = dfdx->get_val(inst); - _jac_vals[idx] = res[idx]; - idx++; - } } } } @@ -1248,10 +1390,11 @@ void Model::fill_in_var_linearity(Ipopt::TNLP::LinearityType* param_types){ void Model::fill_in_cstr_linearity(Ipopt::TNLP::LinearityType* const_types){ Constraint* c = nullptr; bool lin = false; - size_t cid = 0; + size_t cid = 0, id = 0; for(auto& c_p :_cons) { c = c_p.second.get(); + id = 0; if (c->is_linear() || c->is_constant()) { lin = true; } @@ -1260,12 +1403,14 @@ void Model::fill_in_cstr_linearity(Ipopt::TNLP::LinearityType* const_types){ } auto nb_ins = c->_nb_instances; for (int i = 0; i< nb_ins; i++){ - cid = c->_id +i; - if (lin) { - const_types[cid]=Ipopt::TNLP::LINEAR; - } - else { - const_types[cid] = Ipopt::TNLP::NON_LINEAR; + if (!*c->_all_lazy || !c->_lazy[i]) { + cid = c->_id+id++; + if (lin) { + const_types[cid]=Ipopt::TNLP::LINEAR; + } + else { + const_types[cid] = Ipopt::TNLP::NON_LINEAR; + } } } } @@ -1300,33 +1445,36 @@ void Model::fill_in_hess_nnz(int* iRow , int* jCol){ auto d2f = f_pair.second; size_t nb_inst = f->_nb_instances; for (unsigned inst = 0; inst_is_matrix) { - for (unsigned i = 0; i < d2f->_dim[0]; i++) { - for (unsigned j = i; j < d2f->_dim[1]; j++) { + if (!(f->_is_constraint && *c->_all_lazy && c->_lazy[inst])) { + + if (d2f->_is_matrix) { + for (unsigned i = 0; i < d2f->_dim[0]; i++) { + for (unsigned j = i; j < d2f->_dim[1]; j++) { + idx_all++; + iRow[idx] = vid + vi->get_id_inst(i); + jCol[idx] = vjd + vj->get_id_inst(j); + idx++; + } + } + } + else if(d2f->_is_vector){ + // if (d2f->_dim[0] != d2f->_nb_instances) { + // throw invalid_argument("error"); + // } + for (unsigned j = 0; j < d2f->_dim[0]; j++) { idx_all++; - iRow[idx] = vid + vi->get_id_inst(i); + iRow[idx] = vid + vi->get_id_inst(j); jCol[idx] = vjd + vj->get_id_inst(j); idx++; } } - } - else if(d2f->_is_vector){ -// if (d2f->_dim[0] != d2f->_nb_instances) { -// throw invalid_argument("error"); -// } - for (unsigned j = 0; j < d2f->_dim[0]; j++) { + else { idx_all++; - iRow[idx] = vid + vi->get_id_inst(j); - jCol[idx] = vjd + vj->get_id_inst(j); + iRow[idx] = vid + vi->get_id_inst(inst); + jCol[idx] = vjd + vj->get_id_inst(inst); idx++; } } - else { - idx_all++; - iRow[idx] = vid + vi->get_id_inst(inst); - jCol[idx] = vjd + vj->get_id_inst(inst); - idx++; - } } } } @@ -1360,38 +1508,46 @@ void Model::fill_in_hess(const double* x , double obj_factor, const double* lamb } auto d2f = f_pair.second; size_t nb_inst = f->_nb_instances; + unsigned id_inst = 0; for (unsigned inst = 0; inst_is_constraint) { - c_inst = c->get_id_inst(inst); - if (d2f->_is_matrix) { - for (unsigned i = 0; i < d2f->_dim[0]; i++) { - for (unsigned j = i; j < d2f->_dim[1]; j++) { - hess = d2f->get_val(i,j); + if (!*c->_all_lazy || !c->_lazy[inst]) { + idx_inc = false; + c_inst = c->get_id_inst(id_inst++); + + if (d2f->_is_matrix) { + for (unsigned i = 0; i < d2f->_dim[0]; i++) { + for (unsigned j = i; j < d2f->_dim[1]; j++) { + hess = d2f->get_val(i,j); + _hess_vals[idx_in++] = hess; + res[idx++] += lambda[c->_id + c_inst] * hess; + idx_inc = true; + } + } + } + else if(d2f->_is_vector){ + for (unsigned j = 0; j < d2f->_dim[0]; j++) { + hess = d2f->get_val(j); _hess_vals[idx_in++] = hess; - res[idx++] += lambda[c->_id + c_inst] * hess; + res[idx] += lambda[c->_id + c_inst] * hess; + idx++; idx_inc = true; + // } } } - } - else if(d2f->_is_vector){ - for (unsigned j = 0; j < d2f->_dim[0]; j++) { - hess = d2f->get_val(j); + else { + if (d2f->is_number()) { + hess = d2f->_val->at(0); + } + else { + hess = d2f->_val->at(inst); + } _hess_vals[idx_in++] = hess; res[idx] += lambda[c->_id + c_inst] * hess; - idx++; - idx_inc = true; - // } } } else { - if (d2f->is_number()) { - hess = d2f->_val->at(0); - } - else { - hess = d2f->_val->at(inst); - } - _hess_vals[idx_in++] = hess; - res[idx] += lambda[c->_id + c_inst] * hess; + idx_inc = true; } } else { @@ -1420,6 +1576,7 @@ void Model::fill_in_hess(const double* x , double obj_factor, const double* lamb return; } if ((_type==lin_m || _type==quad_m)) { /* No need to recompute Hessian for quadratic models, used stored values */ + unsigned id_inst = 0; for (auto &pairs: _hess_link) { idx_pair = idx; for (auto &f_pair:pairs.second) { @@ -1431,27 +1588,34 @@ void Model::fill_in_hess(const double* x , double obj_factor, const double* lamb } auto d2f = f_pair.second; size_t nb_inst = f->_nb_instances; + id_inst = 0; for (unsigned inst = 0; inst_is_constraint) { - c_inst = c->get_id_inst(inst); - if (d2f->_is_matrix) { - for (unsigned i = 0; i < d2f->_dim[0]; i++) { - for (unsigned j = i; j < d2f->_dim[1]; j++) { + if (!*c->_all_lazy || !c->_lazy[inst]) { + idx_inc = false; + c_inst = c->get_id_inst(id_inst++); + if (d2f->_is_matrix) { + for (unsigned i = 0; i < d2f->_dim[0]; i++) { + for (unsigned j = i; j < d2f->_dim[1]; j++) { + res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; + idx++; + idx_inc = true; + } + } + } + else if(d2f->_is_vector){ + for (unsigned j = 0; j < d2f->_dim[0]; j++) { res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; idx++; idx_inc = true; } } - } - else if(d2f->_is_vector){ - for (unsigned j = 0; j < d2f->_dim[0]; j++) { + else { res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; - idx++; - idx_inc = true; } } else { - res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; + idx_inc = true; } } else { @@ -1485,57 +1649,64 @@ void Model::fill_in_hess(const double* x , double obj_factor, const double* lamb } auto d2f = f_pair.second; size_t nb_inst = f->_nb_instances; + unsigned id_inst = 0; for (unsigned inst = 0; inst_is_constraint) { - c_inst = c->get_id_inst(inst); - if (c->is_quadratic()) { - if (d2f->_is_matrix) { + if (!*c->_all_lazy || !c->_lazy[inst]) { + idx_inc = false; + c_inst = c->get_id_inst(id_inst++); + if (c->is_quadratic()) { + if (d2f->_is_matrix) { + for (unsigned i = 0; i < d2f->_dim[0]; i++) { + for (unsigned j = i; j < d2f->_dim[1]; j++) { + res[idx++] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; + idx_inc = true; + } + } + } + else if(d2f->_is_vector){ + for (unsigned j = 0; j < d2f->_dim[0]; j++) { + res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; + idx++; + idx_inc = true; + } + } + else { + res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; + } + } + else if (d2f->_is_matrix) { for (unsigned i = 0; i < d2f->_dim[0]; i++) { for (unsigned j = i; j < d2f->_dim[1]; j++) { - res[idx++] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; + hess = d2f->get_val(i,j); + _hess_vals[idx_in++] = hess; + res[idx++] += lambda[c->_id + c_inst] * hess; idx_inc = true; } } } else if(d2f->_is_vector){ for (unsigned j = 0; j < d2f->_dim[0]; j++) { - res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; + hess = d2f->get_val(j); + _hess_vals[idx_in++] = hess; + res[idx] += lambda[c->_id + c_inst] * hess; idx++; idx_inc = true; } } else { - res[idx] += lambda[c->_id + c_inst] * _hess_vals[idx_in++]; - } - } - else if (d2f->_is_matrix) { - for (unsigned i = 0; i < d2f->_dim[0]; i++) { - for (unsigned j = i; j < d2f->_dim[1]; j++) { - hess = d2f->get_val(i,j); - _hess_vals[idx_in++] = hess; - res[idx++] += lambda[c->_id + c_inst] * hess; - idx_inc = true; + if (d2f->is_number()) { + hess = d2f->_val->at(0); + } + else { + hess = d2f->_val->at(inst); } - } - } - else if(d2f->_is_vector){ - for (unsigned j = 0; j < d2f->_dim[0]; j++) { - hess = d2f->get_val(j); _hess_vals[idx_in++] = hess; res[idx] += lambda[c->_id + c_inst] * hess; - idx++; - idx_inc = true; } } else { - if (d2f->is_number()) { - hess = d2f->_val->at(0); - } - else { - hess = d2f->_val->at(inst); - } - _hess_vals[idx_in++] = hess; - res[idx] += lambda[c->_id + c_inst] * hess; + idx_inc = true; } } else { @@ -1679,26 +1850,28 @@ void Model::fill_in_maps() { void Model::fill_in_duals(double* lambda, double* z_L, double* z_U){ for (auto &cp: _cons) { - if(cp.second->_dual.size()==0){ - cp.second->_dual.resize(cp.second->_nb_instances); - } + unsigned idx = 0; for (unsigned inst = 0; inst < cp.second->_nb_instances; inst++) { - lambda[cp.second->_id + inst] = 0; + if (!*cp.second->_all_lazy || !cp.second->_lazy[inst]) { + lambda[cp.second->_id + idx++] = 0; + } } + idx = 0; for (unsigned inst = 0; inst < cp.second->_dual.size(); inst++) { - lambda[cp.second->_id + inst] = cp.second->_dual[inst]; + if (!*cp.second->_all_lazy || !cp.second->_lazy[inst]) { + auto index = cp.second->_id + idx++; + lambda[index] = cp.second->_dual[inst]; + } } } for (auto &vp: _vars) { - if(vp.second->_l_dual.size()==0 || vp.second->_u_dual.size()==0 ){ - vp.second->_l_dual.resize(vp.second->get_nb_instances()); - vp.second->_u_dual.resize(vp.second->get_nb_instances()); - } auto nb_inst = vp.second->get_nb_instances(); + auto vid = vp.second->get_id(); for (unsigned inst = 0; inst < nb_inst; inst++) { - z_L[vp.second->get_id() + vp.second->get_id_inst(inst)] = vp.second->_l_dual[inst]; - z_U[vp.second->get_id() + vp.second->get_id_inst(inst)] = vp.second->_u_dual[inst]; + auto id_inst = vp.second->get_id_inst(inst); + z_L[vid + id_inst] = vp.second->_l_dual[inst]; + z_U[vid + id_inst] = vp.second->_u_dual[inst]; // z_L[vp.second->get_id() + vp.second->get_id_inst(inst)] = 0; // z_U[vp.second->get_id() + vp.second->get_id_inst(inst)] = 0; } @@ -1777,28 +1950,37 @@ void Model::fill_in_cstr_bounds(double* g_l ,double* g_u) { switch (c->get_type()) { case eq:{ auto nb_ins = c->_nb_instances; + unsigned idx= 0; for (int i = 0; i< nb_ins; i++){ - cid = c->_id+i; - g_l[cid] = c->_rhs; - g_u[cid] = c->_rhs; + if (!*c->_all_lazy || !c->_lazy[i]) { + cid = c->_id+idx++; + g_l[cid] = c->_rhs; + g_u[cid] = c->_rhs; + } } break; } case leq:{ auto nb_ins = c->_nb_instances; + unsigned idx= 0; for (int i = 0; i< nb_ins; i++){ - cid = c->_id+i; - g_l[cid] = numeric_limits::lowest(); - g_u[cid] = c->_rhs; + if (!*c->_all_lazy || !c->_lazy[i]) { + cid = c->_id+idx++; + g_l[cid] = numeric_limits::lowest(); + g_u[cid] = c->_rhs; + } } break; } case geq:{ auto nb_ins = c->_nb_instances; + unsigned idx= 0; for (int i = 0; i< nb_ins; i++){ - cid = c->_id+i; - g_l[cid] = c->_rhs; - g_u[cid] = numeric_limits::max(); + if (!*c->_all_lazy || !c->_lazy[i]) { + cid = c->_id+idx++; + g_l[cid] = c->_rhs; + g_u[cid] = numeric_limits::max(); + } } break; } diff --git a/src/solver.cpp b/src/solver.cpp index 001a8ea9b..36c68fec1 100755 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -132,10 +132,17 @@ int solver::run(int print_level, bool relax, double tol, const string& lin_solve iapp->Options()->SetNumericValue("tol", tol); iapp->Options()->SetIntegerValue("print_level", print_level); + /** Bonmin options */ +// iapp->Options()->SetStringValue("mu_strategy", "adaptive"); +// iapp->Options()->SetStringValue("mu_oracle", "probing"); +// iapp->Options()->SetNumericValue("gamma_phi", 1e-8); +// iapp->Options()->SetNumericValue("gamma_theta", 1e-4); + + /** Hot start if already solved */ if (!_model->_first_run) { // if (false) { - DebugOff("Using Hot Start!\n"); + DebugOn("Using Hot Start!\n"); iapp->Options()->SetNumericValue("mu_init", mu_init); iapp->Options()->SetStringValue("warm_start_init_point", "yes"); } @@ -171,7 +178,6 @@ int solver::run(int print_level, bool relax, double tol, const string& lin_solve SmartPtr tmp = new IpoptProgram(_model); status = iapp->OptimizeTNLP(tmp); -// violated_constraints = _model->has_violated_constraints(tol); if (status == Solve_Succeeded) { // Retrieve some statistics about the solve @@ -292,7 +298,11 @@ int solver::run(int print_level, bool relax, double tol, const string& lin_solve bonminNotAvailable(); #endif } - violated_constraints = false; + violated_constraints = _model->has_violated_constraints(tol); + if (violated_constraints) { + _model->reindex(); + } +// violated_constraints = false; } return return_status; }