Skip to content

Commit

Permalink
Merge pull request #27 from Allinsights/SDP_cuts
Browse files Browse the repository at this point in the history 
Sdp cuts, fixing bug in param dim
  • Loading branch information
@hhijazi
hhijazi authored Feb 19, 2018
2 parents 6f4e034 + f87a112 commit 30f5afc
Show file tree
Hide file tree
Showing 8 changed files with 518 additions and 157 deletions.
135 changes: 4 additions & 131 deletions examples/MILP/StableSet/Stable_Set_main.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,35 +1,27 @@
//
//
// Stable_Set.cpp
// Gravity
//
// Created by Hijazi, Hassan (Data61, Canberra City) on 6/12/17.
// Created by Hijazi, Hassan on 6/12/17.
//
//
#include <stdio.h>
#include <iostream>
#include <string>
#include <stdio.h>
#include <cstring>
#include <fstream>
#include <gravity/Net.h>
#include <gravity/model.h>
#include <gravity/solver.h>
#include <stdio.h>
#include <stdlib.h>
#include <gravity/solver.h>

using namespace std;
using namespace gravity;


int main (int argc, const char * argv[])
{
// Start Timers
std::cout << "WELCOME, THIS IS AN IMPLEMENTATION OF THE STABLE SET PROBLEM AND SOME OF ITS RELAXATIONS IN GRAVITY\n";

unsigned i, j;

Net graph;
// const char* fname = "../../data_sets/stable_set/toy.txt";
const char* fname = "../../data_sets/stable_set/p.3n150.txt";
graph.read_adjacency_matrix(fname);

Expand All @@ -53,129 +45,10 @@ int main (int argc, const char * argv[])
c <= 1;
c.print();
model.add_constraint(c.in(graph.arcs));

/** Solver **/
SolverType stype = gurobi;
solver s(model,stype);
double wall0 = get_wall_time();
double cpu0 = get_cpu_time();
cout << "Running the IP model\n";
s.run();
double wall1 = get_wall_time();
double cpu1 = get_cpu_time();
cout << "Done running the IP model\n";
cout << "Wall clock computing time = " << wall1 - wall0 << "\n";
cout << "CPU computing time = " << cpu1 - cpu0 << "\n";

/* Schriver's SDP relaxation for the stable set problem */
Model SDP;
/* Variable declaration */
var<double> Xii("Xii", 0, 1);
var<double> Xij("Xij", 0, 1);
SDP.add_var(Xii.in(R(n))); /*< Diagonal entries of the matrix */
SDP.add_var(Xij.in(R(n*(n-1)/2))); /*< Lower left triangular part of the matrix excluding the diagonal*/

/* Constraints declaration */
ordered_pairs indices(1, n);
Constraint SOCP("SOCP");
SOCP = power(Xij, 2) - Xii.from()*Xii.to() ;
SDP.add_constraint(SOCP.in(indices._keys) <= 0);
Constraint diag("diag");
diag = sum(Xii);
SDP.add_constraint(diag == 1); // diagonal sum is 1

Constraint zeros("zeros");
zeros = Xij.in(graph.arcs);
SDP.add_constraint(zeros == 0); // zero elements

auto Xij_ = Xij.pairs_in(graph._bags, 3);
auto Xii_ = Xii.in(graph._bags, 3);
Constraint SDP3("SDP_3D");
SDP3 = -2*Xij_[0]*Xij_[1]*Xij_[2];
SDP3 -= Xii_[0]*Xii_[1]*Xii_[2];
SDP3 += power(Xij_[0],2)*Xii_[2];
SDP3 += power(Xij_[2],2)*Xii_[1];
SDP3 += power(Xij_[1],2)*Xii_[0];
// SDP3.print();
// SDP.add_constraint(SDP3);
// set<tuple<int,int,int>> ids;
// for (i = 0; i < complement_graph._bags.size(); i++){
//// for (i = 0; i < 1; i++){
// auto bag = complement_graph._bags.at(i);
// if (bag.size() < 3) {
// continue;
// }
//// for (int j = 0; j < bag.size()-2; j++){
// for (int j = 0; j < 1; j++){
// i1 = bag[j]->ID;
// i2 = bag[j+1]->ID;
// i3 = bag[j+2]->ID;
// assert(i2>i1 && i3>i2);
// if(ids.count(make_tuple(i1, i2, i3))==0){
// ids.insert(make_tuple(i1, i2, i3));
// }
// else {
// continue;
// }
// Constraint SDP3("SDP3("+to_string(i1)+","+to_string(i2)+","+to_string(i3)+")");
// SDP3 = -2*Xij(i1,i2)*Xij(i2,i3)*Xij(i1,i3);
// SDP3 -= Xii(i1)*Xii(i2)*Xii(i3);
// SDP3 += power(Xij(i1,i2),2)*Xii(i3);
// SDP3 += power(Xij(i1,i3),2)*Xii(i2);
// SDP3 += power(Xij(i2,i3),2)*Xii(i1);
// SDP3.print();
// SDP.add_constraint(SDP3);
// }
// }

/* Objective declaration */
SDP.max(2*sum(Xij) + sum(Xii));


// solver s1(SDP,ipopt);
solver s1(SDP,gurobi);

wall0 = get_wall_time();
cpu0 = get_cpu_time();
cout << "Running the SDP relaxation\n";
s1.run();
wall1 = get_wall_time();
cpu1 = get_cpu_time();
cout << "Done running the SDP relaxation\n";
cout << "\nWall clock computing time = " << wall1 - wall0 << "\n";
cout << "CPU computing time = " << cpu1 - cpu0 << "\n";
return 0;

/* Outer-Approximation of SOCP and 3d-SDP cuts in Shriver's SDP relaxation for the stable set problem */
Model OA;
OA.add_var(Xii);
OA.add_var(Xij);
for (auto &cs_p: SDP._cons) {
if (!cs_p.second->is_linear() && cs_p.second->is_active()){ //Active nonlinear constraint
DebugOff("Active constraint:" << cs_p.second->to_str() << endl);
Constraint oa("OA_"+cs_p.second->get_name());
oa = cs_p.second->get_outer_app();
// oa.print();
OA.add_constraint(oa);
}
}
// OA.add_constraint(diag=1);
for(auto a: graph.arcs){
i = (a->_src)->_id;
j = (a->_dest)->_id;
Constraint zeros("zeros("+to_string(i)+","+to_string(j)+")");
zeros = Xij(i,j);
OA.add_constraint(zeros==0);
}
OA.max(2*sum(Xij) + sum(Xii));
solver s2(OA,cplex);
wall0 = get_wall_time();
cpu0 = get_cpu_time();
cout << "Running the OA-SDP relaxation\n";
s2.run();
wall1 = get_wall_time();
cpu1 = get_cpu_time();
cout << "Done running the OA-SDP relaxation\n";
cout << "\nWall clock computing time = " << wall1 - wall0 << "\n";
cout << "CPU computing time = " << cpu1 - cpu0 << "\n";
return 0;
};
113 changes: 102 additions & 11 deletions examples/MINLP/Power/SDPOPF/Bag.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
#include <gravity/param.h>
#include "Bag.h"

#define FLAT

Bag::Bag(int id, const PowerNet& grid, vector<Node*> nodes):_id(id),_grid((PowerNet*)&grid),_nodes(nodes) {
Arc* aij = NULL;
_wstarp.set_name("wstar");
Expand Down Expand Up @@ -45,6 +47,16 @@ Bag::Bag(int id, const PowerNet& grid, vector<Node*> nodes):_id(id),_grid((Power
_wmax.set_val(namewr,_grid->wr_max(namepair).eval());
_wmax.set_val(namewi,_grid->wi_max(namepair).eval());

_Wmin.set_val(to_string(i)+","+to_string(j),_grid->wr_min(namepair).eval());
_Wmin.set_val(to_string(i+_nodes.size())+","+to_string(j+_nodes.size()),_grid->wr_min(namepair).eval());
_Wmin.set_val(to_string(i)+","+to_string(j+_nodes.size()),_grid->wi_min(namepair).eval());
_Wmin.set_val(to_string(j)+","+to_string(i+_nodes.size()),_grid->wi_min(namepair).eval());

_Wmax.set_val(to_string(i)+","+to_string(j),_grid->wr_max(namepair).eval());
_Wmax.set_val(to_string(i+_nodes.size())+","+to_string(j+_nodes.size()),_grid->wr_max(namepair).eval());
_Wmax.set_val(to_string(i)+","+to_string(j+_nodes.size()),_grid->wi_max(namepair).eval());
_Wmax.set_val(to_string(j)+","+to_string(i+_nodes.size()),_grid->wi_max(namepair).eval());

// _wstarp.set_val(namewr,((Line*)aij)->wr);
// if(!reversed) _wstarp.set_val(namewi,((Line*)aij)->wi);
// else _wstarp.set_val(namewi,-((Line*)aij)->wi);
Expand All @@ -61,6 +73,15 @@ Bag::Bag(int id, const PowerNet& grid, vector<Node*> nodes):_id(id),_grid((Power
_indices.push_back(index_(namew));
_wmin.set_val(namew,_grid->w_min(_nodes[i]->_name).eval());
_wmax.set_val(namew,_grid->w_max(_nodes[i]->_name).eval());

_Wmin.set_val(to_string(i)+","+to_string(i),_grid->w_min(_nodes[i]->_name).eval());
_Wmin.set_val(to_string(i+_nodes.size())+","+to_string(i+_nodes.size()),_grid->w_min(_nodes[i]->_name).eval());
_Wmin.set_val(to_string(i)+","+to_string(i+_nodes.size()),0);

_Wmax.set_val(to_string(i)+","+to_string(i),_grid->w_max(_nodes[i]->_name).eval());
_Wmax.set_val(to_string(i+_nodes.size())+","+to_string(i+_nodes.size()),_grid->w_max(_nodes[i]->_name).eval());
_Wmax.set_val(to_string(i)+","+to_string(i+_nodes.size()),0);

// _wstarp.set_val(namew,((Bus*)_nodes[i])->w);
// cout << "\n" << namew << " = " << ((Bus*)_nodes[i])->w;

Expand Down Expand Up @@ -108,8 +129,13 @@ param<double> Bag::fill_wstar(){

_W_star.set_val(to_string(i)+","+to_string(j),aij->wr);
_W_star.set_val(to_string(i+_nodes.size())+","+to_string(j+_nodes.size()),aij->wr);
_W_star.set_val(to_string(i)+","+to_string(j+_nodes.size()),aij->wi);
_W_star.set_val(to_string(j)+","+to_string(i+_nodes.size()),aij->wi);
if(!reversed) {
_W_star.set_val(to_string(i) + "," + to_string(j + _nodes.size()), aij->wi);
_W_star.set_val(to_string(j) + "," + to_string(i + _nodes.size()), -aij->wi);
}else{
_W_star.set_val(to_string(i) + "," + to_string(j + _nodes.size()), -aij->wi);
_W_star.set_val(to_string(j) + "," + to_string(i + _nodes.size()), aij->wi);
}
}
}
for(int i = 0; i < _nodes.size(); i++){
Expand All @@ -120,7 +146,7 @@ param<double> Bag::fill_wstar(){
_W_star.set_val(to_string(i+_nodes.size())+","+to_string(i+_nodes.size()),((Bus*)_nodes[i])->w);
_W_star.set_val(to_string(i)+","+to_string(i+_nodes.size()),0.0);
}
return W_star;//todo: use this
return W_star;
}

bool Bag::add_lines(){
Expand Down Expand Up @@ -334,6 +360,16 @@ void Bag::update_PSD(){
_is_psd = is_PSD();
}

vector<index_> triang_indices(int n) {
vector<index_> res;
for(int i = 0; i < n; i++) {
for(int j = i; j < n; j++) {
res.push_back(index_(to_string(i)+","+to_string(j)));
}
}
return res;
}

param<double> Bag::nfp(){
param<double> what;
fill_wstar();
Expand All @@ -346,16 +382,24 @@ param<double> Bag::nfp(){
int n = _nodes.size();
DebugOff("\nn = " << n);

// _wstarp.print(true);

#ifndef FLAT
sdpvar<double> W("W");
NPP.add_var(W^(2*n));

// var<double> W("W");
// W._psd = true;
// NPP.add_var(W ^ (n*(2*n+1)));
#else
var<double> W("W");
W._psd = true;
NPP.add_var(W.in(R(2*n,2*n))); // note: number of instances in apper triangle is n*(2*n+1)
#endif

var<double> z("z");
z.in_q_cone();
#ifndef FLAT
auto Rn = R(n*n+1);
#else
auto Rn = R(2*n*2*n+1);
#endif
NPP.add_var(z.in(Rn));

var<double> w("w", _wmin, _wmax);
Expand All @@ -377,14 +421,40 @@ param<double> Bag::nfp(){

/* w*-w = z */
Constraint svec("svec");

#ifndef FLAT
svec = _wstarp - w - z;
NPP.add_constraint(svec.in(_indices)==0);

//todo: flatten(W*-W) = z;

#else
auto idxs = triang_indices(2*n);
svec = _W_star - W - z;
NPP.add_constraint(svec.in(idxs)==0);
// svec.print_expanded();

// vector<index_> zero_idxs;
// for(int i = 0; i < n; i++) {
// zero_idxs.push_back(index_(to_string(i)+","+to_string(i+n)));
// }
// Constraint zeros("zeros");
// zeros = W.in(zero_idxs);
// NPP.add_constraint(zeros==0);
// zeros.print_expanded();

Constraint W_lb("W_lb");
W_lb = _Wmin - W;
NPP.add_constraint(W_lb.in(idxs) <= 0);
// W_lb.print_expanded();

Constraint W_ub("W_ub");
W_ub = _Wmax - W;
NPP.add_constraint(W_ub.in(idxs) >= 0);
// W_ub.print_expanded();

#endif

string namew, namewr, namewi;

#ifndef FLAT
/* matrix structure */
for(int i = 0; i < n; i++){
for(int j = i; j < n; j++){
Expand Down Expand Up @@ -428,14 +498,15 @@ param<double> Bag::nfp(){
}
}
}
#endif

solver s(NPP,Mosek);
s.run(0,0);

// z.print(); cout << "\n";
// W.print(true); cout << "\n";
// w.print(); cout << "\n";

#ifndef FLAT
for(int i = 0; i < n; i++){
for(int j = i; j < n; j++){
if(i==j){
Expand All @@ -453,7 +524,27 @@ param<double> Bag::nfp(){
}
}
}
#else
for(int i = 0; i < n; i++){
for(int j = i; j < n; j++){
if(i==j){
namew = "w(" + _nodes[i]->_name + ")";
what.set_val(namew,W(to_string(i)+","+to_string(j)).eval());
}
else {
namewr = "wr(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")";
namewi = "wi(" + _nodes[i]->_name + "," + _nodes[j]->_name + ")";
what.set_val(namewr, W(to_string(i)+","+to_string(j)).eval());
if(_grid->get_directed_arc(_nodes[i]->_name,_nodes[j]->_name)!=nullptr)
what.set_val(namewi, W(to_string(i)+","+to_string(j+_nodes.size())).eval());
else
what.set_val(namewi, -W(to_string(i)+","+to_string(j+_nodes.size())).eval());
}
}
}
#endif
what.set_name("w_hat");
// what.print(true);
// exit(0);
return what;
}
3 changes: 3 additions & 0 deletions examples/MINLP/Power/SDPOPF/Bag.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,9 @@ class Bag{
param<double> _wmin;
param<double> _wmax;

param<double> _Wmin;
param<double> _Wmax;

// gravity::node_pairs _bus_pairs;

gravity::Model _model;
Expand Down
Loading

0 comments on commit 30f5afc

Please sign in to comment.