Skip to content

Commit

Permalink
NZ model does not work
Browse files Browse the repository at this point in the history
  • Loading branch information
@hhijazi
hhijazi committed Dec 16, 2023
1 parent 6335bbd commit af516bf
Show file tree
Hide file tree
Showing 4 changed files with 108 additions and 20 deletions.
83 changes: 73 additions & 10 deletions examples/Optimization/MINLP/LABS/LABS.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ int main(int argc, char * argv[]){
y.set_lb("0",1);

// M_obj.add(c.in(c_ids));
indices pairs("pairs"), quad_terms("quad_terms"), multi_terms("multi_terms"), multi_quad_terms("multi_quad_terms");
indices pairs("pairs"), quad_terms("quad_terms"), multi_terms("multi_terms"), multi_lin_terms("multi_lin_terms"), multi_quad_terms("multi_quad_terms");
func<> obj;
for (int k = 1; k<=n-1; k++) {
func<> cterm;
Expand All @@ -91,6 +91,7 @@ int main(int argc, char * argv[]){
}
quad_terms = pairs.subset();
multi_quad_terms = pairs.subset();
multi_lin_terms = s_ids.subset();
string pair_idx;
int nb_row = 0;
for(auto p: *obj._pterms)
Expand All @@ -101,6 +102,9 @@ int main(int argc, char * argv[]){
bool existing_pair = false;
auto idx1 = it->first->_indices->_ids->front().at(0);
mult_idx += to_string(idx1);
multi_lin_terms.add_in_row(nb_row, to_string(idx1));
if(next(it)!= p.second._l->end())
mult_idx += ",";
if(tabu.count(idx1)>0)
continue;
for(auto it2 = next(it); it2 != p.second._l->end(); it2++) {
Expand Down Expand Up @@ -133,19 +137,78 @@ int main(int argc, char * argv[]){
multi_p1._ids->at(0).push_back(multi_quad_terms._ids->at(i).at(0));
multi_p2._ids->at(0).push_back(multi_quad_terms._ids->at(i).at(1));
}
var<int> p("p", -1, 1);
var<int> pp("pp", -1, 1);
var<> p("p", -1, 1);
var<> pp("pp", -1, 1);
var<> nz("nz", pos_);

M_obj.add(nz.in(multi_terms));
M_obj.add(p.in(pairs));
M_obj.add(pp.in(multi_terms));

Constraint<> p_def("p_def");
p_def = p - (2*y.from(pairs) - 1)*(2*y.to(pairs) - 1);
M_obj.add(p_def.in(pairs) == 0);

Constraint<> pp_def("pp_def");
pp_def = pp - p.in(multi_p1)*p.in(multi_p2);
M_obj.add(pp_def.in(multi_terms) == 0);
// Constraint<> pp_def("pp_def");
// pp_def = pp - p.in(multi_p1)*p.in(multi_p2);
// M_obj.add(pp_def.in(multi_terms) == 0);


Constraint<> nz_def("nz_def");
nz_def = nz - (4 - y.in(multi_lin_terms));
M_obj.add(nz_def.in(multi_terms) == 0);

var<int> nz0("nz0", 0,1);
var<int> nz1("nz1", 0,1);
var<int> nz2("nz2", 0,1);
var<int> nz3("nz3", 0,1);
var<int> nz4("nz4", 0,1);
M_obj.add(nz0.in(multi_terms), nz1.in(multi_terms), nz2.in(multi_terms), nz3.in(multi_terms), nz4.in(multi_terms));

Constraint<> nz_one("nz_one");
nz_one = nz0 + nz1 + nz2 + nz3 + nz4 - 1;
M_obj.add(nz_one.in(multi_terms) == 0);

Constraint<> nz0_def("nz0_def");
nz0_def = nz - 0;
M_obj.add_on_off(nz0_def.in(multi_terms) == 0, nz0, true);

Constraint<> nz1_def("nz1_def");
nz1_def = nz - 1;
M_obj.add_on_off(nz1_def.in(multi_terms) == 0, nz1, true);

Constraint<> nz2_def("nz2_def");
nz2_def = nz - 2;
M_obj.add_on_off(nz2_def.in(multi_terms) == 0, nz2, true);

Constraint<> nz3_def("nz3_def");
nz3_def = nz - 3;
M_obj.add_on_off(nz3_def.in(multi_terms) == 0, nz3, true);

Constraint<> nz4_def("nz4_def");
nz4_def = nz - 4;
M_obj.add_on_off(nz4_def.in(multi_terms) == 0, nz4, true);


Constraint<> pp_nz0("pp_nz0");
pp_nz0 = pp - 1;
M_obj.add_on_off(pp_nz0.in(multi_terms) == 0, nz0, true);

Constraint<> pp_nz1("pp_nz1");
pp_nz1 = pp + 1;
M_obj.add_on_off(pp_nz1.in(multi_terms) == 0, nz1, true);

Constraint<> pp_nz2("pp_nz2");
pp_nz2 = pp - 1;
M_obj.add_on_off(pp_nz2.in(multi_terms) == 0, nz2, true);

Constraint<> pp_nz3("pp_nz3");
pp_nz3 = pp + 1;
M_obj.add_on_off(pp_nz3.in(multi_terms) == 0, nz3, true);

Constraint<> pp_nz4("pp_nz4");
pp_nz4 = pp - 1;
M_obj.add_on_off(pp_nz4.in(multi_terms) == 0, nz4, true);

// var<> obj_var("obj_var", pos_);
// M_obj.add(obj_var.in(R(1)));
Expand All @@ -166,9 +229,9 @@ int main(int argc, char * argv[]){
// s.initialize_all(1);
solver<> g_sol(M_obj,gurobi);
g_sol.run();
M_obj.print_solution();
M_obj.round_solution();
M_obj.print_solution();
// M_obj.print_solution();
// M_obj.round_solution();
// M_obj.print_solution();
opt_obj = round(M_obj.get_obj_val());
}
else{
Expand Down Expand Up @@ -240,7 +303,7 @@ int main(int argc, char * argv[]){
auto f = ones.tr()*c*c;
M.min(ones.tr()*c*c);
// M.min(sum(cs));
// M.print();
M.print();
solver<> mip_solver(M,gurobi);
// solver<> nlp_solver(M,ipopt);
mip_solver.run(1e-6, 7200);
Expand Down
8 changes: 6 additions & 2 deletions include/gravity/constraint.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ class Constraint_{
vector<bool> _violated;
param<double> _onCoef; /**< Coefficient vector for on in on/off constraints */
param<double> _offCoef; /**< Coefficient vector for off in on/off constraints */


shared_ptr<param_> _on_off_bin = nullptr; /**< Pointer to the binary variable for on/off constraints */
bool _on_off = false; /**< if true, the constraints should be enforced if _on_off_bin = 1. If false, the constraints should be enforced if _on_off_bin = 0.*/


/* Accessors */
Expand Down Expand Up @@ -183,6 +183,8 @@ class Constraint: public Constraint_, public func<type>{
_all_satisfied = c._all_satisfied;
_violated = c._violated;
_relaxed = c._relaxed;
_on_off_bin = c._on_off_bin;
_on_off = c._on_off;
this->_name = c._name;
this->_is_constraint = true;
_onCoef = c._onCoef.deep_copy();
Expand All @@ -202,6 +204,8 @@ class Constraint: public Constraint_, public func<type>{
_all_satisfied = c._all_satisfied;
_violated = c._violated;
_relaxed = c._relaxed;
_on_off_bin = c._on_off_bin;
_on_off = c._on_off;
this->func<type>::operator=(c);
this->_name = c._name;
this->_is_constraint = true;
Expand Down
9 changes: 9 additions & 0 deletions include/gravity/model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6311,6 +6311,15 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str

}

void add_on_off(Constraint<type>& c, const var<int>& on, bool on_status){
if (c.get_ftype() != lin_) {
throw invalid_argument("Calling add_on_off with a nonlinear constraint.\n");
}
c._on_off_bin = make_shared<var<int>>(on);
c._on_off = on_status;
add_constraint(c);
}

//This function adds on-off version of a given linear constraint and the binary variables to activate. The boolean option handles all the facet definining inequalities of the convex hull (if true), else it only adds the Big_M version of the constraint
//INPUT: linear constraint to be activated, corresponding binary variables to form a disjunctive union, big_M version of the constraint or the whole convex hull
//OUTPUT: disjunctive union of the constraints provided by "c" and linked by "on"
Expand Down
28 changes: 20 additions & 8 deletions src/GurobiProgram.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -507,17 +507,18 @@ bool GurobiProgram::solve(bool relax, double mipgap, double time_limit){
// grb_mod->set(GRB_DoubleParam_IntFeasTol, 1e-9);
// grb_mod->set(GRB_IntParam_NumericFocus,3);
// grb_mod->set(GRB_IntParam_PreCrush,0);
// grb_mod->set(GRB_IntParam_MIPFocus,1);
grb_mod->set(GRB_IntParam_MIPFocus,1);
// grb_mod->set(GRB_IntParam_IntegralityFocus,1);
// grb_mod->set(GRB_IntParam_MIPFocus,1);
// grb_mod->set(GRB_DoubleParam_Heuristics,0.1);
// grb_mod->set(GRB_IntParam_PumpPasses,50);
// grb_mod->set(GRB_IntParam_RINS,1);
grb_mod->set(GRB_IntParam_Cuts,0);
// grb_mod->set(GRB_IntParam_Cuts,1);
// grb_mod->set(GRB_IntParam_PreQLinearize,1);

grb_mod->set(GRB_DoubleParam_TimeLimit,time_limit);
grb_mod->set(GRB_IntParam_OutputFlag,1);
// grb_mod->set(GRB_IntParam_StartNodeLimit, -3);
// grb_mod->set(GRB_DoubleParam_Cutoff,0);
// grb_mod->set(GRB_IntParam_MinRelNodes,0);
// grb_mod->set(GRB_DoubleParam_Heuristics, 0);
Expand Down Expand Up @@ -966,13 +967,15 @@ void GurobiProgram::create_grb_constraints(){
GRBVar gvar1, gvar2;
double coeff;
for(auto& p: _model->_cons){
bool lazy = false;
auto c = p.second;
if (!c->_new && c->_all_satisfied) {
continue;
}
if(c->_callback)
if(!c->is_linear() && c->is_ineq() && c->_callback)
{
DebugOn(c->_name<<" added as callback"<<endl);
DebugOn(c->_name<<" will be added in callback."<<endl);
lazy = true;
continue;
}
c->_new = false;
Expand Down Expand Up @@ -1020,10 +1023,19 @@ void GurobiProgram::create_grb_constraints(){
linlhs += lterm;
}
linlhs += c->eval(c->get_cst(), i);
if(c->_indices)
grb_mod->addConstr(linlhs,sense,0,c->get_name()+"("+c->_indices->_keys->at(i)+")");
else
grb_mod->addConstr(linlhs,sense,0,c->get_name());
string cname = c->get_name();
if(c->_indices){
cname += "("+c->_indices->_keys->at(i)+")";
}
if(c->_on_off_bin){
auto gvar = _grb_vars[c->_on_off_bin->get_id() + c->_on_off_bin->get_id_inst(i)];
grb_mod->addGenConstrIndicator(gvar, c->_on_off, linlhs, sense, 0, cname);
}
else{
auto gc = grb_mod->addConstr(linlhs,sense,0,cname);
if(c->_callback)
gc.set(GRB_IntAttr_Lazy, 3);
}
}

// }
Expand Down

0 comments on commit af516bf

Please sign in to comment.