GUROBISolver.cc 28.5 KB
Newer Older
David Bommes's avatar
David Bommes committed
1
2
3
4
5
6
7
8
9
10
11
12
//=============================================================================
//
//  CLASS GUROBISolver - IMPLEMENTATION
//
//=============================================================================

//== INCLUDES =================================================================

//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_GUROBI_AVAILABLE
//=============================================================================
13
 #include <gurobi_c++.h>
David Bommes's avatar
David Bommes committed
14
15

#include "GUROBISolver.hh"
16
17
18
19
#if (COMISO_QT_AVAILABLE)
#include "GurobiHelper.hh"
#endif//COMISO_QT_AVAILABLE
#include <CoMISo/Utils/CoMISoError.hh>
Max Lyon's avatar
Max Lyon committed
20
#include <CoMISo/Utils/StopWatch.hh>
21

David Bommes's avatar
David Bommes committed
22

23
#include <Base/Debug/DebTime.hh>
24
25
#include <stdexcept>

David Bommes's avatar
David Bommes committed
26
27
28
29
30
31
//== NAMESPACES ===============================================================

namespace COMISO {

//== IMPLEMENTATION ========================================================== 

32
33
34
35
36
37
38
39
40
41
42
double* P(std::vector<double>& _v)
{
  if( !_v.empty())
    return ((double*)&_v[0]);
  else
    return 0;
}


//-----------------------------------------------------------------------------

David Bommes's avatar
David Bommes committed
43

44
void add_constraint_to_model(COMISO::NConstraintInterface* _constraint, 
45
  GRBModel& _model, const std::vector<GRBVar>& _vars, const double* _x, int _lazy = 0)
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
{
  DEB_enter_func;
  DEB_warning_if(!_constraint->is_linear(), 1,
     "GUROBISolver received a problem with non-linear constraints!!!\n");

  GRBLinExpr lin_expr;
  NConstraintInterface::SVectorNC gc;
  _constraint->eval_gradient(_x, gc);

  NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
  for(; v_it; ++v_it)
    lin_expr += _vars[v_it.index()]*v_it.value();

  double b = _constraint->eval_constraint(_x);

61
62
  GRBConstr constr;
  switch(_constraint->constraint_type())
63
  {
64
65
66
    case NConstraintInterface::NC_EQUAL         : constr = _model.addConstr(lin_expr + b == 0); break;
    case NConstraintInterface::NC_LESS_EQUAL    : constr = _model.addConstr(lin_expr + b <= 0); break;
    case NConstraintInterface::NC_GREATER_EQUAL : constr = _model.addConstr(lin_expr + b >= 0); break;
67
  }
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
  if (_lazy > 0)
    constr.set(GRB_IntAttr_Lazy, _lazy);
}


//-----------------------------------------------------------------------------

std::vector<GRBVar> allocate_variables(NProblemInterface* _problem, const std::vector<PairIndexVtype>& _discrete_constraints, GRBModel& _model)
{
  // determine variable types: 0->real, 1->integer, 2->bool
  std::vector<char> vtypes(_problem->n_unknowns(),0);
  for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
    switch(_discrete_constraints[i].second)
    {
      case Integer: vtypes[_discrete_constraints[i].first] = 1; break;
      case Binary : vtypes[_discrete_constraints[i].first] = 2; break;
      default     : break;
    }

  // GUROBI variables
  std::vector<GRBVar> vars;
  // first all
  for( int i=0; i<_problem->n_unknowns(); ++i)
    switch(vtypes[i])
    {
      case 0 : vars.push_back( _model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_CONTINUOUS) ); break;
      case 1 : vars.push_back( _model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_INTEGER   ) ); break;
      case 2 : vars.push_back( _model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_BINARY    ) ); break;
    }

  // Integrate new variables
  _model.update();

  return vars;
}


//-----------------------------------------------------------------------------


void set_start(NProblemInterface* _problem, GRBModel& _model, std::vector<GRBVar>& _vars, const std::string& _start_solution_output_path)
{
  // set start
  std::vector<double> start(_vars.size());
  _problem->initial_x(start.data());

  for (int i = 0; i < _problem->n_unknowns(); ++i)
    _vars[i].set(GRB_DoubleAttr_Start, start[i]);

  _model.update();

  if (!_start_solution_output_path.empty())
120
  {
121
122
    std::cout << "Writing problem's start solution into file \"" << _start_solution_output_path << "\"." << std::endl;
    _model.write(_start_solution_output_path);
123
124
125
126
127
128
  }
}


//-----------------------------------------------------------------------------

David Bommes's avatar
David Bommes committed
129

130
void setup_constraints(NProblemInterface* _problem, const std::vector<NConstraintInterface *> &_constraints, GRBModel& _model,  const std::vector<GRBVar>& _vars)
131
132
133
{
  DEB_enter_func;

134
135
136
137
  // get zero vector
  std::vector<double> x(_problem->n_unknowns(), 0.0);

  for(unsigned int i=0; i<_constraints.size();  ++i)
138
  {
139
140
141
142
    DEB_warning_if(!_constraints[i]->is_linear(), 1,
      "GUROBISolver received a problem with non-linear constraints!!!");

    add_constraint_to_model(_constraints[i], _model, _vars, P(x));
143
  }
144
  _model.update();
145
}
146

147
148
149
150
151

//-----------------------------------------------------------------------------


void setup_energy(NProblemInterface* _problem, GRBModel& _model,  const std::vector<GRBVar>& _vars)
David Bommes's avatar
David Bommes committed
152
{
153
  DEB_enter_func;
David Bommes's avatar
David Bommes committed
154

155
156
  DEB_warning_if(!_problem->constant_hessian(), 1,
    "GUROBISolver received a problem with non-constant hessian!!!");
David Bommes's avatar
David Bommes committed
157

158
  GRBQuadExpr objective;
David Bommes's avatar
David Bommes committed
159

160
161
  // get zero vector
  std::vector<double> x(_problem->n_unknowns(), 0.0);
David Bommes's avatar
David Bommes committed
162

163
164
165
166
167
168
169
170
  // add quadratic part
  NProblemInterface::SMatrixNP H;
  _problem->eval_hessian(P(x), H);
  for( int i=0; i<H.outerSize(); ++i)
  {
    for (NProblemInterface::SMatrixNP::InnerIterator it(H,i); it; ++it)
      objective += 0.5*it.value()*_vars[it.row()]*_vars[it.col()];
  }
David Bommes's avatar
David Bommes committed
171

172
173
174
175
176
  // add linear part
  std::vector<double> g(_problem->n_unknowns());
  _problem->eval_gradient(P(x), P(g));
  for(unsigned int i=0; i<g.size(); ++i)
    objective += g[i]*_vars[i];
David Bommes's avatar
David Bommes committed
177

178
179
  // add constant part
  objective += _problem->eval_f(P(x));
Max Lyon's avatar
Max Lyon committed
180

181
182
183
184
  _model.set(GRB_IntAttr_ModelSense, 1);
  _model.setObjective(objective);
  _model.update();
}
David Bommes's avatar
David Bommes committed
185

186
//-----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
187
188


189
190
191
192
double solve_problem_two_phase(GRBModel& _model, double _time_limit0, double _gap0, double _time_limit1, double _gap1,
                             const std::string& _solution_input_path, const std::string& _problem_env_output_path, const std::string& _problem_output_path)
{
  DEB_enter_func;
David Bommes's avatar
David Bommes committed
193

194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
  double final_gap = -1;
  if (_solution_input_path.empty())
  {
    if (!_problem_env_output_path.empty())
    {
      DEB_line(5, "Writing problem's environment into file \"" << _problem_env_output_path << "\".");
      _model.getEnv().writeParams(_problem_env_output_path);
    }
#if (COMISO_QT_AVAILABLE)
    if (!_problem_output_path.empty())
    {
      DEB_line(5, "Writing problem into file \"" << _problem_output_path << "\".");
      GurobiHelper::outputModelToMpsGz(_model, _problem_output_path);
    }
#endif//COMISO_QT_AVAILABLE
David Bommes's avatar
David Bommes committed
209

210
211
212
213
214
    // optimize
    _model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit0);
    _model.getEnv().set(GRB_DoubleParam_MIPGap, _gap0);
    _model.optimize();
    final_gap = _model.get(GRB_DoubleAttr_MIPGap);
David Bommes's avatar
David Bommes committed
215

216
217
218
219
220
221
222
    // jump into phase 2?
    if(_model.get(GRB_DoubleAttr_MIPGap) > _gap1 && _time_limit1 > 0)
    {
      _model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit1);
      _model.getEnv().set(GRB_DoubleParam_MIPGap, _gap1);
      _model.optimize();
      final_gap = _model.get(GRB_DoubleAttr_MIPGap);
David Bommes's avatar
David Bommes committed
223
    }
224
225
226
227
228
  }
  else
  {
      DEB_line(5, "Reading solution from file \"" << _solution_input_path << "\".");
  }
David Bommes's avatar
David Bommes committed
229

230
231
  return final_gap;
}
David Bommes's avatar
David Bommes committed
232
233


234
//-----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
235
236


237
238
239
240
241
double solve_problem(GRBModel& _model, double _time_limit, double _gap, const std::string& _solution_input_path, const std::string& _problem_env_output_path, const std::string& _problem_output_path)
{
  return solve_problem_two_phase(_model, _time_limit, _gap, 0, 1000,
                                 _solution_input_path, _problem_env_output_path, _problem_output_path);
}
David Bommes's avatar
David Bommes committed
242
243
244



245
//-----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
246

David Bommes's avatar
David Bommes committed
247

248
249
250
void store_result(NProblemInterface* _problem, GRBModel& _model, const std::vector<GRBVar>& _vars, const std::string& _solution_input_path)
{
  DEB_enter_func;
David Bommes's avatar
David Bommes committed
251

252
253
  // get zero vector
  std::vector<double> x(_problem->n_unknowns(), 0.0);
254

255
256
257
258
259
260
261
262
  if (_solution_input_path.empty())
  {
    // store computed result
    for(unsigned int i=0; i<_vars.size(); ++i)
      x[i] = _vars[i].get(GRB_DoubleAttr_X);
  }
  else
  {
263
#if (COMISO_QT_AVAILABLE)
264
265
266
267
268
269
270
271
      DEB_line(5, "Loading stored solution from \"" << _solution_input_path << "\".");
      // store loaded result
      const size_t oldSize = x.size();
      x.clear();
      GurobiHelper::readSolutionVectorFromSOL(x, _solution_input_path);
      COMISO_THROW_TODO_if(oldSize != x.size(),
                           "oldSize != x.size() <=> " << oldSize << " != " << x.size() <<
                           "\nLoaded solution vector doesn't have expected dimension.");
272
#endif//COMISO_QT_AVAILABLE
273
  }
David Bommes's avatar
David Bommes committed
274

275
  _problem->store_result(P(x));
David Bommes's avatar
David Bommes committed
276

277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
  // ObjVal is only available if the optimize was called.
  DEB_out_if(_solution_input_path.empty(), 2,
      "GUROBI Objective: " << _model.get(GRB_DoubleAttr_ObjVal) << "\n");
}


//-----------------------------------------------------------------------------


static void process_gurobi_exception(const GRBException& _exc)
{
  DEB_enter_func;
  DEB_error("Gurobi exception error code = " << _exc.getErrorCode() << 
    " and message: [" << _exc.getMessage() << "]");
    
  // NOTE: we could propagate e.getMessage() either using std::exception, 
  // or a specialized Reform exception type

  // The GRB_ error codes are defined in gurobi_c.h Gurobi header. 
  switch (_exc.getErrorCode())
  {
  case GRB_ERROR_NO_LICENSE: 
    COMISO_THROW(GUROBI_LICENCE_ABSENT); 
  case GRB_ERROR_SIZE_LIMIT_EXCEEDED: 
    COMISO_THROW(GUROBI_LICENCE_MODEL_TOO_LARGE); break;
  default: 
    COMISO_THROW(UNSPECIFIED_GUROBI_EXCEPTION);
  }
}

bool
GUROBISolver::
solve(NProblemInterface*                  _problem,
      const std::vector<NConstraintInterface *> &_constraints,
      const std::vector<PairIndexVtype>&  _discrete_constraints,
      const double                        _time_limit,
      const double                        _gap)
{
  DEB_enter_func;
  try
  {
    GRBEnv   env   = GRBEnv();
    GRBModel model = GRBModel(env);
David Bommes's avatar
David Bommes committed
320

321
322
323
324
325
326
    auto vars = allocate_variables(_problem, _discrete_constraints, model);
    set_start(_problem, model, vars, start_solution_output_path_);
    setup_constraints(_problem, _constraints, model, vars);
    setup_energy(_problem, model, vars);
    solve_problem(model, _time_limit, _gap, solution_input_path_, problem_env_output_path_, problem_output_path_);
    store_result(_problem, model, vars, solution_input_path_);
David Bommes's avatar
David Bommes committed
327
328
329

    return true;
  }
Jan Möbius's avatar
Jan Möbius committed
330
  catch(GRBException& e)
David Bommes's avatar
David Bommes committed
331
  {
332
    //process_gurobi_exception(e);
333
334
    DEB_error("Error code = " << e.getErrorCode());
    DEB_error(e.getMessage());
David Bommes's avatar
David Bommes committed
335
336
337
338
    return false;
  }
  catch(...)
  {
339
    DEB_error("Exception during optimization");
David Bommes's avatar
David Bommes committed
340
341
342
343
344
    return false;
  }
}


David Bommes's avatar
David Bommes committed
345
346
//-----------------------------------------------------------------------------

Marcel Campen's avatar
Marcel Campen committed
347
348
349
bool
GUROBISolver::
solve_two_phase(NProblemInterface*                  _problem,                // problem instance
350
351
            const std::vector<NConstraintInterface*>& _constraints,            // linear constraints
            const std::vector<PairIndexVtype>&        _discrete_constraints,   // discrete constraints
Marcel Campen's avatar
Marcel Campen committed
352
353
354
355
356
357
358
359
360
            const double                        _time_limit0, // time limit phase 1 in seconds
            const double                        _gap0,     // MIP gap phase 1
            const double                        _time_limit1, // time limit phase 2 in seconds
            const double                        _gap1)       // MIP gap phase 2
{
  double dummy;
  return solve_two_phase(_problem, _constraints, _discrete_constraints, _time_limit0, _gap0, _time_limit1, _gap1, dummy);
}

David Bommes's avatar
David Bommes committed
361

362
363
364
bool
GUROBISolver::
solve_two_phase(NProblemInterface*                  _problem,                // problem instance
365
366
           const std::vector<NConstraintInterface *> &_constraints,            // linear constraints
           const std::vector<PairIndexVtype> &_discrete_constraints,   // discrete constraints
367
368
369
           const double                        _time_limit0, // time limit phase 1 in seconds
           const double                        _gap0,     // MIP gap phase 1
           const double                        _time_limit1, // time limit phase 2 in seconds
Marcel Campen's avatar
Marcel Campen committed
370
371
           const double                        _gap1,       // MIP gap phase 2
           double&                             _final_gap)  //return final gap
372
{
Max Lyon's avatar
Max Lyon committed
373
374
  DEB_enter_func;

375
376
377
378
379
  try
  {
    GRBEnv   env   = GRBEnv();
    GRBModel model = GRBModel(env);

380
381
382
383
384
385
386
    auto vars = allocate_variables(_problem, _discrete_constraints, model);
    set_start(_problem, model, vars, start_solution_output_path_);
    setup_constraints(_problem, _constraints, model, vars);
    setup_energy(_problem, model, vars);
    _final_gap = solve_problem_two_phase(model, _time_limit0, _gap0, _time_limit1, _gap1,
                                         solution_input_path_, problem_env_output_path_, problem_output_path_);
    store_result(_problem, model, vars, solution_input_path_);
387
388
389
390
391

    return true;
  }
  catch(GRBException& e)
  {
392
    //process_gurobi_exception(e);
393
394
    DEB_error("Error code = " << e.getErrorCode());
    DEB_error(e.getMessage());
395
396
397
398
    return false;
  }
  catch(...)
  {
399
    DEB_error("Exception during optimization");
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
    return false;
  }

  return false;
}


//-----------------------------------------------------------------------------


bool
GUROBISolver::
solve(NProblemInterface*                        _problem,
      const std::vector<NConstraintInterface*>& _constraints,
      const std::vector<NConstraintInterface*>& _lazy_constraints,
      std::vector<PairIndexVtype>&              _discrete_constraints,   // discrete constraints
      const double                              _almost_infeasible,
      const int                                 _max_passes,
      const double                              _time_limit,
      const bool                                _silent)
{
421
  DEB_time_func_def;
422
423
424
425
426
427
428
//  // hack! solve with all constraints
//  std::vector<NConstraintInterface*> all_constraints;
//  std::copy(_constraints.begin(),_constraints.end(),std::back_inserter(all_constraints));
//  std::copy(_lazy_constraints.begin(),_lazy_constraints.end(),std::back_inserter(all_constraints));
//
//  return solve(_problem, all_constraints, _discrete_constraints, _time_limit);

Max Lyon's avatar
Max Lyon committed
429
430
  StopWatch sw; sw.start();

431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
  bool feasible_point_found = false;
  int  cur_pass = 0;
  double acceptable_tolerance = 0.01; // hack: read out from ipopt!!!
  std::vector<bool> lazy_added(_lazy_constraints.size(),false);

  // cache statistics of all iterations
  std::vector<int> n_inf;
  std::vector<int> n_almost_inf;

  try
  {
    //----------------------------------------------
    // 0. set up environment
    //----------------------------------------------

    GRBEnv   env   = GRBEnv();
    GRBModel model = GRBModel(env);

449
450
451
452
    auto vars = allocate_variables(_problem, _discrete_constraints, model);
    set_start(_problem, model, vars, start_solution_output_path_);
    setup_constraints(_problem, _constraints, model, vars);
    setup_energy(_problem, model, vars);
453
454

    //----------------------------------------------
455
    // 4. iteratively solve problem
456
457
458
459
    //----------------------------------------------

    // get zero vector
    std::vector<double> x(_problem->n_unknowns(), 0.0);
460
    model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
    bool solution_found = false;

    while(!feasible_point_found && cur_pass <(_max_passes-1))
    {
      ++cur_pass;
      //----------------------------------------------------------------------------
      // 1. optimize current Model and get result
      //----------------------------------------------------------------------------

      model.update();
      model.optimize();
      // store computed result
      if(model.get(GRB_IntAttr_Status) == GRB_OPTIMAL)
        for(unsigned int i=0; i<vars.size(); ++i)
          x[i] = vars[i].get(GRB_DoubleAttr_X);

      //----------------------------------------------------------------------------
      // 2. Check lazy constraints
      //----------------------------------------------------------------------------
      // check lazy constraints
      n_inf.push_back(0);
      n_almost_inf.push_back(0);
      feasible_point_found = true;
      for(unsigned int i=0; i<_lazy_constraints.size(); ++i)
        if(!lazy_added[i])
        {
          NConstraintInterface* lc = _lazy_constraints[i];

          double v = lc->eval_constraint(P(x));

          bool inf        = false;
          bool almost_inf = false;

          if(lc->constraint_type() == NConstraintInterface::NC_EQUAL)
          {
            v = std::abs(v);
            if(v>acceptable_tolerance)
              inf = true;
            else
              if(v>_almost_infeasible)
                almost_inf = true;
          }
503
504
505
506
          else if(lc->constraint_type() == NConstraintInterface::NC_GREATER_EQUAL)
          {
            if(v<-acceptable_tolerance)
              inf = true;
507
            else
508
509
510
511
512
513
514
515
516
517
518
              if(v<_almost_infeasible)
                almost_inf = true;
          }
          else if(lc->constraint_type() == NConstraintInterface::NC_LESS_EQUAL)
          {
            if(v>acceptable_tolerance)
              inf = true;
            else
              if(v>-_almost_infeasible)
                almost_inf = true;
          }
519
520

          // infeasible?
521
522
523
524
525
526
527
528
529
530
531
532
533
          if(inf)
          {
            add_constraint_to_model( lc, model, vars, P(x));
            lazy_added[i] = true;
            feasible_point_found = false;
            ++n_inf.back();
          }
          else if(almost_inf) // almost violated or violated? -> add to constraints
          {
            add_constraint_to_model( lc, model, vars, P(x));
            lazy_added[i] = true;
            ++n_almost_inf.back();
          }
534
535
536
537
538
539
540
541
        }
    }

    // no termination after max number of passes?
    if(!feasible_point_found)
    {
      ++cur_pass;

542
543
      DEB_out(1, "*************** could not find feasible point after " 
        << _max_passes-1 << " -> solving with all lazy constraints...\n");
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
      for(unsigned int i=0; i<_lazy_constraints.size(); ++i)
        if(!lazy_added[i])
        {
          add_constraint_to_model( _lazy_constraints[i], model, vars, P(x));
        }

      model.update();
      model.optimize();

      // store computed result
      if(model.get(GRB_IntAttr_Status) == GRB_OPTIMAL)
        for(unsigned int i=0; i<vars.size(); ++i)
          x[i] = vars[i].get(GRB_DoubleAttr_X);
    }

    const double overall_time = sw.stop()/1000.0;

    //----------------------------------------------------------------------------
    // 4. output statistics
    //----------------------------------------------------------------------------
564
565
566
567
    // Make this code DEB_only
    DEB_line(2,
      "############# GUROBI with lazy constraints statistics ###############");
    DEB_line(2, "#passes     : " << cur_pass << "( of " << _max_passes << ")");
568
    for(unsigned int i=0; i<n_inf.size(); ++i)
569
570
571
572
573
    {
      DEB_line(2, "pass " << i << " induced " << n_inf[i]
        << " infeasible and " << n_almost_inf[i] << " almost infeasible");
    }
    DEB_line(2, "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal));
574
575
576
577

    //----------------------------------------------
    // 5. store result
    //----------------------------------------------
578
579
580
    //COMISO_THROW_TODO_if(model.get(GRB_IntAttr_Status) != GRB_OPTIMAL, 
    //  "Gurobi solution not optimal");

581
582
583
584
585
586
587
588
589
590
    if(model.get(GRB_IntAttr_Status) != GRB_OPTIMAL)
      return false;
    else
    {
      _problem->store_result(P(x));
      return true;
    }
  }
  catch(GRBException& e)
  {
591
    //process_gurobi_exception(e);
592
593
    DEB_error("Error code = " << e.getErrorCode());
    DEB_error(e.getMessage());
594
595
596
597
    return false;
  }
  catch(...)
  {
598
    DEB_error("Exception during optimization");
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
    return false;
  }

  return false;

//    //----------------------------------------------
//    // 4. iteratively solve problem
//    //----------------------------------------------
//    IloBool solution_found = IloFalse;
//
//    while(!feasible_point_found && cur_pass <(_max_passes-1))
//    {
//      ++cur_pass;
//      //----------------------------------------------------------------------------
//      // 1. Create an instance of current Model and optimize
//      //----------------------------------------------------------------------------
//
//      if(!_silent)
//        std::cerr << "cplex -> setup IloCPlex...\n";
//      IloCplex cplex(model);
//      cplex.setParam(IloCplex::TiLim, _time_limit);
//      // silent mode?
//      if(_silent)
//        cplex.setOut(env_.getNullStream());
//
//      if(!_silent)
//        std::cerr << "cplex -> optimize...\n";
//      solution_found = cplex.solve();
//
//      if(solution_found != IloFalse)
//      {
//        for(unsigned int i=0; i<vars.size(); ++i)
//          x[i] = cplex.getValue(vars[i]);
//
//        _problem->store_result(P(x));
//      }
//      else continue;
//
//      //----------------------------------------------------------------------------
//      // 2. Check lazy constraints
//      //----------------------------------------------------------------------------
//      // check lazy constraints
//      n_inf.push_back(0);
//      n_almost_inf.push_back(0);
//      feasible_point_found = true;
//      for(unsigned int i=0; i<_lazy_constraints.size(); ++i)
//        if(!lazy_added[i])
//        {
//          NConstraintInterface* lc = _lazy_constraints[i];
//
//          double v = lc->eval_constraint(P(x));
//
//          bool inf        = false;
//          bool almost_inf = false;
//
//          if(lc->constraint_type() == NConstraintInterface::NC_EQUAL)
//          {
//            v = std::abs(v);
//            if(v>acceptable_tolerance)
//              inf = true;
//            else
//              if(v>_almost_infeasible)
//                almost_inf = true;
//          }
//          else
//            if(lc->constraint_type() == NConstraintInterface::NC_GREATER_EQUAL)
//            {
//              if(v<-acceptable_tolerance)
//                inf = true;
//              else
//                if(v<_almost_infeasible)
//                  almost_inf = true;
//            }
//            else
//              if(lc->constraint_type() == NConstraintInterface::NC_LESS_EQUAL)
//              {
//                if(v>acceptable_tolerance)
//                  inf = true;
//                else
//                  if(v>-_almost_infeasible)
//                    almost_inf = true;
//              }
//
//          // infeasible?
//          if(inf)
//          {
//            add_constraint_to_model( lc, vars, model, env_, P(x));
//            lazy_added[i] = true;
//            feasible_point_found = false;
//            ++n_inf.back();
//          }
//          else // almost violated or violated? -> add to constraints
//            if(almost_inf)
//            {
//              add_constraint_to_model( lc, vars, model, env_, P(x));
//              lazy_added[i] = true;
//              ++n_almost_inf.back();
//            }
//        }
//    }
//
//    // no termination after max number of passes?
//    if(!feasible_point_found)
//    {
//      ++cur_pass;
//
//      std::cerr << "*************** could not find feasible point after " << _max_passes-1 << " -> solving with all lazy constraints..." << std::endl;
//      for(unsigned int i=0; i<_lazy_constraints.size(); ++i)
//        if(!lazy_added[i])
//        {
//          add_constraint_to_model( _lazy_constraints[i], vars, model, env_, P(x));
//        }
//
//      //----------------------------------------------------------------------------
//      // 1. Create an instance of current Model and optimize
//      //----------------------------------------------------------------------------
//
//      IloCplex cplex(model);
//      cplex.setParam(IloCplex::TiLim, _time_limit);
//      // silent mode?
//      if(_silent)
//        cplex.setOut(env_.getNullStream());
//      solution_found = cplex.solve();
//
//      if(solution_found != IloFalse)
//      {
//        for(unsigned int i=0; i<vars.size(); ++i)
//          x[i] = cplex.getValue(vars[i]);
//
//        _problem->store_result(P(x));
//      }
//    }
//
//    const double overall_time = sw.stop()/1000.0;
//
//    //----------------------------------------------------------------------------
//    // 4. output statistics
//    //----------------------------------------------------------------------------
////    if (solution_found != IloFalse)
////    {
////      // Retrieve some statistics about the solve
////      Ipopt::Index iter_count = app_->Statistics()->IterationCount();
////      printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count);
////
////      Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
////      printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
////    }
//
//    std::cerr <<"############# CPLEX with lazy constraints statistics ###############" << std::endl;
//    std::cerr << "overall time: " << overall_time << "s" << std::endl;
//    std::cerr << "#passes     : " << cur_pass << "( of " << _max_passes << ")" << std::endl;
//    for(unsigned int i=0; i<n_inf.size(); ++i)
//      std::cerr << "pass " << i << " induced " << n_inf[i] << " infeasible and " << n_almost_inf[i] << " almost infeasible" << std::endl;
//
753
754
755
756
757
758
759
760
761
762
  //    return (solution_found != IloFalse);
}

bool
GUROBISolver::
solve(NProblemInterface* _problem,
      const std::vector<NConstraintInterface*>& _constraints,
      const std::vector<NConstraintInterface*>& _lazy_constraints,
      const std::vector<PairIndexVtype>& _discrete_constraints,
      const double _time_limit,
763
764
      const double _gap,
      const int _lazy_level)
765
766
767
768
769
770
771
772
773
774
775
776
777
778
{
  DEB_enter_func;
  try
  {
    //----------------------------------------------
    // 0. set up environment
    //----------------------------------------------

    GRBEnv   env   = GRBEnv();
    GRBModel model = GRBModel(env);

    //----------------------------------------------
    // 1. allocate variables
    //----------------------------------------------
779
780
    auto vars = allocate_variables(_problem, _discrete_constraints, model);
    set_start(_problem, model, vars, start_solution_output_path_);
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796

    //----------------------------------------------
    // 2. setup constraints
    //----------------------------------------------

    // get zero vector
    std::vector<double> x(_problem->n_unknowns(), 0.0);

    for(unsigned int i=0; i<_constraints.size();  ++i)
    {
      add_constraint_to_model(_constraints[i], model, vars, P(x));
    }
    model.update();

    for(unsigned int i=0; i<_lazy_constraints.size();  ++i)
    {
797
      add_constraint_to_model(_lazy_constraints[i], model, vars, P(x), _lazy_level);
798
799
800
801
802
803
    }
    model.update();

    //----------------------------------------------
    // 3. setup energy
    //----------------------------------------------
804
    setup_energy(_problem, model, vars);
805
806
807
808

    //----------------------------------------------
    // 4. solve problem
    //----------------------------------------------
809
    solve_problem(model, _time_limit, _gap, solution_input_path_, problem_env_output_path_, problem_output_path_);
810
811
812
813

    //----------------------------------------------
    // 5. store result
    //----------------------------------------------
814
    store_result(_problem, model, vars, solution_input_path_);
815
816
817
818
819
820

    return true;
  }
  catch(GRBException& e)
  {
    //process_gurobi_exception(e);
821
822
    DEB_error("Error code = " << e.getErrorCode());
    DEB_error(e.getMessage());
823
824
825
826
    return false;
  }
  catch(...)
  {
827
    DEB_error("Exception during optimization");
828
829
830
831
    return false;
  }

  return false;
832
833
834
835
836
837
}


//-----------------------------------------------------------------------------


David Bommes's avatar
David Bommes committed
838
839
840
841
842
843
844
845
846
847
848
void
GUROBISolver::
set_problem_output_path( const std::string &_problem_output_path)
{
  problem_output_path_ = _problem_output_path;
}


//-----------------------------------------------------------------------------


Max Lyon's avatar
Max Lyon committed
849
850
851
852
853
854
855
856
857
858
859
860
861
void
GUROBISolver::
set_start_solution_output_path(const std::string& _start_solution_output_path)
{
  std::stringstream ss;
  ss << _start_solution_output_path << ".mst";
  start_solution_output_path_ = ss.str();
}


//-----------------------------------------------------------------------------


David Bommes's avatar
David Bommes committed
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
void
GUROBISolver::
set_problem_env_output_path( const std::string &_problem_env_output_path)
{
  problem_env_output_path_ = _problem_env_output_path;
}


//-----------------------------------------------------------------------------


void
GUROBISolver::
set_solution_input_path(const std::string &_solution_input_path)
{
  solution_input_path_ = _solution_input_path;
}


David Bommes's avatar
David Bommes committed
881
882
883
884
885
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_GUROBI_AVAILABLE
//=============================================================================