Merge pull request #25 from Allinsights/SDP_cuts

Dropping the ^ operator, introducing R(n)

examples/Gravity101_main.cpp

@@ -417,7 +417,7 @@ int main (int argc, const char * argv[])
 
     /* Declare model */
     Model SOCP("Second-Order Cone Model");
-    SOCP.add_var(x^10); /* Will add a vector of size 10 representing variables named x */
+    SOCP.add_var(x.in(R(10))); /* Will add a vector of size 10 representing variables named x */
     SOCP.add_constraint(SOC.in(bus_pairs) <= 0); /* Will add second-order constraints indexed by bus_pairs (see previous Code Block) */
     SOCP.min(x+2*y); /* Declaring the objective function */
     return 0;

examples/MILP/MinKPart/Minkmodel.cpp

@@ -73,7 +73,8 @@ void Minkmodel::build() {
 void Minkmodel::reset() {};
 void Minkmodel::add_vars_origin() {
     var<bool> zij("zij");
-    _model.add_var(zij^(_graph->nodes.size()*(_graph->nodes.size()-1)/2));
+    auto Rn = R(_graph->nodes.size()*(_graph->nodes.size()-1)/2);
+    _model.add_var(zij.in(Rn));
 
     func_ obj_MIP;
     int i=0, j=0;
@@ -90,7 +91,8 @@ void Minkmodel::add_vars_origin() {
 
 void Minkmodel::add_vars_origin_tree() {
     var<bool> zij("zij");
-    _model.add_var(zij^((_chordal_extension)->arcs.size()));
+    auto Rn = R(((_chordal_extension)->arcs.size()));
+    _model.add_var(zij.in(Rn));
 
     func_ obj_MIP;
     int i=0, j=0;
@@ -503,10 +505,12 @@ void Minkmodel::add_bicycle() {}
 //}
 void Minkmodel::node_edge_formulation(){
      var<bool> x("x");
-     _model.add_var(x^((int)(_K*_graph->nodes.size())));
+    auto Rk = R(_K*_graph->nodes.size());
+     _model.add_var(x.in(Rk));
      var<bool> y("y");
     // the number of arcs in the chordal extension
-    _model.add_var(y^(_graph->arcs.size()));
+    auto Rm = R(_graph->arcs.size());
+    _model.add_var(y.in(Rm));
 
     func_ obj_node_edge;
     int i=0, j=0;

examples/MILP/StableSet/Stable_Set_main.cpp

@@ -39,9 +39,10 @@ int main (int argc, const char * argv[])
     /**  IP model for the stable set problem. **/
 
     /**  Variables **/
+    /* Declaring the n-dimensional Real space */
+    auto Rn = R(n);
     var<bool> x("x");
-//    var<> x("x", 0,1);
-    model.add_var(x^n);
+    model.add_var(x.in(Rn));
 
     /**  Objective **/
     model.max(sum(x));
@@ -70,8 +71,8 @@ int main (int argc, const char * argv[])
     /* Variable declaration */
     var<double> Xii("Xii", 0, 1);
     var<double> Xij("Xij", 0, 1);
-    SDP.add_var(Xii^n); /*< Diagonal entries of the matrix */
-    SDP.add_var(Xij^(n*(n-1)/2)); /*< Lower left triangular part of the matrix excluding the diagonal*/
+    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);

examples/MINLP/Power/ACUC/ACUC_main.cpp

@@ -71,18 +71,20 @@ int main (int argc, const char * argv[])
     // POWER GENERATION
     var<Real> Pg("Pg", grid->pg_min.in(grid->gens, T), grid->pg_max.in(grid->gens, T));
     var<Real> Qg ("Qg", grid->qg_min.in(grid->gens, T), grid->qg_max.in(grid->gens, T));
-    ACUC.add_var(Pg^(T*nb_gen));
-    ACUC.add_var(Qg^(T*nb_gen));
+    auto R_Tng = R(T*nb_gen);
+    ACUC.add_var(Pg.in(R_Tng));
+    ACUC.add_var(Qg.in(R_Tng));
 
      //power flow
     var<Real> Pf_from("Pf_from", grid->S_max.in(grid->arcs, T));
     var<Real> Qf_from("Qf_from", grid->S_max.in(grid->arcs, T));
     var<Real> Pf_to("Pf_to", grid->S_max.in(grid->arcs, T));
     var<Real> Qf_to("Qf_to", grid->S_max.in(grid->arcs, T));
-    ACUC.add_var(Pf_from^(T*nb_lines));
-    ACUC.add_var(Qf_from^(T*nb_lines));
-    ACUC.add_var(Pf_to^(T*nb_lines));
-    ACUC.add_var(Qf_to^(T*nb_lines));
+    auto R_Tnl = R(T*nb_lines);
+    ACUC.add_var(Pf_from.in(R_Tnl));
+    ACUC.add_var(Qf_from.in(R_Tnl));
+    ACUC.add_var(Pf_to.in(R_Tnl));
+    ACUC.add_var(Qf_to.in(R_Tnl));
 
      //Lifted variables.
     var<Real>  R_Wij("R_Wij", grid->wr_min.in(bus_pairs, T), grid->wr_max.in(bus_pairs, T)); // real part of Wij
@@ -91,19 +93,21 @@ int main (int argc, const char * argv[])
     R_Wij.print(true);
     Im_Wij.print(true);
 
-    ACUC.add_var(Wii^(T*nb_buses));
-    ACUC.add_var(R_Wij^(T*nb_bus_pairs));
-    ACUC.add_var(Im_Wij^(T*nb_bus_pairs));
+    auto R_Tnb = R(T*nb_buses);
+    auto R_Tnbp = R(T*nb_bus_pairs);
+    ACUC.add_var(Wii.in(R_Tnb));
+    ACUC.add_var(R_Wij.in(R_Tnbp));
+    ACUC.add_var(Im_Wij.in(R_Tnbp));
     R_Wij.initialize_all(1.0);
     Wii.initialize_all(1.001);
 
     // Commitment variables
     var<Real>  On_off("On_off");
     var<Real>  Start_up("Start_up", 0, 1);
     var<Real>  Shut_down("Shut_down", 0, 1);
-    ACUC.add_var(On_off^(T*nb_gen));
-    ACUC.add_var(Start_up^(T*nb_gen));
-    ACUC.add_var(Shut_down^(T*nb_gen));
+    ACUC.add_var(On_off.in(R_Tng));
+    ACUC.add_var(Start_up.in(R_Tng));
+    ACUC.add_var(Shut_down.in(R_Tng));
 
     /* Construct the objective function*/
     func_ obj;

examples/MINLP/Power/PowerNet.cpp

@@ -22,6 +22,7 @@
 #include <math.h>
 #include <queue>
 #include <time.h>
+#include <gravity/solver.h>
 //#define USEDEBUG
 #ifdef USEDEBUG
 #define Debug(x) cout << x
@@ -745,3 +746,200 @@ void PowerNet::update_net(){
     }
 }
 
+double PowerNet::solve_acopf(PowerModelType pmt, int output, double tol){
+
+    bool polar = (pmt==ACPOL);
+    if (polar) {
+        DebugOn("Using polar model\n");
+    }
+    else {
+        DebugOn("Using rectangular model\n");
+    }
+    Model ACOPF("AC-OPF Model");
+    /** Variables */
+    /* Power generation variables */
+    var<Real> Pg("Pg", pg_min, pg_max);
+    var<Real> Qg ("Qg", qg_min, qg_max);
+    ACOPF.add_var(Pg.in(gens));
+    ACOPF.add_var(Qg.in(gens));
+
+    /* Power flow variables */
+    var<Real> Pf_from("Pf_from", S_max);
+    var<Real> Qf_from("Qf_from", S_max);
+    var<Real> Pf_to("Pf_to", S_max);
+    var<Real> Qf_to("Qf_to", S_max);
+
+    ACOPF.add_var(Pf_from.in(arcs));
+    ACOPF.add_var(Qf_from.in(arcs));
+    ACOPF.add_var(Pf_to.in(arcs));
+    ACOPF.add_var(Qf_to.in(arcs));
+
+    /** Voltage related variables */
+    var<Real> theta("theta");
+    var<Real> v("|V|", v_min, v_max);
+    //    var<Real> vr("vr");
+    //    var<Real> vi("vi");
+    var<Real> vr("vr", v_max);
+    var<Real> vi("vi", v_max);
+
+    if (polar) {
+        ACOPF.add_var(v.in(nodes));
+        ACOPF.add_var(theta.in(nodes));
+        v.initialize_all(1);
+    }
+    else {
+        ACOPF.add_var(vr.in(nodes));
+        ACOPF.add_var(vi.in(nodes));
+        vr.initialize_all(1.0);
+    }
+
+    /** Construct the objective function */
+    func_ obj = product(c1, Pg) + product(c2, power(Pg,2)) + sum(c0);
+    ACOPF.min(obj.in(gens));
+
+    /** Define constraints */
+
+    /* REF BUS */
+    Constraint Ref_Bus("Ref_Bus");
+    if (polar) {
+        Ref_Bus = theta(get_ref_bus());
+    }
+    else {
+        Ref_Bus = vi(get_ref_bus());
+    }
+    ACOPF.add_constraint(Ref_Bus == 0);
+
+    /** KCL Flow conservation */
+    Constraint KCL_P("KCL_P");
+    Constraint KCL_Q("KCL_Q");
+    KCL_P  = sum(Pf_from.out_arcs()) + sum(Pf_to.in_arcs()) + pl - sum(Pg.in_gens());
+    KCL_Q  = sum(Qf_from.out_arcs()) + sum(Qf_to.in_arcs()) + ql - sum(Qg.in_gens());
+    /* Shunts */
+    if (polar) {
+        KCL_P +=  gs*power(v,2);
+        KCL_Q -=  bs*power(v,2);
+    }
+    else {
+        KCL_P +=  gs*(power(vr,2)+power(vi,2));
+        KCL_Q -=  bs*(power(vr,2)+power(vi,2));
+    }
+    ACOPF.add_constraint(KCL_P.in(nodes) == 0);
+    ACOPF.add_constraint(KCL_Q.in(nodes) == 0);
+
+    /** AC Power Flows */
+    /** TODO write the constraints in Complex form */
+    Constraint Flow_P_From("Flow_P_From");
+    Flow_P_From += Pf_from;
+    if (polar) {
+        Flow_P_From -= g/power(tr,2)*power(v.from(),2);
+        Flow_P_From += g/tr*(v.from()*v.to()*cos(theta.from() - theta.to() - as));
+        Flow_P_From += b/tr*(v.from()*v.to()*sin(theta.from() - theta.to() - as));
+    }
+    else {
+        Flow_P_From -= g_ff*(power(vr.from(), 2) + power(vi.from(), 2));
+        Flow_P_From -= g_ft*(vr.from()*vr.to() + vi.from()*vi.to());
+        Flow_P_From -= b_ft*(vi.from()*vr.to() - vr.from()*vi.to());
+    }
+    ACOPF.add_constraint(Flow_P_From.in(arcs)==0);
+
+    Constraint Flow_P_To("Flow_P_To");
+    Flow_P_To += Pf_to;
+    if (polar) {
+        Flow_P_To -= g*power(v.to(), 2);
+        Flow_P_To += g/tr*(v.from()*v.to()*cos(theta.to() - theta.from() + as));
+        Flow_P_To += b/tr*(v.from()*v.to()*sin(theta.to() - theta.from() + as));
+    }
+    else {
+        Flow_P_To -= g_tt*(power(vr.to(), 2) + power(vi.to(), 2));
+        Flow_P_To -= g_tf*(vr.from()*vr.to() + vi.from()*vi.to());
+        Flow_P_To -= b_tf*(vi.to()*vr.from() - vr.to()*vi.from());
+    }
+    ACOPF.add_constraint(Flow_P_To.in(arcs)==0);
+
+    Constraint Flow_Q_From("Flow_Q_From");
+    Flow_Q_From += Qf_from;
+    if (polar) {
+        Flow_Q_From += (0.5*ch+b)/power(tr,2)*power(v.from(),2);
+        Flow_Q_From -= b/tr*(v.from()*v.to()*cos(theta.from() - theta.to() - as));
+        Flow_Q_From += g/tr*(v.from()*v.to()*sin(theta.from() - theta.to() - as));
+    }
+    else {
+        Flow_Q_From += b_ff*(power(vr.from(), 2) + power(vi.from(), 2));
+        Flow_Q_From += b_ft*(vr.from()*vr.to() + vi.from()*vi.to());
+        Flow_Q_From -= g_ft*(vi.from()*vr.to() - vr.from()*vi.to());
+    }
+    ACOPF.add_constraint(Flow_Q_From.in(arcs)==0);
+    Constraint Flow_Q_To("Flow_Q_To");
+    Flow_Q_To += Qf_to;
+    if (polar) {
+        Flow_Q_To += (0.5*ch+b)*power(v.to(),2);
+        Flow_Q_To -= b/tr*(v.from()*v.to()*cos(theta.to() - theta.from() + as));
+        Flow_Q_To += g/tr*(v.from()*v.to()*sin(theta.to() - theta.from() + as));
+    }
+    else {
+        Flow_Q_To += b_tt*(power(vr.to(), 2) + power(vi.to(), 2));
+        Flow_Q_To += b_tf*(vr.from()*vr.to() + vi.from()*vi.to());
+        Flow_Q_To -= g_tf*(vi.to()*vr.from() - vr.to()*vi.from());
+    }
+    ACOPF.add_constraint(Flow_Q_To.in(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(v_max, 2);
+        ACOPF.add_constraint(Vol_limit_UB.in(nodes) <= 0);
+
+        Constraint Vol_limit_LB("Vol_limit_LB");
+        Vol_limit_LB = power(vr, 2) + power(vi, 2);
+        Vol_limit_LB -= power(v_min,2);
+        ACOPF.add_constraint(Vol_limit_LB.in(nodes) >= 0);
+        DebugOff(v_min.to_str(true) << endl);
+        DebugOff(v_max.to_str(true) << endl);
+    }
+
+
+    /* Phase Angle Bounds constraints */
+    Constraint PAD_UB("PAD_UB");
+    Constraint PAD_LB("PAD_LB");
+    auto bus_pairs = get_bus_pairs();
+    if (polar) {
+        PAD_UB = theta.from() - theta.to();
+        PAD_UB -= th_max;
+        PAD_LB = theta.from() - theta.to();
+        PAD_LB -= th_min;
+        DebugOff(th_min.to_str(true) << endl);
+        DebugOff(th_max.to_str(true) << endl);
+    }
+    else {
+        DebugOff("Number of bus_pairs = " << bus_pairs.size() << endl);
+        PAD_UB = vi.from()*vr.to() - vr.from()*vi.to();
+        PAD_UB -= tan_th_max*(vr.from()*vr.to() + vi.from()*vi.to());
+
+        PAD_LB = vi.from()*vr.to() - vr.from()*vi.to();
+        PAD_LB -= tan_th_min*(vr.from()*vr.to() + vi.from()*vi.to());
+        DebugOff(th_min.to_str(true) << endl);
+        DebugOff(th_max.to_str(true) << endl);
+    }
+    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(S_max, 2);
+    ACOPF.add_constraint(Thermal_Limit_from.in(arcs) <= 0);
+
+    Constraint Thermal_Limit_to("Thermal_Limit_to");
+    Thermal_Limit_to += power(Pf_to, 2) + power(Qf_to, 2);
+    Thermal_Limit_to -= power(S_max,2);
+    ACOPF.add_constraint(Thermal_Limit_to.in(arcs) <= 0);
+    DebugOff(S_max.in(arcs).to_str(true) << endl);
+    bool relax;
+    solver OPF(ACOPF,ipopt);
+    OPF.run(output, relax = false, tol = 1e-6, "ma27");
+    return ACOPF._obj_val;
+}
+
+

examples/MINLP/Power/PowerNet.h

@@ -62,5 +62,8 @@ class PowerNet: public Net {
     unsigned get_nb_active_nodes() const;
     void time_expand(unsigned T); /* < Time expansion of the grid parameters */
     void update_net();
+
+    /** Power Models */
+    double solve_acopf(PowerModelType model=ACPOL, int output=0, double tol=1e-6);
 };
 #endif

examples/MINLP/Power/SDPOPF/Bag.cpp

@@ -330,6 +330,10 @@ bool Bag::is_PSD(){
     }
 }
 
+void Bag::update_PSD(){
+    _is_psd = is_PSD();
+}
+
 param<double> Bag::nfp(){
     param<double> what;
     fill_wstar();
@@ -340,18 +344,19 @@ param<double> Bag::nfp(){
 
     Model NPP("NPP model");
     int n = _nodes.size();
+    DebugOff("\nn = " << n);
 
     sdpvar<double> W("W");
     NPP.add_var(W^(2*n));
 
 //    var<double> W("W");
 //    W._psd = true;
-//    W._is_matrix = true;
-//    NPP.add_var(W ^ (n*(2*n-1)));
+//    NPP.add_var(W ^ (n*(2*n+1)));
 
     var<double> z("z");
     z.in_q_cone();
-    NPP.add_var(z^(n*n+1));
+    auto Rn = R(n*n+1);
+    NPP.add_var(z.in(Rn));
 
     var<double> w("w", _wmin, _wmax);
     NPP.add_var(w.in(_indices));

examples/MINLP/Power/SDPOPF/Bag.h

@@ -16,6 +16,7 @@
 class Bag{
 public:
     int _id;
+    bool _is_psd;
     PowerNet* _grid;
     std::vector<Node*> _nodes;
 //    bool _all_lines = true;
@@ -43,6 +44,8 @@ class Bag{
     bool add_lines();
 
     param<double> fill_wstar();
+
+    void update_PSD();
 };
 
 #endif //GRAVITY_BAG_H