source: trunk/ADOL-C/examples/additional_examples/ipopt/MittelmannDistCntrlNeumA/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: 10.1 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 a distributed control problem with homogeneous
21 *  Neumann boundary conditions, as formulated by Hans Mittelmann as
22 *  Examples 4-6 in "Optimization Techniques for Solving Elliptic
23 *  Control Problems with Control and State Constraints. Part 2:
24 *  Distributed Control" ignoring sparsity.
25 *
26 *  no exploitation of sparsity !!
27 *
28 */
29
30#include <cassert>
31
32#include "ADOL-C_NLP.hpp"
33
34#define  alpha_ 0.01
35#define  b_0j_  1.
36#define  b_1j_  1.
37#define  b_i0_  1.
38#define  b_i1_  1.
39
40double *y_d_;
41int  N_;
42double  h_;
43double  hh_;
44
45int y_index(int i, int j) 
46  {
47    return j + (N_+2)*i;
48  }
49
50 int u_index(int i, int j)
51  {
52    return (N_+2)*(N_+2) + (j-1) + (N_)*(i-1);
53  }
54
55int pde_index(int i, int j) 
56  {
57    return (j-1) + N_*(i-1);
58  }
59
60double y_d_cont(double x1, double x2)
61  {
62    return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.));
63  }
64
65adouble fint_cont(double x1, double x2, adouble y, adouble u)
66  {
67    adouble diff_y = y-y_d_cont(x1,x2);
68    return 0.5*(diff_y*diff_y + alpha_*u*u);
69  }
70
71adouble d_cont(double x1, double x2, adouble y, adouble u)
72  {
73    return -exp(y) - u;
74  }
75
76double fint_cont(double x1, double x2, double y, double u)
77  {
78    double diff_y = y-y_d_cont(x1,x2);
79    return 0.5*(diff_y*diff_y + alpha_*u*u);
80  }
81
82double d_cont(double x1, double x2, double y, double u)
83  {
84    return -exp(y) - u;
85  }
86
87
88using namespace Ipopt;
89
90/* Constructor. */
91MyADOLC_NLP::MyADOLC_NLP()
92{}
93
94MyADOLC_NLP::~MyADOLC_NLP(){}
95
96bool MyADOLC_NLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
97                         Index& nnz_h_lag, IndexStyleEnum& index_style)
98{
99  N_ = 10;
100  h_ = (1.0/(N_+1));
101  hh_= (h_*h_);
102
103  y_d_ = new double[(N_+2)*(N_+2)];
104  for (int j=0; j<= N_+1; j++) {
105    for (int i=0; i<= N_+1; i++) {
106      y_d_[y_index(i,j)] = y_d_cont(h_*1.*i,h_*1.*j);
107    }
108  }
109
110  n = (N_+2)*(N_+2) + N_*N_;
111
112  m = N_*N_ + 4*N_;
113
114  // in this example the jacobian is dense. Hence, it contains n*m nonzeros
115  nnz_jac_g = n*m;
116
117  // the hessian is also dense and has n*n total nonzeros, but we
118  // only need the lower left corner (since it is symmetric)
119  nnz_h_lag = n*(n-1)/2+n;
120
121  generate_tapes(n, m);
122
123  // use the C style indexing (0-based)
124  index_style = C_STYLE;
125
126  return true;
127}
128
129bool MyADOLC_NLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
130                            Index m, Number* g_l, Number* g_u)
131{
132  // none of the variables have bounds
133  for (Index i=0; i<n; i++) {
134    x_l[i] = -1e20;
135    x_u[i] =  1e20;
136  }
137
138  // Set the bounds for the constraints
139  for (Index i=0; i<m; i++) {
140    g_l[i] = 0;
141    g_u[i] = 0;
142  }
143
144  return true;
145}
146
147bool MyADOLC_NLP::get_starting_point(Index n, bool init_x, Number* x,
148                               bool init_z, Number* z_L, Number* z_U,
149                               Index m, bool init_lambda,
150                               Number* lambda)
151{
152  // Here, we assume we only have starting values for x, if you code
153  // your own NLP, you can provide starting values for the others if
154  // you wish.
155  assert(init_x == true);
156  assert(init_z == false);
157  assert(init_lambda == false);
158
159  // set all y's to the perfect match with y_d
160  for (Index i=0; i<= N_+1; i++) {
161    for (Index j=0; j<= N_+1; j++) {
162      x[y_index(i,j)] = y_d_[y_index(i,j)];
163      //x[y_index(i,j)] += h_*x1_grid(i) + 2*h_*x2_grid(j);
164    }
165  }
166
167  // Set the initial (constant) value for the u's
168  for (Index i=1; i<= N_; i++) {
169    for (Index j=1; j<= N_; j++) {
170      x[u_index(i,j)] = 0;
171    }
172  }
173
174  return true;
175}
176
177template<class T> bool  MyADOLC_NLP::eval_obj(Index n, const T *x, T& obj_value)
178{
179  // return the value of the objective function
180  obj_value = 0.;
181  for (int i=1; i<=N_; i++) {
182    for (int j=1; j<= N_; j++) {
183      int iy = y_index(i,j);
184      int iu = u_index(i,j);
185      obj_value += fint_cont(h_*1.*i, h_*1.*j, x[iy], x[iu]);
186    }
187  }
188  obj_value *= hh_;
189
190  return true;
191}
192
193template<class T> bool  MyADOLC_NLP::eval_constraints(Index n, const T *x, Index m, T* g)
194{
195  T val;
196  // compute the discretized PDE for each interior grid point
197  for (int i=1; i<=N_; i++) {
198    for (int j=1; j<=N_; j++) {
199
200      // Start with the discretized Laplacian operator
201      val = 4.* x[y_index(i,j)]
202            - x[y_index(i-1,j)] - x[y_index(i+1,j)]
203            - x[y_index(i,j-1)] - x[y_index(i,j+1)];
204
205      // Add the forcing term (including the step size here)
206      val += hh_*d_cont(h_*1.*i, h_*1.*j,x[y_index(i,j)], x[u_index(i,j)]);
207      g[pde_index(i,j)] = val;
208    }
209  }
210
211  int ig = N_*N_;
212  // set up the Neumann boundary conditions
213  for (int i=1; i<= N_; i++) {
214    g[ig] = (1.+h_*b_i0_)*x[y_index(i,0)] - x[y_index(i,1)];
215    ig++;
216  }
217  for (int i=1; i<= N_; i++) {
218    g[ig] = (1.+h_*b_i1_)*x[y_index(i,N_+1)] - x[y_index(i,N_)];
219    ig++;
220  }
221  for (int j=1; j<= N_; j++) {
222    g[ig] = (1.+h_*b_0j_)*x[y_index(0,j)] - x[y_index(1,j)];
223    ig++;
224  }
225  for (int j=1; j<= N_; j++) {
226    g[ig] = (1.+h_*b_1j_)*x[y_index(N_+1,j)] - x[y_index(N_,j)];
227    ig++;
228  } 
229   
230
231  return true;
232}
233
234//*************************************************************************
235//
236//
237//         Nothing has to be changed below this point !!
238//
239//
240//*************************************************************************
241
242
243bool MyADOLC_NLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
244{
245  eval_obj(n,x,obj_value);
246
247  return true;
248}
249
250bool MyADOLC_NLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
251{
252
253  gradient(tag_f,n,x,grad_f);
254
255  return true;
256}
257
258bool MyADOLC_NLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
259{
260
261  eval_constraints(n,x,m,g);
262
263  return true;
264}
265
266bool MyADOLC_NLP::eval_jac_g(Index n, const Number* x, bool new_x,
267                       Index m, Index nele_jac, Index* iRow, Index *jCol,
268                       Number* values)
269{
270  if (values == NULL) {
271    // return the structure of the jacobian,
272    // assuming that the Jacobian is dense
273
274    Index idx = 0;
275    for(Index i=0; i<m; i++)
276      for(Index j=0; j<n; j++)
277        {
278          iRow[idx] = i;
279          jCol[idx++] = j;
280        }
281 }
282  else {
283    // return the values of the jacobian of the constraints
284
285    jacobian(tag_g,m,n,x,Jac);
286
287    Index idx = 0;
288    for(Index i=0; i<m; i++)
289      for(Index j=0; j<n; j++)
290          values[idx++] = Jac[i][j];
291
292  }
293
294  return true;
295}
296
297bool MyADOLC_NLP::eval_h(Index n, const Number* x, bool new_x,
298                   Number obj_factor, Index m, const Number* lambda,
299                   bool new_lambda, Index nele_hess, Index* iRow,
300                   Index* jCol, Number* values)
301{
302  if (values == NULL) {
303    // return the structure. This is a symmetric matrix, fill the lower left
304    // triangle only.
305
306    // the hessian for this problem is actually dense
307    Index idx=0;
308    for (Index row = 0; row < n; row++) {
309      for (Index col = 0; col <= row; col++) {
310        iRow[idx] = row;
311        jCol[idx] = col;
312        idx++;
313      }
314    }
315
316    assert(idx == nele_hess);
317  }
318  else {
319    // return the values. This is a symmetric matrix, fill the lower left
320    // triangle only
321
322    for(Index i = 0; i<n ; i++)
323      x_lam[i] = x[i];
324    for(Index i = 0; i<m ; i++)
325      x_lam[n+i] = lambda[i];
326    x_lam[n+m] = obj_factor;
327
328    hessian(tag_L,n+m+1,x_lam,Hess);
329
330    Index idx = 0;
331
332    for(Index i = 0; i<n ; i++)
333      {
334        for(Index j = 0; j<=i ; j++)
335          {
336            values[idx++] = Hess[i][j];
337          }
338      }
339  }
340
341  return true;
342}
343
344void MyADOLC_NLP::finalize_solution(SolverReturn status,
345                              Index n, const Number* x, const Number* z_L, const Number* z_U,
346                              Index m, const Number* g, const Number* lambda,
347                              Number obj_value,
348                              const IpoptData* ip_data,
349                              IpoptCalculatedQuantities* ip_cq)
350{
351 
352  printf("\n\nObjective value\n");
353  printf("f(x*) = %e\n", obj_value);
354
355// Memory deallocation for ADOL-C variables
356
357  delete[] x_lam;
358
359  for(Index i=0;i<m;i++)
360    delete[] Jac[i];
361  delete[] Jac;
362
363  for(Index i=0;i<n+m+1;i++)
364    delete[] Hess[i];
365  delete[] Hess;
366
367}
368
369
370//***************    ADOL-C part ***********************************
371
372void MyADOLC_NLP::generate_tapes(Index n, Index m)
373{
374  Number *xp    = new double[n];
375  Number *lamp  = new double[m];
376  Number *zl    = new double[m];
377  Number *zu    = new double[m];
378
379  adouble *xa   = new adouble[n];
380  adouble *g    = new adouble[m];
381  adouble *lam  = new adouble[m];
382  adouble sig;
383  adouble obj_value;
384 
385  double dummy;
386
387  Jac = new double*[m];
388  for(Index i=0;i<m;i++)
389    Jac[i] = new double[n];
390
391  x_lam   = new double[n+m+1];
392
393  Hess = new double*[n+m+1];
394  for(Index i=0;i<n+m+1;i++)
395    Hess[i] = new double[i+1];
396
397  get_starting_point(n, 1, xp, 0, zl, zu, m, 0, lamp);
398
399  trace_on(tag_f);
400   
401    for(Index i=0;i<n;i++)
402      xa[i] <<= xp[i];
403
404    eval_obj(n,xa,obj_value);
405
406    obj_value >>= dummy;
407
408  trace_off();
409 
410  trace_on(tag_g);
411   
412    for(Index i=0;i<n;i++)
413      xa[i] <<= xp[i];
414
415    eval_constraints(n,xa,m,g);
416
417
418    for(Index i=0;i<m;i++)
419      g[i] >>= dummy;
420
421  trace_off();
422
423   trace_on(tag_L);
424   
425    for(Index i=0;i<n;i++)
426      xa[i] <<= xp[i];
427    for(Index i=0;i<m;i++)
428      lam[i] <<= 1.0;
429    sig <<= 1.0;
430
431    eval_obj(n,xa,obj_value);
432
433    obj_value *= sig;
434    eval_constraints(n,xa,m,g);
435 
436    for(Index i=0;i<m;i++)
437      obj_value += g[i]*lam[i];
438
439    obj_value >>= dummy;
440
441  trace_off();
442
443  delete[] xa;
444  delete[] xp;
445  delete[] g;
446  delete[] lam;
447  delete[] lamp;
448  delete[] zu;
449  delete[] zl;
450}
Note: See TracBrowser for help on using the repository browser.