IPOPTSolverLean.cc 15.7 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
39


// smart pointer to IpoptApplication to set options etc.
40
class IPOPTSolverLean::Impl
41
42
{// Create an instance of the IpoptApplication
public:
43
  Impl()
44
    : app_(IpoptApplicationFactory()), max_iter_(200), alm_infsb_thrsh_(0.5),
45
    incr_lazy_cnstr_max_iter_nmbr_(5), enbl_all_lzy_cnstr_(true)
46
  {}
47
48
49

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

51
  int max_iter_;
52
  double alm_infsb_thrsh_;
53
54
  int incr_lazy_cnstr_max_iter_nmbr_;
  bool enbl_all_lzy_cnstr_;
55
56
57
58
59
60
61
};

// Constructor
IPOPTSolverLean::IPOPTSolverLean()
  : impl_(new Impl)
{

62
  // Switch to HSL if available
63
64
65
66
67
68
#if COMISO_HSL_AVAILABLE
  impl_->app_->Options()->SetStringValue("linear_solver", "ma57");
#else
  impl_->app_->Options()->SetStringValue("linear_solver", "mumps");
#endif

69
70
#ifdef DEB_ON
  if (!Debug::Config::query().console())
71
#endif
72
  {// Block any output on cout and cerr from Ipopt.
73
74
    impl_->app_->Options()->SetStringValue("suppress_all_output", "yes");
  }
75

76
77
78
79
80
81
82
83
84
85
86
87
88
89
#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
  impl_->app_->Options()->SetIntegerValue("mumps_mem_percent", 5);
#endif

  // set default parameters
  impl_->app_->Options()->SetIntegerValue("max_iter", 100);
  //  app->Options()->SetStringValue("derivative_test", "second-order");
  //  app->Options()->SetIntegerValue("print_level", 0);
  //  app->Options()->SetStringValue("expect_infeasible_problem", "yes");
}

90
91
IPOPTSolverLean::
~IPOPTSolverLean()
92
93
{ delete impl_; }

94
95
96
double
IPOPTSolverLean::
energy()
97
98
99
100
101
102
{
  return impl_->app_->Statistics()->FinalObjective();
}

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

103
104
105
106
void
IPOPTSolverLean::
set_max_iterations
(const int _max_iterations)
107
{
108
    impl_->max_iter_ = _max_iterations;
109
110
}

111
112
113
int
IPOPTSolverLean::
max_iterations() const
114
{
115
    return impl_->max_iter_;
116
117
}

118
119
120
double
IPOPTSolverLean::
almost_infeasible_threshold() const
121
122
123
124
{
  return impl_->alm_infsb_thrsh_;
}

125
126
127
128
void
IPOPTSolverLean::
set_almost_infeasible_threshold
(const double _alm_infsb_thrsh)
129
130
131
132
{
  impl_->alm_infsb_thrsh_ = _alm_infsb_thrsh;
}

133
134
135
int
IPOPTSolverLean::
incremental_lazy_constraint_max_iteration_number() const
136
{
137
  return impl_->incr_lazy_cnstr_max_iter_nmbr_;
138
139
}

140
141
142
143
void
IPOPTSolverLean::
set_incremental_lazy_constraint_max_iteration_number
(const int _incr_lazy_cnstr_max_iter_nmbr)
144
{
145
  impl_->incr_lazy_cnstr_max_iter_nmbr_ = _incr_lazy_cnstr_max_iter_nmbr;
146
147
}

148
149
150
bool
IPOPTSolverLean::
enable_all_lazy_contraints() const
151
{
152
  return impl_->enbl_all_lzy_cnstr_;
153
154
}

155
156
157
158
void
IPOPTSolverLean::
set_enable_all_lazy_contraints
(const bool _enbl_all_lzy_cnstr)
159
{
160
  impl_->enbl_all_lzy_cnstr_ = _enbl_all_lzy_cnstr;
161
162
163
}

//-----------------------------------------------------------------------------
164

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

209
210
211
static void
check_ipopt_status
(Ipopt::ApplicationReturnStatus const _stat)
212
{
213
214
215
  if (_stat != Ipopt::Solve_Succeeded &&
      _stat != Ipopt::Solved_To_Acceptable_Level)
      throw_ipopt_solve_failure(_stat);
216
217
}

218
219
220
221
222
void
IPOPTSolverLean::
solve
(NProblemInterface* _problem,
 const std::vector<NConstraintInterface*>& _constraints)
223
224
225
226
227
228
229
230
231
232
233
234
235
{
  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: ");
236
  if (np2->hessian_constant())
237
238
239
240
241
  {
    DEB_out(2,"*constant hessian* ");
    impl_->app_->Options()->SetStringValue("hessian_constant", "yes");
  }

242
  if (np2->jac_c_constant())
243
244
245
246
247
  {
    DEB_out(2, "*constant jacobian of equality constraints* ");
    impl_->app_->Options()->SetStringValue("jac_c_constant", "yes");
  }

248
  if (np2->jac_d_constant())
249
250
251
252
253
254
255
256
257
258
259
260
  {
    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();
261
  if (status != Ipopt::Solve_Succeeded)
262
263
264
265
266
267
268
    COMISO_THROW(IPOPT_INITIALIZATION_FAILED);

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

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
269
  check_ipopt_status(status);
270

271
272
  // Retrieve some statistics about the solve
  Ipopt::Index iter_count = impl_->app_->Statistics()->IterationCount();
273
  DEB_out(1,"\n*** IPOPT: The problem solved in "
274
275
276
277
278
279
280
281
282
283
284
    << 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");
}


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


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

307
308
309
310
311
312
313
314
315
  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;

316
317
  // set max iterations
  impl_->app_->Options()->SetIntegerValue("max_iter", impl_->max_iter_);
318
319

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

330
    //--------------------------------------------------------------------------
331
    // 2. exploit special characteristics of problem
332
    //--------------------------------------------------------------------------
333
334

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

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

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

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

362
363
    check_ipopt_status(status);

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

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

376
377
      bool inf = false;
      bool almost_inf = false;
378

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

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

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

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

449
    //--------------------------------------------------------------------------
450
    // 2. exploit special characteristics of problem
451
    //--------------------------------------------------------------------------
452
453

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

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

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

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

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
482
  check_ipopt_status(status);
483
484
485

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

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


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

504
505
506
void
IPOPTSolverLean::
solve(NProblemInterface* _problem)
507
508
509
510
511
512
513
514
{
  std::vector<NConstraintInterface*> constraints;
  solve(_problem, constraints);
}

//=============================================================================
} // namespace COMISO
//=============================================================================
515
#endif // COMISO_IPOPT_AVAILABLE
516
//=============================================================================