From 6c04e270a63085cd74f9ee22b782a2301d6f1bbf Mon Sep 17 00:00:00 2001 From: kbestuzheva Date: Wed, 21 Feb 2018 16:53:24 +1100 Subject: [PATCH 01/10] Using eigenvalues to compute W_hat, in progress --- examples/MINLP/Power/SDPOPF/Bag.cpp | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index e19b01fba..45f4d41fb 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -338,8 +338,9 @@ bool Bag::is_PSD(){ } } // A.print(); - arma::cx_mat R; - arma::vec v = arma::eig_sym(A); + arma::cx_mat P; + arma::vec v;// = arma::eig_sym(A); + arma::eig_sym(v,P,A); DebugOff("\n"); double min_eig = 0, max_eig = -1; for(auto eig: v) { @@ -351,7 +352,26 @@ bool Bag::is_PSD(){ return true; } else { - DebugOff("\nBag is not PSD"); + double pos_tol = -0.00001; + DebugOn("\nBag is not PSD"); +// for(int i = 0; i < n; i++) { +// if(v[i] < pos_tol) v[i] = 0; +// } + arma::cx_mat D(n,n); + for(int i = 0; i < n; i++) { + for(int j = 0; j < n; j++) { + if(i==j) D[i,j] = v[i]; + else D[i,j] = 0; + } + } + arma::cx_mat W_hat = P*D*P.t(); //todo: the values of W_hat don't make sense... + for(int i = 0; i < n; i++) { + for(int j = i; j < n; j++) { + cout << "\nW_hat(" << i << "," << j << ") = " << W_hat[i,j].real() << " + i" << W_hat[i,j].imag(); + } + cout << "\n"; + } + exit(0); return false; } } From 43a14617e2cce7ca9d7f9664782134663154fc61 Mon Sep 17 00:00:00 2001 From: Hassan Date: Thu, 22 Feb 2018 00:15:28 +0100 Subject: [PATCH 02/10] minor touch-ups --- examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 2 +- examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp | 2 +- src/func.cpp | 3 +++ 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index bedc600af..48c7ffd30 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -186,7 +186,7 @@ int main_new (int argc, char * argv[]) { auto I_sgn = signs(grid, grid._bags); Constraint Wii_LL("Wii_LL"); Wii_LL = Wii_[0]; - SDP.add(Wii_LL = 0); + SDP.add(Wii_LL == 0); /* Second-order cone constraints */ diff --git a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp index b4979c2a7..41f3cb5d9 100644 --- a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp +++ b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp @@ -124,7 +124,7 @@ 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); + 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; diff --git a/src/func.cpp b/src/func.cpp index ced753537..18db56850 100755 --- a/src/func.cpp +++ b/src/func.cpp @@ -7156,6 +7156,9 @@ namespace gravity{ size_t n = 0; for (auto &p: *_vars) { if (p.second.first->_is_vector) { + if (p.second.first->_dim.size()>1) { + throw invalid_argument("get_nb_vars() should be called with an index, i.e. get_nb_vars(int)"); + } n += p.second.first->get_dim(); } else { From 0a6606d40256efc956a37a2972820be3e5916a9e Mon Sep 17 00:00:00 2001 From: Hassan Date: Thu, 22 Feb 2018 00:53:12 +0100 Subject: [PATCH 03/10] fixing armadillo issue --- examples/MINLP/Power/PowerNet.h | 2 +- examples/MINLP/Power/SDPOPF/Bag.cpp | 33 +++++++++++++++++------------ 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/examples/MINLP/Power/PowerNet.h b/examples/MINLP/Power/PowerNet.h index 092dcfa1b..a0c7449b5 100644 --- a/examples/MINLP/Power/PowerNet.h +++ b/examples/MINLP/Power/PowerNet.h @@ -26,7 +26,7 @@ typedef enum { ACPOL, ACRECT, QC, QC_SDP, OTS, DF, SOCP, SDP, DC, QC_OTS_L, QC_O class PowerNet: public Net { public: - bool add_3d_nlin = true; + bool add_3d_nlin = false; string ref_bus; double bMVA; /**< Base MVA */ double bV; /**< Base Voltage */ diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index 45f4d41fb..1a4528616 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -353,24 +353,29 @@ bool Bag::is_PSD(){ } else { double pos_tol = -0.00001; - DebugOn("\nBag is not PSD"); + DebugOn("\nBag is not PSD\n"); // for(int i = 0; i < n; i++) { // if(v[i] < pos_tol) v[i] = 0; // } - arma::cx_mat D(n,n); - for(int i = 0; i < n; i++) { - for(int j = 0; j < n; j++) { - if(i==j) D[i,j] = v[i]; - else D[i,j] = 0; - } - } + A.print(); + arma::mat D(n,n); + D.zeros(); + D.diag() = v; +// for(int i = 0; i < n; i++) { +// for(int j = 0; j < n; j++) { +// if(i==j) D[i,j] = v[i]; +// else D[i,j] = 0; +// } +// } arma::cx_mat W_hat = P*D*P.t(); //todo: the values of W_hat don't make sense... - for(int i = 0; i < n; i++) { - for(int j = i; j < n; j++) { - cout << "\nW_hat(" << i << "," << j << ") = " << W_hat[i,j].real() << " + i" << W_hat[i,j].imag(); - } - cout << "\n"; - } + cout << endl; + W_hat.print(); +// for(int i = 0; i < n; i++) { +// for(int j = i; j < n; j++) { +// cout << "\nW_hat(" << i << "," << j << ") = " << W_hat[i,j].real() << " + i" << W_hat[i,j].imag(); +// } +// cout << "\n"; +// } exit(0); return false; } From 6a6c3ae25a87cf9acaf38ff0f5b9825e410d883c Mon Sep 17 00:00:00 2001 From: kbestuzheva Date: Thu, 22 Feb 2018 17:03:49 +1100 Subject: [PATCH 04/10] Trying to use eigenvalues to compute nfp --- examples/MINLP/Power/SDPOPF/Bag.cpp | 147 ++++++++++++++++---- examples/MINLP/Power/SDPOPF/Bag.h | 1 + examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 31 ++--- 3 files changed, 135 insertions(+), 44 deletions(-) diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index 1a4528616..210054925 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -5,7 +5,7 @@ #include #include "Bag.h" -#define FLAT +//#define FLAT Bag::Bag(int id, const PowerNet& grid, vector nodes):_id(id),_grid((PowerNet*)&grid),_nodes(nodes) { Arc* aij = NULL; @@ -152,6 +152,16 @@ param Bag::fill_wstar(){ bool Bag::add_lines(){ // this is only called for 3d bags with 1 missing line // if(_nodes.size() != 3 || _all_lines) return; + + int free_lines = 0; + for (int i = 0; i < _nodes.size() - 1; i++) { + for (int j = i + 1; j < _nodes.size(); j++) { + Arc *aij = _grid->get_arc(_nodes[i]->_name, _nodes[j]->_name); + if (aij->_free) free_lines++; + } + } + if (free_lines !=1) return true; + Line *a12, *a13, *a23; Bus *n1, *n2, *n3; double tol = 0.00001; @@ -162,8 +172,6 @@ bool Bag::add_lines(){ a12 = (Line*)_grid->get_arc(n1,n2); a13 = (Line*)_grid->get_arc(n1,n3); a23 = (Line*)_grid->get_arc(n2,n3); -// if((a12->_free && a13->_free) || (a23->_free && a13->_free) || (a12->_free && a23->_free)) //more than 2 unassigned lines -// return; double wr12, wi12, wr13, wi13, wr23, wi23; double w1 = n1->w; double w2 = n2->w; double w3 = n3->w; @@ -353,30 +361,18 @@ bool Bag::is_PSD(){ } else { double pos_tol = -0.00001; - DebugOn("\nBag is not PSD\n"); -// for(int i = 0; i < n; i++) { -// if(v[i] < pos_tol) v[i] = 0; -// } - A.print(); + DebugOff("\nBag is not PSD\n"); + for(int i = 0; i < n; i++) { + if(v[i] < pos_tol) { + v[i] = 0; + P.col(i).zeros(); + } + } arma::mat D(n,n); D.zeros(); D.diag() = v; -// for(int i = 0; i < n; i++) { -// for(int j = 0; j < n; j++) { -// if(i==j) D[i,j] = v[i]; -// else D[i,j] = 0; -// } -// } - arma::cx_mat W_hat = P*D*P.t(); //todo: the values of W_hat don't make sense... - cout << endl; - W_hat.print(); -// for(int i = 0; i < n; i++) { -// for(int j = i; j < n; j++) { -// cout << "\nW_hat(" << i << "," << j << ") = " << W_hat[i,j].real() << " + i" << W_hat[i,j].imag(); -// } -// cout << "\n"; -// } - exit(0); + arma::cx_mat W_hat = P*D*P.t(); + return false; } } @@ -427,7 +423,7 @@ param Bag::nfp(){ #endif NPP.add_var(z.in(Rn)); - var w("w", _wmin, _wmax); + var w("w");//, _wmin, _wmax); NPP.add_var(w.in(_indices)); func_ obj; @@ -569,7 +565,104 @@ param Bag::nfp(){ } #endif what.set_name("w_hat"); -// what.print(true); -// exit(0); + what.print(true); + return what; +} + +param Bag::nfp1(){ + param what; + int n = _nodes.size(); + DebugOff("\nn = " << n); + + if(n == 2) { + DebugOff("Returning empty w_hat."); + return what; + } + + int free_lines = 0; + for (int i = 0; i < _nodes.size() - 1; i++) { + for (int j = i + 1; j < _nodes.size(); j++) { + Arc *aij = _grid->get_arc(_nodes[i]->_name, _nodes[j]->_name); + if (aij->_free) free_lines++; + } + } + + if (free_lines > 1) { + DebugOff("Returning empty w_hat."); + return what; + } + + //the bag has all lines or has > 3 nodes + + double tol = 0.0001; + arma::cx_mat A(n,n); + + for(int i = 0; i < n; i++){ + for(int j = 0; j <= i; j++) { + if(i==j) { + A(i,j) = arma::cx_double(((Bus*)_grid->get_node(_nodes[i]->_name))->w,0); + }else{ + double AijI; + if(_grid->has_directed_arc(_nodes[j],_nodes[i])) { + AijI = ((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wi; + } + else { + AijI = -((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wi; + } + A(i,j) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,AijI); + A(j,i) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,-AijI); + } + } + } + arma::cx_mat Eigval; arma::cx_mat W_hat; + arma::vec eigvec; + arma::eig_sym(eigvec,Eigval,A); + DebugOff("\n"); + double min_eig = 0, max_eig = -1; + for(auto eig: eigvec) { + if(eig < min_eig) min_eig = eig; + if(eig > max_eig) max_eig = eig; + } + + if(min_eig/max_eig > -tol) { + DebugOn("\nBag is PSD"); + return what; + } else { + DebugOn("\nBag is not PSD\n"); + for(int i = 0; i < n; i++) { + if(eigvec[i] < 0) { + eigvec[i] = 0; +// Eigval.col(i).zeros(); + } + } + arma::mat Eigvec(n,n); + Eigvec.zeros(); + Eigvec.diag() = eigvec; + cout << "\nEigval:" << endl; Eigval.print(); cout << "\nD:" << endl; Eigvec.print(); + W_hat = Eigval*Eigvec*Eigval.t(); + } + + string namew, namewr, namewi; + + 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_hat[i,i].real()); + } + else { + namewr = "wr(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; + namewi = "wi(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; + what.set_val(namewr, W_hat[i,j].real()); + if(_grid->get_directed_arc(_nodes[i]->_name,_nodes[j]->_name)!=nullptr) + what.set_val(namewi, W_hat[i,j].imag()); + else + what.set_val(namewi, -W_hat[i,j].imag()); + } + } + } + + what.set_name("w_hat"); + what.print(true); return what; } diff --git a/examples/MINLP/Power/SDPOPF/Bag.h b/examples/MINLP/Power/SDPOPF/Bag.h index 928e76281..a7e26359c 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.h +++ b/examples/MINLP/Power/SDPOPF/Bag.h @@ -43,6 +43,7 @@ class Bag{ /* find nearest feasible point */ param nfp(); + param nfp1(); // uses eigenvalues, disregards bounds bool add_lines(); diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index 48c7ffd30..53ff383d6 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -22,7 +22,8 @@ using namespace gravity; void nearest_point(int i, int j, vector>& w_hat, const vector& bags){ for (unsigned index = i; index 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; + unsigned nr_threads = 1; 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; + +// for(auto& b: bags) { +// if(b._nodes.size()==3) b.add_lines(); +// } /* Update _is_psd */ vector violated_bags; @@ -785,10 +782,10 @@ int main (int argc, char * argv[]) { 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()); @@ -804,12 +801,12 @@ int main (int argc, char * argv[]) { 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(); + 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(); @@ -819,7 +816,7 @@ int main (int argc, char * argv[]) { 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); From 37e2d31225b456e22db882da7ab358b2d8faf4e3 Mon Sep 17 00:00:00 2001 From: Hassan Date: Thu, 22 Feb 2018 10:28:27 +0100 Subject: [PATCH 05/10] Sign issue fixed, still getting different cuts, numerical issues? --- examples/MINLP/Power/SDPOPF/Bag.cpp | 31 +++++++++++++++++------------ 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index 210054925..baa9b3755 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -363,9 +363,9 @@ bool Bag::is_PSD(){ double pos_tol = -0.00001; DebugOff("\nBag is not PSD\n"); for(int i = 0; i < n; i++) { - if(v[i] < pos_tol) { + if(v[i] < 0) { v[i] = 0; - P.col(i).zeros(); +// P.col(i).zeros(); } } arma::mat D(n,n); @@ -602,22 +602,25 @@ param Bag::nfp1(){ if(i==j) { A(i,j) = arma::cx_double(((Bus*)_grid->get_node(_nodes[i]->_name))->w,0); }else{ - double AijI; + double AijI = ((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wi; if(_grid->has_directed_arc(_nodes[j],_nodes[i])) { - AijI = ((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wi; + A(i,j) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,AijI); + A(j,i) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,-AijI); } else { - AijI = -((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wi; + A(i,j) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,-AijI); + A(j,i) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,AijI); } - A(i,j) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,AijI); - A(j,i) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,-AijI); + } } } arma::cx_mat Eigval; arma::cx_mat W_hat; arma::vec eigvec; + DebugOn("x_star = "); + A.print(); arma::eig_sym(eigvec,Eigval,A); - DebugOff("\n"); + DebugOn("\n"); double min_eig = 0, max_eig = -1; for(auto eig: eigvec) { if(eig < min_eig) min_eig = eig; @@ -638,8 +641,10 @@ param Bag::nfp1(){ arma::mat Eigvec(n,n); Eigvec.zeros(); Eigvec.diag() = eigvec; - cout << "\nEigval:" << endl; Eigval.print(); cout << "\nD:" << endl; Eigvec.print(); + cout << "\nEigvec:" << endl; Eigval.print(); cout << "\nD:" << endl; Eigvec.print(); W_hat = Eigval*Eigvec*Eigval.t(); + DebugOn("x_hat = "); + W_hat.print(); } string namew, namewr, namewi; @@ -648,16 +653,16 @@ param Bag::nfp1(){ for(int j = i; j < n; j++){ if(i==j){ namew = "w(" + _nodes[i]->_name + ")"; - what.set_val(namew,W_hat[i,i].real()); + what.set_val(namew,W_hat(i,i).real()); } else { namewr = "wr(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; namewi = "wi(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; - what.set_val(namewr, W_hat[i,j].real()); + what.set_val(namewr, W_hat(i,j).real()); if(_grid->get_directed_arc(_nodes[i]->_name,_nodes[j]->_name)!=nullptr) - what.set_val(namewi, W_hat[i,j].imag()); + what.set_val(namewi, -1*W_hat(i,j).imag()); else - what.set_val(namewi, -W_hat[i,j].imag()); + what.set_val(namewi, W_hat(i,j).imag()); } } } From ff400e748bdb3a302a7910a52dca445c80e1e715 Mon Sep 17 00:00:00 2001 From: Hassan Date: Thu, 22 Feb 2018 23:39:18 +0100 Subject: [PATCH 06/10] Switching back to nfp, nfp1 seems to be numerically unreliable, needs investigation --- examples/MINLP/Power/PowerNet.h | 2 +- examples/MINLP/Power/SDPOPF/Bag.cpp | 69 ++++++++++++--------- examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 54 +++++++++------- examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp | 36 +++++------ 4 files changed, 89 insertions(+), 72 deletions(-) diff --git a/examples/MINLP/Power/PowerNet.h b/examples/MINLP/Power/PowerNet.h index a0c7449b5..092dcfa1b 100644 --- a/examples/MINLP/Power/PowerNet.h +++ b/examples/MINLP/Power/PowerNet.h @@ -26,7 +26,7 @@ typedef enum { ACPOL, ACRECT, QC, QC_SDP, OTS, DF, SOCP, SDP, DC, QC_OTS_L, QC_O class PowerNet: public Net { public: - bool add_3d_nlin = false; + bool add_3d_nlin = true; string ref_bus; double bMVA; /**< Base MVA */ double bV; /**< Base Voltage */ diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index baa9b3755..99d9fa1a5 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -564,12 +564,13 @@ param Bag::nfp(){ } } #endif - what.set_name("w_hat"); - what.print(true); +// what.set_name("w_hat"); +// what.print(true); return what; } param Bag::nfp1(){ +// fill_wstar(); param what; int n = _nodes.size(); DebugOff("\nn = " << n); @@ -588,13 +589,13 @@ param Bag::nfp1(){ } if (free_lines > 1) { - DebugOff("Returning empty w_hat."); + DebugOn("Returning empty w_hat."); return what; } //the bag has all lines or has > 3 nodes - double tol = 0.0001; + double tol = 1e-6; arma::cx_mat A(n,n); for(int i = 0; i < n; i++){ @@ -602,48 +603,60 @@ param Bag::nfp1(){ if(i==j) { A(i,j) = arma::cx_double(((Bus*)_grid->get_node(_nodes[i]->_name))->w,0); }else{ - double AijI = ((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wi; - if(_grid->has_directed_arc(_nodes[j],_nodes[i])) { - A(i,j) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,AijI); - A(j,i) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,-AijI); + + if(_grid->has_directed_arc(_nodes[i],_nodes[j])) { + double AijI = ((Line*)_grid->get_arc(_nodes[i], _nodes[j]))->wi; + double AijR = ((Line*)_grid->get_arc(_nodes[i], _nodes[j]))->wr; + A(i,j) = arma::cx_double(AijR,AijI); + A(j,i) = arma::cx_double(AijR,-AijI); } else { - A(i,j) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,-AijI); - A(j,i) = arma::cx_double(((Line*)_grid->get_arc(_nodes[j],_nodes[i]))->wr,AijI); + double AijI = ((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wi; + double AijR = ((Line*)_grid->get_arc(_nodes[j], _nodes[i]))->wr; + A(i,j) = arma::cx_double(AijR,-AijI); + A(j,i) = arma::cx_double(AijR,AijI); } } } } - arma::cx_mat Eigval; arma::cx_mat W_hat; - arma::vec eigvec; + arma::cx_mat Eigvec; arma::cx_mat W_hat; + arma::vec eigval; DebugOn("x_star = "); A.print(); - arma::eig_sym(eigvec,Eigval,A); + arma::cx_mat B = (A+A.t())/2; + bool success = arma::eig_sym(eigval,Eigvec,B); + if (!success) { + throw invalid_argument("Matrix could not be decomposed"); + } DebugOn("\n"); +// DebugOn("w_star = "); +// _wstarp.print(true); + double min_eig = 0, max_eig = -1; - for(auto eig: eigvec) { + for(auto eig: eigval) { if(eig < min_eig) min_eig = eig; if(eig > max_eig) max_eig = eig; } if(min_eig/max_eig > -tol) { - DebugOn("\nBag is PSD"); + DebugOff("\nBag is PSD"); return what; } else { - DebugOn("\nBag is not PSD\n"); + DebugOff("\nBag is not PSD\n"); for(int i = 0; i < n; i++) { - if(eigvec[i] < 0) { - eigvec[i] = 0; -// Eigval.col(i).zeros(); + if(eigval(i) < 0) { + eigval.at(i) = 0; +// Eigvec.col(i).zeros(); +// Eigvec.row(i).zeros(); } } - arma::mat Eigvec(n,n); - Eigvec.zeros(); - Eigvec.diag() = eigvec; - cout << "\nEigvec:" << endl; Eigval.print(); cout << "\nD:" << endl; Eigvec.print(); - W_hat = Eigval*Eigvec*Eigval.t(); - DebugOn("x_hat = "); + + arma::mat D(n,n); + D.zeros(); + D.diag() = eigval; + W_hat = Eigvec*D*Eigvec.t(); + DebugOn("x_hat = \n"); W_hat.print(); } @@ -659,10 +672,10 @@ param Bag::nfp1(){ namewr = "wr(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; namewi = "wi(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")"; what.set_val(namewr, W_hat(i,j).real()); - if(_grid->get_directed_arc(_nodes[i]->_name,_nodes[j]->_name)!=nullptr) - what.set_val(namewi, -1*W_hat(i,j).imag()); - else + if(_grid->has_directed_arc(_nodes[i],_nodes[j])) what.set_val(namewi, W_hat(i,j).imag()); + else + what.set_val(namewi, -1*W_hat(i,j).imag()); } } } diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index 53ff383d6..21e31acc7 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -22,8 +22,8 @@ using namespace gravity; void nearest_point(int i, int j, vector>& w_hat, const vector& bags){ for (unsigned index = i; index ("I_star"+to_string(numcuts)); - W_star.resize(numcuts+1); - W_star[numcuts] = param<>("W_star"+to_string(numcuts)); +// 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); @@ -369,7 +369,7 @@ int main_new (int argc, char * argv[]) { 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); +// W_star[numcuts].set_val(ni->_name,((Bus *) ni)->w); } else { aij = grid.get_arc(b._nodes[i]->_name, b._nodes[j]->_name); @@ -386,13 +386,15 @@ int main_new (int argc, char * argv[]) { 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); +// 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++; - + DebugOn("W_hat = \n"); + W_hat[numcuts].print(true); + cout << "\n"; 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))); @@ -698,7 +700,7 @@ int main (int argc, char * argv[]) { node_pairs bus_pairs_sdp; DebugOff("\nNumber of bags = " << bags.size()); - vector> R_star, I_star, W_star; +// vector> R_star, I_star, W_star; vector> R_diff, I_diff, W_diff; vector> R_hat, I_hat, W_hat; @@ -743,12 +745,12 @@ int main (int argc, char * argv[]) { 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_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); @@ -773,7 +775,7 @@ int main (int argc, char * argv[]) { 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); +// W_star[numcuts].set_val(ni->_name,((Bus *) ni)->w); } else { aij = grid.get_arc(b._nodes[i]->_name, b._nodes[j]->_name); @@ -790,18 +792,22 @@ int main (int argc, char * argv[]) { 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); +// 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++; +// cout << "W_HAT= \n"; +// W_hat[numcuts].print(true); +// cout << "WI_HAT= \n"; +// I_hat[numcuts].print(true); 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(); +// lin.print_expanded(); SDP.add_constraint(lin <= 0); cuts_added++; numcuts++; diff --git a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp index 41f3cb5d9..592f213a7 100644 --- a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp +++ b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp @@ -151,41 +151,39 @@ int main (int argc, char * argv[]) Constraint PAD_UB("PAD_UB"); PAD_UB = Im_Wij; PAD_UB <= grid.tan_th_max*R_Wij; - SOCP.add_lazy(PAD_UB.in(bus_pairs)); + SOCP.add(PAD_UB.in(bus_pairs)); Constraint PAD_LB("PAD_LB"); PAD_LB = Im_Wij; PAD_LB >= grid.tan_th_min*R_Wij; - SOCP.add_lazy(PAD_LB.in(bus_pairs)); + SOCP.add(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_lazy(Thermal_Limit_from.in(grid.arcs)); -// SOCP.add(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_lazy(Thermal_Limit_to.in(grid.arcs)); -// SOCP.add(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_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_lazy(LNC2.in(bus_pairs) >= 0); +// 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(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(LNC2.in(bus_pairs) >= 0); /* Solver selection */ From 3591882ad5c02ca593af79bde13fe5ef7b957bf7 Mon Sep 17 00:00:00 2001 From: kbestuzheva Date: Tue, 27 Feb 2018 17:14:02 +1100 Subject: [PATCH 07/10] Checking the feasibility of the solution from Marecek et al --- examples/MINLP/Power/ACOPF/ACOPF_main.cpp | 36 +++++++++++++++++++---- 1 file changed, 30 insertions(+), 6 deletions(-) diff --git a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp index f6bc0116d..7057ef64c 100644 --- a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp +++ b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp @@ -46,6 +46,7 @@ int main (int argc, char * argv[]) fname = opt["f"]; mtype = opt["m"]; output = op::str2int(opt["l"]); + output = 5; bool has_help = op::str2bool(opt["h"]); /** show help */ if(has_help) { @@ -95,8 +96,31 @@ int main (int argc, char * argv[]) ACOPF.add_var(Pf_to.in(grid.arcs)); ACOPF.add_var(Qf_to.in(grid.arcs)); + vector rev, imv; + ifstream file("/home/kbestuzheva/gravity/nesta_case30_fsr__api.log", std::ifstream::in); + if(!file.is_open()) throw invalid_argument("Could not open file\n"); + string word; + for(int i = 0; i < 30; i++){ + file >> word; + rev.push_back(atof(word.c_str())); + cout << "\nrev[" << i << "] = " << atof(word.c_str()); + } + for(int i = 0; i < 30; i++){ + file >> word; + imv.push_back(atof(word.c_str())); + cout << "\nimv[" << i << "] = " << atof(word.c_str()); + } + + param t_min, t_max; + for(int i = 0; i < 30; i++){ + grid.v_min.set_val(grid.nodes[i]->_name, sqrt(rev[i]*rev[i] + imv[i]*imv[i])); + grid.v_max.set_val(grid.nodes[i]->_name, sqrt(rev[i]*rev[i] + imv[i]*imv[i])); + t_min.set_val(grid.nodes[i]->_name, atan(imv[i]/rev[i])); + t_max.set_val(grid.nodes[i]->_name, atan(imv[i]/rev[i])); + } + /** Voltage related variables */ - var theta("theta"); + var theta("theta", t_min, t_max); var v("|V|", grid.v_min, grid.v_max); // var vr("vr"); // var vi("vi"); @@ -242,25 +266,25 @@ 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); double solver_time_start = get_wall_time(); - OPF.run(output, relax = false, tol = 1e-6, "ma27", mehrotra = "no"); + OPF.run(output, relax = false, tol = 1e-4, "ma27", mehrotra = "no"); double solver_time_end = get_wall_time(); double total_time_end = get_wall_time(); auto solve_time = solver_time_end - solver_time_start; From 3f0c70f04ebf37a1347d630d842b18ba8819d75d Mon Sep 17 00:00:00 2001 From: kbestuzheva Date: Tue, 27 Feb 2018 17:21:40 +1100 Subject: [PATCH 08/10] Fixed the filename --- examples/MINLP/Power/ACOPF/ACOPF_main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp index 7057ef64c..2637fd2d7 100644 --- a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp +++ b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp @@ -97,7 +97,7 @@ int main (int argc, char * argv[]) ACOPF.add_var(Qf_to.in(grid.arcs)); vector rev, imv; - ifstream file("/home/kbestuzheva/gravity/nesta_case30_fsr__api.log", std::ifstream::in); + ifstream file("../nesta_case30_fsr__api.log", std::ifstream::in); if(!file.is_open()) throw invalid_argument("Could not open file\n"); string word; for(int i = 0; i < 30; i++){ From 976867268abf83f0cf9a1192d640b2486ad60e1c Mon Sep 17 00:00:00 2001 From: Hassan Date: Tue, 13 Mar 2018 17:03:05 -0600 Subject: [PATCH 09/10] fixing reverse parallel arcs --- examples/MINLP/Power/PowerNet.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/examples/MINLP/Power/PowerNet.cpp b/examples/MINLP/Power/PowerNet.cpp index d3fcbe199..53503f724 100755 --- a/examples/MINLP/Power/PowerNet.cpp +++ b/examples/MINLP/Power/PowerNet.cpp @@ -375,13 +375,16 @@ int PowerNet::readgrid(const char* fname) { string src,dest,key; file >> word; index = 0; + bool reversed = false; while(word.compare("];")) { src = word; file >> dest; key = dest+","+src;//Taking care of reversed direction arcs + reversed = false; if(arcID.find(key)!=arcID.end()) {//Reverse arc direction DebugOn("Adding arc linking " +src+" and "+dest); DebugOn(" with reversed direction, reversing source and destination.\n"); + reversed = true; key = src; src = dest; dest = key; @@ -438,6 +441,13 @@ int PowerNet::readgrid(const char* fname) { arc->tbound.max = 60*pi/180; } + if (reversed) { + arc->tr /= 1; + arc->as *= -1; + auto temp = arc->tbound.max; + arc->tbound.max = -1*arc->tbound.min; + arc->tbound.min = -1*temp; + } // arc->tbound.max = 30*pi/180; m_theta_ub += arc->tbound.max; From 51f6dfdc58e5b8a857e2927ab1a577882ce1c706 Mon Sep 17 00:00:00 2001 From: Hassan Date: Tue, 13 Mar 2018 17:31:14 -0600 Subject: [PATCH 10/10] supporting add(var) and add(constr) and cleaning up acopf and socopf, to merge with master --- examples/MINLP/Power/ACOPF/ACOPF_main.cpp | 81 +++++++-------------- examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp | 26 +++---- include/gravity/model.h | 10 +++ src/model.cpp | 1 + 4 files changed, 49 insertions(+), 69 deletions(-) diff --git a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp index 2637fd2d7..7644a37a4 100644 --- a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp +++ b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp @@ -82,8 +82,8 @@ int main (int argc, char * argv[]) /* Power generation variables */ var Pg("Pg", grid.pg_min, grid.pg_max); var Qg ("Qg", grid.qg_min, grid.qg_max); - ACOPF.add_var(Pg.in(grid.gens)); - ACOPF.add_var(Qg.in(grid.gens)); + ACOPF.add(Pg.in(grid.gens)); + ACOPF.add(Qg.in(grid.gens)); /* Power flow variables */ var Pf_from("Pf_from", grid.S_max); @@ -91,50 +91,26 @@ int main (int argc, char * argv[]) var Pf_to("Pf_to", grid.S_max); var Qf_to("Qf_to", grid.S_max); - ACOPF.add_var(Pf_from.in(grid.arcs)); - ACOPF.add_var(Qf_from.in(grid.arcs)); - ACOPF.add_var(Pf_to.in(grid.arcs)); - ACOPF.add_var(Qf_to.in(grid.arcs)); + ACOPF.add(Pf_from.in(grid.arcs)); + ACOPF.add(Qf_from.in(grid.arcs)); + ACOPF.add(Pf_to.in(grid.arcs)); + ACOPF.add(Qf_to.in(grid.arcs)); - vector rev, imv; - ifstream file("../nesta_case30_fsr__api.log", std::ifstream::in); - if(!file.is_open()) throw invalid_argument("Could not open file\n"); - string word; - for(int i = 0; i < 30; i++){ - file >> word; - rev.push_back(atof(word.c_str())); - cout << "\nrev[" << i << "] = " << atof(word.c_str()); - } - for(int i = 0; i < 30; i++){ - file >> word; - imv.push_back(atof(word.c_str())); - cout << "\nimv[" << i << "] = " << atof(word.c_str()); - } - - param t_min, t_max; - for(int i = 0; i < 30; i++){ - grid.v_min.set_val(grid.nodes[i]->_name, sqrt(rev[i]*rev[i] + imv[i]*imv[i])); - grid.v_max.set_val(grid.nodes[i]->_name, sqrt(rev[i]*rev[i] + imv[i]*imv[i])); - t_min.set_val(grid.nodes[i]->_name, atan(imv[i]/rev[i])); - t_max.set_val(grid.nodes[i]->_name, atan(imv[i]/rev[i])); - } /** Voltage related variables */ - var theta("theta", t_min, t_max); + var theta("theta"); var v("|V|", grid.v_min, grid.v_max); -// var vr("vr"); -// var vi("vi"); var vr("vr", grid.v_max); var vi("vi", grid.v_max); if (polar) { - ACOPF.add_var(v.in(grid.nodes)); - ACOPF.add_var(theta.in(grid.nodes)); + ACOPF.add(v.in(grid.nodes)); + ACOPF.add(theta.in(grid.nodes)); v.initialize_all(1); } else { - ACOPF.add_var(vr.in(grid.nodes)); - ACOPF.add_var(vi.in(grid.nodes)); + ACOPF.add(vr.in(grid.nodes)); + ACOPF.add(vi.in(grid.nodes)); vr.initialize_all(1.0); } @@ -152,7 +128,7 @@ int main (int argc, char * argv[]) else { Ref_Bus = vi(grid.get_ref_bus()); } - ACOPF.add_constraint(Ref_Bus == 0); + ACOPF.add(Ref_Bus == 0); /** KCL Flow conservation */ Constraint KCL_P("KCL_P"); @@ -168,8 +144,8 @@ int main (int argc, char * argv[]) KCL_P += grid.gs*(power(vr,2)+power(vi,2)); KCL_Q -= grid.bs*(power(vr,2)+power(vi,2)); } - ACOPF.add_constraint(KCL_P.in(grid.nodes) == 0); - ACOPF.add_constraint(KCL_Q.in(grid.nodes) == 0); + ACOPF.add(KCL_P.in(grid.nodes) == 0); + ACOPF.add(KCL_Q.in(grid.nodes) == 0); /** AC Power Flows */ /** TODO write the constraints in Complex form */ @@ -185,7 +161,7 @@ int main (int argc, char * argv[]) Flow_P_From -= grid.g_ft*(vr.from()*vr.to() + vi.from()*vi.to()); Flow_P_From -= grid.b_ft*(vi.from()*vr.to() - vr.from()*vi.to()); } - ACOPF.add_constraint(Flow_P_From.in(grid.arcs)==0); + ACOPF.add(Flow_P_From.in(grid.arcs)==0); Constraint Flow_P_To("Flow_P_To"); Flow_P_To += Pf_to; @@ -199,7 +175,7 @@ int main (int argc, char * argv[]) Flow_P_To -= grid.g_tf*(vr.from()*vr.to() + vi.from()*vi.to()); Flow_P_To -= grid.b_tf*(vi.to()*vr.from() - vr.to()*vi.from()); } - ACOPF.add_constraint(Flow_P_To.in(grid.arcs)==0); + ACOPF.add(Flow_P_To.in(grid.arcs)==0); Constraint Flow_Q_From("Flow_Q_From"); Flow_Q_From += Qf_from; @@ -213,7 +189,7 @@ int main (int argc, char * argv[]) Flow_Q_From += grid.b_ft*(vr.from()*vr.to() + vi.from()*vi.to()); Flow_Q_From -= grid.g_ft*(vi.from()*vr.to() - vr.from()*vi.to()); } - ACOPF.add_constraint(Flow_Q_From.in(grid.arcs)==0); + ACOPF.add(Flow_Q_From.in(grid.arcs)==0); Constraint Flow_Q_To("Flow_Q_To"); Flow_Q_To += Qf_to; if (polar) { @@ -226,21 +202,19 @@ int main (int argc, char * argv[]) Flow_Q_To += grid.b_tf*(vr.from()*vr.to() + vi.from()*vi.to()); Flow_Q_To -= grid.g_tf*(vi.to()*vr.from() - vr.to()*vi.from()); } - ACOPF.add_constraint(Flow_Q_To.in(grid.arcs)==0); + ACOPF.add(Flow_Q_To.in(grid.arcs)==0); /** AC voltage limit constraints. */ if (!polar) { Constraint Vol_limit_UB("Vol_limit_UB"); Vol_limit_UB = power(vr, 2) + power(vi, 2); Vol_limit_UB -= power(grid.v_max, 2); - ACOPF.add_constraint(Vol_limit_UB.in(grid.nodes) <= 0); + ACOPF.add(Vol_limit_UB.in(grid.nodes) <= 0); Constraint Vol_limit_LB("Vol_limit_LB"); Vol_limit_LB = power(vr, 2) + power(vi, 2); Vol_limit_LB -= power(grid.v_min,2); - ACOPF.add_constraint(Vol_limit_LB.in(grid.nodes) >= 0); - DebugOff(grid.v_min.to_str(true) << endl); - DebugOff(grid.v_max.to_str(true) << endl); + ACOPF.add(Vol_limit_LB.in(grid.nodes) >= 0); } @@ -253,8 +227,6 @@ int main (int argc, char * argv[]) PAD_UB -= grid.th_max; PAD_LB = theta.from() - theta.to(); PAD_LB -= grid.th_min; - DebugOff(grid.th_min.to_str(true) << endl); - DebugOff(grid.th_max.to_str(true) << endl); } else { DebugOff("Number of bus_pairs = " << bus_pairs.size() << endl); @@ -263,28 +235,25 @@ int main (int argc, char * argv[]) PAD_LB = vi.from()*vr.to() - vr.from()*vi.to(); PAD_LB -= grid.tan_th_min*(vr.from()*vr.to() + vi.from()*vi.to()); - 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(PAD_UB.in(bus_pairs) <= 0); + ACOPF.add(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(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); - DebugOff(grid.S_max.in(grid.arcs).to_str(true) << endl); + ACOPF.add(Thermal_Limit_to.in(grid.arcs) <= 0); solver OPF(ACOPF,ipopt); double solver_time_start = get_wall_time(); - OPF.run(output, relax = false, tol = 1e-4, "ma27", mehrotra = "no"); + OPF.run(output, relax = false, tol = 1e-6, "ma27", mehrotra = "no"); double solver_time_end = get_wall_time(); double total_time_end = get_wall_time(); auto solve_time = solver_time_end - solver_time_start; diff --git a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp index 592f213a7..9eb034c49 100644 --- a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp +++ b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp @@ -171,19 +171,19 @@ int main (int argc, char * argv[]) 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(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(LNC2.in(bus_pairs) >= 0); + 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(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(LNC2.in(bus_pairs) >= 0); /* Solver selection */ diff --git a/include/gravity/model.h b/include/gravity/model.h index 33b6c268d..0ee8b9b95 100755 --- a/include/gravity/model.h +++ b/include/gravity/model.h @@ -158,6 +158,11 @@ namespace gravity { } }; + template + void add(var& v){//Add variables by copy + add_var(v); + } + template void add_var(var&& v){//Add variables by copy @@ -200,6 +205,11 @@ namespace gravity { } }; + template + void add(var&& v){ + add_var(move(v)); + } + void del_var(const param_& v); diff --git a/src/model.cpp b/src/model.cpp index 1dd33b26e..34560bf88 100755 --- a/src/model.cpp +++ b/src/model.cpp @@ -232,6 +232,7 @@ void Model::del_param(const param_& v){ } }; + void Model::add(const Constraint& c){ add_constraint(c); }