source: trunk/ADOL-C/examples/additional_examples/ipopt/MittelmannDistCntrlNeumA_sparse/ADOL-C_sparseNLP.cpp @ 332

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

update of ipopt examples

File size: 15.7 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  jacval=NULL;
436  sparse_jac(tag_g, m, n, 0, xp, &nnz_jac, &rind_g, &cind_g, &jacval, options_g); 
437
438  options_g[2] = 0;
439  nnz_jac_g = nnz_jac;
440
441  unsigned int  **JP_f=NULL;                /* compressed block row storage */
442  unsigned int  **JP_g=NULL;                /* compressed block row storage */
443  unsigned int  **HP_f=NULL;                /* compressed block row storage */
444  unsigned int  **HP_g=NULL;                /* compressed block row storage */
445  unsigned int  *HP_length=NULL;            /* length of arrays */
446  unsigned int  *temp=NULL;                 /* help array */
447
448  int ctrl_H;
449
450  JP_f = (unsigned int **) malloc(sizeof(unsigned int*));
451  JP_g = (unsigned int **) malloc(m*sizeof(unsigned int*));
452  HP_f = (unsigned int **) malloc(n*sizeof(unsigned int*)); 
453  HP_g = (unsigned int **) malloc(n*sizeof(unsigned int*));
454  HP_t = (unsigned int **) malloc((n+m+1)*sizeof(unsigned int*));
455  HP_length = (unsigned int *) malloc((n)*sizeof(unsigned int));
456  ctrl_H = 0;
457
458  hess_pat(tag_f, n, xp, HP_f, ctrl_H);
459
460  indopro_forward_safe(tag_f, 1, n, xp, JP_f);
461  indopro_forward_safe(tag_g, m, n, xp, JP_g);
462  nonl_ind_forward_safe(tag_g, m, n, xp, HP_g);
463
464  for (i=0;i<n;i++) 
465    {
466      if (HP_f[i][0]+HP_g[i][0]!=0)
467        {
468          if (HP_f[i][0]==0)
469            {
470              HP_t[i] = (unsigned int *) malloc((HP_g[i][0]+HPOFF)*sizeof(unsigned int));
471              for(j=0;j<=(int) HP_g[i][0];j++)
472                {
473                  HP_t[i][j] = HP_g[i][j];
474                }
475              HP_length[i] = HP_g[i][0]+HPOFF;
476            }
477          else
478            {
479              if (HP_g[i][0]==0)
480                {
481                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HPOFF)*sizeof(unsigned int));
482                  for(j=0;j<=(int) HP_f[i][0];j++)
483                    {
484                      HP_t[i][j] = HP_f[i][j];
485                    }
486                  HP_length[i] = HP_f[i][0]+HPOFF;
487                }
488              else
489                {
490                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HP_g[i][0]+HPOFF)*sizeof(unsigned int));
491                  k = l = j = 1;
492                  while ((k<=(int) HP_f[i][0]) && (l <= (int) HP_g[i][0]))
493                    {
494                      if (HP_f[i][k] < HP_g[i][l])
495                        {
496                          HP_t[i][j]=HP_f[i][k];
497                          j++; k++;
498                        }
499                      else
500                        {
501                          if (HP_f[i][k] == HP_g[i][l])
502                            {
503                              HP_t[i][j]=HP_f[i][k];
504                              l++;j++;k++;
505                            }
506                          else
507                            {
508                              HP_t[i][j]=HP_g[i][l];
509                              j++;l++;               
510                            }
511                        }
512                    } // end while
513                  for(ii=k;ii<=(int) HP_f[i][0];ii++)
514                    {
515                      HP_t[i][j] = HP_f[i][ii];
516                      j++;
517                    }
518                  for(ii=l;ii<=(int) HP_g[i][0];ii++)
519                    {
520                      HP_t[i][j] = HP_g[i][ii];
521                      j++;
522                    }
523                 
524                }
525            }
526          HP_t[i][0]=j-1;
527          HP_length[i] = HP_f[i][0]+HP_g[i][0]+HPOFF;
528        }
529      else
530        {
531          HP_t[i] = (unsigned int *) malloc((HPOFF+1)*sizeof(unsigned int));
532          HP_t[i][0]=0;
533          HP_length[i]=HPOFF;
534        }
535    }   
536
537  for (i=0;i<m;i++) 
538    {
539      HP_t[n+i] = (unsigned int *) malloc((JP_g[i][0]+1)*sizeof(unsigned int));
540      HP_t[n+i][0]=JP_g[i][0];
541      for(j=1;j<= (int) JP_g[i][0];j++)
542        {
543          HP_t[n+i][j]=JP_g[i][j];
544          if (HP_length[JP_g[i][j]] < HP_t[JP_g[i][j]][0]+1)
545            {
546              temp = (unsigned int *) malloc((HP_t[JP_g[i][j]][0])*sizeof(unsigned int));
547              for(l=0;l<=(int)HP_t[JP_g[i][j]][0];l++)
548                temp[l] = HP_t[JP_g[i][j]][l];
549              free(HP_t[JP_g[i][j]]);
550              HP_t[JP_g[i][j]] = (unsigned int *) malloc(2*HP_length[JP_g[i][j]]*sizeof(unsigned int));
551              HP_length[JP_g[i][j]] = 2*HP_length[JP_g[i][j]];
552              for(l=0;l<=(int)temp[0];l++)
553                HP_t[JP_g[i][j]][l] =temp[l];
554              free(temp);
555            }
556          HP_t[JP_g[i][j]][0] = HP_t[JP_g[i][j]][0]+1;
557          HP_t[JP_g[i][j]][HP_t[JP_g[i][j]][0]] = i+n;
558        }
559    }
560
561
562  for(j=1;j<= (int) JP_f[0][0];j++)
563    {
564      if (HP_length[JP_f[0][j]] < HP_t[JP_f[0][j]][0]+1)
565        {
566          temp = (unsigned int *) malloc((HP_t[JP_f[0][j]][0])*sizeof(unsigned int));
567          for(l=0;l<=(int)HP_t[JP_f[0][j]][0];l++)
568            temp[l] = HP_t[JP_f[0][j]][l];
569          free(HP_t[JP_f[0][j]]);
570          HP_t[JP_f[0][j]] = (unsigned int *) malloc(2*HP_length[JP_f[0][j]]*sizeof(unsigned int));
571          HP_length[JP_f[0][j]] = 2*HP_length[JP_f[0][j]];
572          for(l=0;l<=(int)temp[0];l++)
573            HP_t[JP_f[0][j]][l] =temp[l];
574          free(temp);
575        }
576      HP_t[JP_f[0][j]][0] = HP_t[JP_f[0][j]][0]+1;
577      HP_t[JP_f[0][j]][HP_t[JP_f[0][j]][0]] = n+m;
578    }
579
580
581  HP_t[n+m] = (unsigned int *) malloc((JP_f[0][0]+2)*sizeof(unsigned int));
582  HP_t[n+m][0]=JP_f[0][0]+1;
583  for(j=1;j<= (int) JP_f[0][0];j++)
584    HP_t[n+m][j]=JP_f[0][j];
585  HP_t[n+m][JP_f[0][0]+1]=n+m;
586
587
588
589  set_HP(tag_L,n+m+1,HP_t);
590
591  nnz_h_lag = 0;
592   for (i=0;i<n;i++) {
593    for (j=1;j<=(int) HP_t[i][0];j++)
594      if ((int) HP_t[i][j] <= i)
595        nnz_h_lag++;
596     free(HP_f[i]);
597     free(HP_g[i]);
598   }
599  nnz_L = nnz_h_lag;
600
601  options_L[0] = 0;         
602  options_L[1] = 1;       
603
604  sparse_hess(tag_L, n+m+1, -1, xp, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);
605
606  rind_L = new unsigned int[nnz_L];
607  cind_L = new unsigned int[nnz_L];
608  rind_L_total = new unsigned int[nnz_L_total];
609  cind_L_total = new unsigned int[nnz_L_total];
610
611  unsigned int ind = 0;
612
613  for (int i=0;i<n;i++) 
614    for (unsigned int j=1;j<=HP_t[i][0];j++)
615      {
616        if (((int) HP_t[i][j]>=i) &&((int) HP_t[i][j]<n)) 
617          {
618            rind_L[ind] = i;
619            cind_L[ind++] = HP_t[i][j];
620          }
621      }
622
623   ind = 0;
624   for (int i=0;i<n+m+1;i++) 
625     for (unsigned int j=1;j<=HP_t[i][0];j++)
626       {
627        if ((int) HP_t[i][j]>=i) 
628          {
629            rind_L_total[ind] = i;
630            cind_L_total[ind++] = HP_t[i][j];
631          }
632       }
633
634  for (i=0;i<m;i++) {
635     free(JP_g[i]);
636   }
637
638  free(JP_f[0]);
639  free(JP_f);
640  free(JP_g);
641  free(HP_f);
642  free(HP_g);
643  free(HP_length);
644
645  delete[] lam;
646  delete[] g;
647  delete[] xa;
648  delete[] zu;
649  delete[] zl;
650  delete[] lamp;
651  delete[] xp;
652}
Note: See TracBrowser for help on using the repository browser.