From 90d38ea45f1dc656bb314dd78dde590a76c57af4 Mon Sep 17 00:00:00 2001 From: Ksenia Bestuzheva Date: Thu, 25 Jan 2018 11:12:32 -0700 Subject: [PATCH 1/5] Making sure bags are checked in the correct order --- examples/MINLP/Power/SDPOPF/Bag.cpp | 8 +++++-- examples/MINLP/Power/SDPOPF/Bag.h | 3 +++ examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 26 +++++++++++++-------- 3 files changed, 25 insertions(+), 12 deletions(-) diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index 5e1c3cbf7..2b90bbce8 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -330,6 +330,10 @@ bool Bag::is_PSD(){ } } +void Bag::update_PSD(){ + _is_psd = is_PSD(); +} + param Bag::nfp(){ param what; fill_wstar(); @@ -340,14 +344,14 @@ param Bag::nfp(){ Model NPP("NPP model"); int n = _nodes.size(); + DebugOff("\nn = " << n); sdpvar W("W"); NPP.add_var(W^(2*n)); // var 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 z("z"); z.in_q_cone(); diff --git a/examples/MINLP/Power/SDPOPF/Bag.h b/examples/MINLP/Power/SDPOPF/Bag.h index 9f088d436..b6c810708 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.h +++ b/examples/MINLP/Power/SDPOPF/Bag.h @@ -16,6 +16,7 @@ class Bag{ public: int _id; + bool _is_psd; PowerNet* _grid; std::vector _nodes; // bool _all_lines = true; @@ -43,6 +44,8 @@ class Bag{ bool add_lines(); param fill_wstar(); + + void update_PSD(); }; #endif //GRAVITY_BAG_H diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index a80eff0df..4e77dd848 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 ("psd"); DebugOff("Bag " << to_string(index) << " is PSD"< 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)); // SDP.add_var(R_Wij.in(bus_pairs)); @@ -334,9 +336,8 @@ int main (int argc, char * argv[]) { // int unchanged = 0; int iter = 0, hdim_cuts = 0, cuts_added = 1; - // while(unchanged < 3) { - unsigned nr_threads = 1, nb_bags = bags.size(); + unsigned nr_threads = 4, nb_bags = bags.size(); while(cuts_added > 0) { cuts_added = 0; vector> w_hat_vec; @@ -344,8 +345,13 @@ int main (int argc, char * argv[]) { iter++; // while((SDP._obj_val - prev_opt)/SDP._obj_val > fp_tol) { // prev_opt = SDP._obj_val; - - + + /* Update _is_psd */ + for(auto& b: bags){ + b.update_PSD(); + } + } + vector threads; /* Split subproblems into nr_threads parts */ vector limits = bounds(nr_threads, nb_bags); @@ -383,8 +389,8 @@ int main (int argc, char * argv[]) { I_hat[numcuts] = param<>("I_hat"+to_string(numcuts)); W_hat.resize(numcuts+1); W_hat[numcuts] = param<>("W_hat"+to_string(numcuts)); - - + +// what = b.nfp(); node_pairs b_pairs("node_pairs"); Constraint sdpcut("sdpcut_" + to_string(numcuts)); Node *ni; @@ -419,7 +425,7 @@ int main (int argc, char * argv[]) { } } 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))); From 8f4ab48372671c2e2e01dafd934fbe559d268b4b Mon Sep 17 00:00:00 2001 From: Ksenia Bestuzheva Date: Thu, 25 Jan 2018 13:17:21 -0700 Subject: [PATCH 2/5] Fixing typo --- examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index 4e77dd848..a6301568d 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -350,7 +350,6 @@ int main (int argc, char * argv[]) { for(auto& b: bags){ b.update_PSD(); } - } vector threads; /* Split subproblems into nr_threads parts */ From 6b6eb4957e8ef11cf2214be092ddc9b1044a7ef3 Mon Sep 17 00:00:00 2001 From: Hassan Date: Thu, 25 Jan 2018 14:57:27 -0700 Subject: [PATCH 3/5] Fixing bug in Constraint::is_active() --- examples/MINLP/Power/PowerNet.h | 2 +- examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 7 ++++--- src/constraint.cpp | 3 ++- src/solver.cpp | 2 +- 4 files changed, 8 insertions(+), 6 deletions(-) diff --git a/examples/MINLP/Power/PowerNet.h b/examples/MINLP/Power/PowerNet.h index 8b5283d32..a1ca745b6 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/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index a6301568d..d13247cc3 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -95,8 +95,8 @@ int main (int argc, char * argv[]) { SolverType solv_type = ipopt; double tol = 1e-6; string mehrotra = "no"; -// string fname = "../data_sets/Power/nesta_case3_lmbd.m"; - string fname = "../nesta-0.7.0/opf/api/nesta_case24_ieee_rts__api.m"; + string fname = "../data_sets/Power/nesta_case3_lmbd.m"; +// string fname = "../nesta-0.7.0/opf/api/nesta_case24_ieee_rts__api.m"; // fname = "/Users/hlh/Dropbox/Work/Dev/power_models/data/nesta_api/nesta_case30_fsr__api.m"; // create a OptionParser with options @@ -337,7 +337,7 @@ int main (int argc, char * argv[]) { int iter = 0, hdim_cuts = 0, cuts_added = 1; // while(unchanged < 3) { - unsigned nr_threads = 4, nb_bags = bags.size(); + unsigned nr_threads = 10, nb_bags = bags.size(); while(cuts_added > 0) { cuts_added = 0; vector> w_hat_vec; @@ -479,6 +479,7 @@ int main (int argc, char * argv[]) { 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++; } } diff --git a/src/constraint.cpp b/src/constraint.cpp index 8cac9a535..61ef5d303 100755 --- a/src/constraint.cpp +++ b/src/constraint.cpp @@ -55,7 +55,8 @@ int Constraint::get_type() const{ }; bool Constraint::is_active(unsigned inst) const{ - return fabs(_dual[inst]) > EPS; + return fabs(get_val(inst) - _rhs) < EPS; +// return fabs(_dual[inst]) > EPS; } /* Operators */ diff --git a/src/solver.cpp b/src/solver.cpp index d95dbf5fc..001a8ea9b 100755 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -171,7 +171,7 @@ int solver::run(int print_level, bool relax, double tol, const string& lin_solve SmartPtr tmp = new IpoptProgram(_model); status = iapp->OptimizeTNLP(tmp); - violated_constraints = _model->has_violated_constraints(tol); +// violated_constraints = _model->has_violated_constraints(tol); if (status == Solve_Succeeded) { // Retrieve some statistics about the solve From 5a74beba86280fd742aefd5c2ab2d1a3a8c3a64d Mon Sep 17 00:00:00 2001 From: Hassan Date: Thu, 25 Jan 2018 22:05:27 -0700 Subject: [PATCH 4/5] adding upperbound into sdpopf and removing psd bags from bag vector --- examples/MINLP/Power/PowerNet.cpp | 198 ++++++++++++++++++++ examples/MINLP/Power/PowerNet.h | 5 +- examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp | 55 +++--- src/CplexProgram.cpp | 9 +- 4 files changed, 230 insertions(+), 37 deletions(-) diff --git a/examples/MINLP/Power/PowerNet.cpp b/examples/MINLP/Power/PowerNet.cpp index ab0417de5..7afe99431 100755 --- a/examples/MINLP/Power/PowerNet.cpp +++ b/examples/MINLP/Power/PowerNet.cpp @@ -22,6 +22,7 @@ #include #include #include +#include //#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 Pg("Pg", pg_min, pg_max); + var Qg ("Qg", qg_min, qg_max); + ACOPF.add_var(Pg.in(gens)); + ACOPF.add_var(Qg.in(gens)); + + /* Power flow variables */ + var Pf_from("Pf_from", S_max); + var Qf_from("Qf_from", S_max); + var Pf_to("Pf_to", S_max); + var 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 theta("theta"); + var v("|V|", v_min, v_max); + // var vr("vr"); + // var vi("vi"); + var vr("vr", v_max); + var 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; +} + + diff --git a/examples/MINLP/Power/PowerNet.h b/examples/MINLP/Power/PowerNet.h index a1ca745b6..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 */ @@ -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 diff --git a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp index d13247cc3..47b6c2144 100644 --- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp +++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp @@ -22,14 +22,7 @@ using namespace gravity; void nearest_point(int i, int j, vector>& w_hat, const vector& bags){ for (unsigned index = i; index ("psd"); - DebugOff("Bag " << to_string(index) << " is PSD"< 0) { + 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; - w_hat_vec.resize(nb_bags); + 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(bags))); + 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 = 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])); - if (w_hat._name.compare("psd")==0) { - continue; - } R_star.resize(numcuts+1); R_star[numcuts] = param<>("R_star"+to_string(numcuts)); I_star.resize(numcuts+1); @@ -389,7 +389,6 @@ int main (int argc, char * argv[]) { W_hat.resize(numcuts+1); W_hat[numcuts] = param<>("W_hat"+to_string(numcuts)); -// what = b.nfp(); node_pairs b_pairs("node_pairs"); Constraint sdpcut("sdpcut_" + to_string(numcuts)); Node *ni; @@ -429,10 +428,6 @@ 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 = product(R_star.in(b_pairs._keys) - R_hat.in(b_pairs._keys),(R_Wij.in(b_pairs._keys) - R_hat.in(b_pairs._keys))); - // lin += product(I_star.in(b_pairs._keys) - I_hat.in(b_pairs._keys),(Im_Wij.in(b_pairs._keys) - I_hat.in(b_pairs._keys))); - // lin += product(W_star.in(b._nodes) - W_hat.in(b._nodes),(Wii.in(b._nodes) - W_hat.in(b._nodes))); - // lin.print(); // lin.print_expanded(); SDP.add_constraint(lin <= 0); @@ -440,17 +435,7 @@ int main (int argc, char * argv[]) { } if(!bus_pairs_sdp._keys.empty()) { -// cout << "Adding SOC indices: {"; -// for(auto& bp: bus_pairs_sdp._keys) { -// cout << bp->_name << ", "; -// } -// cout << "}" << endl; - -// Constraint SOC_im("SOC_im"+to_string(numcuts)); -// SOC_im = power(R_Wij, 2) + power(Im_Wij, 2) - Wii.from() * Wii.to(); -// SDP.add_constraint(SOC_im.in(bus_pairs_sdp) <= 0); SDP.add_indices("SOC",bus_pairs_sdp); -// SDP.get_constraint("SOC")->print_expanded(); bus_pairs_sdp.clear(); } @@ -460,6 +445,7 @@ int main (int argc, char * argv[]) { } 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); @@ -494,6 +480,9 @@ int main (int argc, char * argv[]) { 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 <_obj_val = cplex.getObjValue(); for (auto i = 0; i < _cplex_vars.size(); i++) { - IloNumArray vals(*_cplex_env); - cplex.getValues(vals,_cplex_vars[i]); for (auto j = 0; j < _model->_vars[i]->get_dim(); j++) { - poly_set_val(j, vals[j], _model->_vars[i]); + if(cplex.isExtracted(_cplex_vars[i][j])){ + poly_set_val(j, cplex.getValue(_cplex_vars[i][j]), _model->_vars[i]); + } + else { + poly_set_val(j, 0, _model->_vars[i]); + } } } } From 751f415a2b7873726a2747486fe8558cd63b6830 Mon Sep 17 00:00:00 2001 From: Hassan Date: Wed, 31 Jan 2018 10:17:35 -0500 Subject: [PATCH 5/5] adding Real sets, dropping the ^ operator for dimension setting --- examples/Gravity101_main.cpp | 2 +- examples/MILP/MinKPart/Minkmodel.cpp | 12 ++-- examples/MILP/StableSet/Stable_Set_main.cpp | 9 +-- examples/MINLP/Power/ACUC/ACUC_main.cpp | 28 +++++---- examples/MINLP/Power/SDPOPF/Bag.cpp | 3 +- examples/NLP/Ising/Ising_main.cpp | 5 +- include/gravity/model.h | 4 +- include/gravity/param.h | 31 ++++++++-- include/gravity/types.h | 2 + include/gravity/var.h | 66 ++++++++++++++++++++- src/var.cpp | 5 ++ 11 files changed, 135 insertions(+), 32 deletions(-) diff --git a/examples/Gravity101_main.cpp b/examples/Gravity101_main.cpp index 208d3f9b6..466f88c36 100755 --- a/examples/Gravity101_main.cpp +++ b/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; diff --git a/examples/MILP/MinKPart/Minkmodel.cpp b/examples/MILP/MinKPart/Minkmodel.cpp index 09f636bf6..87fae883d 100755 --- a/examples/MILP/MinKPart/Minkmodel.cpp +++ b/examples/MILP/MinKPart/Minkmodel.cpp @@ -73,7 +73,8 @@ void Minkmodel::build() { void Minkmodel::reset() {}; void Minkmodel::add_vars_origin() { var 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 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 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 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; diff --git a/examples/MILP/StableSet/Stable_Set_main.cpp b/examples/MILP/StableSet/Stable_Set_main.cpp index 90cb29a06..48f3c40ef 100755 --- a/examples/MILP/StableSet/Stable_Set_main.cpp +++ b/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 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 Xii("Xii", 0, 1); var 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); diff --git a/examples/MINLP/Power/ACUC/ACUC_main.cpp b/examples/MINLP/Power/ACUC/ACUC_main.cpp index 8a0af109e..f499287c5 100644 --- a/examples/MINLP/Power/ACUC/ACUC_main.cpp +++ b/examples/MINLP/Power/ACUC/ACUC_main.cpp @@ -71,18 +71,20 @@ int main (int argc, const char * argv[]) // POWER GENERATION var Pg("Pg", grid->pg_min.in(grid->gens, T), grid->pg_max.in(grid->gens, T)); var 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 Pf_from("Pf_from", grid->S_max.in(grid->arcs, T)); var Qf_from("Qf_from", grid->S_max.in(grid->arcs, T)); var Pf_to("Pf_to", grid->S_max.in(grid->arcs, T)); var 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 R_Wij("R_Wij", grid->wr_min.in(bus_pairs, T), grid->wr_max.in(bus_pairs, T)); // real part of Wij @@ -91,9 +93,11 @@ 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); @@ -101,9 +105,9 @@ int main (int argc, const char * argv[]) var On_off("On_off"); var Start_up("Start_up", 0, 1); var 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; diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp index 2b90bbce8..9ca9da9de 100644 --- a/examples/MINLP/Power/SDPOPF/Bag.cpp +++ b/examples/MINLP/Power/SDPOPF/Bag.cpp @@ -355,7 +355,8 @@ param Bag::nfp(){ var 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 w("w", _wmin, _wmax); NPP.add_var(w.in(_indices)); diff --git a/examples/NLP/Ising/Ising_main.cpp b/examples/NLP/Ising/Ising_main.cpp index 47166ae09..6f46a69ca 100644 --- a/examples/NLP/Ising/Ising_main.cpp +++ b/examples/NLP/Ising/Ising_main.cpp @@ -116,8 +116,9 @@ void solve_spin(unsigned spin1, unsigned spin2, int log_lev=0, bool relax=false, Model Ising("Ising Model"); /** Variables */ var x("x"), z("z", pos_), obj("obj"); - Ising.add_var(x^nb_spins); - Ising.add_var(z^nb_spins); + auto Rn = R(nb_spins); + Ising.add_var(x.in(Rn)); + Ising.add_var(z.in(Rn)); Ising.add_var(obj); Ising.min(obj); diff --git a/include/gravity/model.h b/include/gravity/model.h index 61636044f..003d202a0 100755 --- a/include/gravity/model.h +++ b/include/gravity/model.h @@ -144,7 +144,7 @@ namespace gravity { v.set_vec_id(_vars.size()); param_* newv; if (v._dim[0]==0) { - newv = (param_*)copy(v^1); + newv = (param_*)copy(v.in(R(1))); } else { newv = (param_*)copy(v); @@ -186,7 +186,7 @@ namespace gravity { v.set_vec_id(_vars.size()); param_* newv; if (v._dim[0]==0) { - newv = (param_*)copy(v^1); + newv = (param_*)copy(v.in(R(1))); } else { newv = (param_*)copy(v); diff --git a/include/gravity/param.h b/include/gravity/param.h index ed36a7a7b..868a29e8f 100755 --- a/include/gravity/param.h +++ b/include/gravity/param.h @@ -508,6 +508,29 @@ namespace gravity { /* Modifiers */ + void set_size(vector dims){ + if (dims.size()==1) { + set_size(dims[0]); + } + else if (dims.size()==2){ + set_size(dims[0],dims[1]); + } + else { + _dim = dims; + //TODO allocate val + } + } + + void set_size(size_t s1, size_t s2, type val = 0) { + _is_matrix = true; + _dim.resize(2); + _dim[0] = s1; + _dim[1] = s2; + auto index = _dim[1]*s1-1+s2-1; + _val->resize(index+1); + _val->at(index) = val; + }; + void set_size(size_t s, type val = 0) { _val->resize(s, val); _dim[0] = s; @@ -641,10 +664,10 @@ namespace gravity { //return (get_name()==p.get_name() && _type==p._type && _intype==p._intype && _dim==p._dim && _indices==p._indices && _val==p._val); } - param& operator^(size_t d) { - set_size(d); - return *this; - } +// param& operator^(size_t d) { +// set_size(d); +// return *this; +// } void initialize_all(type v) { for (int i = 0; i<_val->size(); i++) { diff --git a/include/gravity/types.h b/include/gravity/types.h index 9ca6afbe0..ab7e79cf6 100755 --- a/include/gravity/types.h +++ b/include/gravity/types.h @@ -35,6 +35,8 @@ namespace gravity{ typedef enum { minimize, maximize } ObjectiveType; typedef enum { id_, number_, plus_, minus_, product_, div_, power_, cos_, sin_, sqrt_, exp_, log_} OperatorType; /* Operation type in the expression tree */ + typedef enum { R_, R_p_, C_, C_p_} SpaceType; /* Operation type in the expression tree */ + typedef enum { ordered_pairs_, unordered_ } SetType; // typedef enum { vec_=0, in_ordered_pairs_=1, from_ordered_pairs_=2, to_ordered_pairs_=3, in_arcs_=4, from_arcs_=5, to_arcs_=6, in_nodes_=7, in_set_=8, mask_=9, in_bags_=10, time_expand_ = 11, in_set_at_} IndexType; /* Index type */ diff --git a/include/gravity/var.h b/include/gravity/var.h index 3cc1441fb..a31a70be2 100644 --- a/include/gravity/var.h +++ b/include/gravity/var.h @@ -22,7 +22,62 @@ using namespace std; namespace gravity { + + class func_; + class space{ + public: + SpaceType _type; + vector _dim; + }; + + class R: public space{ + public: + R(){}; + template + R(size_t t1, Args&&... args) { + _type = R_; + list dims = {forward(args)...}; + dims.push_front(t1); + size_t size = dims.size(); + _dim.resize(size); + auto it = dims.begin(); + size_t index = 0; + while (it!=dims.end()) { + _dim[index++] = *it++; + } + } + + R operator^(size_t n){return R(n);}; + + }; + + class R_plus: public space{ + public: + R_plus(){}; + template + R_plus(size_t t1, Args&&... args) { + _type = R_p_; + list dims = {forward(args)...}; + dims.push_front(t1); + size_t size = dims.size(); + _dim.resize(size); + auto it = dims.begin(); + size_t index = 0; + while (it!=dims.end()) { + _dim[index++] = *it++; + } + } + + R_plus operator^(size_t n){return R_plus(n);}; + + }; + + class C: public space{ + public: + /* TODO */ + }; + /** Backbone class for parameter */ class var_ { public: @@ -356,6 +411,7 @@ class var: public param, public var_ { Sign get_sign(int idx = 0) const; /* Modifiers */ + void set_size(vector); void set_size(size_t s, type val = 0); void add_bounds(type lb, type ub); @@ -383,10 +439,16 @@ class var: public param, public var_ { bool operator==(const var& v) const; bool operator!=(const var& v) const; - var& operator^(size_t d) { - set_size(d); + + var& in(const space& s){ + set_size(s._dim); return *this; } + +// var& operator^(size_t d) { +// set_size(d); +// return *this; +// } var tr() const { auto v = var(*this); diff --git a/src/var.cpp b/src/var.cpp index c6884d100..019a69cc8 100755 --- a/src/var.cpp +++ b/src/var.cpp @@ -139,6 +139,11 @@ template var& var::operator=(var&& v) { }; /* Modifiers */ + +template void var::set_size(vector dims) { + param::set_size(dims); +}; + template void var::set_size(size_t s, type val) { param::set_size(s,val); };