Skip to content

Commit

Permalink
working on solution conversion from nl to Gurobi .sol
Browse files Browse the repository at this point in the history
  • Loading branch information
@hhijazi
hhijazi committed Jul 14, 2024
1 parent 0e0a59e commit e6b1ace
Show file tree
Hide file tree
Showing 9 changed files with 538 additions and 222 deletions.
16 changes: 9 additions & 7 deletions examples/Gravity_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,18 +23,20 @@ using namespace std;
using namespace gravity;

TEST_CASE("testing Gurobi's compound NL") {
double pi = 3.123456789;
cout << to_string_with_precision(pi, 4) << endl;
cout << to_string_with_precision(pi, 8) << endl;
cout << to_string_with_precision(pi, 10) << endl;
Model<> M;
M.set_name("M");
var<> x("x", -2,2);
var<> y("y", 0,0);
var<> z("z", pi/2.,pi/2.);
M.set_name("M0");
var<> x("x", -1,1);
var<> y("y", -1,1);
var<> resvar("resvar", -10,10);
M.add(resvar);
M.add(x);
M.add(y);
M.add(z);
Constraint<> C("C");
C = cos(cos(y) - sin(z + cos(y))) - x - 2 - resvar;
C = 6*sin(x)/cos(y) - resvar;
M.add(C==0);
M.max(resvar);
solver<> GRB(M,gurobi);
Expand All @@ -45,7 +47,7 @@ TEST_CASE("testing Gurobi's compound NL") {
#ifdef USE_MP
TEST_CASE("testing readNL() function") {
Model<> M;
string NL_file = string(prj_dir)+"/data_sets/NL/camshape100.nl";
string NL_file = "/Users/hassan.hijazi/Downloads/polygon25.nl.txt";
int status = M.readNL(NL_file);
solver<> GRB(M,gurobi);
double solver_time_start = get_wall_time();
Expand Down
11 changes: 11 additions & 0 deletions examples/Optimization/NonLinear/Power/ACOPF/ACOPF_main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,17 @@ int main (int argc, char * argv[])
mtype=argv[2];
}
#endif
Model<> M;
int status = M.readNL(fname);
solver<> GRB(M,gurobi);
GRB.run();
M.print();
M.read_solution(fname);
M.write_solution();
// M.reset();
// M.is_feasible(1e-6);
return 0;

double total_time_start = get_wall_time();
PowerNet grid;
grid.readgrid(fname);
Expand Down
5 changes: 4 additions & 1 deletion include/gravity/constant.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,10 @@ namespace gravity {
if(std::numeric_limits<T>::is_specialized && a_value==numeric_limits<T>::max()){
return "+∞";
}
out << std::setprecision(n) << a_value;
out.setf(ios::fixed);
out.precision(n);
out << std::scientific << a_value;
// out << std::setprecision(n) << a_value;
return out.str();
}
/**
Expand Down
272 changes: 142 additions & 130 deletions include/gravity/func.h
View file

Large diffs are not rendered by default.

56 changes: 42 additions & 14 deletions include/gravity/model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2951,9 +2951,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
c->_violated[inst] = false;
diff = std::abs(c->eval(inst));
if(diff > tol) {
DebugOff("Violated equation: ");
// c->print(inst);
DebugOff(", violation = "<< diff << endl);
DebugOn("Violated equation: ");
c->print(inst);
DebugOn(", violation = "<< diff << endl);
nb_viol++;
// violated = true;
if (*c->_all_lazy) {
Expand All @@ -2979,9 +2979,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
c->_violated[inst] = false;
diff = c->eval(inst);
if(diff > tol) {
DebugOff("Violated inequality: ");
// c->print(inst);
DebugOff(", violation = "<< diff << endl);
DebugOn("Violated inequality: ");
c->print(inst);
DebugOn(", violation = "<< diff << endl);
nb_viol++;
// violated = true;
if (*c->_all_lazy) {
Expand Down Expand Up @@ -3013,9 +3013,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
c->_violated[inst] = false;
diff = c->eval(inst);
if(diff < -tol) {
DebugOff("Violated inequality: ");
// c->print(inst);
DebugOff(", violation = "<< diff << endl);
DebugOn("Violated inequality: ");
c->print(inst);
DebugOn(", violation = "<< diff << endl);
nb_viol++;
// violated = true;
if (*c->_all_lazy) {
Expand Down Expand Up @@ -5077,11 +5077,14 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str

for(auto& c_p :_cons)
{
// c_p.second->uneval();
c_p.second->_new = true;
c_p.second->uneval();
c_p.second->eval_all();
c_p.second->_dfdx->clear();
}
_obj->_new = true;
_obj->uneval();
_obj->eval_all();
_obj->_dfdx->clear();
for(auto &v_p: _vars)
{
Expand Down Expand Up @@ -5859,6 +5862,21 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
return size_header;
}

/** Read solution point to file */
void read_solution(const string& fname);

/** Write solution point to file */
void write_solution(int precision = 10){
ofstream myfile;
string fname = _name+".sol";
myfile.open(fname);
std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
std::cout.rdbuf(myfile.rdbuf());
print_solution(precision);
std::cout.rdbuf(coutbuf);
myfile.close();
}

/** Write Model to file */
void write(int precision = 10){
ofstream myfile;
Expand Down Expand Up @@ -5897,6 +5915,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str



bool has_int() const{
return !_int_vars.empty();
}

void replace_integers(){/*< Replace internal type of integer variables so that continuous relaxations can be computed */
bool has_int = false;
Expand Down Expand Up @@ -8195,11 +8216,15 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
}

var<> get_cont_int_var(int index){
auto x = _model->get_var<double>("x");
auto y = _model->get_var<double>("y");
if(index >= y.get_id())/* Integer variable */
return y(index);
if(_model->has_int()){
auto y = _model->get_var<double>("y");
if(y._indices->has_key(to_string(index)))/* Integer variable */
{
return y(index);
}
}
/* Continuous variable */
auto x = _model->get_var<double>("x");
return x(index);
}

Expand Down Expand Up @@ -8271,6 +8296,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
func<> VisitLog(UnaryExpr e) {
return log(Visit(e.arg()));
}
func<> VisitLog10(UnaryExpr e) {
return log(Visit(e.arg()))/std::log(10);
}
// func<> VisitFloor(UnaryExpr e) {
// return IloFloor(Visit(e.arg()));
// }
Expand Down
4 changes: 2 additions & 2 deletions include/gravity/param.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2351,10 +2351,10 @@ namespace gravity {
}
else {
for (size_t i = 0; i < _dim[0]; i++) {
if((!_is_relaxed && !is_integer()) || std::abs(eval(i)) > 1e-4){
// if((!_is_relaxed && !is_integer()) || std::abs(eval(i)) > 1e-4){
str += name+"(" + _indices->_keys->at(i) + ") " + to_string_with_precision(eval(i), prec);
str += " \n";
}
// }
}
}
}
Expand Down
1 change: 1 addition & 0 deletions include/gravity/solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -434,6 +434,7 @@ class solver {
grb_prog->_output = output;
// prog.grb_prog->reset_model();
grb_prog->prepare_model();
return 0;
//DebugOn("calling prep after init"<<endl);
optimal = grb_prog->solve(relax, tol, time_limit);
return_status = optimal ? 0 : -1;
Expand Down
111 changes: 81 additions & 30 deletions src/GurobiProgram.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -430,7 +430,7 @@ GurobiProgram::GurobiProgram(Model<>* m) {
}
}
_model = m;
m->fill_in_maps();
// m->fill_in_maps();
m->compute_funcs();
}

Expand Down Expand Up @@ -814,10 +814,14 @@ bool GurobiProgram::solve(bool relax, double mipgap, double time_limit){
return true;
}

#include <fstream>
#include <iostream>
#include <string>

bool copyWithoutLastLine(const std::string& inputFile) {
std::string tempFile = inputFile + ".temp"; // Temporary file

// First, perform the operation with inputFile as input and tempFile as output
// Open input and output files
std::ifstream in(inputFile);
std::ofstream out(tempFile, std::ios::trunc);
if (!in.is_open() || !out.is_open()) {
Expand All @@ -827,17 +831,22 @@ bool copyWithoutLastLine(const std::string& inputFile) {

std::string line;
std::string prevLine;
while (std::getline(in, line)) {
bool foundGenerals = false;
while (!foundGenerals && std::getline(in, line)) {
if (line.find("Generals") != std::string::npos) {
foundGenerals = true; // Flag that we found the "Generals" line
}
if (!prevLine.empty()) {
out << prevLine << '\n'; // Write previous line, ensuring not to write the last one
out << prevLine << '\n'; // Write previous line
}
prevLine = line; // Keep track of the previous line
prevLine = line; // Update previous line
}


// Close files
in.close();
out.close();

// Now, replace the original file with the temporary file
// Replace original file with temporary file
if (std::remove(inputFile.c_str()) != 0) {
std::cerr << "Error removing original file." << std::endl;
return false;
Expand All @@ -851,8 +860,9 @@ bool copyWithoutLastLine(const std::string& inputFile) {
}



void GurobiProgram::prepare_model(){
_model->fill_in_maps();
// _model->fill_in_maps();
_model->compute_funcs();
fill_in_grb_vmap();
create_grb_constraints();
Expand All @@ -868,7 +878,7 @@ void GurobiProgram::prepare_model(){
write_NLCstr(filename);


exit(0);
return;
// print_constraints();
}
void GurobiProgram::update_model(){
Expand Down Expand Up @@ -1017,40 +1027,81 @@ void GurobiProgram::write_NLCstr(const string &fname){
cerr << "Failed to open file for appending nonlinear constraints.\n";
return;
}
// file << " zero = 0\n";
size_t idx = 0, idx_inst = 0, idx1 = 0, idx2 = 0, idx_inst1 = 0, idx_inst2 = 0, nb_inst = 0, inst = 0;
int pos_i = 0, neg_i = 0;
for(auto& p: _model->_cons){
auto c = p.second;
if (c->is_polynomial() || c->is_nonlinear()) {
nb_inst = c->get_nb_inst();
for (size_t i = 0; i< nb_inst; i++){
switch(c->get_ctype()) {
case geq:
pos_i++;
break;
case leq:
neg_i++;
break;
default:
break;
}
}
}
}
file << " zero(0) = 0\n";
for(int i = 0; i < pos_i; i++){
file << " 0 <= pos(" + to_string(i) + ")\n";
}
for(int i = 0; i < neg_i; i++){
file << " -Inf <= neg(" + to_string(i) + ") <= 0\n";
}
if(_model->has_int()){
auto y = _model->get_var_ptr("y");
file << " Generals\n";
file << " ";
for (auto i = 0; i<y->get_dim(); i++){
file << y->get_name(i) << " " ;
}
file << "\n";
}
file << " General Constraints\n";

char sense;
size_t idx = 0, idx_inst = 0, idx1 = 0, idx2 = 0, idx_inst1 = 0, idx_inst2 = 0, nb_inst = 0, inst = 0;
GRBLinExpr lterm, linlhs;
GRBQuadExpr quadlhs;
GRBVar gvar1, gvar2;
double coeff;
pos_i = 0;
neg_i = 0;
for(auto& p: _model->_cons){
auto c = p.second;
if (c->is_nonlinear() && (!(c->_expr->is_uexpr() && c->get_nb_vars()==2) && !c->_expr->is_mexpr())) {
DebugOff("Founf a nonlinear constraint with more than two variables, writing the expression tree.\n");
switch(c->get_ctype()) {
case geq:
sense = GRB_GREATER_EQUAL;
break;
case leq:
sense = GRB_LESS_EQUAL;
break;
case eq:
sense = GRB_EQUAL;
break;
default:
break;
}
assert(sense==GRB_EQUAL);/* Check what to do with inequalities */
// if(c->get_ctype()==eq)
// DebugOn(c->get_name() << endl);
// if(c->get_name()=="NL_C_eq_73")
// DebugOn(c->get_name() << endl);
if (c->is_polynomial() || c->is_nonlinear()) {
DebugOff("Found a nonlinear constraint with more than two variables, writing the expression tree.\n");
nb_inst = c->get_nb_inst();
for (size_t i = 0; i< nb_inst; i++){
if(c->_indices)
file << " " << c->get_name()+"("+c->_indices->_keys->at(i)+"):\n";
else
file << " " << c->get_name()+"("+to_string(i)+"):\n";
file << " zero(0) = NL :\n";
switch(c->get_ctype()) {
case geq:
sense = GRB_GREATER_EQUAL;
file << " pos("+ to_string(pos_i++) +") = NL :\n";
break;
case leq:
sense = GRB_LESS_EQUAL;
file << " neg("+ to_string(neg_i++) +") = NL :\n";
break;
case eq:
sense = GRB_EQUAL;
file << " zero(0) = NL :\n";
break;
default:
break;
}
int parent_id = -1;
file << c->getNLexpr(-1,parent_id,i);
}
Expand Down Expand Up @@ -1086,7 +1137,7 @@ void GurobiProgram::create_grb_constraints(){
}
c->_new = false;

if (c->is_nonlinear() && (!(c->_expr->is_uexpr() && c->get_nb_vars()==2) && !c->_expr->is_mexpr())) {
if (c->is_polynomial() || c->is_nonlinear()) {
DebugOff("Gurobi cannot handle nonlinear constraints with more than two variables, ignoring it.\n");
continue;
}
Expand Down Expand Up @@ -1228,7 +1279,7 @@ void GurobiProgram::create_grb_constraints(){
}
else {
coeff = c->eval(lt._coef,i);
if (!((coeff==1 && lt._sign) || (coeff==-1 && !lt._sign))) {
if (!((coeff==1 && !lt._sign) || (coeff==-1 && lt._sign))) {
throw invalid_argument("Gurobi does not support this type of nonlinear constraints");
}
gvar1 = _grb_vars[lt._p->get_id() + lt._p->get_id_inst(i)];
Expand Down
Loading

0 comments on commit e6b1ace

Please sign in to comment.