Skip to content

Commit

Permalink
Merge pull request #26 from Allinsights/SDP_cuts
Browse files Browse the repository at this point in the history
implementing lazy constraint generation
  • Loading branch information
hhijazi authored Feb 16, 2018
2 parents c0e1866 + 71492ae commit 6f4e034
Show file tree
Hide file tree
Showing 15 changed files with 596 additions and 354 deletions.
4 changes: 2 additions & 2 deletions examples/Classification/NeuralNet/NeuralNet_main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
8 changes: 4 additions & 4 deletions examples/MINLP/Power/ACOPF/ACOPF_main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
3 changes: 3 additions & 0 deletions examples/MINLP/Power/PowerNet.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand Down
15 changes: 8 additions & 7 deletions examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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();
}

Expand Down Expand Up @@ -257,39 +258,39 @@ 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");
LNC1 = (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);
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<Bag> bags;
int n3 = 0;
Expand Down
114 changes: 58 additions & 56 deletions examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -81,29 +81,29 @@ int main (int argc, char * argv[])

/** Variables */
/* power generation variables */
var<Real> Pg("Pg", grid->pg_min, grid->pg_max);
var<Real> Qg ("Qg", grid->qg_min, grid->qg_max);
SOCP.add_var(Pg.in(grid->gens));
SOCP.add_var(Qg.in(grid->gens));
var<Real> Pg("Pg", grid.pg_min, grid.pg_max);
var<Real> 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<Real> Pf_from("Pf_from", grid->S_max);
var<Real> Qf_from("Qf_from", grid->S_max);
var<Real> Pf_to("Pf_to", grid->S_max);
var<Real> 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<Real> Pf_from("Pf_from", grid.S_max);
var<Real> Qf_from("Qf_from", grid.S_max);
var<Real> Pf_to("Pf_to", grid.S_max);
var<Real> 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<Real> R_Wij("R_Wij", grid->wr_min, grid->wr_max);
var<Real> R_Wij("R_Wij", grid.wr_min, grid.wr_max);
/* Imaginary part of Wij = ViVj */
var<Real> Im_Wij("Im_Wij", grid->wi_min, grid->wi_max);
var<Real> Im_Wij("Im_Wij", grid.wi_min, grid.wi_max);
/* Magnitude of Wii = Vi^2 */
var<Real> Wii("Wii", grid->w_min, grid->w_max);
SOCP.add_var(Wii.in(grid->nodes));
var<Real> 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));

Expand All @@ -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 */
Expand All @@ -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 */
Expand Down Expand Up @@ -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<double>::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<double>::infinity()) + ", " + to_string(solve_time) + ", LocalOptimal, " + to_string(total_time);
DebugOn(out <<endl);
return 0;
}
14 changes: 14 additions & 0 deletions include/gravity/constraint.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,13 @@ namespace gravity {
ConstraintType _ctype = leq; /**< Constraint type: leq, geq or eq */
double _rhs = 0;
vector<double> _dual ; /**< Lagrange multipliers at a KKT point */
bool _all_active = true;
vector<bool> _active;
shared_ptr<bool> _all_lazy;
vector<bool> _lazy;
bool _all_satisfied = true;
vector<bool> _violated;


/** Constructor */
//@{
Expand Down Expand Up @@ -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;
Expand Down
2 changes: 0 additions & 2 deletions include/gravity/func.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool> _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;
Expand Down
9 changes: 7 additions & 2 deletions include/gravity/model.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 */

Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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<Constraint> add_constraint(const Constraint& c);

shared_ptr<func_> embed(shared_ptr<func_> 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. */
Expand All @@ -216,6 +220,7 @@ namespace gravity {
void set_objective(pair<func_*, ObjectiveType> 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();
Expand Down
2 changes: 1 addition & 1 deletion include/gravity/solver.h
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ namespace gravity {
public:
Model* _model = nullptr;
Program* _prog = nullptr;
SolverType _stype;
SolverType _stype;
double _tol = 1e-6; /*<< Solver tolerance. */

/** Constructor */
Expand Down
Loading

0 comments on commit 6f4e034

Please sign in to comment.