IPOPTSolverLean.cc 17 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
//=============================================================================
//
//  CLASS IPOPTSolverLean - IMPLEMENTATION
//
//=============================================================================

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

//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
11
#if COMISO_IPOPT_AVAILABLE
12
13
14
15
16
17
//=============================================================================


#include "IPOPTSolverLean.hh"
#include "NProblemGmmInterface.hh"
#include "NProblemInterface.hh"
18
#include "NProblemIPOPT.hh"
19
20
21
22
#include "NConstraintInterface.hh"
#include "BoundConstraint.hh"
#include "CoMISo/Utils/CoMISoError.hh"

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

Max Lyon's avatar
Max Lyon committed
26
#include <CoMISo/Utils/gmm.hh>
27
28
29
30
31
32
33
34
35

#include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>

//== NAMESPACES ===============================================================

namespace COMISO {

36
//== IMPLEMENTATION ===========================================================
37

38
class IPOPTSolverLean::Impl
39
{
40
public:
41
  Impl()
42
43
44
45
46
47
48
49
    : app_(IpoptApplicationFactory()), // Create an instance of IpoptApplication
      alm_infsb_thrsh_(0.5),
      incr_lazy_cnstr_max_iter_nmbr_(5),
      enbl_all_lzy_cnstr_(true)
  {
    setup_ipopt_defaults();
  }

50
51
52
53
54
55
56
  void set_ipopt_option(std::string, int);
  void set_ipopt_option(std::string, double);
  void set_ipopt_option(std::string, std::string);

  template <typename T>
  T get_ipopt_option(std::string);

57
58
private:
  void setup_ipopt_defaults();
59
60
61

public:
  Ipopt::SmartPtr<Ipopt::IpoptApplication> app_;
62

63
  double alm_infsb_thrsh_;
64
65
  int incr_lazy_cnstr_max_iter_nmbr_;
  bool enbl_all_lzy_cnstr_;
66
67
68
69
70
71

private:
  // default ipopt options
  static const std::string ipopt_default_hsl_solver;
  static const int ipopt_default_max_iter;
  static const int ipopt_default_mumps_mem_percent;
72
73
};

74
75
76
77
const std::string IPOPTSolverLean::Impl::ipopt_default_hsl_solver = "ma57";
const int IPOPTSolverLean::Impl::ipopt_default_max_iter = 200;
const int IPOPTSolverLean::Impl::ipopt_default_mumps_mem_percent = 5;

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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
void
IPOPTSolverLean::Impl::
set_ipopt_option
(std::string option, int value)
{
  app_->Options()->SetIntegerValue(option, value);
}

void
IPOPTSolverLean::Impl::
set_ipopt_option
(std::string option, double value)
{
  app_->Options()->SetNumericValue(option, value);
}

void
IPOPTSolverLean::Impl::
set_ipopt_option
(std::string option, std::string value)
{
  app_->Options()->SetStringValue(option, value);
}

template <typename T> T
IPOPTSolverLean::Impl::
get_ipopt_option(std::string)
{
  // @TODO print warning about unsupported option type!
}

template <> int
IPOPTSolverLean::Impl::
get_ipopt_option<int>(std::string option)
{
  int value;
  app_->Options()->GetIntegerValue(option, value, "");
  return value;
}

template <> double
IPOPTSolverLean::Impl::
get_ipopt_option<double>(std::string option)
{
  double value;
  app_->Options()->GetNumericValue(option, value, "");
  return value;
}

template <> std::string
IPOPTSolverLean::Impl::
get_ipopt_option<std::string>(std::string option)
{
  std::string value;
  app_->Options()->GetStringValue(option, value, "");
  return value;
}
135

136
137
void IPOPTSolverLean::Impl::setup_ipopt_defaults()
{
138
  // Switch to HSL if available
139
#if COMISO_HSL_AVAILABLE
140
  set_ipopt_option("linear_solver", ipopt_default_hsl_solver);
141
#else
142
  set_ipopt_option("linear_solver", "mumps");
143
144
#endif

145
146
#ifdef DEB_ON
  if (!Debug::Config::query().console())
147
#endif
148
  {// Block any output on cout and cerr from Ipopt.
149
    set_ipopt_option("suppress_all_output", "yes");
150
  }
151

152
153
154
155
#ifdef WIN32
  // Restrict memory to be able to run larger problems on windows
  // with the default mumps solver
  // TODO: find out what this does and whether it makes sense to do it
156
157
  set_ipopt_option("mumps_mem_percent",
                   ipopt_default_mumps_mem_percent);
158
159
#endif

160
161
  // set maximum solver iterations
  set_ipopt_option("max_iter", ipopt_default_max_iter);
162
163
}

164
165
//-----------------------------------------------------------------------------

166
167
168
IPOPTSolverLean::IPOPTSolverLean()
  : impl_(new Impl)
{
169
170
}

171
172
IPOPTSolverLean::
~IPOPTSolverLean()
173
174
175
{
  delete impl_;
}
176

177
template <typename T>
178
179
void
IPOPTSolverLean::
180
181
set_ipopt_option
(std::string option, const T& value)
182
{
183
  impl_->set_ipopt_option(option, value);
184
185
}

186
template <typename T> T
187
IPOPTSolverLean::
188
get_ipopt_option(std::string option)
189
{
190
  return impl_->get_ipopt_option<T>(option);
191
192
}

193
void
194
IPOPTSolverLean::
195
196
set_max_iterations
(const int _max_iterations)
197
{
198
199
200
201
202
203
204
205
  impl_->set_ipopt_option("max_iter", _max_iterations);
}

int
IPOPTSolverLean::
get_max_iterations() const
{
  return impl_->get_ipopt_option<int>("max_iter");
206
207
}

208
209
210
211
void
IPOPTSolverLean::
set_almost_infeasible_threshold
(const double _alm_infsb_thrsh)
212
213
214
215
{
  impl_->alm_infsb_thrsh_ = _alm_infsb_thrsh;
}

216
double
217
IPOPTSolverLean::
218
get_almost_infeasible_threshold() const
219
{
220
  return impl_->alm_infsb_thrsh_;
221
222
}

223
224
225
226
void
IPOPTSolverLean::
set_incremental_lazy_constraint_max_iteration_number
(const int _incr_lazy_cnstr_max_iter_nmbr)
227
{
228
  impl_->incr_lazy_cnstr_max_iter_nmbr_ = _incr_lazy_cnstr_max_iter_nmbr;
229
230
}

231
int
232
IPOPTSolverLean::
233
get_incremental_lazy_constraint_max_iteration_number() const
234
{
235
  return impl_->incr_lazy_cnstr_max_iter_nmbr_;
236
237
}

238
239
240
241
void
IPOPTSolverLean::
set_enable_all_lazy_contraints
(const bool _enbl_all_lzy_cnstr)
242
{
243
  impl_->enbl_all_lzy_cnstr_ = _enbl_all_lzy_cnstr;
244
245
}

246
247
248
249
250
251
252
bool
IPOPTSolverLean::
get_enable_all_lazy_contraints() const
{
  return impl_->enbl_all_lzy_cnstr_;
}

253
254
255
static void
throw_ipopt_solve_failure
(Ipopt::ApplicationReturnStatus const status)
256
257
{
  DEB_enter_func
Max Lyon's avatar
Max Lyon committed
258
  DEB_error(" IPOPT solve failure code is " << status)
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
  // TODO: we could translate these return codes, but will not do it for now
  //  enum ApplicationReturnStatus
  //    {
  //      Solve_Succeeded=0,
  //      Solved_To_Acceptable_Level=1,
  //      Infeasible_Problem_Detected=2,
  //      Search_Direction_Becomes_Too_Small=3,
  //      Diverging_Iterates=4,
  //      User_Requested_Stop=5,
  //      Feasible_Point_Found=6,
  //
  //      Maximum_Iterations_Exceeded=-1,
  //      Restoration_Failed=-2,
  //      Error_In_Step_Computation=-3,
  //      Maximum_CpuTime_Exceeded=-4,
  //      Not_Enough_Degrees_Of_Freedom=-10,
  //      Invalid_Problem_Definition=-11,
  //      Invalid_Option=-12,
  //      Invalid_Number_Detected=-13,
  //
  //      Unrecoverable_Exception=-100,
  //      NonIpopt_Exception_Thrown=-101,
  //      Insufficient_Memory=-102,
  //      Internal_Error=-199
  //    };
  //------------------------------------------------------
285
  switch (status)
286
  {
287
  case Ipopt::Maximum_Iterations_Exceeded:
288
    COMISO_THROW(IPOPT_MAXIMUM_ITERATIONS_EXCEEDED);
289
  case Ipopt::NonIpopt_Exception_Thrown:
290
291
    // this could be due to a thrown PROGRESS_ABORTED exception, ...
    PROGRESS_RESUME_ABORT; // ... so check if we need to resume it
292
293
  default:
    COMISO_THROW(IPOPT_OPTIMIZATION_FAILED);
294
  }
295
296
}

297
298
299
static void
check_ipopt_status
(Ipopt::ApplicationReturnStatus const _stat)
300
{
301
302
303
  if (_stat != Ipopt::Solve_Succeeded &&
      _stat != Ipopt::Solved_To_Acceptable_Level)
      throw_ipopt_solve_failure(_stat);
304
305
}

306
307
308
309
310
void
IPOPTSolverLean::
solve
(NProblemInterface* _problem,
 const std::vector<NConstraintInterface*>& _constraints)
311
312
313
314
315
316
317
318
319
320
321
322
323
{
  DEB_time_func_def;
  //----------------------------------------------------------------------------
  // 1. Create an instance of IPOPT NLP
  //----------------------------------------------------------------------------
  Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, _constraints);
  NProblemIPOPT* np2 = dynamic_cast<NProblemIPOPT*> (Ipopt::GetRawPtr(np));

  //----------------------------------------------------------------------------
  // 2. exploit special characteristics of problem
  //----------------------------------------------------------------------------

  DEB_out(2,"exploit detected special properties: ");
324
  if (np2->hessian_constant())
325
326
327
328
329
  {
    DEB_out(2,"*constant hessian* ");
    impl_->app_->Options()->SetStringValue("hessian_constant", "yes");
  }

330
  if (np2->jac_c_constant())
331
332
333
334
335
  {
    DEB_out(2, "*constant jacobian of equality constraints* ");
    impl_->app_->Options()->SetStringValue("jac_c_constant", "yes");
  }

336
  if (np2->jac_d_constant())
337
338
339
340
341
342
343
344
345
346
347
348
  {
    DEB_out(2, "*constant jacobian of in-equality constraints*");
    impl_->app_->Options()->SetStringValue("jac_d_constant", "yes");
  }
  DEB_out(2,"\n");

  //----------------------------------------------------------------------------
  // 3. solve problem
  //----------------------------------------------------------------------------

  // Initialize the IpoptApplication and process the options
  Ipopt::ApplicationReturnStatus status = impl_->app_->Initialize();
349
  if (status != Ipopt::Solve_Succeeded)
350
351
352
353
354
355
356
    COMISO_THROW(IPOPT_INITIALIZATION_FAILED);

  status = impl_->app_->OptimizeTNLP( np);

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
357
  check_ipopt_status(status);
358

359
360
  // Retrieve some statistics about the solve
  Ipopt::Index iter_count = impl_->app_->Statistics()->IterationCount();
361
  DEB_out(1,"\n*** IPOPT: The problem solved in "
362
363
364
365
366
367
368
    << iter_count << " iterations!\n");

  Ipopt::Number final_obj = impl_->app_->Statistics()->FinalObjective();
  DEB_out(1,"\n*** IPOPT: The final value of the objective function is "
    << final_obj << "\n");
}

369
370
371
372
373
374
void
IPOPTSolverLean::
solve
(NProblemInterface*                        _problem,
 const std::vector<NConstraintInterface*>& _constraints,
 const std::vector<NConstraintInterface*>& _lazy_constraints)
375
376
377
378
379
380
381
382
383
384
385
386
387
{
  DEB_time_func_def;
  //----------------------------------------------------------------------------
  // 0. Initialize IPOPT Application
  //----------------------------------------------------------------------------

  // Initialize the IpoptApplication and process the options
  Ipopt::ApplicationReturnStatus status;
  status = impl_->app_->Initialize();
  if (status != Ipopt::Solve_Succeeded)
    COMISO_THROW(IPOPT_INITIALIZATION_FAILED);

  bool feasible_point_found = false;
388
389
  int  cur_pass = impl_->enbl_all_lzy_cnstr_ ? 1 : 0;
  const int max_passes = impl_->incr_lazy_cnstr_max_iter_nmbr_;
390

391
392
393
394
395
396
397
398
399
  double acceptable_tolerance = 0.01; // hack: read out from ipopt!!!
  // copy default constraints
  std::vector<NConstraintInterface*> constraints = _constraints;
  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;

400
  while(!feasible_point_found && cur_pass < max_passes)
401
402
  {
    ++cur_pass;
403
    //--------------------------------------------------------------------------
404
    // 1. Create an instance of current IPOPT NLP
405
    //--------------------------------------------------------------------------
406
407
408
409
410
    Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, constraints);
    NProblemIPOPT* np2 = dynamic_cast<NProblemIPOPT*> (Ipopt::GetRawPtr(np));
    // enable caching of solution
    np2->store_solution() = true;

411
    //--------------------------------------------------------------------------
412
    // 2. exploit special characteristics of problem
413
    //--------------------------------------------------------------------------
414
415

    DEB_out(2, "detected special properties which will be exploit: ");
416
    if (np2->hessian_constant())
417
418
419
420
421
    {
      DEB_out(2, "*constant hessian* ");
      impl_->app_->Options()->SetStringValue("hessian_constant", "yes");
    }

422
    if (np2->jac_c_constant())
423
424
425
426
427
    {
      DEB_out(2, "*constant jacobian of equality constraints* ");
      impl_->app_->Options()->SetStringValue("jac_c_constant", "yes");
    }

428
    if (np2->jac_d_constant())
429
430
431
432
433
434
    {
      DEB_out(2, "*constant jacobian of in-equality constraints*");
      impl_->app_->Options()->SetStringValue("jac_d_constant", "yes");
    }
    DEB_out(2, "\n");

435
    //--------------------------------------------------------------------------
436
    // 3. solve problem
437
    //--------------------------------------------------------------------------
438
439
440
441
    {
      DEB_time_session_def("IPOPT App OptimizeTNLP(np)");
      status = impl_->app_->OptimizeTNLP(np);
    }
442

443
444
    check_ipopt_status(status);

445
446
447
448
    // check lazy constraints
    n_inf.push_back(0);
    n_almost_inf.push_back(0);
    feasible_point_found = true;
449
450
451
452
453
    for (unsigned int i = 0; i < _lazy_constraints.size(); ++i)
    {
      if (lazy_added[i])
        continue;
      NConstraintInterface* lc = _lazy_constraints[i];
454

455
      double v = lc->eval_constraint(&(np2->solution()[0]));
456

457
458
      bool inf = false;
      bool almost_inf = false;
459

460
461
462
463
464
465
466
467
468
469
470
      if (lc->constraint_type() == NConstraintInterface::NC_EQUAL)
      {
        v = std::abs(v);
        if (v > acceptable_tolerance)
          inf = true;
        else
          if (v > impl_->alm_infsb_thrsh_)
            almost_inf = true;
      }
      else
        if (lc->constraint_type() == NConstraintInterface::NC_GREATER_EQUAL)
471
        {
472
          if (v < -acceptable_tolerance)
473
474
            inf = true;
          else
475
            if (v < impl_->alm_infsb_thrsh_)
476
477
478
              almost_inf = true;
        }
        else
479
          if (lc->constraint_type() == NConstraintInterface::NC_LESS_EQUAL)
480
          {
481
            if (v > acceptable_tolerance)
482
483
              inf = true;
            else
484
              if (v > -impl_->alm_infsb_thrsh_)
485
486
487
                almost_inf = true;
          }

488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
      // infeasible?
      if (inf)
      {
        constraints.push_back(lc);
        lazy_added[i] = true;
        feasible_point_found = false;
        ++n_inf.back();
      }

      // almost violated or violated? -> add to constraints
      if (almost_inf)
      {
        constraints.push_back(lc);
        lazy_added[i] = true;
        ++n_almost_inf.back();
503
      }
504
    }
505
506
507
  }

  // no termination after max number of passes?
508
  if (!feasible_point_found)
509
  {
510
511
    DEB_warning(2, "Could not find a feasible point after " << max_passes - 1
                << " incremental lazy constraint iterations");
512
    if (!impl_->enbl_all_lzy_cnstr_)
513
514
      throw_ipopt_solve_failure(Ipopt::Maximum_Iterations_Exceeded);

515
    DEB_line(2, "Solving with ALL lazy constraints...");
516
    ++cur_pass;
517
518
519
    for (unsigned int i = 0; i < _lazy_constraints.size(); ++i)
    {
      if (!lazy_added[i])
520
        constraints.push_back(_lazy_constraints[i]);
521
    }
522
    //--------------------------------------------------------------------------
523
    // 1. Create an instance of current IPOPT NLP
524
    //--------------------------------------------------------------------------
525
526
527
528
529
    Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, constraints);
    NProblemIPOPT* np2 = dynamic_cast<NProblemIPOPT*> (Ipopt::GetRawPtr(np));
    // enable caching of solution
    np2->store_solution() = true;

530
    //--------------------------------------------------------------------------
531
    // 2. exploit special characteristics of problem
532
    //--------------------------------------------------------------------------
533
534

    DEB_out(2, "exploit detected special properties: ");
535
    if (np2->hessian_constant())
536
537
538
539
540
    {
      DEB_out(2, "*constant hessian* ");
      impl_->app_->Options()->SetStringValue("hessian_constant", "yes");
    }

541
    if (np2->jac_c_constant())
542
543
544
545
546
    {
      DEB_out(2, "*constant jacobian of equality constraints* ");
      impl_->app_->Options()->SetStringValue("jac_c_constant", "yes");
    }

547
    if (np2->jac_d_constant())
548
549
550
551
552
553
    {
      DEB_out(2, "*constant jacobian of in-equality constraints*");
      impl_->app_->Options()->SetStringValue("jac_d_constant", "yes");
    }
    std::cerr << std::endl;

554
    //--------------------------------------------------------------------------
555
    // 3. solve problem
556
    //--------------------------------------------------------------------------
557
558
559
560
561
562
    status = impl_->app_->OptimizeTNLP( np);
  }

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
563
  check_ipopt_status(status);
564
565
566

  // Retrieve some statistics about the solve
  Ipopt::Index iter_count = impl_->app_->Statistics()->IterationCount();
567
  DEB_out(1, "\n*** IPOPT: The problem solved in "
568
569
570
571
572
573
    << iter_count << " iterations!\n");

  Ipopt::Number final_obj = impl_->app_->Statistics()->FinalObjective();
  DEB_out(1, "\n*** IPOPT: The final value of the objective function is "
    << final_obj << "\n");

574
575
  DEB_out(2, "############# IPOPT with "
          "lazy constraints statistics ###############\n");
576
  DEB_out(2, "#passes     : " << cur_pass << "( of " << max_passes << ")\n");
577
  for(unsigned int i=0; i<n_inf.size(); ++i)
578
    DEB_out(3, "pass " << i << " induced " << n_inf[i]
579
580
581
      << " infeasible and " << n_almost_inf[i] << " almost infeasible\n")
}

582
583
584
585
586
587
588
double
IPOPTSolverLean::
energy()
{
  return impl_->app_->Statistics()->FinalObjective();
}

589
590
591
//=============================================================================
} // namespace COMISO
//=============================================================================
592
#endif // COMISO_IPOPT_AVAILABLE
593
//=============================================================================