IPOPTSolverLean.cc 15.6 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
50
51
    : 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();
  }

private:
  void setup_ipopt_defaults();
52
53
54

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

56
  double alm_infsb_thrsh_;
57
58
  int incr_lazy_cnstr_max_iter_nmbr_;
  bool enbl_all_lzy_cnstr_;
59
60
61
62
63
64

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;
65
66
};

67
68
69
70
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;

71

72
73
void IPOPTSolverLean::Impl::setup_ipopt_defaults()
{
74
  // Switch to HSL if available
75
#if COMISO_HSL_AVAILABLE
76
  app_->Options()->SetStringValue("linear_solver", ipopt_default_hsl_solver);
77
#else
78
  app_->Options()->SetStringValue("linear_solver", "mumps");
79
80
#endif

81
82
#ifdef DEB_ON
  if (!Debug::Config::query().console())
83
#endif
84
  {// Block any output on cout and cerr from Ipopt.
85
    app_->Options()->SetStringValue("suppress_all_output", "yes");
86
  }
87

88
89
90
91
#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
92
93
  app_->Options()->SetIntegerValue("mumps_mem_percent",
                                   ipopt_default_mumps_mem_percent);
94
95
#endif

96
97
98
99
100
101
102
103
  // set max iterations
  app_->Options()->SetIntegerValue("max_iter", ipopt_default_max_iter);
}

// Constructor
IPOPTSolverLean::IPOPTSolverLean()
  : impl_(new Impl)
{
104
105
}

106
107
IPOPTSolverLean::
~IPOPTSolverLean()
108
109
110
111
{ delete impl_; }

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

112
113
114
115
void
IPOPTSolverLean::
set_max_iterations
(const int _max_iterations)
116
{
117
  impl_->app_->Options()->SetIntegerValue("max_iter", _max_iterations);
118
119
}

120
121
122
int
IPOPTSolverLean::
max_iterations() const
123
{
124
125
126
  int max_iter;
  impl_->app_->Options()->GetIntegerValue("max_iter", max_iter, "");
  return max_iter;
127
128
}

129
130
131
double
IPOPTSolverLean::
almost_infeasible_threshold() const
132
133
134
135
{
  return impl_->alm_infsb_thrsh_;
}

136
137
138
139
void
IPOPTSolverLean::
set_almost_infeasible_threshold
(const double _alm_infsb_thrsh)
140
141
142
143
{
  impl_->alm_infsb_thrsh_ = _alm_infsb_thrsh;
}

144
145
146
int
IPOPTSolverLean::
incremental_lazy_constraint_max_iteration_number() const
147
{
148
  return impl_->incr_lazy_cnstr_max_iter_nmbr_;
149
150
}

151
152
153
154
void
IPOPTSolverLean::
set_incremental_lazy_constraint_max_iteration_number
(const int _incr_lazy_cnstr_max_iter_nmbr)
155
{
156
  impl_->incr_lazy_cnstr_max_iter_nmbr_ = _incr_lazy_cnstr_max_iter_nmbr;
157
158
}

159
160
161
bool
IPOPTSolverLean::
enable_all_lazy_contraints() const
162
{
163
  return impl_->enbl_all_lzy_cnstr_;
164
165
}

166
167
168
169
void
IPOPTSolverLean::
set_enable_all_lazy_contraints
(const bool _enbl_all_lzy_cnstr)
170
{
171
  impl_->enbl_all_lzy_cnstr_ = _enbl_all_lzy_cnstr;
172
173
}

174
175
176
static void
throw_ipopt_solve_failure
(Ipopt::ApplicationReturnStatus const status)
177
178
{
  DEB_enter_func
Max Lyon's avatar
Max Lyon committed
179
  DEB_error(" IPOPT solve failure code is " << status)
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
  // 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
  //    };
  //------------------------------------------------------
206
  switch (status)
207
  {
208
  case Ipopt::Maximum_Iterations_Exceeded:
209
    COMISO_THROW(IPOPT_MAXIMUM_ITERATIONS_EXCEEDED);
210
  case Ipopt::NonIpopt_Exception_Thrown:
211
212
    // this could be due to a thrown PROGRESS_ABORTED exception, ...
    PROGRESS_RESUME_ABORT; // ... so check if we need to resume it
213
214
  default:
    COMISO_THROW(IPOPT_OPTIMIZATION_FAILED);
215
  }
216
217
}

218
219
220
static void
check_ipopt_status
(Ipopt::ApplicationReturnStatus const _stat)
221
{
222
223
224
  if (_stat != Ipopt::Solve_Succeeded &&
      _stat != Ipopt::Solved_To_Acceptable_Level)
      throw_ipopt_solve_failure(_stat);
225
226
}

227
228
229
230
231
void
IPOPTSolverLean::
solve
(NProblemInterface* _problem,
 const std::vector<NConstraintInterface*>& _constraints)
232
233
234
235
236
237
238
239
240
241
242
243
244
{
  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: ");
245
  if (np2->hessian_constant())
246
247
248
249
250
  {
    DEB_out(2,"*constant hessian* ");
    impl_->app_->Options()->SetStringValue("hessian_constant", "yes");
  }

251
  if (np2->jac_c_constant())
252
253
254
255
256
  {
    DEB_out(2, "*constant jacobian of equality constraints* ");
    impl_->app_->Options()->SetStringValue("jac_c_constant", "yes");
  }

257
  if (np2->jac_d_constant())
258
259
260
261
262
263
264
265
266
267
268
269
  {
    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();
270
  if (status != Ipopt::Solve_Succeeded)
271
272
273
274
275
276
277
    COMISO_THROW(IPOPT_INITIALIZATION_FAILED);

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

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
278
  check_ipopt_status(status);
279

280
281
  // Retrieve some statistics about the solve
  Ipopt::Index iter_count = impl_->app_->Statistics()->IterationCount();
282
  DEB_out(1,"\n*** IPOPT: The problem solved in "
283
284
285
286
287
288
289
    << 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");
}

290
291
292
293
294
295
void
IPOPTSolverLean::
solve
(NProblemInterface*                        _problem,
 const std::vector<NConstraintInterface*>& _constraints,
 const std::vector<NConstraintInterface*>& _lazy_constraints)
296
297
298
299
300
301
302
303
304
305
306
307
308
{
  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;
309
310
  int  cur_pass = impl_->enbl_all_lzy_cnstr_ ? 1 : 0;
  const int max_passes = impl_->incr_lazy_cnstr_max_iter_nmbr_;
311

312
313
314
315
316
317
318
319
320
  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;

321
  while(!feasible_point_found && cur_pass < max_passes)
322
323
  {
    ++cur_pass;
324
    //--------------------------------------------------------------------------
325
    // 1. Create an instance of current IPOPT NLP
326
    //--------------------------------------------------------------------------
327
328
329
330
331
    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;

332
    //--------------------------------------------------------------------------
333
    // 2. exploit special characteristics of problem
334
    //--------------------------------------------------------------------------
335
336

    DEB_out(2, "detected special properties which will be exploit: ");
337
    if (np2->hessian_constant())
338
339
340
341
342
    {
      DEB_out(2, "*constant hessian* ");
      impl_->app_->Options()->SetStringValue("hessian_constant", "yes");
    }

343
    if (np2->jac_c_constant())
344
345
346
347
348
    {
      DEB_out(2, "*constant jacobian of equality constraints* ");
      impl_->app_->Options()->SetStringValue("jac_c_constant", "yes");
    }

349
    if (np2->jac_d_constant())
350
351
352
353
354
355
    {
      DEB_out(2, "*constant jacobian of in-equality constraints*");
      impl_->app_->Options()->SetStringValue("jac_d_constant", "yes");
    }
    DEB_out(2, "\n");

356
    //--------------------------------------------------------------------------
357
    // 3. solve problem
358
    //--------------------------------------------------------------------------
359
360
361
362
    {
      DEB_time_session_def("IPOPT App OptimizeTNLP(np)");
      status = impl_->app_->OptimizeTNLP(np);
    }
363

364
365
    check_ipopt_status(status);

366
367
368
369
    // check lazy constraints
    n_inf.push_back(0);
    n_almost_inf.push_back(0);
    feasible_point_found = true;
370
371
372
373
374
    for (unsigned int i = 0; i < _lazy_constraints.size(); ++i)
    {
      if (lazy_added[i])
        continue;
      NConstraintInterface* lc = _lazy_constraints[i];
375

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

378
379
      bool inf = false;
      bool almost_inf = false;
380

381
382
383
384
385
386
387
388
389
390
391
      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)
392
        {
393
          if (v < -acceptable_tolerance)
394
395
            inf = true;
          else
396
            if (v < impl_->alm_infsb_thrsh_)
397
398
399
              almost_inf = true;
        }
        else
400
          if (lc->constraint_type() == NConstraintInterface::NC_LESS_EQUAL)
401
          {
402
            if (v > acceptable_tolerance)
403
404
              inf = true;
            else
405
              if (v > -impl_->alm_infsb_thrsh_)
406
407
408
                almost_inf = true;
          }

409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
      // 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();
424
      }
425
    }
426
427
428
  }

  // no termination after max number of passes?
429
  if (!feasible_point_found)
430
  {
431
432
    DEB_warning(2, "Could not find a feasible point after " << max_passes - 1
                << " incremental lazy constraint iterations");
433
    if (!impl_->enbl_all_lzy_cnstr_)
434
435
      throw_ipopt_solve_failure(Ipopt::Maximum_Iterations_Exceeded);

436
    DEB_line(2, "Solving with ALL lazy constraints...");
437
    ++cur_pass;
438
439
440
    for (unsigned int i = 0; i < _lazy_constraints.size(); ++i)
    {
      if (!lazy_added[i])
441
        constraints.push_back(_lazy_constraints[i]);
442
    }
443
    //--------------------------------------------------------------------------
444
    // 1. Create an instance of current IPOPT NLP
445
    //--------------------------------------------------------------------------
446
447
448
449
450
    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;

451
    //--------------------------------------------------------------------------
452
    // 2. exploit special characteristics of problem
453
    //--------------------------------------------------------------------------
454
455

    DEB_out(2, "exploit detected special properties: ");
456
    if (np2->hessian_constant())
457
458
459
460
461
    {
      DEB_out(2, "*constant hessian* ");
      impl_->app_->Options()->SetStringValue("hessian_constant", "yes");
    }

462
    if (np2->jac_c_constant())
463
464
465
466
467
    {
      DEB_out(2, "*constant jacobian of equality constraints* ");
      impl_->app_->Options()->SetStringValue("jac_c_constant", "yes");
    }

468
    if (np2->jac_d_constant())
469
470
471
472
473
474
    {
      DEB_out(2, "*constant jacobian of in-equality constraints*");
      impl_->app_->Options()->SetStringValue("jac_d_constant", "yes");
    }
    std::cerr << std::endl;

475
    //--------------------------------------------------------------------------
476
    // 3. solve problem
477
    //--------------------------------------------------------------------------
478
479
480
481
482
483
    status = impl_->app_->OptimizeTNLP( np);
  }

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
484
  check_ipopt_status(status);
485
486
487

  // Retrieve some statistics about the solve
  Ipopt::Index iter_count = impl_->app_->Statistics()->IterationCount();
488
  DEB_out(1, "\n*** IPOPT: The problem solved in "
489
490
491
492
493
494
    << 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");

495
496
  DEB_out(2, "############# IPOPT with "
          "lazy constraints statistics ###############\n");
497
  DEB_out(2, "#passes     : " << cur_pass << "( of " << max_passes << ")\n");
498
  for(unsigned int i=0; i<n_inf.size(); ++i)
499
    DEB_out(3, "pass " << i << " induced " << n_inf[i]
500
501
502
      << " infeasible and " << n_almost_inf[i] << " almost infeasible\n")
}

503
504
505
506
507
508
509
double
IPOPTSolverLean::
energy()
{
  return impl_->app_->Statistics()->FinalObjective();
}

510
511
512
//=============================================================================
} // namespace COMISO
//=============================================================================
513
#endif // COMISO_IPOPT_AVAILABLE
514
//=============================================================================