IPOPTSolver.cc 15 KB
Newer Older
David Bommes's avatar
David Bommes committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
//=============================================================================
//
//  CLASS IPOPTSolver - IMPLEMENTATION
//
//=============================================================================

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

//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_IPOPT_AVAILABLE
//=============================================================================


#include "IPOPTSolver.hh"

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

namespace COMISO {

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


David Bommes's avatar
David Bommes committed
24
25
26
27
28
29
30
// Constructor
IPOPTSolver::
IPOPTSolver()
{
  // Create an instance of the IpoptApplication
  app_ = IpoptApplicationFactory();

31
32
33
34
35
36
37
38
39
40
41
  // Switch to HSL if available in Comiso
  #if COMISO_HSL_AVAILABLE
    app_->Options()->SetStringValue("linear_solver", "ma57");
  #endif

  // Restrict memory to be able to run larger problems on windows
  // with the default mumps solver
  #ifdef WIN32
    app_->Options()->SetIntegerValue("mumps_mem_percent", 5);
  #endif

David Bommes's avatar
David Bommes committed
42
43
44
45
46
47
  // set default parameters
  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");

David Bommes's avatar
David Bommes committed
48
  print_level_ = 5;
David Bommes's avatar
David Bommes committed
49
50
51
52
53
54
}


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


David Bommes's avatar
David Bommes committed
55
56
57

int
IPOPTSolver::
58
solve(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints)
David Bommes's avatar
David Bommes committed
59
{
David Bommes's avatar
David Bommes committed
60
61
62
63
64
65
66
67
68
69
70
71
72
  //----------------------------------------------------------------------------
  // 0. Check whether hessian_approximation is active
  //----------------------------------------------------------------------------
  bool hessian_approximation = false;
  std::string ha, p;
  app().Options()->GetStringValue("hessian_approximation", ha, p);
  if(ha != "exact")
  {
    if(print_level_>=2)
      std::cerr << "Hessian approximation is enabled" << std::endl;
    hessian_approximation = true;
  }

David Bommes's avatar
David Bommes committed
73
74
75
  //----------------------------------------------------------------------------
  // 1. Create an instance of IPOPT NLP
  //----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
76
  Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, _constraints, hessian_approximation);
77
  NProblemIPOPT* np2 = dynamic_cast<NProblemIPOPT*> (Ipopt::GetRawPtr(np));
David Bommes's avatar
David Bommes committed
78
79

  //----------------------------------------------------------------------------
80
81
82
  // 2. exploit special characteristics of problem
  //----------------------------------------------------------------------------

David Bommes's avatar
David Bommes committed
83
84
  if(print_level_>=2)
    std::cerr << "exploit detected special properties: ";
85
86
  if(np2->hessian_constant())
  {
David Bommes's avatar
David Bommes committed
87
88
    if(print_level_>=2)
      std::cerr << "*constant hessian* ";
89
90
91
92
93
    app().Options()->SetStringValue("hessian_constant", "yes");
  }

  if(np2->jac_c_constant())
  {
David Bommes's avatar
David Bommes committed
94
95
    if(print_level_>=2)
      std::cerr << "*constant jacobian of equality constraints* ";
96
97
98
99
100
    app().Options()->SetStringValue("jac_c_constant", "yes");
  }

  if(np2->jac_d_constant())
  {
David Bommes's avatar
David Bommes committed
101
102
    if(print_level_>=2)
      std::cerr << "*constant jacobian of in-equality constraints*";
103
104
    app().Options()->SetStringValue("jac_d_constant", "yes");
  }
David Bommes's avatar
David Bommes committed
105
106
107

  if(print_level_>=2)
    std::cerr << std::endl;
108
109
110

  //----------------------------------------------------------------------------
  // 3. solve problem
David Bommes's avatar
David Bommes committed
111
112
  //----------------------------------------------------------------------------

David Bommes's avatar
David Bommes committed
113
114
115
116
117
  // Initialize the IpoptApplication and process the options
  Ipopt::ApplicationReturnStatus status;
  status = app_->Initialize();
  if (status != Ipopt::Solve_Succeeded)
  {
David Bommes's avatar
David Bommes committed
118
119
    if(print_level_>=2)
      printf("\n\n*** Error IPOPT during initialization!\n");
David Bommes's avatar
David Bommes committed
120
121
  }

122
123
  status = app_->OptimizeTNLP( np);

David Bommes's avatar
David Bommes committed
124
  //----------------------------------------------------------------------------
125
  // 4. output statistics
David Bommes's avatar
David Bommes committed
126
  //----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
127
128
129
130
131
132
  if(print_level_>=2)
    if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
    {
      // 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);
133

David Bommes's avatar
David Bommes committed
134
135
136
      Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
      printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
    }
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153

  return status;
}


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



int
IPOPTSolver::
solve(NProblemInterface*                        _problem,
      const std::vector<NConstraintInterface*>& _constraints,
      const std::vector<NConstraintInterface*>& _lazy_constraints,
      const double                              _almost_infeasible,
      const int                                 _max_passes        )
{
David Bommes's avatar
David Bommes committed
154
155
156
157
158
159
160
161
  //----------------------------------------------------------------------------
  // 0. Check whether hessian_approximation is active
  //----------------------------------------------------------------------------
  bool hessian_approximation = false;
  std::string ha, p;
  app().Options()->GetStringValue("hessian_approximation", ha, p);
  if(ha != "exact")
  {
Marcel Campen's avatar
Marcel Campen committed
162
163
    if(print_level_>=2)
      std::cerr << "Hessian approximation is enabled" << std::endl;
David Bommes's avatar
David Bommes committed
164
165
166
    hessian_approximation = true;
  }

167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
  //----------------------------------------------------------------------------
  // 0. Initialize IPOPT Applicaiton
  //----------------------------------------------------------------------------

  StopWatch sw; sw.start();

  // Initialize the IpoptApplication and process the options
  Ipopt::ApplicationReturnStatus status;
  status = app_->Initialize();
  if (status != Ipopt::Solve_Succeeded)
  {
    printf("\n\n*** Error IPOPT during initialization!\n");
  }

  bool feasible_point_found = false;
  int  cur_pass = 0;
  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;

  while(!feasible_point_found && cur_pass <(_max_passes-1))
  {
Marcel Campen's avatar
Marcel Campen committed
194

195
196
197
198
    ++cur_pass;
    //----------------------------------------------------------------------------
    // 1. Create an instance of current IPOPT NLP
    //----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
199
    Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, constraints, hessian_approximation);
200
201
202
203
204
205
206
207
    NProblemIPOPT* np2 = dynamic_cast<NProblemIPOPT*> (Ipopt::GetRawPtr(np));
    // enable caching of solution
    np2->store_solution() = true;

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

Marcel Campen's avatar
Marcel Campen committed
208
    if(print_level_>=2) std::cerr << "detected special properties which will be exploit: ";
209
210
    if(np2->hessian_constant())
    {
Marcel Campen's avatar
Marcel Campen committed
211
      if(print_level_>=2) std::cerr << "*constant hessian* ";
212
213
214
215
216
      app().Options()->SetStringValue("hessian_constant", "yes");
    }

    if(np2->jac_c_constant())
    {
Marcel Campen's avatar
Marcel Campen committed
217
      if(print_level_>=2) std::cerr << "*constant jacobian of equality constraints* ";
218
219
220
221
222
      app().Options()->SetStringValue("jac_c_constant", "yes");
    }

    if(np2->jac_d_constant())
    {
Marcel Campen's avatar
Marcel Campen committed
223
      if(print_level_>=2) std::cerr << "*constant jacobian of in-equality constraints*";
224
225
      app().Options()->SetStringValue("jac_d_constant", "yes");
    }
Marcel Campen's avatar
Marcel Campen committed
226
    if(print_level_>=2) std::cerr << std::endl;
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
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
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

    //----------------------------------------------------------------------------
    // 3. solve problem
    //----------------------------------------------------------------------------
    status = app_->OptimizeTNLP( np);

    // 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(&(np2->solution()[0]));
        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)
        {
          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();
        }
      }
  }

  // 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])
        constraints.push_back(_lazy_constraints[i]);

    //----------------------------------------------------------------------------
    // 1. Create an instance of current IPOPT NLP
    //----------------------------------------------------------------------------
    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;

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

Marcel Campen's avatar
Marcel Campen committed
314
    if(print_level_>=2) std::cerr << "exploit detected special properties: ";
315
316
    if(np2->hessian_constant())
    {
Marcel Campen's avatar
Marcel Campen committed
317
      if(print_level_>=2) std::cerr << "*constant hessian* ";
318
319
320
321
322
      app().Options()->SetStringValue("hessian_constant", "yes");
    }

    if(np2->jac_c_constant())
    {
Marcel Campen's avatar
Marcel Campen committed
323
      if(print_level_>=2) std::cerr << "*constant jacobian of equality constraints* ";
324
325
326
327
328
      app().Options()->SetStringValue("jac_c_constant", "yes");
    }

    if(np2->jac_d_constant())
    {
Marcel Campen's avatar
Marcel Campen committed
329
      if(print_level_>=2) std::cerr << "*constant jacobian of in-equality constraints*";
330
331
      app().Options()->SetStringValue("jac_d_constant", "yes");
    }
Marcel Campen's avatar
Marcel Campen committed
332
    if(print_level_>=2) std::cerr << std::endl;
333
334
335
336
337
338
339
340

    //----------------------------------------------------------------------------
    // 3. solve problem
    //----------------------------------------------------------------------------
    status = app_->OptimizeTNLP( np);
  }

  const double overall_time = sw.stop()/1000.0;
David Bommes's avatar
David Bommes committed
341
342
343
344
345
346

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
  if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
  {
Marcel Campen's avatar
Marcel Campen committed
347
348
349
350
351
    if(print_level_>=2)
    {
      // 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);
David Bommes's avatar
David Bommes committed
352

Marcel Campen's avatar
Marcel Campen committed
353
354
355
      Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
      printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
    }
David Bommes's avatar
David Bommes committed
356
357
  }

Marcel Campen's avatar
Marcel Campen committed
358
359
360
361
362
363
364
365
  if(print_level_>=2)
  {
    std::cerr <<"############# IPOPT 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;
  }
366

David Bommes's avatar
David Bommes committed
367
368
  return status;
}
David Bommes's avatar
David Bommes committed
369

David Bommes's avatar
David Bommes committed
370
371
372
373

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


374
375
376
377
378
379
380
381
382
383
384
385
int
IPOPTSolver::
solve(NProblemInterface*    _problem)
{
  std::vector<NConstraintInterface*> constraints;
  return this->solve(_problem, constraints);
}


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


David Bommes's avatar
David Bommes committed
386
387
388
389
390
391
392
393
394
395
396
397
398
399
int
IPOPTSolver::
solve(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints)
{
  std::cerr << "****** Warning: NProblemGmmInterface is deprecated!!! -> use NProblemInterface *******\n";

  //----------------------------------------------------------------------------
  // 1. Create an instance of IPOPT NLP
  //----------------------------------------------------------------------------
  Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemGmmIPOPT(_problem, _constraints);

  //----------------------------------------------------------------------------
  // 2. solve problem
  //----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
400

David Bommes's avatar
David Bommes committed
401
402
  // Initialize the IpoptApplication and process the options
  Ipopt::ApplicationReturnStatus status;
David Bommes's avatar
David Bommes committed
403
  status = app_->Initialize();
David Bommes's avatar
David Bommes committed
404
405
406
407
408
409
410
411
  if (status != Ipopt::Solve_Succeeded)
  {
    printf("\n\n*** Error IPOPT during initialization!\n");
  }

  //----------------------------------------------------------------------------
  // 3. solve problem
  //----------------------------------------------------------------------------
David Bommes's avatar
David Bommes committed
412
  status = app_->OptimizeTNLP(np);
David Bommes's avatar
David Bommes committed
413
414
415
416
417
418
419

  //----------------------------------------------------------------------------
  // 4. output statistics
  //----------------------------------------------------------------------------
  if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
  {
    // Retrieve some statistics about the solve
David Bommes's avatar
David Bommes committed
420
    Ipopt::Index iter_count = app_->Statistics()->IterationCount();
David Bommes's avatar
David Bommes committed
421
422
    printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count);

David Bommes's avatar
David Bommes committed
423
    Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
David Bommes's avatar
David Bommes committed
424
425
426
427
428
429
430
431
432
433
434
    printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
  }

  return status;
}


//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_IPOPT_AVAILABLE
David Bommes's avatar
David Bommes committed
435
//=============================================================================