diff --git a/examples/MINLP/Power/ACOPF/ACOPF_main.cpp b/examples/MINLP/Power/ACOPF/ACOPF_main.cpp
index f6bc0116d..7644a37a4 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) {
@@ -81,8 +82,8 @@ 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);
@@ -90,27 +91,26 @@ int main (int argc, char * argv[])
     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);
     }
 
@@ -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");
@@ -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 */
@@ -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;
@@ -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;
@@ -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) {
@@ -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);
     }
 
     
@@ -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);
@@ -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();
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;
 
diff --git a/examples/MINLP/Power/SDPOPF/Bag.cpp b/examples/MINLP/Power/SDPOPF/Bag.cpp
index e19b01fba..99d9fa1a5 100644
--- a/examples/MINLP/Power/SDPOPF/Bag.cpp
+++ b/examples/MINLP/Power/SDPOPF/Bag.cpp
@@ -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;
@@ -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;
@@ -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;
@@ -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) {
@@ -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;
     }
 }
@@ -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;
@@ -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;
 }
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<double> nfp();
+    param<double> 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 bedc600af..21e31acc7 100644
--- a/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp
+++ b/examples/MINLP/Power/SDPOPF/SDPOPF_main.cpp
@@ -23,6 +23,7 @@ void nearest_point(int i, int j, vector<param<>>& w_hat, const vector<Bag>& bags
     for (unsigned index = i; index <j; index++) {
         auto b = bags[index];
         w_hat[index] = b.nfp();
+//        w_hat[index] = b.nfp1();
     }
 }
 
@@ -186,7 +187,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 */
@@ -293,7 +294,7 @@ int main_new (int argc, char * argv[]) {
     //    double fp_tol = 0.0001;
     DebugOff("\nNumber of bags = " << bags.size());
     
-    vector<param<double>> R_star, I_star, W_star;
+//    vector<param<double>> R_star, I_star, W_star;
     vector<param<double>> R_diff, I_diff, W_diff;
     vector<param<double>> R_hat, I_hat, W_hat;
     
@@ -338,12 +339,12 @@ int main_new (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);
@@ -368,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);
@@ -385,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)));
@@ -692,33 +695,29 @@ int main (int argc, char * argv[]) {
     }
     for(auto& node: grid.nodes) ((Bus*)node)->w = Wii(node->_name).eval();
 
-
-//    param<double> 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<param<double>> R_star, I_star, W_star;
+//    vector<param<double>> R_star, I_star, W_star;
     vector<param<double>> R_diff, I_diff, W_diff;
     vector<param<double>> 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) << "%."<<endl);
         cuts_added = 0;
         vector<param<>> 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<Bag> violated_bags;
@@ -746,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);
@@ -776,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);
@@ -785,31 +784,35 @@ 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());
                         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++;
         }
-        
+
         if(!bus_pairs_sdp._keys.empty()) {
             SDP.add_indices("SOC",bus_pairs_sdp);
             bus_pairs_sdp.clear();
@@ -819,7 +822,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);
diff --git a/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp b/examples/MINLP/Power/SOCOPF/SOCOPF_main.cpp
index b4979c2a7..9eb034c49 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;
@@ -151,26 +151,24 @@ 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");
@@ -178,14 +176,14 @@ int main (int argc, char * argv[])
     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);
+    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_lazy(LNC2.in(bus_pairs) >= 0);
+    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 <typename type>
+        void add(var<type>& v){//Add variables by copy
+            add_var(v);
+        }
+        
         
         template <typename type>
         void add_var(var<type>&& v){//Add variables by copy
@@ -200,6 +205,11 @@ namespace gravity {
             }
         };
         
+        template <typename type>
+        void add(var<type>&& v){
+            add_var(move(v));
+        }
+        
         void del_var(const param_& v);
         
         
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 {
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);
 }