From f1476f5dcd56c759588d3bd037cd7480ae6ce258 Mon Sep 17 00:00:00 2001 From: kbestuzheva Date: Mon, 19 Feb 2018 13:31:23 +1100 Subject: [PATCH 1/3] More compact formulation of the nearest point problem --- examples/MINLP/Power/SDPOPF/Bag.cpp | 115 +++++++++++++++++++++++++--- examples/MINLP/Power/SDPOPF/Bag.h | 3 + include/gravity/param.h | 2 +- src/MosekProgram.cpp | 30 ++++++-- 4 files changed, 131 insertions(+), 19 deletions(-) diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index 9ca9da9de..351b7eaf2 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -5,6 +5,8 @@ #include #include "Bag.h" +#define FLAT + Bag::Bag(int id, const PowerNet& grid, vector nodes):_id(id),_grid((PowerNet*)&grid),_nodes(nodes) { Arc* aij = NULL; _wstarp.set_name("wstar"); @@ -45,6 +47,16 @@ Bag::Bag(int id, const PowerNet& grid, vector nodes):_id(id),_grid((Power _wmax.set_val(namewr,_grid->wr_max(namepair).eval()); _wmax.set_val(namewi,_grid->wi_max(namepair).eval()); + _Wmin.set_val(to_string(i)+","+to_string(j),_grid->wr_min(namepair).eval()); + _Wmin.set_val(to_string(i+_nodes.size())+","+to_string(j+_nodes.size()),_grid->wr_min(namepair).eval()); + _Wmin.set_val(to_string(i)+","+to_string(j+_nodes.size()),_grid->wi_min(namepair).eval()); + _Wmin.set_val(to_string(j)+","+to_string(i+_nodes.size()),_grid->wi_min(namepair).eval()); + + _Wmax.set_val(to_string(i)+","+to_string(j),_grid->wr_max(namepair).eval()); + _Wmax.set_val(to_string(i+_nodes.size())+","+to_string(j+_nodes.size()),_grid->wr_max(namepair).eval()); + _Wmax.set_val(to_string(i)+","+to_string(j+_nodes.size()),_grid->wi_max(namepair).eval()); + _Wmax.set_val(to_string(j)+","+to_string(i+_nodes.size()),_grid->wi_max(namepair).eval()); + // _wstarp.set_val(namewr,((Line*)aij)->wr); // if(!reversed) _wstarp.set_val(namewi,((Line*)aij)->wi); // else _wstarp.set_val(namewi,-((Line*)aij)->wi); @@ -61,6 +73,15 @@ Bag::Bag(int id, const PowerNet& grid, vector nodes):_id(id),_grid((Power _indices.push_back(index_(namew)); _wmin.set_val(namew,_grid->w_min(_nodes[i]->_name).eval()); _wmax.set_val(namew,_grid->w_max(_nodes[i]->_name).eval()); + + _Wmin.set_val(to_string(i)+","+to_string(i),_grid->w_min(_nodes[i]->_name).eval()); + _Wmin.set_val(to_string(i+_nodes.size())+","+to_string(i+_nodes.size()),_grid->w_min(_nodes[i]->_name).eval()); + _Wmin.set_val(to_string(i)+","+to_string(i+_nodes.size()),0); + + _Wmax.set_val(to_string(i)+","+to_string(i),_grid->w_max(_nodes[i]->_name).eval()); + _Wmax.set_val(to_string(i+_nodes.size())+","+to_string(i+_nodes.size()),_grid->w_max(_nodes[i]->_name).eval()); + _Wmax.set_val(to_string(i)+","+to_string(i+_nodes.size()),0); + // _wstarp.set_val(namew,((Bus*)_nodes[i])->w); // cout << "\n" << namew << " = " << ((Bus*)_nodes[i])->w; @@ -108,8 +129,13 @@ param Bag::fill_wstar(){ _W_star.set_val(to_string(i)+","+to_string(j),aij->wr); _W_star.set_val(to_string(i+_nodes.size())+","+to_string(j+_nodes.size()),aij->wr); - _W_star.set_val(to_string(i)+","+to_string(j+_nodes.size()),aij->wi); - _W_star.set_val(to_string(j)+","+to_string(i+_nodes.size()),aij->wi); + if(!reversed) { + _W_star.set_val(to_string(i) + "," + to_string(j + _nodes.size()), aij->wi); + _W_star.set_val(to_string(j) + "," + to_string(i + _nodes.size()), -aij->wi); + }else{ + _W_star.set_val(to_string(i) + "," + to_string(j + _nodes.size()), -aij->wi); + _W_star.set_val(to_string(j) + "," + to_string(i + _nodes.size()), aij->wi); + } } } for(int i = 0; i < _nodes.size(); i++){ @@ -120,7 +146,7 @@ param Bag::fill_wstar(){ _W_star.set_val(to_string(i+_nodes.size())+","+to_string(i+_nodes.size()),((Bus*)_nodes[i])->w); _W_star.set_val(to_string(i)+","+to_string(i+_nodes.size()),0.0); } - return W_star;//todo: use this + return W_star; } bool Bag::add_lines(){ @@ -334,6 +360,16 @@ void Bag::update_PSD(){ _is_psd = is_PSD(); } +vector triang_indices(int n) { + vector res; + for(int i = 0; i < n; i++) { + for(int j = i; j < n; j++) { + res.push_back(index_(to_string(i)+","+to_string(j))); + } + } + return res; +} + param Bag::nfp(){ param what; fill_wstar(); @@ -346,16 +382,24 @@ param Bag::nfp(){ int n = _nodes.size(); DebugOff("\nn = " << n); + _wstarp.print(true); + +#ifndef FLAT sdpvar W("W"); NPP.add_var(W^(2*n)); - -// var W("W"); -// W._psd = true; -// NPP.add_var(W ^ (n*(2*n+1))); +#else + var W("W"); + W._psd = true; + NPP.add_var(W.in(R(2*n,2*n))); // note: number of instances in apper triangle is n*(2*n+1) +#endif var z("z"); z.in_q_cone(); +#ifndef FLAT auto Rn = R(n*n+1); +#else + auto Rn = R(2*n*2*n+1); +#endif NPP.add_var(z.in(Rn)); var w("w", _wmin, _wmax); @@ -377,14 +421,40 @@ param Bag::nfp(){ /* w*-w = z */ Constraint svec("svec"); + +#ifndef FLAT svec = _wstarp - w - z; NPP.add_constraint(svec.in(_indices)==0); - - //todo: flatten(W*-W) = z; - +#else + auto idxs = triang_indices(2*n); + svec = _W_star - W - z; + NPP.add_constraint(svec.in(idxs)==0); + svec.print_expanded(); + +// vector zero_idxs; +// for(int i = 0; i < n; i++) { +// zero_idxs.push_back(index_(to_string(i)+","+to_string(i+n))); +// } +// Constraint zeros("zeros"); +// zeros = W.in(zero_idxs); +// NPP.add_constraint(zeros==0); +// zeros.print_expanded(); + + Constraint W_lb("W_lb"); + W_lb = _Wmin - W; + NPP.add_constraint(W_lb.in(idxs) <= 0); + W_lb.print_expanded(); + + Constraint W_ub("W_ub"); + W_ub = _Wmax - W; + NPP.add_constraint(W_ub.in(idxs) >= 0); + W_ub.print_expanded(); + +#endif string namew, namewr, namewi; +#ifndef FLAT /* matrix structure */ for(int i = 0; i < n; i++){ for(int j = i; j < n; j++){ @@ -428,6 +498,7 @@ param Bag::nfp(){ } } } +#endif solver s(NPP,Mosek); s.run(0,0); @@ -435,7 +506,7 @@ param Bag::nfp(){ // z.print(); cout << "\n"; // W.print(true); cout << "\n"; // w.print(); cout << "\n"; - +#ifndef FLAT for(int i = 0; i < n; i++){ for(int j = i; j < n; j++){ if(i==j){ @@ -453,7 +524,27 @@ param Bag::nfp(){ } } } +#else + for(int i = 0; i < n; i++){ + for(int j = i; j < n; j++){ + if(i==j){ + namew = "w(" + _nodes[i]->_name + ")"; + what.set_val(namew,W(to_string(i)+","+to_string(j)).eval()); + } + else { + namewr = "wr(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; + namewi = "wi(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; + what.set_val(namewr, W(to_string(i)+","+to_string(j)).eval()); + if(_grid->get_directed_arc(_nodes[i]->_name,_nodes[j]->_name)!=nullptr) + what.set_val(namewi, W(to_string(i)+","+to_string(j+_nodes.size())).eval()); + else + what.set_val(namewi, -W(to_string(i)+","+to_string(j+_nodes.size())).eval()); + } + } + } +#endif what.set_name("w_hat"); -// what.print(true); + what.print(true); + exit(0); return what; } diff --git a/examples/MINLP/Power/SDPOPF/Bag.h b/examples/MINLP/Power/SDPOPF/Bag.h index b6c810708..928e76281 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.h +++ b/examples/MINLP/Power/SDPOPF/Bag.h @@ -25,6 +25,9 @@ class Bag{ param _wmin; param _wmax; + param _Wmin; + param _Wmax; + // gravity::node_pairs _bus_pairs; gravity::Model _model; diff --git a/include/gravity/param.h b/include/gravity/param.h index 868a29e8f..bc863908e 100755 --- a/include/gravity/param.h +++ b/include/gravity/param.h @@ -99,7 +99,7 @@ namespace gravity { if(it->second==idx) break; } string key = it->first; - cout << "\nkey found: " << it->first.at(0) << "," << it->first.at(2); + DebugOff("\nkey found: " << it->first.at(0) << "," << it->first.at(2)); pair res; res = make_pair(std::stoi(it->first.substr(0,1)),std::stoi(it->first.substr(2,1))); return res; diff --git a/src/MosekProgram.cpp b/src/MosekProgram.cpp index da9d1a2fc..4a3c00807 100644 --- a/src/MosekProgram.cpp +++ b/src/MosekProgram.cpp @@ -77,9 +77,21 @@ bool MosekProgram::solve(bool relax) { poly_set_val(j, val , _model->_vars[i]); } } + else if(!_model->_vars[i]->_is_matrix){ + for (auto j = 0; j < _model->_vars[i]->get_nb_instances(); j++) { + poly_set_val(j, (*sol)[j], _model->_vars[i]); + } + } else { - for (auto j = 0; j < _model->_vars[i]->get_nb_instances(); j++) - poly_set_val(j,(*sol)[j] , _model->_vars[i]); + int n = _model->_vars[i]->_dim[0]; + for(auto j1 = 0; j1 < n; j1++) { + for(auto j2 = j1; j2 < n; j2++) { + string key = to_string(j1)+","+to_string(j2); + ((param*)_model->_vars[i])->set_val(key, (*sol)[j1*n+j2]); + } + } + + } } return 0; @@ -109,7 +121,8 @@ void MosekProgram::fill_in_mosek_vars() { if(!real_var->_in_q_cone && !real_var->_psd) _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), real_var->get_nb_instances(), fusion::Domain::inRange(lb,ub))); else if(real_var->_in_q_cone) _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inQCone(real_var->get_nb_instances()))); else { - _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inPSDCone((1 + sqrt(1+8*real_var->get_nb_instances()))/2))); +// _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inPSDCone((1 + sqrt(1+8*real_var->get_nb_instances()))/2))); + _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inPSDCone(real_var->_dim[0]))); } } else { @@ -154,7 +167,8 @@ void MosekProgram::fill_in_mosek_vars() { if(!real_var->_in_q_cone && !real_var->_psd) _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), real_var->get_nb_instances(), fusion::Domain::inRange(lb,ub))); else if(real_var->_in_q_cone) _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inQCone(real_var->get_nb_instances()))); else { - _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inPSDCone((1 + sqrt(1+8*real_var->get_nb_instances()))/2.))); +// _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inPSDCone((1 + sqrt(1+8*real_var->get_nb_instances()))/2.))); + _mosek_vars.push_back(_mosek_model->variable(real_var->get_name(), fusion::Domain::inPSDCone(real_var->_dim[0]))); } } else { @@ -294,7 +308,7 @@ void MosekProgram::create_mosek_constraints() { } else { c_idx_inst = inst; } - cout << "\nindex: " << pair.first << ", " << pair.second; + DebugOff("\nindex: " << pair.first << ", " << pair.second); auto lterm = fusion::Expr::mul(poly_eval(it1.second._coef, c_idx_inst), _mosek_vars[idx]->index(pair.first,pair.second)); @@ -357,7 +371,7 @@ void MosekProgram::create_mosek_constraints() { else { c_idx_inst = inst; } - auto lterm = fusion::Expr::constTerm(poly_eval(it1.second._coef, c_idx_inst) * poly_eval(it1.second._p, idx_inst)); + auto lterm = fusion::Expr::constTerm(poly_eval(it1.second._coef, c_idx_inst) * poly_eval(it1.second._p, i)); // was idx_inst if (!it1.second._sign) { lterm = fusion::Expr::mul(-1, lterm); } @@ -365,13 +379,17 @@ void MosekProgram::create_mosek_constraints() { } } } + DebugOff("\nconstr in Mosek: " << expr->toString()); if(c->get_type()==geq) { + DebugOff("\n >= " << c->get_rhs()); _mosek_model->constraint(c->get_name()+to_string(i), expr, fusion::Domain::greaterThan(c->get_rhs())); } else if(c->get_type()==leq) { + DebugOff("\n <= " << c->get_rhs()); _mosek_model->constraint(c->get_name()+to_string(i), expr, fusion::Domain::lessThan(c->get_rhs())); } else if(c->get_type()==eq) { + DebugOff("\n = " << c->get_rhs()); _mosek_model->constraint(c->get_name()+to_string(i), expr, fusion::Domain::equalsTo(c->get_rhs())); } inst++; From cd1584764b1f98d0a3674b9de44ed29a6ef5c3a5 Mon Sep 17 00:00:00 2001 From: kbestuzheva Date: Mon, 19 Feb 2018 15:58:10 +1100 Subject: [PATCH 2/3] Minor changes in Bag --- examples/MINLP/Power/SDPOPF/Bag.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index 351b7eaf2..e19b01fba 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -382,7 +382,7 @@ param Bag::nfp(){ int n = _nodes.size(); DebugOff("\nn = " << n); - _wstarp.print(true); +// _wstarp.print(true); #ifndef FLAT sdpvar W("W"); @@ -429,7 +429,7 @@ param Bag::nfp(){ auto idxs = triang_indices(2*n); svec = _W_star - W - z; NPP.add_constraint(svec.in(idxs)==0); - svec.print_expanded(); +// svec.print_expanded(); // vector zero_idxs; // for(int i = 0; i < n; i++) { @@ -443,12 +443,12 @@ param Bag::nfp(){ Constraint W_lb("W_lb"); W_lb = _Wmin - W; NPP.add_constraint(W_lb.in(idxs) <= 0); - W_lb.print_expanded(); +// W_lb.print_expanded(); Constraint W_ub("W_ub"); W_ub = _Wmax - W; NPP.add_constraint(W_ub.in(idxs) >= 0); - W_ub.print_expanded(); +// W_ub.print_expanded(); #endif @@ -544,7 +544,7 @@ param Bag::nfp(){ } #endif what.set_name("w_hat"); - what.print(true); - exit(0); +// what.print(true); +// exit(0); return what; } From 2520b91f651cda7916fa42affe52b1bb66940ea4 Mon Sep 17 00:00:00 2001 From: Hassan Date: Mon, 19 Feb 2018 08:03:49 +0100 Subject: [PATCH 3/3] working on cholesky sdp and fixing bug in pl dim --- examples/MILP/StableSet/Stable_Set_main.cpp | 135 +------ examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 377 +++++++++++++++++++- include/gravity/param.h | 6 +- include/gravity/var.h | 8 +- src/Arc.cpp | 1 + 5 files changed, 388 insertions(+), 139 deletions(-) diff --git a/examples/MILP/StableSet/Stable_Set_main.cpp b/examples/MILP/StableSet/Stable_Set_main.cpp index 48f3c40ef..69d4b4b77 100755 --- a/examples/MILP/StableSet/Stable_Set_main.cpp +++ b/examples/MILP/StableSet/Stable_Set_main.cpp @@ -1,35 +1,27 @@ - // +// // Stable_Set.cpp // Gravity // -// Created by Hijazi, Hassan (Data61, Canberra City) on 6/12/17. +// Created by Hijazi, Hassan on 6/12/17. // // #include #include #include -#include #include #include -#include -#include -#include -#include #include +#include using namespace std; using namespace gravity; - int main (int argc, const char * argv[]) { // Start Timers std::cout << "WELCOME, THIS IS AN IMPLEMENTATION OF THE STABLE SET PROBLEM AND SOME OF ITS RELAXATIONS IN GRAVITY\n"; - unsigned i, j; - Net graph; -// const char* fname = "../../data_sets/stable_set/toy.txt"; const char* fname = "../../data_sets/stable_set/p.3n150.txt"; graph.read_adjacency_matrix(fname); @@ -53,129 +45,10 @@ int main (int argc, const char * argv[]) c <= 1; c.print(); model.add_constraint(c.in(graph.arcs)); + /** Solver **/ SolverType stype = gurobi; solver s(model,stype); - double wall0 = get_wall_time(); - double cpu0 = get_cpu_time(); - cout << "Running the IP model\n"; s.run(); - double wall1 = get_wall_time(); - double cpu1 = get_cpu_time(); - cout << "Done running the IP model\n"; - cout << "Wall clock computing time = " << wall1 - wall0 << "\n"; - cout << "CPU computing time = " << cpu1 - cpu0 << "\n"; - - /* Schriver's SDP relaxation for the stable set problem */ - Model SDP; - /* Variable declaration */ - var Xii("Xii", 0, 1); - var Xij("Xij", 0, 1); - SDP.add_var(Xii.in(R(n))); /*< Diagonal entries of the matrix */ - SDP.add_var(Xij.in(R(n*(n-1)/2))); /*< Lower left triangular part of the matrix excluding the diagonal*/ - - /* Constraints declaration */ - ordered_pairs indices(1, n); - Constraint SOCP("SOCP"); - SOCP = power(Xij, 2) - Xii.from()*Xii.to() ; - SDP.add_constraint(SOCP.in(indices._keys) <= 0); - Constraint diag("diag"); - diag = sum(Xii); - SDP.add_constraint(diag == 1); // diagonal sum is 1 - - Constraint zeros("zeros"); - zeros = Xij.in(graph.arcs); - SDP.add_constraint(zeros == 0); // zero elements - - auto Xij_ = Xij.pairs_in(graph._bags, 3); - auto Xii_ = Xii.in(graph._bags, 3); - Constraint SDP3("SDP_3D"); - SDP3 = -2*Xij_[0]*Xij_[1]*Xij_[2]; - SDP3 -= Xii_[0]*Xii_[1]*Xii_[2]; - SDP3 += power(Xij_[0],2)*Xii_[2]; - SDP3 += power(Xij_[2],2)*Xii_[1]; - SDP3 += power(Xij_[1],2)*Xii_[0]; -// SDP3.print(); -// SDP.add_constraint(SDP3); -// set> ids; -// for (i = 0; i < complement_graph._bags.size(); i++){ -//// for (i = 0; i < 1; i++){ -// auto bag = complement_graph._bags.at(i); -// if (bag.size() < 3) { -// continue; -// } -//// for (int j = 0; j < bag.size()-2; j++){ -// for (int j = 0; j < 1; j++){ -// i1 = bag[j]->ID; -// i2 = bag[j+1]->ID; -// i3 = bag[j+2]->ID; -// assert(i2>i1 && i3>i2); -// if(ids.count(make_tuple(i1, i2, i3))==0){ -// ids.insert(make_tuple(i1, i2, i3)); -// } -// else { -// continue; -// } -// Constraint SDP3("SDP3("+to_string(i1)+","+to_string(i2)+","+to_string(i3)+")"); -// SDP3 = -2*Xij(i1,i2)*Xij(i2,i3)*Xij(i1,i3); -// SDP3 -= Xii(i1)*Xii(i2)*Xii(i3); -// SDP3 += power(Xij(i1,i2),2)*Xii(i3); -// SDP3 += power(Xij(i1,i3),2)*Xii(i2); -// SDP3 += power(Xij(i2,i3),2)*Xii(i1); -// SDP3.print(); -// SDP.add_constraint(SDP3); -// } -// } - - /* Objective declaration */ - SDP.max(2*sum(Xij) + sum(Xii)); - - -// solver s1(SDP,ipopt); - solver s1(SDP,gurobi); - - wall0 = get_wall_time(); - cpu0 = get_cpu_time(); - cout << "Running the SDP relaxation\n"; - s1.run(); - wall1 = get_wall_time(); - cpu1 = get_cpu_time(); - cout << "Done running the SDP relaxation\n"; - cout << "\nWall clock computing time = " << wall1 - wall0 << "\n"; - cout << "CPU computing time = " << cpu1 - cpu0 << "\n"; - return 0; - - /* Outer-Approximation of SOCP and 3d-SDP cuts in Shriver's SDP relaxation for the stable set problem */ - Model OA; - OA.add_var(Xii); - OA.add_var(Xij); - for (auto &cs_p: SDP._cons) { - if (!cs_p.second->is_linear() && cs_p.second->is_active()){ //Active nonlinear constraint - DebugOff("Active constraint:" << cs_p.second->to_str() << endl); - Constraint oa("OA_"+cs_p.second->get_name()); - oa = cs_p.second->get_outer_app(); - // oa.print(); - OA.add_constraint(oa); - } - } -// OA.add_constraint(diag=1); - for(auto a: graph.arcs){ - i = (a->_src)->_id; - j = (a->_dest)->_id; - Constraint zeros("zeros("+to_string(i)+","+to_string(j)+")"); - zeros = Xij(i,j); - OA.add_constraint(zeros==0); - } - OA.max(2*sum(Xij) + sum(Xii)); - solver s2(OA,cplex); - wall0 = get_wall_time(); - cpu0 = get_cpu_time(); - cout << "Running the OA-SDP relaxation\n"; - s2.run(); - wall1 = get_wall_time(); - cpu1 = get_cpu_time(); - cout << "Done running the OA-SDP relaxation\n"; - cout << "\nWall clock computing time = " << wall1 - wall0 << "\n"; - cout << "CPU computing time = " << cpu1 - cpu0 << "\n"; return 0; }; diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index 55ac8acf1..bedc600af 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -81,6 +81,381 @@ vector> signs(Net& net, const std::vector>& bags) { return res; } +int main_new (int argc, char * argv[]) { + int output = 0; + bool relax = false; + string solver_str = "ipopt"; + SolverType solv_type = ipopt; + double tol = 1e-6; + string mehrotra = "no"; + string fname = "../data_sets/Power/nesta_case5_pjm.m"; + + // create a OptionParser with options + op::OptionParser opt; + opt.add_option("h", "help", + "shows option help"); // no default value means boolean options, which default value is false + opt.add_option("f", "file", "Input file name", fname); + opt.add_option("s", "solver", "Solvers: ipopt/cplex/gurobi, default = ipopt", solver_str); + /* parse the options and verify that all went well. If not, errors and help will be shown */ + bool correct_parsing = opt.parse_options(argc, argv); + + if (!correct_parsing) { + return EXIT_FAILURE; + } + + fname = opt["f"]; + bool has_help = op::str2bool(opt["h"]); + if (has_help) { + opt.show_help(); + exit(0); + } + solver_str = opt["s"]; + if (solver_str.compare("gurobi")==0) { + solv_type = gurobi; + } + else if(solver_str.compare("cplex")==0) { + solv_type = cplex; + } + + double total_time_start = get_wall_time(); + PowerNet grid; + grid.readgrid(fname.c_str()); + + grid.get_tree_decomp_bags(); + grid.update_net(); + + /* Grid Parameters */ + auto bus_pairs = grid.get_bus_pairs(); + auto bus_pairs_chord = grid.get_bus_pairs_chord(); + auto nb_lines = grid.get_nb_active_arcs(); + auto nb_buses = grid.get_nb_active_nodes(); + DebugOff("nb gens = " << nb_gen << endl); + DebugOff("nb lines = " << nb_lines << endl); + DebugOff("nb buses = " << nb_buses << endl); + DebugOff("nb bus_pairs = " << nb_bus_pairs_chord << endl); + + + double upper_bound = grid.solve_acopf(); + + + /** Build model */ + Model SDP("SDP Model"); + + /** Variables */ + /* power generation variables */ + var Pg("Pg", grid.pg_min, grid.pg_max); + var Qg ("Qg", grid.qg_min, grid.qg_max); + SDP.add_var(Pg.in(grid.gens)); + SDP.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); + SDP.add_var(Pf_from.in(grid.arcs)); + SDP.add_var(Qf_from.in(grid.arcs)); + SDP.add_var(Pf_to.in(grid.arcs)); + SDP.add_var(Qf_to.in(grid.arcs)); + + /* Real part of Wij = ViVj */ + 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); + /* Magnitude of Wii = Vi^2 */ + var Wii("Wii", grid.w_min, grid.w_max); + SDP.add_var(Wii.in(grid.nodes)); + SDP.add_var(R_Wij.in(bus_pairs_chord)); + SDP.add_var(Im_Wij.in(bus_pairs_chord)); + + /* Initialize variables */ + R_Wij.initialize_all(1.0); + Wii.initialize_all(1.001); + + /** Objective */ + auto obj = product(grid.c1, Pg) + product(grid.c2, power(Pg,2)) + sum(grid.c0); + SDP.min(obj.in(grid.gens)); + + /** Constraints */ + + DebugOn("Adding Cholesky 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); + auto I_sgn = signs(grid, grid._bags); + Constraint Wii_LL("Wii_LL"); + Wii_LL = Wii_[0]; + SDP.add(Wii_LL = 0); + + + /* Second-order cone constraints */ + Constraint SOC("SOC"); + SOC = power(R_Wij, 2) + power(Im_Wij, 2) - Wii.from()*Wii.to(); + SDP.add_constraint(SOC.in(bus_pairs) <= 0); + + /* 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; + SDP.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; + SDP.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()); + SDP.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()); + SDP.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()); + SDP.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()); + SDP.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; + 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_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_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_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_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_lazy(LNC2.in(bus_pairs) >= 0); + + vector bags; + int n3 = 0; + int bagid = 0; + for(auto& b: grid._bags){ + bags.push_back(Bag(bagid,grid,b)); + if(b.size()==3) n3++; + bagid++; + } + + DebugOn("\nNum of 3d bags = " << n3); + + + /* Solver selection */ + solver SDPOPF(SDP,solv_type); + double solver_time_start = get_wall_time(); + SDPOPF.run(output = 5, relax = false); + // SDPOPF.run(output = 0, relax = false, tol = 1e-6, "ma27", mehrotra = "no"); + + for(auto& arc: grid.arcs){ + ((Line*)arc)->wr = R_Wij(arc->_src->_name+","+arc->_dest->_name).eval(); + ((Line*)arc)->wi = Im_Wij(arc->_src->_name+","+arc->_dest->_name).eval(); + if(grid.add_3d_nlin && !arc->_free) arc->_active = true; + } + for(auto& node: grid.nodes) ((Bus*)node)->w = Wii(node->_name).eval(); + + + // param w_hat; + string namew, namewr, namewi; + int numcuts = 0; + node_pairs bus_pairs_sdp; + // double prev_opt = 0; + // double fp_tol = 0.0001; + DebugOff("\nNumber of bags = " << bags.size()); + + vector> R_star, I_star, W_star; + vector> R_diff, I_diff, W_diff; + vector> R_hat, I_hat, W_hat; + + // int unchanged = 0; + int iter = 0, hdim_cuts = 0, cuts_added = 1; + + // while(unchanged < 3) { + unsigned nr_threads = 10; + double gap = 100*(upper_bound - SDP._obj_val)/upper_bound; + while(cuts_added > 0 && gap > 1) { + DebugOn("Current Gap = " << to_string(gap) << "%."<> w_hat_vec; + + iter++; + // while((SDP._obj_val - prev_opt)/SDP._obj_val > fp_tol) { + // prev_opt = SDP._obj_val; + + /* Update _is_psd */ + vector violated_bags; + for(auto& b: bags){ + b.update_PSD(); + if (!b._is_psd) { + violated_bags.push_back(b); + } + } + auto nb_bags = violated_bags.size(); + w_hat_vec.resize(nb_bags); + vector threads; + /* Split subproblems into nr_threads parts */ + vector limits = bounds(nr_threads, nb_bags); + /* Launch all threads in parallel */ + for (int i = 0; i < nr_threads; ++i) { + threads.push_back(thread(nearest_point, limits[i], limits[i+1], ref(w_hat_vec), ref(violated_bags))); + } + /* Join the threads with the main thread */ + for(auto &t : threads){ + t.join(); + } + for (unsigned index = 0; index < nb_bags; index++) { + auto b = violated_bags[index]; + DebugOff("\nBag number " << b._id); + DebugOff("\nNodes: n1 = " << b._nodes[0]->_name << ", n2 = " << b._nodes[1]->_name << ", n3 = " << b._nodes[2]->_name); + auto w_hat(move(w_hat_vec[index])); + R_star.resize(numcuts+1); + R_star[numcuts] = param<>("R_star"+to_string(numcuts)); + I_star.resize(numcuts+1); + I_star[numcuts] = param<>("I_star"+to_string(numcuts)); + W_star.resize(numcuts+1); + W_star[numcuts] = param<>("W_star"+to_string(numcuts)); + R_diff.resize(numcuts+1); + R_diff[numcuts] = param<>("R_diff"+to_string(numcuts)); + I_diff.resize(numcuts+1); + I_diff[numcuts] = param<>("I_diff"+to_string(numcuts)); + W_diff.resize(numcuts+1); + W_diff[numcuts] = param<>("W_diff"+to_string(numcuts)); + R_hat.resize(numcuts+1); + R_hat[numcuts] = param<>("R_hat"+to_string(numcuts)); + I_hat.resize(numcuts+1); + I_hat[numcuts] = param<>("I_hat"+to_string(numcuts)); + W_hat.resize(numcuts+1); + W_hat[numcuts] = param<>("W_hat"+to_string(numcuts)); + + node_pairs b_pairs("node_pairs"); + Constraint sdpcut("sdpcut_" + to_string(numcuts)); + Node *ni; + Arc *aij; + for (int i = 0; i < b._nodes.size(); i++) { + for (int j = i; j < b._nodes.size(); j++) { + if (i == j) { + namew = "w(" + b._nodes[i]->_name + ")"; + ni = b._nodes[i]; + W_diff[numcuts].set_val(ni->_name,((Bus *) ni)->w - w_hat(namew).eval()); + W_hat[numcuts].set_val(ni->_name,w_hat(namew).eval()); + W_star[numcuts].set_val(ni->_name,((Bus *) ni)->w); + } + else { + aij = grid.get_arc(b._nodes[i]->_name, b._nodes[j]->_name); + if(aij->_imaginary && !aij->_active) { + bus_pairs_sdp._keys.push_back(new index_pair(index_(aij->_src->_name), index_(aij->_dest->_name))); + aij->_active = true; + } + b_pairs._keys.push_back(new index_pair(index_(aij->_src->_name), index_(aij->_dest->_name))); + + namewr = "wr(" + b._nodes[i]->_name + "," + b._nodes[j]->_name + ")"; + namewi = "wi(" + b._nodes[i]->_name + "," + b._nodes[j]->_name + ")"; + + R_diff[numcuts].set_val(aij->_src->_name + "," + aij->_dest->_name,((Line *) aij)->wr - w_hat(namewr).eval()); + I_diff[numcuts].set_val(aij->_src->_name + "," + aij->_dest->_name,((Line *) aij)->wi - w_hat(namewi).eval()); + R_hat[numcuts].set_val(aij->_src->_name + "," + aij->_dest->_name,w_hat(namewr).eval()); + I_hat[numcuts].set_val(aij->_src->_name + "," + aij->_dest->_name,w_hat(namewi).eval()); + R_star[numcuts].set_val(aij->_src->_name + "," + aij->_dest->_name,((Line *) aij)->wr); + I_star[numcuts].set_val(aij->_src->_name + "," + aij->_dest->_name,((Line *) aij)->wi); + } + } + } + if(b._nodes.size()>3) hdim_cuts++; + + Constraint lin("lin"+to_string(numcuts)); + lin = product(R_diff[numcuts].in(b_pairs),(R_Wij.in(b_pairs) - R_hat[numcuts].in(b_pairs))); + lin += product(I_diff[numcuts].in(b_pairs),(Im_Wij.in(b_pairs) - I_hat[numcuts].in(b_pairs))); + lin += product(W_diff[numcuts].in(b._nodes),(Wii.in(b._nodes) - W_hat[numcuts].in(b._nodes))); + // lin.print_expanded(); + SDP.add_constraint(lin <= 0); + + cuts_added++; numcuts++; + } + + if(!bus_pairs_sdp._keys.empty()) { + SDP.add_indices("SOC",bus_pairs_sdp); + bus_pairs_sdp.clear(); + } + + + for(auto& a: grid.arcs){ + if(a->_imaginary && !a->_active) a->_free = true; + } + + SDPOPF.run(output = 5, relax = false, tol = 1e-6, "ma27", mehrotra = "no"); + gap = 100*(upper_bound - SDP._obj_val)/upper_bound; + DebugOff("\nPrev = " << prev_opt << ", cur = " << SDP._obj_val); + DebugOff("\n(opt - prev)/opt = " << (SDP._obj_val - prev_opt)/SDP._obj_val << endl); + + /* Update values for w_star */ + for(auto& arc: grid.arcs){ + ((Line*)arc)->wr = R_Wij(arc->_src->_name+","+arc->_dest->_name).eval(); + ((Line*)arc)->wi = Im_Wij(arc->_src->_name+","+arc->_dest->_name).eval(); + } + for(auto& node: grid.nodes) ((Bus*)node)->w = Wii(node->_name).eval(); + + + cout << "\nNum of iterations = " << iter << ", number of cuts = " << numcuts << endl; + // if((SDP._obj_val - prev_opt)/SDP._obj_val < fp_tol) unchanged++; + } + string lin = "lin"; + int num_act_cuts = 0; + for (auto &cp: SDP._cons) { + if(cp.second->get_name().find(lin) != std::string::npos) { + for (unsigned inst = 0; inst < cp.second->_nb_instances; inst++) { + cp.second->eval(inst); + if(cp.second->is_active(inst)) num_act_cuts++; + } + } + } + + cout << "\nEnd: num of iterations = " << iter << ", number of cuts = " << numcuts << ", number of higher dim cuts = " << hdim_cuts; + cout << "\nNumber of active cuts = " << num_act_cuts; + + double solver_time_end = get_wall_time(); + double total_time_end = get_wall_time(); + auto solve_time = solver_time_end - solver_time_start; + auto total_time = total_time_end - total_time_start; + string out = "\nDATA_OPF, " + grid._name + ", " + to_string(nb_buses) + ", " + to_string(nb_lines) +", " + to_string(SDP._obj_val) + ", " + to_string(-numeric_limits::infinity()) + ", " + to_string(solve_time) + ", LocalOptimal, " + to_string(total_time); + DebugOn(out <_name; // in_cycle = false; // parallel = false; // connect();