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