source: trunk/ADOL-C/examples/additional_examples/ipopt/LuksanVlcek1/ADOL-C_NLP.cpp @ 332

Last change on this file since 332 was 332, checked in by awalther, 7 years ago

update of ipopt examples

File size: 7.6 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     ADOL-C_NLP.cpp
4 Revision: $$
5 Contents: class myADOLC_NPL for interfacing with Ipopt
6 
7 Copyright (c) Andrea Walther
8   
9 This file is part of ADOL-C. This software is provided as open source.
10 Any use, reproduction, or distribution of the software constitutes
11 recipient's acceptance of the terms of the accompanying license file.
12 
13 This code is based on the file  MyNLP.cpp contained in the Ipopt package
14 with the authors:  Carl Laird, Andreas Waechter   
15----------------------------------------------------------------------------*/
16
17/** C++ Example NLP for interfacing a problem with IPOPT and ADOL-C.
18 *  MyADOL-C_NLP implements a C++ example showing how to interface
19 *  with IPOPT and ADOL-C through the TNLP interface. This class
20 *  implements the Example 5.1 from "Sparse and Parially Separable
21 *  Test Problems for Unconstrained and Equality Constrained
22 *  Optimization" by L. Luksan and J. Vlcek ignoring sparsity.
23 *
24 *  no exploitation of sparsity !!
25 *
26 */
27#include <cassert>
28
29#include "ADOL-C_NLP.hpp"
30
31using namespace Ipopt;
32
33/* Constructor. */
34MyADOLC_NLP::MyADOLC_NLP()
35{}
36
37MyADOLC_NLP::~MyADOLC_NLP(){}
38
39bool MyADOLC_NLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
40                         Index& nnz_h_lag, IndexStyleEnum& index_style)
41{
42  n = 20;
43
44  m = n-2;
45
46  // in this example the jacobian is dense. Hence, it contains n*m nonzeros
47  nnz_jac_g = n*m;
48
49  // the hessian is also dense and has n*n total nonzeros, but we
50  // only need the lower left corner (since it is symmetric)
51  nnz_h_lag = n*(n-1)/2+n;
52
53  generate_tapes(n, m);
54
55  // use the C style indexing (0-based)
56  index_style = C_STYLE;
57
58  return true;
59}
60
61bool MyADOLC_NLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
62                            Index m, Number* g_l, Number* g_u)
63{
64  // none of the variables have bounds
65  for (Index i=0; i<n; i++) {
66    x_l[i] = -1e20;
67    x_u[i] =  1e20;
68  }
69
70  // Set the bounds for the constraints
71  for (Index i=0; i<m; i++) {
72    g_l[i] = 0;
73    g_u[i] = 0;
74  }
75
76  return true;
77}
78
79bool MyADOLC_NLP::get_starting_point(Index n, bool init_x, Number* x,
80                               bool init_z, Number* z_L, Number* z_U,
81                               Index m, bool init_lambda,
82                               Number* lambda)
83{
84  // Here, we assume we only have starting values for x, if you code
85  // your own NLP, you can provide starting values for the others if
86  // you wish.
87  assert(init_x == true);
88  assert(init_z == false);
89  assert(init_lambda == false);
90
91  // set the starting point
92  for (Index i=0; i<n/2; i++) {
93    x[2*i] = -1.2;
94    x[2*i+1] = 1.;
95  }
96  if (n != 2*(n/2)) {
97    x[n-1] = -1.2;
98  }
99
100  return true;
101}
102
103template<class T> bool  MyADOLC_NLP::eval_obj(Index n, const T *x, T& obj_value)
104{
105  T a1, a2;
106  obj_value = 0.;
107  for (Index i=0; i<n-1; i++) {
108    a1 = x[i]*x[i]-x[i+1];
109    a2 = x[i] - 1.;
110    obj_value += 100.*a1*a1 + a2*a2;
111  }
112
113  return true;
114}
115
116template<class T> bool  MyADOLC_NLP::eval_constraints(Index n, const T *x, Index m, T* g)
117{
118  for (Index i=0; i<m; i++) {
119    g[i] = 3.*pow(x[i+1],3.) + 2.*x[i+2] - 5.
120           + sin(x[i+1]-x[i+2])*sin(x[i+1]+x[i+2]) + 4.*x[i+1]
121           - x[i]*exp(x[i]-x[i+1]) - 3.;
122  }
123
124  return true;
125}
126
127//*************************************************************************
128//
129//
130//         Nothing has to be changed below this point !!
131//
132//
133//*************************************************************************
134
135
136bool MyADOLC_NLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
137{
138  eval_obj(n,x,obj_value);
139
140  return true;
141}
142
143bool MyADOLC_NLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
144{
145
146  gradient(tag_f,n,x,grad_f);
147
148  return true;
149}
150
151bool MyADOLC_NLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
152{
153
154  eval_constraints(n,x,m,g);
155
156  return true;
157}
158
159bool MyADOLC_NLP::eval_jac_g(Index n, const Number* x, bool new_x,
160                       Index m, Index nele_jac, Index* iRow, Index *jCol,
161                       Number* values)
162{
163  if (values == NULL) {
164    // return the structure of the jacobian,
165    // assuming that the Jacobian is dense
166
167    Index idx = 0;
168    for(Index i=0; i<m; i++)
169      for(Index j=0; j<n; j++)
170        {
171          iRow[idx] = i;
172          jCol[idx++] = j;
173        }
174 }
175  else {
176    // return the values of the jacobian of the constraints
177
178    jacobian(tag_g,m,n,x,Jac);
179
180    Index idx = 0;
181    for(Index i=0; i<m; i++)
182      for(Index j=0; j<n; j++)
183          values[idx++] = Jac[i][j];
184
185  }
186
187  return true;
188}
189
190bool MyADOLC_NLP::eval_h(Index n, const Number* x, bool new_x,
191                   Number obj_factor, Index m, const Number* lambda,
192                   bool new_lambda, Index nele_hess, Index* iRow,
193                   Index* jCol, Number* values)
194{
195  if (values == NULL) {
196    // return the structure. This is a symmetric matrix, fill the lower left
197    // triangle only.
198
199    // the hessian for this problem is actually dense
200    Index idx=0;
201    for (Index row = 0; row < n; row++) {
202      for (Index col = 0; col <= row; col++) {
203        iRow[idx] = row;
204        jCol[idx] = col;
205        idx++;
206      }
207    }
208
209    assert(idx == nele_hess);
210  }
211  else {
212    // return the values. This is a symmetric matrix, fill the lower left
213    // triangle only
214
215    for(Index i = 0; i<n ; i++)
216      x_lam[i] = x[i];
217    for(Index i = 0; i<m ; i++)
218      x_lam[n+i] = lambda[i];
219    x_lam[n+m] = obj_factor;
220
221    hessian(tag_L,n+m+1,x_lam,Hess);
222
223    Index idx = 0;
224
225    for(Index i = 0; i<n ; i++)
226      {
227        for(Index j = 0; j<=i ; j++)
228          {
229            values[idx++] = Hess[i][j];
230          }
231      }
232  }
233
234  return true;
235}
236
237void MyADOLC_NLP::finalize_solution(SolverReturn status,
238                              Index n, const Number* x, const Number* z_L, const Number* z_U,
239                              Index m, const Number* g, const Number* lambda,
240                              Number obj_value,
241                              const IpoptData* ip_data,
242                              IpoptCalculatedQuantities* ip_cq)
243{
244
245  printf("\n\nObjective value\n");
246  printf("f(x*) = %e\n", obj_value);
247
248// Memory deallocation for ADOL-C variables
249
250  delete[] x_lam;
251
252  for(Index i=0;i<m;i++)
253    delete[] Jac[i];
254  delete[] Jac;
255
256  for(Index i=0;i<n+m+1;i++)
257    delete[] Hess[i];
258  delete[] Hess;
259}
260
261
262//***************    ADOL-C part ***********************************
263
264void MyADOLC_NLP::generate_tapes(Index n, Index m)
265{
266  Number *xp    = new double[n];
267  Number *lamp  = new double[m];
268  Number *zl    = new double[m];
269  Number *zu    = new double[m];
270
271  adouble *xa   = new adouble[n];
272  adouble *g    = new adouble[m];
273  adouble *lam  = new adouble[m];
274  adouble sig;
275  adouble obj_value;
276 
277  double dummy;
278
279  Jac = new double*[m];
280  for(Index i=0;i<m;i++)
281    Jac[i] = new double[n];
282
283  x_lam   = new double[n+m+1];
284
285  Hess = new double*[n+m+1];
286  for(Index i=0;i<n+m+1;i++)
287    Hess[i] = new double[i+1];
288
289  get_starting_point(n, 1, xp, 0, zl, zu, m, 0, lamp);
290
291  trace_on(tag_f);
292   
293    for(Index i=0;i<n;i++)
294      xa[i] <<= xp[i];
295
296    eval_obj(n,xa,obj_value);
297
298    obj_value >>= dummy;
299
300  trace_off();
301 
302  trace_on(tag_g);
303   
304    for(Index i=0;i<n;i++)
305      xa[i] <<= xp[i];
306
307    eval_constraints(n,xa,m,g);
308
309
310    for(Index i=0;i<m;i++)
311      g[i] >>= dummy;
312
313  trace_off();
314
315   trace_on(tag_L);
316   
317    for(Index i=0;i<n;i++)
318      xa[i] <<= xp[i];
319    for(Index i=0;i<m;i++)
320      lam[i] <<= 1.0;
321    sig <<= 1.0;
322
323    eval_obj(n,xa,obj_value);
324
325    obj_value *= sig;
326    eval_constraints(n,xa,m,g);
327 
328    for(Index i=0;i<m;i++)
329      obj_value += g[i]*lam[i];
330
331    obj_value >>= dummy;
332
333  trace_off();
334
335  delete[] xa;
336  delete[] xp;
337  delete[] g;
338  delete[] lam;
339  delete[] lamp;
340  delete[] zu;
341  delete[] zl;
342
343}
Note: See TracBrowser for help on using the repository browser.