Skip to content

Commit

Permalink
Merge pull request #28 from Allinsights/SDP_cuts
Browse files Browse the repository at this point in the history 
Sdp cuts
  • Loading branch information
@hhijazi
hhijazi authored Mar 13, 2018
2 parents 30f5afc + 51f6dfd commit 26060aa
Show file tree
Hide file tree
Showing 9 changed files with 243 additions and 88 deletions.
57 changes: 25 additions & 32 deletions examples/MINLP/Power/ACOPF/ACOPF_main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -81,36 +82,35 @@ int main (int argc, char * argv[])
/* Power generation variables */
var<Real> Pg("Pg", grid.pg_min, grid.pg_max);
var<Real> 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<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);

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));


/** Voltage related variables */
var<Real> theta("theta");
var<Real> v("|V|", grid.v_min, grid.v_max);
// var<Real> vr("vr");
// var<Real> vi("vi");
var<Real> vr("vr", grid.v_max);
var<Real> 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);
}

Expand All @@ -128,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");
Expand All @@ -144,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 */
Expand All @@ -161,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;
Expand All @@ -175,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;
Expand All @@ -189,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) {
Expand All @@ -202,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);
}


Expand All @@ -229,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);
Expand All @@ -239,24 +235,21 @@ 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();
Expand Down
10 changes: 10 additions & 0 deletions examples/MINLP/Power/PowerNet.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand Down
154 changes: 145 additions & 9 deletions examples/MINLP/Power/SDPOPF/Bag.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
#include <gravity/param.h>
#include "Bag.h"

#define FLAT
//#define FLAT

Bag::Bag(int id, const PowerNet& grid, vector<Node*> nodes):_id(id),_grid((PowerNet*)&grid),_nodes(nodes) {
Arc* aij = NULL;
Expand Down Expand Up @@ -152,6 +152,16 @@ param<double> 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;
Expand All @@ -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;
Expand Down Expand Up @@ -338,8 +346,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) {
Expand All @@ -351,7 +360,19 @@ bool Bag::is_PSD(){
return true;
}
else {
DebugOff("\nBag is not PSD");
double pos_tol = -0.00001;
DebugOff("\nBag is not PSD\n");
for(int i = 0; i < n; i++) {
if(v[i] < 0) {
v[i] = 0;
// P.col(i).zeros();
}
}
arma::mat D(n,n);
D.zeros();
D.diag() = v;
arma::cx_mat W_hat = P*D*P.t();

return false;
}
}
Expand Down Expand Up @@ -402,7 +423,7 @@ param<double> Bag::nfp(){
#endif
NPP.add_var(z.in(Rn));

var<double> w("w", _wmin, _wmax);
var<double> w("w");//, _wmin, _wmax);
NPP.add_var(w.in(_indices));

func_ obj;
Expand Down Expand Up @@ -543,8 +564,123 @@ param<double> Bag::nfp(){
}
}
#endif
what.set_name("w_hat");
// what.set_name("w_hat");
// what.print(true);
// exit(0);
return what;
}

param<double> Bag::nfp1(){
// fill_wstar();
param<double> 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) {
DebugOn("Returning empty w_hat.");
return what;
}

//the bag has all lines or has > 3 nodes

double tol = 1e-6;
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{

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 {
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 Eigvec; arma::cx_mat W_hat;
arma::vec eigval;
DebugOn("x_star = ");
A.print();
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: eigval) {
if(eig < min_eig) min_eig = eig;
if(eig > max_eig) max_eig = eig;
}

if(min_eig/max_eig > -tol) {
DebugOff("\nBag is PSD");
return what;
} else {
DebugOff("\nBag is not PSD\n");
for(int i = 0; i < n; i++) {
if(eigval(i) < 0) {
eigval.at(i) = 0;
// Eigvec.col(i).zeros();
// Eigvec.row(i).zeros();
}
}

arma::mat D(n,n);
D.zeros();
D.diag() = eigval;
W_hat = Eigvec*D*Eigvec.t();
DebugOn("x_hat = \n");
W_hat.print();
}

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->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());
}
}
}

what.set_name("w_hat");
what.print(true);
return what;
}
Loading

0 comments on commit 26060aa

Please sign in to comment.