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

Last change on this file since 71 was 71, checked in by awalther, 10 years ago

add two example for the coupling of Ipopt and ADOL-C

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