# source:branches/parallel/Ipopt/examples/ScalableProblems/LuksanVlcek2.cpp@1729

Last change on this file since 1729 was 1729, checked in by andreasw, 4 years ago

synchronized branches/parallel with trunk rev 1728

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 6.7 KB
Line
3// This code is published under the Common Public License.
4//
5// \$Id\$
6//
7// Authors:  Andreas Waechter              IBM    2005-10-127
8
9#include "LuksanVlcek2.hpp"
10
11#ifdef HAVE_CMATH
12# include <cmath>
13#else
14# ifdef HAVE_MATH_H
15#  include <math.h>
16# else
17#  error "don't have header file for math"
18# endif
19#endif
20
21#include <cstdio>
22
23using namespace Ipopt;
24
25LuksanVlcek2::LuksanVlcek2(Number g_l, Number g_u)
26{
27  g_l_ = g_l;
28  g_u_ = g_u;
29}
30
31bool LuksanVlcek2::InitializeProblem(Index N)
32{
33  N_=N;
34  if (N_<=13 || 2*(N_/2)!=N_) {
35    printf("N needs to be at least 14 and even.\n");
36    return false;
37  }
38  return true;
39}
40
41// returns the size of the problem
42bool LuksanVlcek2::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
43                                Index& nnz_h_lag, IndexStyleEnum& index_style)
44{
45  // The problem described in LuksanVlcek2.hpp has 4 variables, x[0]
46  // through x[3]
47  n = N_+2;
48
49  m = N_-7;
50
51  nnz_jac_g = 25 + (m-5)*8;
52
53  nnz_h_lag = n + N_ + 1;
54
55  // use the C style numbering of matrix indices (starting at 0)
56  index_style = TNLP::C_STYLE;
57
58  return true;
59}
60
61// returns the variable bounds
62bool LuksanVlcek2::get_bounds_info(Index n, Number* x_l, Number* x_u,
63                                   Index m, Number* g_l, Number* g_u)
64{
65  // none of the variables have bounds
66  for (Index i=0; i<n; i++) {
67    x_l[i] = -1e20;
68    x_u[i] =  1e20;
69  }
70
71  // Set the bounds for the constraints
72  for (Index i=0; i<m; i++) {
73    g_l[i] = g_l_;
74    g_u[i] = g_u_;
75  }
76
77  return true;
78}
79
80// returns the initial point for the problem
81bool LuksanVlcek2::get_starting_point(Index n, bool init_x, Number* x,
82                                      bool init_z, Number* z_L, Number* z_U,
83                                      Index m, bool init_lambda,
84                                      Number* lambda)
85{
86  if (!init_x || init_z || init_lambda) {
87    return false;
88  }
89
90  // set the starting point
91  for (Index i=0; i<n/2; i++) {
92    x[2*i] = -2.;
93    x[2*i+1] = 1.;
94  }
95
96  /*
97  // DELETEME
98  for (Index i=0; i<n; i++) {
99    x[i] += 0.1*((Number) i);
100  }
101  */
102
103  return true;
104}
105
106// returns the value of the objective function
107bool LuksanVlcek2::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
108{
109  obj_value = 0.;
110  for (Index i=0; i<N_/2; i++) {
111    Number a1 = x[2*i]*x[2*i] - x[2*i+1];
112    Number a2 = x[2*i] - 1.;
113    Number a3 = x[2*i+2]*x[2*i+2] - x[2*i+3];
114    Number a4 = x[2*i+2] - 1.;
115    Number a5 = x[2*i+1] + x[2*i+3] - 2.;
116    Number a6 = x[2*i+1] - x[2*i+3];
117    obj_value += 100.*a1*a1 + a2*a2 + 90.*a3*a3 + a4*a4 + 10.*a5*a5 + .1*a6*a6;
118  }
119
120  return true;
121}
122
125{
128  for (Index i=0; i<N_/2; i++) {
129    grad_f[2*i] += 400.*x[2*i]*(x[2*i]*x[2*i]-x[2*i+1]) + 2.*(x[2*i]-1.);
131                     + 20*(x[2*i+1] + x[2*i+3] - 2.)
132                     + 0.2*(x[2*i+1] - x[2*i+3]);
133    grad_f[2*i+2] = 360.*x[2*i+2]*(x[2*i+2]*x[2*i+2] - x[2*i+3])
134                    + 2.*(x[2*i+2] -1.);
135    grad_f[2*i+3] = -180.*(x[2*i+2]*x[2*i+2] - x[2*i+3])
136                    + 20.*(x[2*i+1] + x[2*i+3] -2.)
137                    - 0.2*(x[2*i+1] - x[2*i+3]);
138  }
139
140  return true;
141}
142
143// return the value of the constraints: g(x)
144bool LuksanVlcek2::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
145{
146  for (Index i=0; i<N_-7; i++) {
147    g[i] = (2.+5.*x[i+5]*x[i+5])*x[i+5] + 1.;
148    for (Index k=Max(0,i-5); k<=i+1; k++) {
149      g[i] += x[k]*(x[k]+1.);
150    }
151  }
152
153  return true;
154}
155
156// return the structure or values of the jacobian
157bool LuksanVlcek2::eval_jac_g(Index n, const Number* x, bool new_x,
158                              Index m, Index nele_jac, Index* iRow, Index *jCol,
159                              Number* values)
160{
161  if (values == NULL) {
162    // return the structure of the jacobian
163
164    Index ijac=0;
165    for (Index i=0; i<N_-7; i++) {
166      for (Index k=Max(0,i-5); k<=i+1; k++) {
167        iRow[ijac] = i;
168        jCol[ijac] = k;
169        ijac++;
170      }
171      iRow[ijac] = i;
172      jCol[ijac] = i+5;
173      ijac++;
174    }
175    DBG_ASSERT(ijac == nele_jac);
176  }
177  else {
178    // return the values of the jacobian of the constraints
179
180    Index ijac=0;
181    for (Index i=0; i<N_-7; i++) {
182      for (Index k=Max(0,i-5); k<=i+1; k++) {
183        values[ijac] = 2.*x[k] + 1.;
184        ijac++;
185      }
186      values[ijac] = 2. + 15.*x[i+5]*x[i+5];
187      ijac++;
188    }
189  }
190
191  return true;
192}
193
194//return the structure or values of the hessian
195bool LuksanVlcek2::eval_h(Index n, const Number* x, bool new_x,
196                          Number obj_factor, Index m, const Number* lambda,
197                          bool new_lambda, Index nele_hess, Index* iRow,
198                          Index* jCol, Number* values)
199{
200  if (values == NULL) {
201    Index ihes=0;
202    // First the diagonal elements
203    for (Index i=0; i<n; i++) {
204      iRow[ihes] = i;
205      jCol[ihes] = i;
206      ihes++;
207    }
208    // And now the off-diagonal elements
209    for (Index i=0; i<N_/2; i++) {
210      iRow[ihes] = 2*i;
211      jCol[ihes] = 2*i+1;
212      ihes++;
213      iRow[ihes] = 2*i+1;
214      jCol[ihes] = 2*i+3;
215      ihes++;
216    }
217    iRow[ihes] = n-2;
218    jCol[ihes] = n-1;
219    ihes++;
220    DBG_ASSERT(ihes == nele_hess);
221  }
222  else {
223    // First we take care of the diagonal elements coming from the
224    // objective function
225    values[0] = 0.;
226    values[1] = 0.;
227    for (Index i=0; i<N_/2; i++) {
228      values[2*i] += obj_factor*(1200.*x[2*i]*x[2*i] - 400.*x[2*i+1] + 2.);
229      values[2*i+1] += obj_factor*220.2;
230      values[2*i+2] = obj_factor*(1080.*x[2*i+2]*x[2*i+2] - 360*x[2*i+3] + 2.);
231      values[2*i+3] = obj_factor*200.2;
232    }
233    // Now we take care of the off-diagonal elements coming from the
234    // objective function
235    Index ihes = n;
236    values[ihes] = 0.;
237    for (Index i=0; i<N_/2; i++) {
238      values[ihes] += obj_factor*(-400.*x[2*i]);
239      ihes++;
240      values[ihes] = obj_factor*19.8;
241      ihes++;
242      values[ihes] = obj_factor*(-360.*x[2*i+2]);
243    }
244
245    // Ok, now the diagonal elements from the constraints
246    for (Index i=0; i<N_-7; i++) {
247      for (Index k=Max(0,i-5); k<=i+1; k++) {
248        values[k] += lambda[i]*2.;
249      }
250      values[i+5] += lambda[i]*30.*x[i+5];
251    }
252  }
253
254  return true;
255}
256
257void LuksanVlcek2::finalize_solution(SolverReturn status,
258                                     Index n, const Number* x, const Number* z_L, const Number* z_U,
259                                     Index m, const Number* g, const Number* lambda,
260                                     Number obj_value,
261                                     const IpoptData* ip_data,
262                                     IpoptCalculatedQuantities* ip_cq)
263{}
264
Note: See TracBrowser for help on using the repository browser.