source: trunk/ADOL-C/examples/additional_examples/ipopt/LuksanVlcek1_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: 13.3 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 *  MyADOLC_sparseNLP 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 taking sparsity
23 *  into account.
24 *
25 *  exploitation of sparsity !!
26 *
27 */
28
29#include <cassert>
30
31#include "ADOL-C_sparseNLP.hpp"
32
33using namespace Ipopt;
34
35/* Constructor. */
36MyADOLC_sparseNLP::MyADOLC_sparseNLP()
37{}
38
39MyADOLC_sparseNLP::~MyADOLC_sparseNLP()
40{}
41
42bool MyADOLC_sparseNLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
43                         Index& nnz_h_lag, IndexStyleEnum& index_style)
44{
45  n = 4000;
46
47  m = n-2;
48
49  generate_tapes(n, m, nnz_jac_g, nnz_h_lag);
50
51  // use the C style indexing (0-based)
52  index_style = C_STYLE;
53
54  return true;
55}
56
57bool MyADOLC_sparseNLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
58                            Index m, Number* g_l, Number* g_u)
59{
60  // none of the variables have bounds
61  for (Index i=0; i<n; i++) {
62    x_l[i] = -1e20;
63    x_u[i] =  1e20;
64  }
65
66  // Set the bounds for the constraints
67  for (Index i=0; i<m; i++) {
68    g_l[i] = 0;
69    g_u[i] = 0;
70  }
71
72  return true;
73}
74
75bool MyADOLC_sparseNLP::get_starting_point(Index n, bool init_x, Number* x,
76                               bool init_z, Number* z_L, Number* z_U,
77                               Index m, bool init_lambda,
78                               Number* lambda)
79{
80  // Here, we assume we only have starting values for x, if you code
81  // your own NLP, you can provide starting values for the others if
82  // you wish.
83  assert(init_x == true);
84  assert(init_z == false);
85  assert(init_lambda == false);
86
87  // set the starting point
88  for (Index i=0; i<n/2; i++) {
89    x[2*i] = -1.2;
90    x[2*i+1] = 1.;
91  }
92  if (n != 2*(n/2)) {
93    x[n-1] = -1.2;
94  }
95
96  return true;
97}
98
99template<class T> bool  MyADOLC_sparseNLP::eval_obj(Index n, const T *x, T& obj_value)
100{
101  T a1, a2;
102  obj_value = 0.;
103  for (Index i=0; i<n-1; i++) {
104    a1 = x[i]*x[i]-x[i+1];
105    a2 = x[i] - 1.;
106    obj_value += 100.*a1*a1 + a2*a2;
107  }
108
109  return true;
110}
111
112template<class T> bool  MyADOLC_sparseNLP::eval_constraints(Index n, const T *x, Index m, T* g)
113{
114  for (Index i=0; i<m; i++) {
115    g[i] = 3.*pow(x[i+1],3.) + 2.*x[i+2] - 5.
116           + sin(x[i+1]-x[i+2])*sin(x[i+1]+x[i+2]) + 4.*x[i+1]
117           - x[i]*exp(x[i]-x[i+1]) - 3.;
118  }
119
120  return true;
121}
122
123//*************************************************************************
124//
125//
126//         Nothing has to be changed below this point !!
127//
128//
129//*************************************************************************
130
131
132bool MyADOLC_sparseNLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
133{
134  eval_obj(n,x,obj_value);
135
136  return true;
137}
138
139bool MyADOLC_sparseNLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
140{
141
142  gradient(tag_f,n,x,grad_f);
143
144  return true;
145}
146
147bool MyADOLC_sparseNLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
148{
149
150  eval_constraints(n,x,m,g);
151
152  return true;
153}
154
155bool MyADOLC_sparseNLP::eval_jac_g(Index n, const Number* x, bool new_x,
156                       Index m, Index nele_jac, Index* iRow, Index *jCol,
157                       Number* values)
158{
159
160  if (values == NULL) {
161    // return the structure of the jacobian
162
163    for(Index idx=0; idx<nnz_jac; idx++)
164      {
165        iRow[idx] = rind_g[idx];
166        jCol[idx] = cind_g[idx];
167      }
168  }
169  else {
170    // return the values of the jacobian of the constraints
171
172    sparse_jac(tag_g, m, n, 1, x, &nnz_jac, &rind_g, &cind_g, &jacval, options_g); 
173
174    for(Index idx=0; idx<nnz_jac; idx++)
175      {
176        values[idx] = jacval[idx];
177
178      }
179  }
180  return true;
181}
182
183bool MyADOLC_sparseNLP::eval_h(Index n, const Number* x, bool new_x,
184                   Number obj_factor, Index m, const Number* lambda,
185                   bool new_lambda, Index nele_hess, Index* iRow,
186                   Index* jCol, Number* values)
187{
188
189  if (values == NULL) {
190    // return the structure. This is a symmetric matrix, fill the lower left
191    // triangle only.
192
193    for(Index idx=0; idx<nnz_L; idx++)
194      {
195        iRow[idx] = rind_L[idx];
196        jCol[idx] = cind_L[idx];
197      }
198  }
199  else {
200    // return the values. This is a symmetric matrix, fill the lower left
201    // triangle only
202
203    for(Index idx = 0; idx<n ; idx++)
204      x_lam[idx] = x[idx];
205    for(Index idx = 0; idx<m ; idx++)
206      x_lam[n+idx] = lambda[idx];
207    x_lam[n+m] = obj_factor;
208
209    sparse_hess(tag_L, n+m+1, 1, x_lam, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);
210     
211    Index idx = 0;
212    for(Index idx_total = 0; idx_total <nnz_L_total ; idx_total++)
213      {
214        if((rind_L_total[idx_total] < (unsigned int) n) && (cind_L_total[idx_total] < (unsigned int) n))
215          {
216            values[idx] = hessval[idx_total];
217            idx++;
218          }
219      }
220  }
221
222  return true;
223}
224
225void MyADOLC_sparseNLP::finalize_solution(SolverReturn status,
226                              Index n, const Number* x, const Number* z_L, const Number* z_U,
227                              Index m, const Number* g, const Number* lambda,
228                              Number obj_value,
229                              const IpoptData* ip_data,
230                              IpoptCalculatedQuantities* ip_cq)
231{
232
233  printf("\n\nObjective value\n");
234  printf("f(x*) = %e\n", obj_value);
235
236// memory deallocation of ADOL-C variables
237
238  delete x_lam;
239  delete rind_g;
240  delete cind_g;
241  delete rind_L;
242  delete cind_L;
243  delete rind_L_total;
244  delete cind_L_total;
245  delete jacval;
246  delete hessval;
247
248  for (int i=0;i<n+m+1;i++) {
249     free(HP_t[i]);
250   }
251  free(HP_t);
252}
253
254
255//***************    ADOL-C part ***********************************
256
257void MyADOLC_sparseNLP::generate_tapes(Index n, Index m, Index& nnz_jac_g, Index& nnz_h_lag)
258{
259  Number *xp    = new double[n];
260  Number *lamp  = new double[m];
261  Number *zl    = new double[m];
262  Number *zu    = new double[m];
263
264  adouble *xa   = new adouble[n];
265  adouble *g    = new adouble[m];
266  adouble *lam  = new adouble[m];
267  adouble sig;
268  adouble obj_value;
269 
270  double dummy;
271  double *jacval;
272
273  int i,j,k,l,ii;
274
275  x_lam   = new double[n+m+1];
276
277  get_starting_point(n, 1, xp, 0, zl, zu, m, 0, lamp);
278
279  trace_on(tag_f);
280   
281    for(Index idx=0;idx<n;idx++)
282      xa[idx] <<= xp[idx];
283
284    eval_obj(n,xa,obj_value);
285
286    obj_value >>= dummy;
287
288  trace_off();
289 
290  trace_on(tag_g);
291   
292    for(Index idx=0;idx<n;idx++)
293      xa[idx] <<= xp[idx];
294
295    eval_constraints(n,xa,m,g);
296
297
298    for(Index idx=0;idx<m;idx++)
299      g[idx] >>= dummy;
300
301  trace_off();
302
303  trace_on(tag_L);
304   
305    for(Index idx=0;idx<n;idx++)
306      xa[idx] <<= xp[idx];
307    for(Index idx=0;idx<m;idx++)
308      lam[idx] <<= 1.0;
309    sig <<= 1.0;
310
311    eval_obj(n,xa,obj_value);
312
313    obj_value *= sig;
314    eval_constraints(n,xa,m,g);
315 
316    for(Index idx=0;idx<m;idx++)
317      obj_value += g[idx]*lam[idx];
318
319    obj_value >>= dummy;
320
321  trace_off();
322
323  rind_g = NULL; 
324  cind_g = NULL;
325
326  options_g[0] = 0;          /* sparsity pattern by index domains (default) */ 
327  options_g[1] = 0;          /*                         safe mode (default) */ 
328  options_g[2] = -1;         /*                     &jacval is not computed */ 
329  options_g[3] = 0;          /*                column compression (default) */ 
330 
331  jacval=NULL;
332  sparse_jac(tag_g, m, n, 0, xp, &nnz_jac, &rind_g, &cind_g, &jacval, options_g); 
333
334  options_g[2] = 0;
335  nnz_jac_g = nnz_jac;
336
337  unsigned int  **JP_f=NULL;                /* compressed block row storage */
338  unsigned int  **JP_g=NULL;                /* compressed block row storage */
339  unsigned int  **HP_f=NULL;                /* compressed block row storage */
340  unsigned int  **HP_g=NULL;                /* compressed block row storage */
341  unsigned int  *HP_length=NULL;            /* length of arrays */
342  unsigned int  *temp=NULL;                 /* help array */
343
344  int ctrl_H;
345
346  JP_f = (unsigned int **) malloc(sizeof(unsigned int*));
347  JP_g = (unsigned int **) malloc(m*sizeof(unsigned int*));
348  HP_f = (unsigned int **) malloc(n*sizeof(unsigned int*)); 
349  HP_g = (unsigned int **) malloc(n*sizeof(unsigned int*));
350  HP_t = (unsigned int **) malloc((n+m+1)*sizeof(unsigned int*));
351  HP_length = (unsigned int *) malloc((n)*sizeof(unsigned int));
352  ctrl_H = 0;
353
354  hess_pat(tag_f, n, xp, HP_f, ctrl_H);
355
356  indopro_forward_safe(tag_f, 1, n, xp, JP_f);
357  indopro_forward_safe(tag_g, m, n, xp, JP_g);
358  nonl_ind_forward_safe(tag_g, m, n, xp, HP_g);
359
360
361  for (i=0;i<n;i++) 
362    {
363      if (HP_f[i][0]+HP_g[i][0]!=0)
364        {
365          if (HP_f[i][0]==0)
366            {
367              HP_t[i] = (unsigned int *) malloc((HP_g[i][0]+HPOFF)*sizeof(unsigned int));
368              for(j=0;j<=(int) HP_g[i][0];j++)
369                {
370                  HP_t[i][j] = HP_g[i][j];
371                }
372              HP_length[i] = HP_g[i][0]+HPOFF;
373            }
374          else
375            {
376              if (HP_g[i][0]==0)
377                {
378                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HPOFF)*sizeof(unsigned int));
379                  for(j=0;j<=(int) HP_f[i][0];j++)
380                    {
381                      HP_t[i][j] = HP_f[i][j];
382                    }
383                  HP_length[i] = HP_f[i][0]+HPOFF;
384                }
385              else
386                {
387                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HP_g[i][0]+HPOFF)*sizeof(unsigned int));
388                  k = l = j = 1;
389                  while ((k<=(int) HP_f[i][0]) && (l <= (int) HP_g[i][0]))
390                    {
391                      if (HP_f[i][k] < HP_g[i][l])
392                        {
393                          HP_t[i][j]=HP_f[i][k];
394                          j++; k++;
395                        }
396                      else
397                        {
398                          if (HP_f[i][k] == HP_g[i][l])
399                            {
400                              HP_t[i][j]=HP_f[i][k];
401                              l++;j++;k++;
402                            }
403                          else
404                            {
405                              HP_t[i][j]=HP_g[i][l];
406                              j++;l++;               
407                            }
408                        }
409                    } // end while
410                  for(ii=k;ii<=(int) HP_f[i][0];ii++)
411                    {
412                      HP_t[i][j] = HP_f[i][ii];
413                      j++;
414                    }
415                  for(ii=l;ii<=(int) HP_g[i][0];ii++)
416                    {
417                      HP_t[i][j] = HP_g[i][ii];
418                      j++;
419                    }
420                 
421                }
422            }
423          HP_t[i][0]=j-1;
424          HP_length[i] = HP_f[i][0]+HP_g[i][0]+HPOFF;
425        }
426      else
427        {
428          HP_t[i] = (unsigned int *) malloc((HPOFF+1)*sizeof(unsigned int));
429          HP_t[i][0]=0;
430          HP_length[i]=HPOFF;
431        }
432    }   
433
434  for (i=0;i<m;i++) 
435    {
436      HP_t[n+i] = (unsigned int *) malloc((JP_g[i][0]+1)*sizeof(unsigned int));
437      HP_t[n+i][0]=JP_g[i][0];
438      for(j=1;j<= (int) JP_g[i][0];j++)
439        {
440          HP_t[n+i][j]=JP_g[i][j];
441          if (HP_length[JP_g[i][j]] < HP_t[JP_g[i][j]][0]+1)
442            {
443              temp = (unsigned int *) malloc((HP_t[JP_g[i][j]][0])*sizeof(unsigned int));
444              for(l=0;l<=(int)HP_t[JP_g[i][j]][0];l++)
445                temp[l] = HP_t[JP_g[i][j]][l];
446              free(HP_t[JP_g[i][j]]);
447              HP_t[JP_g[i][j]] = (unsigned int *) malloc(2*HP_length[JP_g[i][j]]*sizeof(unsigned int));
448              HP_length[JP_g[i][j]] = 2*HP_length[JP_g[i][j]];
449              for(l=0;l<=(int)temp[0];l++)
450                HP_t[JP_g[i][j]][l] =temp[l];
451              free(temp);
452            }
453          HP_t[JP_g[i][j]][0] = HP_t[JP_g[i][j]][0]+1;
454          HP_t[JP_g[i][j]][HP_t[JP_g[i][j]][0]] = i+n;
455        }
456    }
457
458
459  for(j=1;j<= (int) JP_f[0][0];j++)
460    {
461      if (HP_length[JP_f[0][j]] < HP_t[JP_f[0][j]][0]+1)
462        {
463          temp = (unsigned int *) malloc((HP_t[JP_f[0][j]][0])*sizeof(unsigned int));
464          for(l=0;l<=(int)HP_t[JP_f[0][j]][0];l++)
465            temp[l] = HP_t[JP_f[0][j]][l];
466          free(HP_t[JP_f[0][j]]);
467          HP_t[JP_f[0][j]] = (unsigned int *) malloc(2*HP_length[JP_f[0][j]]*sizeof(unsigned int));
468          HP_length[JP_f[0][j]] = 2*HP_length[JP_f[0][j]];
469          for(l=0;l<=(int)temp[0];l++)
470            HP_t[JP_f[0][j]][l] =temp[l];
471          free(temp);
472        }
473      HP_t[JP_f[0][j]][0] = HP_t[JP_f[0][j]][0]+1;
474      HP_t[JP_f[0][j]][HP_t[JP_f[0][j]][0]] = n+m;
475    }
476
477
478  HP_t[n+m] = (unsigned int *) malloc((JP_f[0][0]+2)*sizeof(unsigned int));
479  HP_t[n+m][0]=JP_f[0][0]+1;
480  for(j=1;j<= (int) JP_f[0][0];j++)
481    HP_t[n+m][j]=JP_f[0][j];
482  HP_t[n+m][JP_f[0][0]+1]=n+m;
483
484
485
486  set_HP(tag_L,n+m+1,HP_t);
487
488  nnz_h_lag = 0;
489   for (i=0;i<n;i++) {
490    for (j=1;j<=(int) HP_t[i][0];j++)
491      if ((int) HP_t[i][j] <= i)
492        nnz_h_lag++;
493     free(HP_f[i]);
494     free(HP_g[i]);
495   }
496  nnz_L = nnz_h_lag;
497
498  options_L[0] = 0;         
499  options_L[1] = 1;       
500
501  sparse_hess(tag_L, n+m+1, -1, xp, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);
502
503  rind_L = new unsigned int[nnz_L];
504  cind_L = new unsigned int[nnz_L];
505  rind_L_total = new unsigned int[nnz_L_total];
506  cind_L_total = new unsigned int[nnz_L_total];
507
508  unsigned int ind = 0;
509
510  for (int i=0;i<n;i++) 
511    for (unsigned int j=1;j<=HP_t[i][0];j++)
512      {
513        if (((int) HP_t[i][j]>=i) &&((int) HP_t[i][j]<n)) 
514          {
515            rind_L[ind] = i;
516            cind_L[ind++] = HP_t[i][j];
517          }
518      }
519
520   ind = 0;
521   for (int i=0;i<n+m+1;i++) 
522     for (unsigned int j=1;j<=HP_t[i][0];j++)
523       {
524        if ((int) HP_t[i][j]>=i) 
525          {
526            rind_L_total[ind] = i;
527            cind_L_total[ind++] = HP_t[i][j];
528          }
529       }
530
531  for (i=0;i<m;i++) {
532     free(JP_g[i]);
533   }
534
535  free(JP_f[0]);
536  free(JP_f);
537  free(JP_g);
538  free(HP_f);
539  free(HP_g);
540  free(HP_length);
541
542  delete[] lam;
543  delete[] g;
544  delete[] xa;
545  delete[] zu;
546  delete[] zl;
547  delete[] lamp;
548  delete[] xp;
549}
Note: See TracBrowser for help on using the repository browser.