# source:branches/devel/Ipopt/examples/ScalableProblems/LuksanVlcek3.cpp@938

Last change on this file since 938 was 938, checked in by andreasw, 7 years ago

provide IpoptData? and IpoptCalculatedQuantities? in finalize_solution methods

• Property svn:eol-style set to `native`
• Property svn:keywords set to `"Author Date Id Revision"`
File size: 8.0 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 "LuksanVlcek3.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
21using namespace Ipopt;
22
23LuksanVlcek3::LuksanVlcek3(Number g_l, Number g_u)
24{
25  g_l_ = g_l;
26  g_u_ = g_u;
27}
28
29bool LuksanVlcek3::InitializeProblem(Index N)
30{
31  N_=N;
32  if (N_<=5 || 4*((N_+2)/4)!=N_+2) {
33    printf("N needs to be at least 6 and N+2 needs to be a multiple of 4.\n");
34    return false;
35  }
36  return true;
37}
38
39// returns the size of the problem
40bool LuksanVlcek3::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
41                                Index& nnz_h_lag, IndexStyleEnum& index_style)
42{
43  // The problem described in LuksanVlcek3.hpp has 4 variables, x[0] through x[3]
44  n = N_+2;
45
46  m = 2;
47
48  nnz_jac_g = 4;
49
50  nnz_h_lag = 5*N_/2 + 3;
51
52  // use the C style numbering of matrix indices (starting at 0)
53  index_style = TNLP::C_STYLE;
54
55  return true;
56}
57
58// returns the variable bounds
59bool LuksanVlcek3::get_bounds_info(Index n, Number* x_l, Number* x_u,
60                                   Index m, Number* g_l, Number* g_u)
61{
62  // none of the variables have bounds
63  for (Index i=0; i<n; i++) {
64    x_l[i] = -1e20;
65    x_u[i] =  1e20;
66  }
67
68  // Set the bounds for the constraints
69  for (Index i=0; i<m; i++) {
70    g_l[i] = g_l_;
71    g_u[i] = g_u_;
72  }
73
74  return true;
75}
76
77// returns the initial point for the problem
78bool LuksanVlcek3::get_starting_point(Index n, bool init_x, Number* x,
79                                      bool init_z, Number* z_L, Number* z_U,
80                                      Index m, bool init_lambda,
81                                      Number* lambda)
82{
83  if (!init_x || init_z || init_lambda) {
84    return false;
85  }
86
87  // set the starting point
88  for (Index i=0; i<n/4; i++) {
89    x[4*i] = 3.;
90    x[4*i+1] = -1.;
91    x[4*i+2] = 0.;
92    x[4*i+3] = 1.;
93  }
94
95  /*
96  // DELETEME
97  for (Index i=0; i<n; i++) {
98    x[i] += 0.1*((Number) i);
99  }
100  */
101
102  return true;
103}
104
105// returns the value of the objective function
106bool LuksanVlcek3::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
107{
108  obj_value = 0.;
109  for (Index i=0; i<N_/2; i++) {
110    Number a1 = x[2*i]+10.*x[2*i+1];
111    Number a2 = x[2*i+2] - x[2*i+3];
112    Number a3 = x[2*i+1] - 2.*x[2*i+2];
113    Number a4 = x[2*i] - x[2*i+3];
114    obj_value += a1*a1 + 5.*a2*a2 + pow(a3,4)+ 10.*pow(a4,4);
115  }
116
117  return true;
118}
119
122{
125  for (Index i=0; i<N_/2; i++) {
126    Number a1 = x[2*i]+10.*x[2*i+1];
127    Number a2 = x[2*i+2] - x[2*i+3];
128    Number a3 = x[2*i+1] - 2.*x[2*i+2];
129    Number a4 = x[2*i] - x[2*i+3];
130    grad_f[2*i] += 2.*a1 + 40.*pow(a4,3);
131    grad_f[2*i+1] += 20.*a1 + 4.*pow(a3,3);
132    grad_f[2*i+2] = 10.*a2 - 8.*pow(a3,3);
133    grad_f[2*i+3] = -10.*a2 - 40.*pow(a4,3);
134  }
135
136  return true;
137}
138
139// return the value of the constraints: g(x)
140bool LuksanVlcek3::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
141{
142  g[0] = 3.*pow(x[0],3) + 2.*x[1] - 5. + sin(x[0]-x[1])*sin(x[0]+x[1]);
143  g[1] = 4.*x[n-3] - x[n-4]*exp(x[n-4]-x[n-3]) - 3;
144  ;
145  return true;
146}
147
148// return the structure or values of the jacobian
149bool LuksanVlcek3::eval_jac_g(Index n, const Number* x, bool new_x,
150                              Index m, Index nele_jac, Index* iRow, Index *jCol,
151                              Number* values)
152{
153  if (values == NULL) {
154    // return the structure of the jacobian
155
156    Index ijac = 0;
157    iRow[ijac] = 0;
158    jCol[ijac] = 0;
159    ijac++;
160    iRow[ijac] = 0;
161    jCol[ijac] = 1;
162    ijac++;
163    iRow[ijac] = 1;
164    jCol[ijac] = n-4;
165    ijac++;
166    iRow[ijac] = 1;
167    jCol[ijac] = n-3;
168    ijac++;
169
170    DBG_ASSERT(ijac == nele_jac);
171  }
172  else {
173    // return the values of the jacobian of the constraints
174
175    Index ijac = 0;
176    values[ijac] = 9.*x[0]*x[0]
177                   + cos(x[0]-x[1])*sin(x[0]+x[1])
178                   + sin(x[0]-x[1])*cos(x[0]+x[1]);
179    ijac++;
180    values[ijac] = 2.
181                   - cos(x[0]-x[1])*sin(x[0]+x[1])
182                   + sin(x[0]-x[1])*cos(x[0]+x[1]);
183    ijac++;
184    values[ijac] = -(1.+x[n-4])*exp(x[n-4]-x[n-3]);
185    ijac++;
186    values[ijac] = 4. + x[n-4]*exp(x[n-4]-x[n-3]);
187    ijac++;
188  }
189
190  return true;
191}
192
193//return the structure or values of the hessian
194bool LuksanVlcek3::eval_h(Index n, const Number* x, bool new_x,
195                          Number obj_factor, Index m, const Number* lambda,
196                          bool new_lambda, Index nele_hess, Index* iRow,
197                          Index* jCol, Number* values)
198{
199  if (values == NULL) {
200    Index ihes=0;
201    for (Index i=0; i<N_/2; i++) {
202      iRow[ihes] = 2*i;
203      jCol[ihes] = 2*i;
204      ihes++;
205      iRow[ihes] = 2*i;
206      jCol[ihes] = 2*i+1;
207      ihes++;
208      iRow[ihes] = 2*i;
209      jCol[ihes] = 2*i+3;
210      ihes++;
211      iRow[ihes] = 2*i+1;
212      jCol[ihes] = 2*i+1;
213      ihes++;
214      iRow[ihes] = 2*i+1;
215      jCol[ihes] = 2*i+2;
216      ihes++;
217    }
218    iRow[ihes] = N_;
219    jCol[ihes] = N_;
220    ihes++;
221    iRow[ihes] = N_;
222    jCol[ihes] = N_+1;
223    ihes++;
224    iRow[ihes] = N_+1;
225    jCol[ihes] = N_+1;
226    ihes++;
227    DBG_ASSERT(ihes == nele_hess);
228  }
229  else {
230    Index ihes=0;
231    values[0] = 0.;
232    values[1] = 0.;
233    values[3] = 0.;
234    for (Index i=0; i<N_/2-1; i++) {
235      Number a3 = x[2*i+1] - 2.*x[2*i+2];
236      Number a4 = x[2*i] - x[2*i+3];
237      // x[2*i] x[2*i]
238      values[ihes] += obj_factor*(2. + 120.*a4*a4);
239      ihes++;
240      // x[2*i] x[2*i+1]
241      values[ihes] += obj_factor*20.;
242      ihes++;
243      // x[2*i] x[2*i+3]
244      values[ihes] = -obj_factor*120.*a4*a4;
245      ihes++;
246      // x[2*i+1] x[2*i+1]
247      values[ihes] += obj_factor*(200. + 12.*a3*a3);
248      ihes++;
249      // x[2*i+1] x[2*i+2]
250      values[ihes] = -obj_factor*24.*a3*a3;
251      ihes++;
252      // x[2*i+2] x[2*i+2]
253      values[ihes] = obj_factor*(10. + 48.*a3*a3);
254      // x[2*i+2] x[2*i+3]
255      values[ihes+1] = -obj_factor*10.;
256      // x[2*i+3] x[2*i+3]
257      values[ihes+3] = obj_factor*(10. + 120.*a4*a4);
258    }
259    {
260      Index i = N_/2-1;
261      Number a3 = x[2*i+1] - 2.*x[2*i+2];
262      Number a4 = x[2*i] - x[2*i+3];
263      // x[2*i] x[2*i]
264      values[ihes] += obj_factor*(2. + 120.*a4*a4);
265      ihes++;
266      // x[2*i] x[2*i+1]
267      values[ihes] += obj_factor*20.;
268      ihes++;
269      // x[2*i] x[2*i+3]
270      values[ihes] = -obj_factor*120.*a4*a4;
271      ihes++;
272      // x[2*i+1] x[2*i+1]
273      values[ihes] += obj_factor*(200. + 12.*a3*a3);
274      ihes++;
275      // x[2*i+1] x[2*i+2]
276      values[ihes] = -obj_factor*24.*a3*a3;
277      ihes++;
278      // x[2*i+2] x[2*i+2]
279      values[ihes] = obj_factor*(10. + 48.*a3*a3);
280      // x[2*i+2] x[2*i+3]
281      values[ihes+1] = -obj_factor*10.;
282      // x[2*i+3] x[2*i+3]
283      values[ihes+2] = obj_factor*(10. + 120.*a4*a4);
284    }
285
286    // Now the constraints
287    ihes = 0;
288    Number d1 = x[0] - x[1];
289    Number d2 = x[0] + x[1];
290    values[ihes] += lambda[0]*(18.*x[0]
291                               -2.*sin(d1)*sin(d2)
292                               +2.*cos(d1)*cos(d2));
293    ihes+=3;
294    values[ihes] += lambda[0]*(-2.*sin(d1)*sin(d2)
295                               -2.*cos(d1)*cos(d2));
296
297    d1 = x[n-4]-x[n-3];
298
299    // x[n-4] x[n-4]
300    ihes = nele_hess - 8;
301    values[ihes] += -lambda[1]*(2.+x[n-4])*exp(d1);
302
303    // x[n-4] x[n-3]
304    ihes++;
305    values[ihes] += lambda[1]*(1.+x[n-4])*exp(d1);
306
307    // x[n-3] x[n-3]
308    ihes += 2;
309    values[ihes] += -lambda[1]*x[n-4]*exp(d1);
310  }
311  return true;
312}
313
314void LuksanVlcek3::finalize_solution(SolverReturn status,
315                                     Index n, const Number* x, const Number* z_L, const Number* z_U,
316                                     Index m, const Number* g, const Number* lambda,
317                                     Number obj_value,
318                                     const IpoptData* ip_data,
319                                     IpoptCalculatedQuantities* ip_cq)
320{}
321
Note: See TracBrowser for help on using the repository browser.