FFSM++  1.1.0
French Forest Sector Model ++
Adolc_debugtest.cpp
Go to the documentation of this file.
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 "Adolc_debugtest.h"
30 
31 using namespace Ipopt;
32 
33 /* Constructor. */
35 {}
36 
38 
39 bool 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 
61 bool 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 
79 bool 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 
103 template<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 
116 template<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 
136 bool 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 
143 bool 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 
151 bool 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 
159 bool 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 
190 bool 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 
237 void 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 
264 void 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 }
#define tag_g
virtual bool eval_f(Index n, const Number *x, bool new_x, Number &obj_value)
virtual bool eval_jac_g(Index n, const Number *x, bool new_x, Index m, Index nele_jac, Index *iRow, Index *jCol, Number *values)
virtual bool get_bounds_info(Index n, Number *x_l, Number *x_u, Index m, Number *g_l, Number *g_u)
#define tag_f
virtual bool eval_g(Index n, const Number *x, bool new_x, Index m, Number *g)
virtual ~MyADOLC_NLP()
#define tag_L
virtual bool get_nlp_info(Index &n, Index &m, Index &nnz_jac_g, Index &nnz_h_lag, IndexStyleEnum &index_style)
virtual bool eval_h(Index n, const Number *x, bool new_x, Number obj_factor, Index m, const Number *lambda, bool new_lambda, Index nele_hess, Index *iRow, Index *jCol, Number *values)
virtual void generate_tapes(Index n, Index m)
bool eval_constraints(Index n, const T *x, Index m, T *g)
bool eval_obj(Index n, const T *x, T &obj_value)
virtual bool eval_grad_f(Index n, const Number *x, bool new_x, Number *grad_f)
virtual void finalize_solution(SolverReturn status, Index n, const Number *x, const Number *z_L, const Number *z_U, Index m, const Number *g, const Number *lambda, Number obj_value, const IpoptData *ip_data, IpoptCalculatedQuantities *ip_cq)
virtual bool get_starting_point(Index n, bool init_x, Number *x, bool init_z, Number *z_L, Number *z_U, Index m, bool init_lambda, Number *lambda)