source: stable/2.4/ADOL-C/examples/additional_examples/ipopt/MittelmannDistCntrlNeumA_sparse/ADOL-C_sparseNLP.cpp

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

set derivative values array pointer to NULL in MittelmannDistCntrlNeumA_sparse

File size: 15.8 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     ADOL-C_sparseNLP.cpp
4 Revision: $$
5 Contents: class myADOLC_sparseNPL 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_sparseNLP 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" taking sparsity into account.
25 *
26 *  exploitation of sparsity !!
27 *
28 */
29
30#include <cassert>
31
32#include "ADOL-C_sparseNLP.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
87using namespace Ipopt;
88
89/* Constructor. */
90MyADOLC_sparseNLP::MyADOLC_sparseNLP()
91{}
92
93MyADOLC_sparseNLP::~MyADOLC_sparseNLP()
94{}
95
96bool MyADOLC_sparseNLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
97                         Index& nnz_h_lag, IndexStyleEnum& index_style)
98{
99  N_ = 15;
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  generate_tapes(n, m, nnz_jac_g, nnz_h_lag);
115
116  // use the C style indexing (0-based)
117  index_style = C_STYLE;
118
119  return true;
120}
121
122bool MyADOLC_sparseNLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
123                            Index m, Number* g_l, Number* g_u)
124{
125  // none of the variables have bounds
126  for (Index i=0; i<n; i++) {
127    x_l[i] = -1e20;
128    x_u[i] =  1e20;
129  }
130
131  // Set the bounds for the constraints
132  for (Index i=0; i<m; i++) {
133    g_l[i] = 0;
134    g_u[i] = 0;
135  }
136
137  return true;
138}
139
140bool MyADOLC_sparseNLP::get_starting_point(Index n, bool init_x, Number* x,
141                               bool init_z, Number* z_L, Number* z_U,
142                               Index m, bool init_lambda,
143                               Number* lambda)
144{
145  // Here, we assume we only have starting values for x, if you code
146  // your own NLP, you can provide starting values for the others if
147  // you wish.
148  assert(init_x == true);
149  assert(init_z == false);
150  assert(init_lambda == false);
151
152  // set all y's to the perfect match with y_d
153  for (Index i=0; i<= N_+1; i++) {
154    for (Index j=0; j<= N_+1; j++) {
155      x[y_index(i,j)] = y_d_[y_index(i,j)];
156      //x[y_index(i,j)] += h_*x1_grid(i) + 2*h_*x2_grid(j);
157    }
158  }
159
160  // Set the initial (constant) value for the u's
161  for (Index i=1; i<= N_; i++) {
162    for (Index j=1; j<= N_; j++) {
163      x[u_index(i,j)] = 0;
164    }
165  }
166
167
168  return true;
169}
170
171template<class T> bool  MyADOLC_sparseNLP::eval_obj(Index n, const T *x, T& obj_value)
172{
173  // return the value of the objective function
174  obj_value = 0.;
175  for (int i=1; i<=N_; i++) {
176    for (int j=1; j<= N_; j++) {
177      int iy = y_index(i,j);
178      int iu = u_index(i,j);
179      obj_value += fint_cont(h_*1.*i, h_*1.*j, x[iy], x[iu]);
180    }
181  }
182  obj_value *= hh_;
183
184  return true;
185}
186
187template<class T> bool  MyADOLC_sparseNLP::eval_constraints(Index n, const T *x, Index m, T* g)
188{
189  T val;
190  // compute the discretized PDE for each interior grid point
191  for (int i=1; i<=N_; i++) {
192    for (int j=1; j<=N_; j++) {
193
194      // Start with the discretized Laplacian operator
195      val = 4.* x[y_index(i,j)]
196            - x[y_index(i-1,j)] - x[y_index(i+1,j)]
197            - x[y_index(i,j-1)] - x[y_index(i,j+1)];
198
199      // Add the forcing term (including the step size here)
200      val += hh_*d_cont(h_*1.*i, h_*1.*j,x[y_index(i,j)], x[u_index(i,j)]);
201      g[pde_index(i,j)] = val;
202    }
203  }
204
205  int ig = N_*N_;
206  // set up the Neumann boundary conditions
207  for (int i=1; i<= N_; i++) {
208    g[ig] = (1.+h_*b_i0_)*x[y_index(i,0)] - x[y_index(i,1)];
209    ig++;
210  }
211  for (int i=1; i<= N_; i++) {
212    g[ig] = (1.+h_*b_i1_)*x[y_index(i,N_+1)] - x[y_index(i,N_)];
213    ig++;
214  }
215  for (int j=1; j<= N_; j++) {
216    g[ig] = (1.+h_*b_0j_)*x[y_index(0,j)] - x[y_index(1,j)];
217    ig++;
218  }
219  for (int j=1; j<= N_; j++) {
220    g[ig] = (1.+h_*b_1j_)*x[y_index(N_+1,j)] - x[y_index(N_,j)];
221    ig++;
222  } 
223   
224  return true;
225}
226
227//*************************************************************************
228//
229//
230//         Nothing has to be changed below this point !!
231//
232//
233//*************************************************************************
234
235
236bool MyADOLC_sparseNLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
237{
238  eval_obj(n,x,obj_value);
239
240  return true;
241}
242
243bool MyADOLC_sparseNLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
244{
245
246  gradient(tag_f,n,x,grad_f);
247
248  return true;
249}
250
251bool MyADOLC_sparseNLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
252{
253
254  eval_constraints(n,x,m,g);
255
256  return true;
257}
258
259bool MyADOLC_sparseNLP::eval_jac_g(Index n, const Number* x, bool new_x,
260                       Index m, Index nele_jac, Index* iRow, Index *jCol,
261                       Number* values)
262{
263
264  if (values == NULL) {
265    // return the structure of the jacobian
266
267    for(Index idx=0; idx<nnz_jac; idx++)
268      {
269        iRow[idx] = rind_g[idx];
270        jCol[idx] = cind_g[idx];
271      }
272  }
273  else {
274    // return the values of the jacobian of the constraints
275
276    sparse_jac(tag_g, m, n, 1, x, &nnz_jac, &rind_g, &cind_g, &jacval, options_g); 
277
278    for(Index idx=0; idx<nnz_jac; idx++)
279      {
280        values[idx] = jacval[idx];
281
282      }
283  }
284  return true;
285}
286
287bool MyADOLC_sparseNLP::eval_h(Index n, const Number* x, bool new_x,
288                   Number obj_factor, Index m, const Number* lambda,
289                   bool new_lambda, Index nele_hess, Index* iRow,
290                   Index* jCol, Number* values)
291{
292
293  if (values == NULL) {
294    // return the structure. This is a symmetric matrix, fill the lower left
295    // triangle only.
296
297    for(Index idx=0; idx<nnz_L; idx++)
298      {
299        iRow[idx] = rind_L[idx];
300        jCol[idx] = cind_L[idx];
301      }
302  }
303  else {
304    // return the values. This is a symmetric matrix, fill the lower left
305    // triangle only
306
307    for(Index idx = 0; idx<n ; idx++)
308      x_lam[idx] = x[idx];
309    for(Index idx = 0; idx<m ; idx++)
310      x_lam[n+idx] = lambda[idx];
311    x_lam[n+m] = obj_factor;
312
313    sparse_hess(tag_L, n+m+1, 1, x_lam, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);
314     
315    Index idx = 0;
316    for(Index idx_total = 0; idx_total <nnz_L_total ; idx_total++)
317      {
318        if((rind_L_total[idx_total] < (unsigned int) n) && (cind_L_total[idx_total] < (unsigned int) n))
319          {
320            values[idx] = hessval[idx_total];
321            idx++;
322          }
323      }
324  }
325
326  return true;
327}
328
329void MyADOLC_sparseNLP::finalize_solution(SolverReturn status,
330                              Index n, const Number* x, const Number* z_L, const Number* z_U,
331                              Index m, const Number* g, const Number* lambda,
332                              Number obj_value,
333                              const IpoptData* ip_data,
334                              IpoptCalculatedQuantities* ip_cq)
335{
336
337  printf("\n\nObjective value\n");
338  printf("f(x*) = %e\n", obj_value);
339
340// memory deallocation of ADOL-C variables
341
342  delete x_lam;
343  delete rind_g;
344  delete cind_g;
345  delete rind_L;
346  delete cind_L;
347  delete rind_L_total;
348  delete cind_L_total;
349  delete jacval;
350  delete hessval;
351
352  for (int i=0;i<n+m+1;i++) {
353     free(HP_t[i]);
354   }
355  free(HP_t);
356}
357
358
359//***************    ADOL-C part ***********************************
360
361void MyADOLC_sparseNLP::generate_tapes(Index n, Index m, Index& nnz_jac_g, Index& nnz_h_lag)
362{
363  Number *xp    = new double[n];
364  Number *lamp  = new double[m];
365  Number *zl    = new double[m];
366  Number *zu    = new double[m];
367
368  adouble *xa   = new adouble[n];
369  adouble *g    = new adouble[m];
370  adouble *lam  = new adouble[m];
371  adouble sig;
372  adouble obj_value;
373 
374  double dummy;
375  double *jacval;
376
377  int i,j,k,l,ii;
378
379  x_lam   = new double[n+m+1];
380
381  get_starting_point(n, 1, xp, 0, zl, zu, m, 0, lamp);
382
383  trace_on(tag_f);
384   
385    for(Index idx=0;idx<n;idx++)
386      xa[idx] <<= xp[idx];
387
388    eval_obj(n,xa,obj_value);
389
390    obj_value >>= dummy;
391
392  trace_off();
393 
394  trace_on(tag_g);
395   
396    for(Index idx=0;idx<n;idx++)
397      xa[idx] <<= xp[idx];
398
399    eval_constraints(n,xa,m,g);
400
401
402    for(Index idx=0;idx<m;idx++)
403      g[idx] >>= dummy;
404
405  trace_off();
406
407  trace_on(tag_L);
408   
409    for(Index idx=0;idx<n;idx++)
410      xa[idx] <<= xp[idx];
411    for(Index idx=0;idx<m;idx++)
412      lam[idx] <<= 1.0;
413    sig <<= 1.0;
414
415    eval_obj(n,xa,obj_value);
416
417    obj_value *= sig;
418    eval_constraints(n,xa,m,g);
419 
420    for(Index idx=0;idx<m;idx++)
421      obj_value += g[idx]*lam[idx];
422
423    obj_value >>= dummy;
424
425  trace_off();
426
427  rind_g = NULL; 
428  cind_g = NULL;
429
430  options_g[0] = 0;          /* sparsity pattern by index domains (default) */ 
431  options_g[1] = 0;          /*                         safe mode (default) */ 
432  options_g[2] = -1;         /*                     &jacval is not computed */ 
433  options_g[3] = 0;          /*                column compression (default) */ 
434 
435  this->jacval=NULL;
436  this->hessval=NULL;
437  sparse_jac(tag_g, m, n, 0, xp, &nnz_jac, &rind_g, &cind_g, &jacval, options_g); 
438
439  options_g[2] = 0;
440  nnz_jac_g = nnz_jac;
441
442  unsigned int  **JP_f=NULL;                /* compressed block row storage */
443  unsigned int  **JP_g=NULL;                /* compressed block row storage */
444  unsigned int  **HP_f=NULL;                /* compressed block row storage */
445  unsigned int  **HP_g=NULL;                /* compressed block row storage */
446  unsigned int  *HP_length=NULL;            /* length of arrays */
447  unsigned int  *temp=NULL;                 /* help array */
448
449  int ctrl_H;
450
451  JP_f = (unsigned int **) malloc(sizeof(unsigned int*));
452  JP_g = (unsigned int **) malloc(m*sizeof(unsigned int*));
453  HP_f = (unsigned int **) malloc(n*sizeof(unsigned int*)); 
454  HP_g = (unsigned int **) malloc(n*sizeof(unsigned int*));
455  HP_t = (unsigned int **) malloc((n+m+1)*sizeof(unsigned int*));
456  HP_length = (unsigned int *) malloc((n)*sizeof(unsigned int));
457  ctrl_H = 0;
458
459  hess_pat(tag_f, n, xp, HP_f, ctrl_H);
460
461  indopro_forward_safe(tag_f, 1, n, xp, JP_f);
462  indopro_forward_safe(tag_g, m, n, xp, JP_g);
463  nonl_ind_forward_safe(tag_g, m, n, xp, HP_g);
464
465  for (i=0;i<n;i++) 
466    {
467      if (HP_f[i][0]+HP_g[i][0]!=0)
468        {
469          if (HP_f[i][0]==0)
470            {
471              HP_t[i] = (unsigned int *) malloc((HP_g[i][0]+HPOFF)*sizeof(unsigned int));
472              for(j=0;j<=(int) HP_g[i][0];j++)
473                {
474                  HP_t[i][j] = HP_g[i][j];
475                }
476              HP_length[i] = HP_g[i][0]+HPOFF;
477            }
478          else
479            {
480              if (HP_g[i][0]==0)
481                {
482                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HPOFF)*sizeof(unsigned int));
483                  for(j=0;j<=(int) HP_f[i][0];j++)
484                    {
485                      HP_t[i][j] = HP_f[i][j];
486                    }
487                  HP_length[i] = HP_f[i][0]+HPOFF;
488                }
489              else
490                {
491                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HP_g[i][0]+HPOFF)*sizeof(unsigned int));
492                  k = l = j = 1;
493                  while ((k<=(int) HP_f[i][0]) && (l <= (int) HP_g[i][0]))
494                    {
495                      if (HP_f[i][k] < HP_g[i][l])
496                        {
497                          HP_t[i][j]=HP_f[i][k];
498                          j++; k++;
499                        }
500                      else
501                        {
502                          if (HP_f[i][k] == HP_g[i][l])
503                            {
504                              HP_t[i][j]=HP_f[i][k];
505                              l++;j++;k++;
506                            }
507                          else
508                            {
509                              HP_t[i][j]=HP_g[i][l];
510                              j++;l++;               
511                            }
512                        }
513                    } // end while
514                  for(ii=k;ii<=(int) HP_f[i][0];ii++)
515                    {
516                      HP_t[i][j] = HP_f[i][ii];
517                      j++;
518                    }
519                  for(ii=l;ii<=(int) HP_g[i][0];ii++)
520                    {
521                      HP_t[i][j] = HP_g[i][ii];
522                      j++;
523                    }
524                 
525                }
526            }
527          HP_t[i][0]=j-1;
528          HP_length[i] = HP_f[i][0]+HP_g[i][0]+HPOFF;
529        }
530      else
531        {
532          HP_t[i] = (unsigned int *) malloc((HPOFF+1)*sizeof(unsigned int));
533          HP_t[i][0]=0;
534          HP_length[i]=HPOFF;
535        }
536    }   
537
538  for (i=0;i<m;i++) 
539    {
540      HP_t[n+i] = (unsigned int *) malloc((JP_g[i][0]+1)*sizeof(unsigned int));
541      HP_t[n+i][0]=JP_g[i][0];
542      for(j=1;j<= (int) JP_g[i][0];j++)
543        {
544          HP_t[n+i][j]=JP_g[i][j];
545          if (HP_length[JP_g[i][j]] < HP_t[JP_g[i][j]][0]+1)
546            {
547              temp = (unsigned int *) malloc((HP_t[JP_g[i][j]][0])*sizeof(unsigned int));
548              for(l=0;l<=(int)HP_t[JP_g[i][j]][0];l++)
549                temp[l] = HP_t[JP_g[i][j]][l];
550              free(HP_t[JP_g[i][j]]);
551              HP_t[JP_g[i][j]] = (unsigned int *) malloc(2*HP_length[JP_g[i][j]]*sizeof(unsigned int));
552              HP_length[JP_g[i][j]] = 2*HP_length[JP_g[i][j]];
553              for(l=0;l<=(int)temp[0];l++)
554                HP_t[JP_g[i][j]][l] =temp[l];
555              free(temp);
556            }
557          HP_t[JP_g[i][j]][0] = HP_t[JP_g[i][j]][0]+1;
558          HP_t[JP_g[i][j]][HP_t[JP_g[i][j]][0]] = i+n;
559        }
560    }
561
562
563  for(j=1;j<= (int) JP_f[0][0];j++)
564    {
565      if (HP_length[JP_f[0][j]] < HP_t[JP_f[0][j]][0]+1)
566        {
567          temp = (unsigned int *) malloc((HP_t[JP_f[0][j]][0])*sizeof(unsigned int));
568          for(l=0;l<=(int)HP_t[JP_f[0][j]][0];l++)
569            temp[l] = HP_t[JP_f[0][j]][l];
570          free(HP_t[JP_f[0][j]]);
571          HP_t[JP_f[0][j]] = (unsigned int *) malloc(2*HP_length[JP_f[0][j]]*sizeof(unsigned int));
572          HP_length[JP_f[0][j]] = 2*HP_length[JP_f[0][j]];
573          for(l=0;l<=(int)temp[0];l++)
574            HP_t[JP_f[0][j]][l] =temp[l];
575          free(temp);
576        }
577      HP_t[JP_f[0][j]][0] = HP_t[JP_f[0][j]][0]+1;
578      HP_t[JP_f[0][j]][HP_t[JP_f[0][j]][0]] = n+m;
579    }
580
581
582  HP_t[n+m] = (unsigned int *) malloc((JP_f[0][0]+2)*sizeof(unsigned int));
583  HP_t[n+m][0]=JP_f[0][0]+1;
584  for(j=1;j<= (int) JP_f[0][0];j++)
585    HP_t[n+m][j]=JP_f[0][j];
586  HP_t[n+m][JP_f[0][0]+1]=n+m;
587
588
589
590  set_HP(tag_L,n+m+1,HP_t);
591
592  nnz_h_lag = 0;
593   for (i=0;i<n;i++) {
594    for (j=1;j<=(int) HP_t[i][0];j++)
595      if ((int) HP_t[i][j] <= i)
596        nnz_h_lag++;
597     free(HP_f[i]);
598     free(HP_g[i]);
599   }
600  nnz_L = nnz_h_lag;
601
602  options_L[0] = 0;         
603  options_L[1] = 1;       
604
605  sparse_hess(tag_L, n+m+1, -1, xp, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);
606
607  rind_L = new unsigned int[nnz_L];
608  cind_L = new unsigned int[nnz_L];
609  rind_L_total = new unsigned int[nnz_L_total];
610  cind_L_total = new unsigned int[nnz_L_total];
611
612  unsigned int ind = 0;
613
614  for (int i=0;i<n;i++) 
615    for (unsigned int j=1;j<=HP_t[i][0];j++)
616      {
617        if (((int) HP_t[i][j]>=i) &&((int) HP_t[i][j]<n)) 
618          {
619            rind_L[ind] = i;
620            cind_L[ind++] = HP_t[i][j];
621          }
622      }
623
624   ind = 0;
625   for (int i=0;i<n+m+1;i++) 
626     for (unsigned int j=1;j<=HP_t[i][0];j++)
627       {
628        if ((int) HP_t[i][j]>=i) 
629          {
630            rind_L_total[ind] = i;
631            cind_L_total[ind++] = HP_t[i][j];
632          }
633       }
634
635  for (i=0;i<m;i++) {
636     free(JP_g[i]);
637   }
638
639  free(JP_f[0]);
640  free(JP_f);
641  free(JP_g);
642  free(HP_f);
643  free(HP_g);
644  free(HP_length);
645
646  delete[] lam;
647  delete[] g;
648  delete[] xa;
649  delete[] zu;
650  delete[] zl;
651  delete[] lamp;
652  delete[] xp;
653}
Note: See TracBrowser for help on using the repository browser.