1 | /////////////////////////////////////////////////////////////////////////////////
|
---|
2 | //
|
---|
3 | // Demonstration driver program for the Levenberg - Marquardt minimization
|
---|
4 | // algorithm
|
---|
5 | // Copyright (C) 2004-05 Manolis Lourakis (lourakis at ics forth gr)
|
---|
6 | // Institute of Computer Science, Foundation for Research & Technology - Hellas
|
---|
7 | // Heraklion, Crete, Greece.
|
---|
8 | //
|
---|
9 | // This program is free software; you can redistribute it and/or modify
|
---|
10 | // it under the terms of the GNU General Public License as published by
|
---|
11 | // the Free Software Foundation; either version 2 of the License, or
|
---|
12 | // (at your option) any later version.
|
---|
13 | //
|
---|
14 | // This program is distributed in the hope that it will be useful,
|
---|
15 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
16 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
17 | // GNU General Public License for more details.
|
---|
18 | //
|
---|
19 | /////////////////////////////////////////////////////////////////////////////////
|
---|
20 |
|
---|
21 | /********************************************************************************
|
---|
22 | * Levenberg-Marquardt minimization demo driver. Only the double precision versions
|
---|
23 | * are tested here. See the Meyer case for an example of verifying the Jacobian
|
---|
24 | ********************************************************************************/
|
---|
25 |
|
---|
26 | #include <stdio.h>
|
---|
27 | #include <stdlib.h>
|
---|
28 | #include <math.h>
|
---|
29 | #include <float.h>
|
---|
30 |
|
---|
31 | #include "levmar.h"
|
---|
32 |
|
---|
33 | #ifndef LM_DBL_PREC
|
---|
34 | #error Demo program assumes that levmar has been compiled with double precision, see LM_DBL_PREC!
|
---|
35 | #endif
|
---|
36 |
|
---|
37 |
|
---|
38 | /* Sample functions to be minimized with LM and their Jacobians.
|
---|
39 | * More test functions at http://www.csit.fsu.edu/~burkardt/f_src/test_nls/test_nls.html
|
---|
40 | * Check also the CUTE problems collection at ftp://ftp.numerical.rl.ac.uk/pub/cute/;
|
---|
41 | * CUTE is searchable through http://numawww.mathematik.tu-darmstadt.de:8081/opti/select.html
|
---|
42 | * CUTE problems can also be solved through the AMPL web interface at http://www.ampl.com/TRYAMPL/startup.html
|
---|
43 | *
|
---|
44 | * Nonlinear optimization models in AMPL can be found at http://www.princeton.edu/~rvdb/ampl/nlmodels/
|
---|
45 | */
|
---|
46 |
|
---|
47 | #define ROSD 105.0
|
---|
48 |
|
---|
49 | /* Rosenbrock function, global minimum at (1, 1) */
|
---|
50 | void ros(double *p, double *x, int m, int n, void *data)
|
---|
51 | {
|
---|
52 | register int i;
|
---|
53 |
|
---|
54 | for(i=0; i<n; ++i)
|
---|
55 | x[i]=((1.0-p[0])*(1.0-p[0]) + ROSD*(p[1]-p[0]*p[0])*(p[1]-p[0]*p[0]));
|
---|
56 | }
|
---|
57 |
|
---|
58 | void jacros(double *p, double *jac, int m, int n, void *data)
|
---|
59 | {
|
---|
60 | register int i, j;
|
---|
61 |
|
---|
62 | for(i=j=0; i<n; ++i){
|
---|
63 | jac[j++]=(-2 + 2*p[0]-4*ROSD*(p[1]-p[0]*p[0])*p[0]);
|
---|
64 | jac[j++]=(2*ROSD*(p[1]-p[0]*p[0]));
|
---|
65 | }
|
---|
66 | }
|
---|
67 |
|
---|
68 |
|
---|
69 | #define MODROSLAM 1E02
|
---|
70 | /* Modified Rosenbrock problem, global minimum at (1, 1) */
|
---|
71 | void modros(double *p, double *x, int m, int n, void *data)
|
---|
72 | {
|
---|
73 | register int i;
|
---|
74 |
|
---|
75 | for(i=0; i<n; i+=3){
|
---|
76 | x[i]=10*(p[1]-p[0]*p[0]);
|
---|
77 | x[i+1]=1.0-p[0];
|
---|
78 | x[i+2]=MODROSLAM;
|
---|
79 | }
|
---|
80 | }
|
---|
81 |
|
---|
82 | void jacmodros(double *p, double *jac, int m, int n, void *data)
|
---|
83 | {
|
---|
84 | register int i, j;
|
---|
85 |
|
---|
86 | for(i=j=0; i<n; i+=3){
|
---|
87 | jac[j++]=-20.0*p[0];
|
---|
88 | jac[j++]=10.0;
|
---|
89 |
|
---|
90 | jac[j++]=-1.0;
|
---|
91 | jac[j++]=0.0;
|
---|
92 |
|
---|
93 | jac[j++]=0.0;
|
---|
94 | jac[j++]=0.0;
|
---|
95 | }
|
---|
96 | }
|
---|
97 |
|
---|
98 |
|
---|
99 | /* Powell's function, minimum at (0, 0) */
|
---|
100 | void powell(double *p, double *x, int m, int n, void *data)
|
---|
101 | {
|
---|
102 | register int i;
|
---|
103 |
|
---|
104 | for(i=0; i<n; i+=2){
|
---|
105 | x[i]=p[0];
|
---|
106 | x[i+1]=10.0*p[0]/(p[0]+0.1) + 2*p[1]*p[1];
|
---|
107 | }
|
---|
108 | }
|
---|
109 |
|
---|
110 | void jacpowell(double *p, double *jac, int m, int n, void *data)
|
---|
111 | {
|
---|
112 | register int i, j;
|
---|
113 |
|
---|
114 | for(i=j=0; i<n; i+=2){
|
---|
115 | jac[j++]=1.0;
|
---|
116 | jac[j++]=0.0;
|
---|
117 |
|
---|
118 | jac[j++]=1.0/((p[0]+0.1)*(p[0]+0.1));
|
---|
119 | jac[j++]=4.0*p[1];
|
---|
120 | }
|
---|
121 | }
|
---|
122 |
|
---|
123 | /* Wood's function, minimum at (1, 1, 1, 1) */
|
---|
124 | void wood(double *p, double *x, int m, int n, void *data)
|
---|
125 | {
|
---|
126 | register int i;
|
---|
127 |
|
---|
128 | for(i=0; i<n; i+=6){
|
---|
129 | x[i]=10.0*(p[1] - p[0]*p[0]);
|
---|
130 | x[i+1]=1.0 - p[0];
|
---|
131 | x[i+2]=sqrt(90.0)*(p[3] - p[2]*p[2]);
|
---|
132 | x[i+3]=1.0 - p[2];
|
---|
133 | x[i+4]=sqrt(10.0)*(p[1]+p[3] - 2.0);
|
---|
134 | x[i+5]=(p[1] - p[3])/sqrt(10.0);
|
---|
135 | }
|
---|
136 | }
|
---|
137 |
|
---|
138 | /* Meyer's (reformulated) problem, minimum at (2.48, 6.18, 3.45) */
|
---|
139 | void meyer(double *p, double *x, int m, int n, void *data)
|
---|
140 | {
|
---|
141 | register int i;
|
---|
142 | double ui;
|
---|
143 |
|
---|
144 | for(i=0; i<n; ++i){
|
---|
145 | ui=0.45+0.05*i;
|
---|
146 | x[i]=p[0]*exp(10.0*p[1]/(ui+p[2]) - 13.0);
|
---|
147 | }
|
---|
148 | }
|
---|
149 |
|
---|
150 | void jacmeyer(double *p, double *jac, int m, int n, void *data)
|
---|
151 | {
|
---|
152 | register int i, j;
|
---|
153 | double ui, tmp;
|
---|
154 |
|
---|
155 | for(i=j=0; i<n; ++i){
|
---|
156 | ui=0.45+0.05*i;
|
---|
157 | tmp=exp(10.0*p[1]/(ui+p[2]) - 13.0);
|
---|
158 |
|
---|
159 | jac[j++]=tmp;
|
---|
160 | jac[j++]=10.0*p[0]*tmp/(ui+p[2]);
|
---|
161 | jac[j++]=-10.0*p[0]*p[1]*tmp/((ui+p[2])*(ui+p[2]));
|
---|
162 | }
|
---|
163 | }
|
---|
164 |
|
---|
165 | /* Osborne's problem, minimum at (0.3754, 1.9358, -1.4647, 0.0129, 0.0221) */
|
---|
166 | void osborne(double *p, double *x, int m, int n, void *data)
|
---|
167 | {
|
---|
168 | register int i;
|
---|
169 | double t;
|
---|
170 |
|
---|
171 | for(i=0; i<n; ++i){
|
---|
172 | t=10*i;
|
---|
173 | x[i]=p[0] + p[1]*exp(-p[3]*t) + p[2]*exp(-p[4]*t);
|
---|
174 | }
|
---|
175 | }
|
---|
176 |
|
---|
177 | void jacosborne(double *p, double *jac, int m, int n, void *data)
|
---|
178 | {
|
---|
179 | register int i, j;
|
---|
180 | double t, tmp1, tmp2;
|
---|
181 |
|
---|
182 | for(i=j=0; i<n; ++i){
|
---|
183 | t=10*i;
|
---|
184 | tmp1=exp(-p[3]*t);
|
---|
185 | tmp2=exp(-p[4]*t);
|
---|
186 |
|
---|
187 | jac[j++]=1.0;
|
---|
188 | jac[j++]=tmp1;
|
---|
189 | jac[j++]=tmp2;
|
---|
190 | jac[j++]=-p[1]*t*tmp1;
|
---|
191 | jac[j++]=-p[2]*t*tmp2;
|
---|
192 | }
|
---|
193 | }
|
---|
194 |
|
---|
195 | /* helical valley function, minimum at (1.0, 0.0, 0.0) */
|
---|
196 | #ifndef M_PI
|
---|
197 | #define M_PI 3.14159265358979323846 /* pi */
|
---|
198 | #endif
|
---|
199 |
|
---|
200 | void helval(double *p, double *x, int m, int n, void *data)
|
---|
201 | {
|
---|
202 | double theta;
|
---|
203 |
|
---|
204 | if(p[0]<0.0)
|
---|
205 | theta=atan(p[1]/p[0])/(2.0*M_PI) + 0.5;
|
---|
206 | else if(0.0<p[0])
|
---|
207 | theta=atan(p[1]/p[0])/(2.0*M_PI);
|
---|
208 | else
|
---|
209 | theta=(p[1]>=0)? 0.25 : -0.25;
|
---|
210 |
|
---|
211 | x[0]=10.0*(p[2] - 10.0*theta);
|
---|
212 | x[1]=10.0*(sqrt(p[0]*p[0] + p[1]*p[1]) - 1.0);
|
---|
213 | x[2]=p[2];
|
---|
214 | }
|
---|
215 |
|
---|
216 | void jachelval(double *p, double *jac, int m, int n, void *data)
|
---|
217 | {
|
---|
218 | register int i=0;
|
---|
219 | double tmp;
|
---|
220 |
|
---|
221 | tmp=p[0]*p[0] + p[1]*p[1];
|
---|
222 |
|
---|
223 | jac[i++]=50.0*p[1]/(M_PI*tmp);
|
---|
224 | jac[i++]=-50.0*p[0]/(M_PI*tmp);
|
---|
225 | jac[i++]=10.0;
|
---|
226 |
|
---|
227 | jac[i++]=10.0*p[0]/sqrt(tmp);
|
---|
228 | jac[i++]=10.0*p[1]/sqrt(tmp);
|
---|
229 | jac[i++]=0.0;
|
---|
230 |
|
---|
231 | jac[i++]=0.0;
|
---|
232 | jac[i++]=0.0;
|
---|
233 | jac[i++]=1.0;
|
---|
234 | }
|
---|
235 |
|
---|
236 | /* Boggs - Tolle problem 3 (linearly constrained), minimum at (-0.76744, 0.25581, 0.62791, -0.11628, 0.25581)
|
---|
237 | * constr1: p[0] + 3*p[1] = 0;
|
---|
238 | * constr2: p[2] + p[3] - 2*p[4] = 0;
|
---|
239 | * constr3: p[1] - p[4] = 0;
|
---|
240 | */
|
---|
241 | void bt3(double *p, double *x, int m, int n, void *data)
|
---|
242 | {
|
---|
243 | register int i;
|
---|
244 | double t1, t2, t3, t4;
|
---|
245 |
|
---|
246 | t1=p[0]-p[1];
|
---|
247 | t2=p[1]+p[2]-2.0;
|
---|
248 | t3=p[3]-1.0;
|
---|
249 | t4=p[4]-1.0;
|
---|
250 |
|
---|
251 | for(i=0; i<n; ++i)
|
---|
252 | x[i]=t1*t1 + t2*t2 + t3*t3 + t4*t4;
|
---|
253 | }
|
---|
254 |
|
---|
255 | void jacbt3(double *p, double *jac, int m, int n, void *data)
|
---|
256 | {
|
---|
257 | register int i, j;
|
---|
258 | double t1, t2, t3, t4;
|
---|
259 |
|
---|
260 | t1=p[0]-p[1];
|
---|
261 | t2=p[1]+p[2]-2.0;
|
---|
262 | t3=p[3]-1.0;
|
---|
263 | t4=p[4]-1.0;
|
---|
264 |
|
---|
265 | for(i=j=0; i<n; ++i){
|
---|
266 | jac[j++]=2.0*t1;
|
---|
267 | jac[j++]=2.0*(t2-t1);
|
---|
268 | jac[j++]=2.0*t2;
|
---|
269 | jac[j++]=2.0*t3;
|
---|
270 | jac[j++]=2.0*t4;
|
---|
271 | }
|
---|
272 | }
|
---|
273 |
|
---|
274 | /* Hock - Schittkowski problem 28 (linearly constrained), minimum at (0.5, -0.5, 0.5)
|
---|
275 | * constr1: p[0] + 2*p[1] + 3*p[2] = 1;
|
---|
276 | */
|
---|
277 | void hs28(double *p, double *x, int m, int n, void *data)
|
---|
278 | {
|
---|
279 | register int i;
|
---|
280 | double t1, t2;
|
---|
281 |
|
---|
282 | t1=p[0]+p[1];
|
---|
283 | t2=p[1]+p[2];
|
---|
284 |
|
---|
285 | for(i=0; i<n; ++i)
|
---|
286 | x[i]=t1*t1 + t2*t2;
|
---|
287 | }
|
---|
288 |
|
---|
289 | void jachs28(double *p, double *jac, int m, int n, void *data)
|
---|
290 | {
|
---|
291 | register int i, j;
|
---|
292 | double t1, t2;
|
---|
293 |
|
---|
294 | t1=p[0]+p[1];
|
---|
295 | t2=p[1]+p[2];
|
---|
296 |
|
---|
297 | for(i=j=0; i<n; ++i){
|
---|
298 | jac[j++]=2.0*t1;
|
---|
299 | jac[j++]=2.0*(t1+t2);
|
---|
300 | jac[j++]=2.0*t2;
|
---|
301 | }
|
---|
302 | }
|
---|
303 |
|
---|
304 | /* Hock - Schittkowski problem 48 (linearly constrained), minimum at (1.0, 1.0, 1.0, 1.0, 1.0)
|
---|
305 | * constr1: sum {i in 0..4} p[i] = 5;
|
---|
306 | * constr2: p[2] - 2*(p[3]+p[4]) = -3;
|
---|
307 | */
|
---|
308 | void hs48(double *p, double *x, int m, int n, void *data)
|
---|
309 | {
|
---|
310 | register int i;
|
---|
311 | double t1, t2, t3;
|
---|
312 |
|
---|
313 | t1=p[0]-1.0;
|
---|
314 | t2=p[1]-p[2];
|
---|
315 | t3=p[3]-p[4];
|
---|
316 |
|
---|
317 | for(i=0; i<n; ++i)
|
---|
318 | x[i]=t1*t1 + t2*t2 + t3*t3;
|
---|
319 | }
|
---|
320 |
|
---|
321 | void jachs48(double *p, double *jac, int m, int n, void *data)
|
---|
322 | {
|
---|
323 | register int i, j;
|
---|
324 | double t1, t2, t3;
|
---|
325 |
|
---|
326 | t1=p[0]-1.0;
|
---|
327 | t2=p[1]-p[2];
|
---|
328 | t3=p[3]-p[4];
|
---|
329 |
|
---|
330 | for(i=j=0; i<n; ++i){
|
---|
331 | jac[j++]=2.0*t1;
|
---|
332 | jac[j++]=2.0*t2;
|
---|
333 | jac[j++]=-2.0*t2;
|
---|
334 | jac[j++]=2.0*t3;
|
---|
335 | jac[j++]=-2.0*t3;
|
---|
336 | }
|
---|
337 | }
|
---|
338 |
|
---|
339 | /* Hock - Schittkowski problem 51 (linearly constrained), minimum at (1.0, 1.0, 1.0, 1.0, 1.0)
|
---|
340 | * constr1: p[0] + 3*p[1] = 4;
|
---|
341 | * constr2: p[2] + p[3] - 2*p[4] = 0;
|
---|
342 | * constr3: p[1] - p[4] = 0;
|
---|
343 | */
|
---|
344 | void hs51(double *p, double *x, int m, int n, void *data)
|
---|
345 | {
|
---|
346 | register int i;
|
---|
347 | double t1, t2, t3, t4;
|
---|
348 |
|
---|
349 | t1=p[0]-p[1];
|
---|
350 | t2=p[1]+p[2]-2.0;
|
---|
351 | t3=p[3]-1.0;
|
---|
352 | t4=p[4]-1.0;
|
---|
353 |
|
---|
354 | for(i=0; i<n; ++i)
|
---|
355 | x[i]=t1*t1 + t2*t2 + t3*t3 + t4*t4;
|
---|
356 | }
|
---|
357 |
|
---|
358 | void jachs51(double *p, double *jac, int m, int n, void *data)
|
---|
359 | {
|
---|
360 | register int i, j;
|
---|
361 | double t1, t2, t3, t4;
|
---|
362 |
|
---|
363 | t1=p[0]-p[1];
|
---|
364 | t2=p[1]+p[2]-2.0;
|
---|
365 | t3=p[3]-1.0;
|
---|
366 | t4=p[4]-1.0;
|
---|
367 |
|
---|
368 | for(i=j=0; i<n; ++i){
|
---|
369 | jac[j++]=2.0*t1;
|
---|
370 | jac[j++]=2.0*(t2-t1);
|
---|
371 | jac[j++]=2.0*t2;
|
---|
372 | jac[j++]=2.0*t3;
|
---|
373 | jac[j++]=2.0*t4;
|
---|
374 | }
|
---|
375 | }
|
---|
376 |
|
---|
377 | /* Hock - Schittkowski problem 01 (box constrained), minimum at (1.0, 1.0)
|
---|
378 | * constr1: p[1]>=-1.5;
|
---|
379 | */
|
---|
380 | void hs01(double *p, double *x, int m, int n, void *data)
|
---|
381 | {
|
---|
382 | double t;
|
---|
383 |
|
---|
384 | t=p[0]*p[0];
|
---|
385 | x[0]=10.0*(p[1]-t);
|
---|
386 | x[1]=1.0-p[0];
|
---|
387 | }
|
---|
388 |
|
---|
389 | void jachs01(double *p, double *jac, int m, int n, void *data)
|
---|
390 | {
|
---|
391 | register int j=0;
|
---|
392 |
|
---|
393 | jac[j++]=-20.0*p[0];
|
---|
394 | jac[j++]=10.0;
|
---|
395 |
|
---|
396 | jac[j++]=-1.0;
|
---|
397 | jac[j++]=0.0;
|
---|
398 | }
|
---|
399 |
|
---|
400 | /* Hock - Schittkowski MODIFIED problem 21 (box constrained), minimum at (2.0, 0.0)
|
---|
401 | * constr1: 2 <= p[0] <=50;
|
---|
402 | * constr2: -50 <= p[1] <=50;
|
---|
403 | *
|
---|
404 | * Original HS21 has the additional constraint 10*p[0] - p[1] >= 10; which is inactive
|
---|
405 | * at the solution, so it is dropped here.
|
---|
406 | */
|
---|
407 | void hs21(double *p, double *x, int m, int n, void *data)
|
---|
408 | {
|
---|
409 | x[0]=p[0]/10.0;
|
---|
410 | x[1]=p[1];
|
---|
411 | }
|
---|
412 |
|
---|
413 | void jachs21(double *p, double *jac, int m, int n, void *data)
|
---|
414 | {
|
---|
415 | register int j=0;
|
---|
416 |
|
---|
417 | jac[j++]=0.1;
|
---|
418 | jac[j++]=0.0;
|
---|
419 |
|
---|
420 | jac[j++]=0.0;
|
---|
421 | jac[j++]=1.0;
|
---|
422 | }
|
---|
423 |
|
---|
424 | /* Problem hatfldb (box constrained), minimum at (0.947214, 0.8, 0.64, 0.4096)
|
---|
425 | * constri: p[i]>=0.0; (i=1..4)
|
---|
426 | * constr5: p[1]<=0.8;
|
---|
427 | */
|
---|
428 | void hatfldb(double *p, double *x, int m, int n, void *data)
|
---|
429 | {
|
---|
430 | register int i;
|
---|
431 |
|
---|
432 | x[0]=p[0]-1.0;
|
---|
433 |
|
---|
434 | for(i=1; i<m; ++i)
|
---|
435 | x[i]=p[i-1]-sqrt(p[i]);
|
---|
436 | }
|
---|
437 |
|
---|
438 | void jachatfldb(double *p, double *jac, int m, int n, void *data)
|
---|
439 | {
|
---|
440 | register int j=0;
|
---|
441 |
|
---|
442 | jac[j++]=1.0;
|
---|
443 | jac[j++]=0.0;
|
---|
444 | jac[j++]=0.0;
|
---|
445 | jac[j++]=0.0;
|
---|
446 |
|
---|
447 | jac[j++]=1.0;
|
---|
448 | jac[j++]=-0.5/sqrt(p[1]);
|
---|
449 | jac[j++]=0.0;
|
---|
450 | jac[j++]=0.0;
|
---|
451 |
|
---|
452 | jac[j++]=0.0;
|
---|
453 | jac[j++]=1.0;
|
---|
454 | jac[j++]=-0.5/sqrt(p[2]);
|
---|
455 | jac[j++]=0.0;
|
---|
456 |
|
---|
457 | jac[j++]=0.0;
|
---|
458 | jac[j++]=0.0;
|
---|
459 | jac[j++]=1.0;
|
---|
460 | jac[j++]=-0.5/sqrt(p[3]);
|
---|
461 | }
|
---|
462 |
|
---|
463 | /* Problem hatfldc (box constrained), minimum at (1.0, 1.0, 1.0, 1.0)
|
---|
464 | * constri: p[i]>=0.0; (i=1..4)
|
---|
465 | * constri+4: p[i]<=10.0; (i=1..4)
|
---|
466 | */
|
---|
467 | void hatfldc(double *p, double *x, int m, int n, void *data)
|
---|
468 | {
|
---|
469 | register int i;
|
---|
470 |
|
---|
471 | x[0]=p[0]-1.0;
|
---|
472 |
|
---|
473 | for(i=1; i<m-1; ++i)
|
---|
474 | x[i]=p[i-1]-sqrt(p[i]);
|
---|
475 |
|
---|
476 | x[m-1]=p[m-1]-1.0;
|
---|
477 | }
|
---|
478 |
|
---|
479 | void jachatfldc(double *p, double *jac, int m, int n, void *data)
|
---|
480 | {
|
---|
481 | register int j=0;
|
---|
482 |
|
---|
483 | jac[j++]=1.0;
|
---|
484 | jac[j++]=0.0;
|
---|
485 | jac[j++]=0.0;
|
---|
486 | jac[j++]=0.0;
|
---|
487 |
|
---|
488 | jac[j++]=1.0;
|
---|
489 | jac[j++]=-0.5/sqrt(p[1]);
|
---|
490 | jac[j++]=0.0;
|
---|
491 | jac[j++]=0.0;
|
---|
492 |
|
---|
493 | jac[j++]=0.0;
|
---|
494 | jac[j++]=1.0;
|
---|
495 | jac[j++]=-0.5/sqrt(p[2]);
|
---|
496 | jac[j++]=0.0;
|
---|
497 |
|
---|
498 | jac[j++]=0.0;
|
---|
499 | jac[j++]=0.0;
|
---|
500 | jac[j++]=0.0;
|
---|
501 | jac[j++]=1.0;
|
---|
502 | }
|
---|
503 |
|
---|
504 | /* Hock - Schittkowski (modified #1) problem 52 (box/linearly constrained), minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)
|
---|
505 | * constr1: p[0] + 3*p[1] = 0;
|
---|
506 | * constr2: p[2] + p[3] - 2*p[4] = 0;
|
---|
507 | * constr3: p[1] - p[4] = 0;
|
---|
508 | *
|
---|
509 | * To the above 3 constraints, we add the following 5:
|
---|
510 | * constr4: -0.09 <= p[0];
|
---|
511 | * constr5: 0.0 <= p[1] <= 0.3;
|
---|
512 | * constr6: p[2] <= 0.25;
|
---|
513 | * constr7: -0.2 <= p[3] <= 0.3;
|
---|
514 | * constr8: 0.0 <= p[4] <= 0.3;
|
---|
515 | *
|
---|
516 | */
|
---|
517 | void mod1hs52(double *p, double *x, int m, int n, void *data)
|
---|
518 | {
|
---|
519 | x[0]=4.0*p[0]-p[1];
|
---|
520 | x[1]=p[1]+p[2]-2.0;
|
---|
521 | x[2]=p[3]-1.0;
|
---|
522 | x[3]=p[4]-1.0;
|
---|
523 | }
|
---|
524 |
|
---|
525 | void jacmod1hs52(double *p, double *jac, int m, int n, void *data)
|
---|
526 | {
|
---|
527 | register int j=0;
|
---|
528 |
|
---|
529 | jac[j++]=4.0;
|
---|
530 | jac[j++]=-1.0;
|
---|
531 | jac[j++]=0.0;
|
---|
532 | jac[j++]=0.0;
|
---|
533 | jac[j++]=0.0;
|
---|
534 |
|
---|
535 | jac[j++]=0.0;
|
---|
536 | jac[j++]=1.0;
|
---|
537 | jac[j++]=1.0;
|
---|
538 | jac[j++]=0.0;
|
---|
539 | jac[j++]=0.0;
|
---|
540 |
|
---|
541 | jac[j++]=0.0;
|
---|
542 | jac[j++]=0.0;
|
---|
543 | jac[j++]=0.0;
|
---|
544 | jac[j++]=1.0;
|
---|
545 | jac[j++]=0.0;
|
---|
546 |
|
---|
547 | jac[j++]=0.0;
|
---|
548 | jac[j++]=0.0;
|
---|
549 | jac[j++]=0.0;
|
---|
550 | jac[j++]=0.0;
|
---|
551 | jac[j++]=1.0;
|
---|
552 | }
|
---|
553 |
|
---|
554 |
|
---|
555 | /* Hock - Schittkowski (modified #2) problem 52 (linear inequality constrained), minimum at (0.5, 2.0, 0.0, 1.0, 1.0)
|
---|
556 | * A fifth term [(p[0]-0.5)^2] is added to the objective function and
|
---|
557 | * the equality contraints are replaced by the following inequalities:
|
---|
558 | * constr1: p[0] + 3*p[1] >= -1.0;
|
---|
559 | * constr2: p[2] + p[3] - 2*p[4] >= -2.0;
|
---|
560 | * constr3: p[1] - p[4] <= 7.0;
|
---|
561 | *
|
---|
562 | *
|
---|
563 | */
|
---|
564 | void mod2hs52(double *p, double *x, int m, int n, void *data)
|
---|
565 | {
|
---|
566 | x[0]=4.0*p[0]-p[1];
|
---|
567 | x[1]=p[1]+p[2]-2.0;
|
---|
568 | x[2]=p[3]-1.0;
|
---|
569 | x[3]=p[4]-1.0;
|
---|
570 | x[4]=p[0]-0.5;
|
---|
571 | }
|
---|
572 |
|
---|
573 | void jacmod2hs52(double *p, double *jac, int m, int n, void *data)
|
---|
574 | {
|
---|
575 | register int j=0;
|
---|
576 |
|
---|
577 | jac[j++]=4.0;
|
---|
578 | jac[j++]=-1.0;
|
---|
579 | jac[j++]=0.0;
|
---|
580 | jac[j++]=0.0;
|
---|
581 | jac[j++]=0.0;
|
---|
582 |
|
---|
583 | jac[j++]=0.0;
|
---|
584 | jac[j++]=1.0;
|
---|
585 | jac[j++]=1.0;
|
---|
586 | jac[j++]=0.0;
|
---|
587 | jac[j++]=0.0;
|
---|
588 |
|
---|
589 | jac[j++]=0.0;
|
---|
590 | jac[j++]=0.0;
|
---|
591 | jac[j++]=0.0;
|
---|
592 | jac[j++]=1.0;
|
---|
593 | jac[j++]=0.0;
|
---|
594 |
|
---|
595 | jac[j++]=0.0;
|
---|
596 | jac[j++]=0.0;
|
---|
597 | jac[j++]=0.0;
|
---|
598 | jac[j++]=0.0;
|
---|
599 | jac[j++]=1.0;
|
---|
600 |
|
---|
601 | jac[j++]=1.0;
|
---|
602 | jac[j++]=0.0;
|
---|
603 | jac[j++]=0.0;
|
---|
604 | jac[j++]=0.0;
|
---|
605 | jac[j++]=0.0;
|
---|
606 | }
|
---|
607 |
|
---|
608 | /* Schittkowski (modified) problem 235 (box/linearly constrained), minimum at (-1.725, 2.9, 0.725)
|
---|
609 | * constr1: p[0] + p[2] = -1.0;
|
---|
610 | *
|
---|
611 | * To the above constraint, we add the following 2:
|
---|
612 | * constr2: p[1] - 4*p[2] = 0;
|
---|
613 | * constr3: 0.1 <= p[1] <= 2.9;
|
---|
614 | * constr4: 0.7 <= p[2];
|
---|
615 | *
|
---|
616 | */
|
---|
617 | void mods235(double *p, double *x, int m, int n, void *data)
|
---|
618 | {
|
---|
619 | x[0]=0.1*(p[0]-1.0);
|
---|
620 | x[1]=p[1]-p[0]*p[0];
|
---|
621 | }
|
---|
622 |
|
---|
623 | void jacmods235(double *p, double *jac, int m, int n, void *data)
|
---|
624 | {
|
---|
625 | register int j=0;
|
---|
626 |
|
---|
627 | jac[j++]=0.1;
|
---|
628 | jac[j++]=0.0;
|
---|
629 | jac[j++]=0.0;
|
---|
630 |
|
---|
631 | jac[j++]=-2.0*p[0];
|
---|
632 | jac[j++]=1.0;
|
---|
633 | jac[j++]=0.0;
|
---|
634 | }
|
---|
635 |
|
---|
636 | /* Boggs and Tolle modified problem 7 (box/linearly constrained), minimum at (0.7, 0.49, 0.19, 1.19, -0.2)
|
---|
637 | * We keep the original objective function & starting point and use the following constraints:
|
---|
638 | *
|
---|
639 | * subject to cons1:
|
---|
640 | * x[1]+x[2] - x[3] = 1.0;
|
---|
641 | * subject to cons2:
|
---|
642 | * x[2] - x[4] + x[1] = 0.0;
|
---|
643 | * subject to cons3:
|
---|
644 | * x[5] + x[1] = 0.5;
|
---|
645 | * subject to cons4:
|
---|
646 | * x[5]>=-0.3;
|
---|
647 | * subject to cons5:
|
---|
648 | * x[1]<=0.7;
|
---|
649 | *
|
---|
650 | */
|
---|
651 | void modbt7(double *p, double *x, int m, int n, void *data)
|
---|
652 | {
|
---|
653 | register int i;
|
---|
654 |
|
---|
655 | for(i=0; i<n; ++i)
|
---|
656 | x[i]=100.0*(p[1]-p[0]*p[0])*(p[1]-p[0]*p[0]) + (p[0]-1.0)*(p[0]-1.0);
|
---|
657 | }
|
---|
658 |
|
---|
659 | void jacmodbt7(double *p, double *jac, int m, int n, void *data)
|
---|
660 | {
|
---|
661 | register int i, j;
|
---|
662 |
|
---|
663 | for(i=j=0; i<m; ++i){
|
---|
664 | jac[j++]=-400.0*(p[1]-p[0]*p[0])*p[0] + 2.0*p[0] - 2.0;
|
---|
665 | jac[j++]=200.0*(p[1]-p[0]*p[0]);
|
---|
666 | jac[j++]=0.0;
|
---|
667 | jac[j++]=0.0;
|
---|
668 | jac[j++]=0.0;
|
---|
669 | }
|
---|
670 | }
|
---|
671 |
|
---|
672 | /* Equilibrium combustion problem, constrained nonlinear equation from the book by Floudas et al.
|
---|
673 | * Minimum at (0.0034, 31.3265, 0.0684, 0.8595, 0.0370)
|
---|
674 | * constri: p[i]>=0.0001; (i=1..5)
|
---|
675 | * constri+5: p[i]<=100.0; (i=1..5)
|
---|
676 | */
|
---|
677 | void combust(double *p, double *x, int m, int n, void *data)
|
---|
678 | {
|
---|
679 | double R, R5, R6, R7, R8, R9, R10;
|
---|
680 |
|
---|
681 | R=10;
|
---|
682 | R5=0.193;
|
---|
683 | R6=4.10622*1e-4;
|
---|
684 | R7=5.45177*1e-4;
|
---|
685 | R8=4.4975*1e-7;
|
---|
686 | R9=3.40735*1e-5;
|
---|
687 | R10=9.615*1e-7;
|
---|
688 |
|
---|
689 | x[0]=p[0]*p[1]+p[0]-3*p[4];
|
---|
690 | x[1]=2*p[0]*p[1]+p[0]+3*R10*p[1]*p[1]+p[1]*p[2]*p[2]+R7*p[1]*p[2]+R9*p[1]*p[3]+R8*p[1]-R*p[4];
|
---|
691 | x[2]=2*p[1]*p[2]*p[2]+R7*p[1]*p[2]+2*R5*p[2]*p[2]+R6*p[2]-8*p[4];
|
---|
692 | x[3]=R9*p[1]*p[3]+2*p[3]*p[3]-4*R*p[4];
|
---|
693 | x[4]=p[0]*p[1]+p[0]+R10*p[1]*p[1]+p[1]*p[2]*p[2]+R7*p[1]*p[2]+R9*p[1]*p[3]+R8*p[1]+R5*p[2]*p[2]+R6*p[2]+p[3]*p[3]-1.0;
|
---|
694 | }
|
---|
695 |
|
---|
696 | void jaccombust(double *p, double *jac, int m, int n, void *data)
|
---|
697 | {
|
---|
698 | register int j=0;
|
---|
699 | double R, R5, R6, R7, R8, R9, R10;
|
---|
700 |
|
---|
701 | R=10;
|
---|
702 | R5=0.193;
|
---|
703 | R6=4.10622*1e-4;
|
---|
704 | R7=5.45177*1e-4;
|
---|
705 | R8=4.4975*1e-7;
|
---|
706 | R9=3.40735*1e-5;
|
---|
707 | R10=9.615*1e-7;
|
---|
708 |
|
---|
709 | for(j=0; j<m*n; ++j) jac[j]=0.0;
|
---|
710 |
|
---|
711 | j=0;
|
---|
712 | jac[j]=p[1]+1;
|
---|
713 | jac[j+1]=p[0];
|
---|
714 | jac[j+4]=-3;
|
---|
715 |
|
---|
716 | j+=m;
|
---|
717 | jac[j]=2*p[1]+1;
|
---|
718 | jac[j+1]=2*p[0]+6*R10*p[1]+p[2]*p[2]+R7*p[2]+R9*p[3]+R8;
|
---|
719 | jac[j+2]=2*p[1]*p[2]+R7*p[1];
|
---|
720 | jac[j+3]=R9*p[1];
|
---|
721 | jac[j+4]=-R;
|
---|
722 |
|
---|
723 | j+=m;
|
---|
724 | jac[j+1]=2*p[2]*p[2]+R7*p[2];
|
---|
725 | jac[j+2]=4*p[1]*p[2]+R7*p[1]+4*R5*p[2]+R6;
|
---|
726 | jac[j+4]=-8;
|
---|
727 |
|
---|
728 | j+=m;
|
---|
729 | jac[j+1]=R9*p[3];
|
---|
730 | jac[j+3]=R9*p[1]+4*p[3];
|
---|
731 | jac[j+4]=-4*R;
|
---|
732 |
|
---|
733 | j+=m;
|
---|
734 | jac[j]=p[1]+1;
|
---|
735 | jac[j+1]=p[0]+2*R10*p[1]+p[2]*p[2]+R7*p[2]+R9*p[3]+R8;
|
---|
736 | jac[j+2]=2*p[1]*p[2]+R7*p[1]+2*R5*p[2]+R6;
|
---|
737 | jac[j+3]=R9*p[1]+2*p[3];
|
---|
738 | }
|
---|
739 |
|
---|
740 | /* Hock - Schittkowski (modified) problem 76 (linear inequalities & equations constrained), minimum at (0.0, 0.00909091, 0.372727, 0.354545)
|
---|
741 | * The non-squared terms in the objective function have been removed, the rhs of constr2 has been changed to 0.4 (from 4)
|
---|
742 | * and constr3 has been changed to an equation.
|
---|
743 | *
|
---|
744 | * constr1: p[0] + 2*p[1] + p[2] + p[3] <= 5;
|
---|
745 | * constr2: 3*p[0] + p[1] + 2*p[2] - p[3] <= 0.4;
|
---|
746 | * constr3: p[1] + 4*p[2] = 1.5;
|
---|
747 | *
|
---|
748 | */
|
---|
749 | void modhs76(double *p, double *x, int m, int n, void *data)
|
---|
750 | {
|
---|
751 | x[0]=p[0];
|
---|
752 | x[1]=sqrt(0.5)*p[1];
|
---|
753 | x[2]=p[2];
|
---|
754 | x[3]=sqrt(0.5)*p[3];
|
---|
755 | }
|
---|
756 |
|
---|
757 | void jacmodhs76(double *p, double *jac, int m, int n, void *data)
|
---|
758 | {
|
---|
759 | register int j=0;
|
---|
760 |
|
---|
761 | jac[j++]=1.0;
|
---|
762 | jac[j++]=0.0;
|
---|
763 | jac[j++]=0.0;
|
---|
764 | jac[j++]=0.0;
|
---|
765 |
|
---|
766 | jac[j++]=0.0;
|
---|
767 | jac[j++]=sqrt(0.5);
|
---|
768 | jac[j++]=0.0;
|
---|
769 | jac[j++]=0.0;
|
---|
770 |
|
---|
771 | jac[j++]=0.0;
|
---|
772 | jac[j++]=0.0;
|
---|
773 | jac[j++]=1.0;
|
---|
774 | jac[j++]=0.0;
|
---|
775 |
|
---|
776 | jac[j++]=0.0;
|
---|
777 | jac[j++]=0.0;
|
---|
778 | jac[j++]=0.0;
|
---|
779 | jac[j++]=sqrt(0.5);
|
---|
780 | }
|
---|
781 |
|
---|
782 |
|
---|
783 |
|
---|
784 | int main()
|
---|
785 | {
|
---|
786 | register int i, j;
|
---|
787 | int problem, ret;
|
---|
788 | double p[5], // 5 is max(2, 3, 5)
|
---|
789 | x[16]; // 16 is max(2, 3, 5, 6, 16)
|
---|
790 | int m, n;
|
---|
791 | double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
|
---|
792 | char *probname[]={
|
---|
793 | "Rosenbrock function",
|
---|
794 | "modified Rosenbrock problem",
|
---|
795 | "Powell's function",
|
---|
796 | "Wood's function",
|
---|
797 | "Meyer's (reformulated) problem",
|
---|
798 | "Osborne's problem",
|
---|
799 | "helical valley function",
|
---|
800 | "Boggs & Tolle's problem #3",
|
---|
801 | "Hock - Schittkowski problem #28",
|
---|
802 | "Hock - Schittkowski problem #48",
|
---|
803 | "Hock - Schittkowski problem #51",
|
---|
804 | "Hock - Schittkowski problem #01",
|
---|
805 | "Hock - Schittkowski modified problem #21",
|
---|
806 | "hatfldb problem",
|
---|
807 | "hatfldc problem",
|
---|
808 | "equilibrium combustion problem",
|
---|
809 | "Hock - Schittkowski modified #1 problem #52",
|
---|
810 | "Schittkowski modified problem #235",
|
---|
811 | "Boggs & Tolle modified problem #7",
|
---|
812 | "Hock - Schittkowski modified #2 problem #52",
|
---|
813 | "Hock - Schittkowski modified problem #76",
|
---|
814 | };
|
---|
815 |
|
---|
816 | opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
|
---|
817 | opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
|
---|
818 | //opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute!
|
---|
819 |
|
---|
820 | /* uncomment the appropriate line below to select a minimization problem */
|
---|
821 | problem=
|
---|
822 | //0; // Rosenbrock function
|
---|
823 | //1; // modified Rosenbrock problem
|
---|
824 | //2; // Powell's function
|
---|
825 | //3; // Wood's function
|
---|
826 | 4; // Meyer's (reformulated) problem
|
---|
827 | //5; // Osborne's problem
|
---|
828 | //6; // helical valley function
|
---|
829 | #ifdef HAVE_LAPACK
|
---|
830 | //7; // Boggs & Tolle's problem 3
|
---|
831 | //8; // Hock - Schittkowski problem 28
|
---|
832 | //9; // Hock - Schittkowski problem 48
|
---|
833 | //10; // Hock - Schittkowski problem 51
|
---|
834 | #else // no LAPACK
|
---|
835 | #ifdef _MSC_VER
|
---|
836 | #pragma message("LAPACK not available, some test problems cannot be used")
|
---|
837 | #else
|
---|
838 | #warning LAPACK not available, some test problems cannot be used
|
---|
839 | #endif // _MSC_VER
|
---|
840 |
|
---|
841 | #endif /* HAVE_LAPACK */
|
---|
842 | //11; // Hock - Schittkowski problem 01
|
---|
843 | //12; // Hock - Schittkowski modified problem 21
|
---|
844 | //13; // hatfldb problem
|
---|
845 | //14; // hatfldc problem
|
---|
846 | //15; // equilibrium combustion problem
|
---|
847 | #ifdef HAVE_LAPACK
|
---|
848 | //16; // Hock - Schittkowski modified #1 problem 52
|
---|
849 | //17; // Schittkowski modified problem 235
|
---|
850 | //18; // Boggs & Tolle modified problem #7
|
---|
851 | //19; // Hock - Schittkowski modified #2 problem 52
|
---|
852 | //20; // Hock - Schittkowski modified problem #76"
|
---|
853 | #endif /* HAVE_LAPACK */
|
---|
854 |
|
---|
855 | switch(problem){
|
---|
856 | default: fprintf(stderr, "unknown problem specified (#%d)! Note that some minimization problems require LAPACK.\n", problem);
|
---|
857 | exit(1);
|
---|
858 | break;
|
---|
859 |
|
---|
860 | case 0:
|
---|
861 | /* Rosenbrock function */
|
---|
862 | m=2; n=2;
|
---|
863 | p[0]=-1.2; p[1]=1.0;
|
---|
864 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
865 | ret=dlevmar_der(ros, jacros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
866 | //ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
|
---|
867 | break;
|
---|
868 |
|
---|
869 | case 1:
|
---|
870 | /* modified Rosenbrock problem */
|
---|
871 | m=2; n=3;
|
---|
872 | p[0]=-1.2; p[1]=1.0;
|
---|
873 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
874 | ret=dlevmar_der(modros, jacmodros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
875 | //ret=dlevmar_dif(modros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
|
---|
876 | break;
|
---|
877 |
|
---|
878 | case 2:
|
---|
879 | /* Powell's function */
|
---|
880 | m=2; n=2;
|
---|
881 | p[0]=3.0; p[1]=1.0;
|
---|
882 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
883 | ret=dlevmar_der(powell, jacpowell, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
884 | //ret=dlevmar_dif(powell, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
|
---|
885 | break;
|
---|
886 |
|
---|
887 | case 3:
|
---|
888 | /* Wood's function */
|
---|
889 | m=4; n=6;
|
---|
890 | p[0]=-3.0; p[1]=-1.0; p[2]=-3.0; p[3]=-1.0;
|
---|
891 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
892 | ret=dlevmar_dif(wood, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
|
---|
893 | break;
|
---|
894 |
|
---|
895 | case 4:
|
---|
896 | /* Meyer's data fitting problem */
|
---|
897 | m=3; n=16;
|
---|
898 | p[0]=8.85; p[1]=4.0; p[2]=2.5;
|
---|
899 | x[0]=34.780; x[1]=28.610; x[2]=23.650; x[3]=19.630;
|
---|
900 | x[4]=16.370; x[5]=13.720; x[6]=11.540; x[7]=9.744;
|
---|
901 | x[8]=8.261; x[9]=7.030; x[10]=6.005; x[11]=5.147;
|
---|
902 | x[12]=4.427; x[13]=3.820; x[14]=3.307; x[15]=2.872;
|
---|
903 | //ret=dlevmar_der(meyer, jacmeyer, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
904 |
|
---|
905 | { double *work, *covar;
|
---|
906 | work=malloc((LM_DIF_WORKSZ(m, n)+m*m)*sizeof(double));
|
---|
907 | if(!work){
|
---|
908 | fprintf(stderr, "memory allocation request failed in main()\n");
|
---|
909 | exit(1);
|
---|
910 | }
|
---|
911 | covar=work+LM_DIF_WORKSZ(m, n);
|
---|
912 |
|
---|
913 | ret=dlevmar_dif(meyer, p, x, m, n, 1000, opts, info, work, covar, NULL); // no Jacobian, caller allocates work memory, covariance estimated
|
---|
914 |
|
---|
915 | printf("Covariance of the fit:\n");
|
---|
916 | for(i=0; i<m; ++i){
|
---|
917 | for(j=0; j<m; ++j)
|
---|
918 | printf("%g ", covar[i*m+j]);
|
---|
919 | printf("\n");
|
---|
920 | }
|
---|
921 | printf("\n");
|
---|
922 |
|
---|
923 | free(work);
|
---|
924 | }
|
---|
925 |
|
---|
926 | /* uncomment the following block to verify Jacobian */
|
---|
927 | /*
|
---|
928 | {
|
---|
929 | double err[16];
|
---|
930 | dlevmar_chkjac(meyer, jacmeyer, p, m, n, NULL, err);
|
---|
931 | for(i=0; i<n; ++i) printf("gradient %d, err %g\n", i, err[i]);
|
---|
932 | }
|
---|
933 | */
|
---|
934 | break;
|
---|
935 |
|
---|
936 | case 5:
|
---|
937 | /* Osborne's data fitting problem */
|
---|
938 | {
|
---|
939 | double x33[]={
|
---|
940 | 8.44E-1, 9.08E-1, 9.32E-1, 9.36E-1, 9.25E-1, 9.08E-1, 8.81E-1,
|
---|
941 | 8.5E-1, 8.18E-1, 7.84E-1, 7.51E-1, 7.18E-1, 6.85E-1, 6.58E-1,
|
---|
942 | 6.28E-1, 6.03E-1, 5.8E-1, 5.58E-1, 5.38E-1, 5.22E-1, 5.06E-1,
|
---|
943 | 4.9E-1, 4.78E-1, 4.67E-1, 4.57E-1, 4.48E-1, 4.38E-1, 4.31E-1,
|
---|
944 | 4.24E-1, 4.2E-1, 4.14E-1, 4.11E-1, 4.06E-1};
|
---|
945 |
|
---|
946 | m=5; n=33;
|
---|
947 | p[0]=0.5; p[1]=1.5; p[2]=-1.0; p[3]=1.0E-2; p[4]=2.0E-2;
|
---|
948 |
|
---|
949 | ret=dlevmar_der(osborne, jacosborne, p, x33, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
950 | //ret=dlevmar_dif(osborne, p, x33, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
|
---|
951 | }
|
---|
952 | break;
|
---|
953 |
|
---|
954 | case 6:
|
---|
955 | /* helical valley function */
|
---|
956 | m=3; n=3;
|
---|
957 | p[0]=-1.0; p[1]=0.0; p[2]=0.0;
|
---|
958 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
959 | ret=dlevmar_der(helval, jachelval, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
960 | //ret=dlevmar_dif(helval, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
|
---|
961 | break;
|
---|
962 |
|
---|
963 | #ifdef HAVE_LAPACK
|
---|
964 | case 7:
|
---|
965 | /* Boggs-Tolle problem 3 */
|
---|
966 | m=5; n=5;
|
---|
967 | p[0]=2.0; p[1]=2.0; p[2]=2.0;
|
---|
968 | p[3]=2.0; p[4]=2.0;
|
---|
969 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
970 |
|
---|
971 | {
|
---|
972 | double A[3*5]={1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -2.0, 0.0, 1.0, 0.0, 0.0, -1.0},
|
---|
973 | b[3]={0.0, 0.0, 0.0};
|
---|
974 |
|
---|
975 | ret=dlevmar_lec_der(bt3, jacbt3, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic Jacobian
|
---|
976 | //ret=dlevmar_lec_dif(bt3, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
|
---|
977 | }
|
---|
978 | break;
|
---|
979 |
|
---|
980 | case 8:
|
---|
981 | /* Hock - Schittkowski problem 28 */
|
---|
982 | m=3; n=3;
|
---|
983 | p[0]=-4.0; p[1]=1.0; p[2]=1.0;
|
---|
984 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
985 |
|
---|
986 | {
|
---|
987 | double A[1*3]={1.0, 2.0, 3.0},
|
---|
988 | b[1]={1.0};
|
---|
989 |
|
---|
990 | ret=dlevmar_lec_der(hs28, jachs28, p, x, m, n, A, b, 1, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic Jacobian
|
---|
991 | //ret=dlevmar_lec_dif(hs28, p, x, m, n, A, b, 1, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
|
---|
992 | }
|
---|
993 | break;
|
---|
994 |
|
---|
995 | case 9:
|
---|
996 | /* Hock - Schittkowski problem 48 */
|
---|
997 | m=5; n=5;
|
---|
998 | p[0]=3.0; p[1]=5.0; p[2]=-3.0;
|
---|
999 | p[3]=2.0; p[4]=-2.0;
|
---|
1000 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1001 |
|
---|
1002 | {
|
---|
1003 | double A[2*5]={1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, -2.0, -2.0},
|
---|
1004 | b[2]={5.0, -3.0};
|
---|
1005 |
|
---|
1006 | ret=dlevmar_lec_der(hs48, jachs48, p, x, m, n, A, b, 2, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic Jacobian
|
---|
1007 | //ret=dlevmar_lec_dif(hs48, p, x, m, n, A, b, 2, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
|
---|
1008 | }
|
---|
1009 | break;
|
---|
1010 |
|
---|
1011 | case 10:
|
---|
1012 | /* Hock - Schittkowski problem 51 */
|
---|
1013 | m=5; n=5;
|
---|
1014 | p[0]=2.5; p[1]=0.5; p[2]=2.0;
|
---|
1015 | p[3]=-1.0; p[4]=0.5;
|
---|
1016 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1017 |
|
---|
1018 | {
|
---|
1019 | double A[3*5]={1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -2.0, 0.0, 1.0, 0.0, 0.0, -1.0},
|
---|
1020 | b[3]={4.0, 0.0, 0.0};
|
---|
1021 |
|
---|
1022 | ret=dlevmar_lec_der(hs51, jachs51, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic Jacobian
|
---|
1023 | //ret=dlevmar_lec_dif(hs51, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
|
---|
1024 | }
|
---|
1025 | break;
|
---|
1026 |
|
---|
1027 | #endif /* HAVE_LAPACK */
|
---|
1028 |
|
---|
1029 | case 11:
|
---|
1030 | /* Hock - Schittkowski problem 01 */
|
---|
1031 | m=2; n=2;
|
---|
1032 | p[0]=-2.0; p[1]=1.0;
|
---|
1033 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1034 | //ret=dlevmar_der(hs01, jachs01, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1035 | {
|
---|
1036 | double lb[2], ub[2];
|
---|
1037 |
|
---|
1038 | lb[0]=-DBL_MAX; lb[1]=-1.5;
|
---|
1039 | ub[0]=ub[1]=DBL_MAX;
|
---|
1040 |
|
---|
1041 | ret=dlevmar_bc_der(hs01, jachs01, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1042 | }
|
---|
1043 | break;
|
---|
1044 |
|
---|
1045 | case 12:
|
---|
1046 | /* Hock - Schittkowski (modified) problem 21 */
|
---|
1047 | m=2; n=2;
|
---|
1048 | p[0]=-1.0; p[1]=-1.0;
|
---|
1049 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1050 | //ret=dlevmar_der(hs21, jachs21, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1051 | {
|
---|
1052 | double lb[2], ub[2];
|
---|
1053 |
|
---|
1054 | lb[0]=2.0; lb[1]=-50.0;
|
---|
1055 | ub[0]=50.0; ub[1]=50.0;
|
---|
1056 |
|
---|
1057 | ret=dlevmar_bc_der(hs21, jachs21, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1058 | }
|
---|
1059 | break;
|
---|
1060 |
|
---|
1061 | case 13:
|
---|
1062 | /* hatfldb problem */
|
---|
1063 | m=4; n=4;
|
---|
1064 | p[0]=p[1]=p[2]=p[3]=0.1;
|
---|
1065 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1066 | //ret=dlevmar_der(hatfldb, jachatfldb, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1067 | {
|
---|
1068 | double lb[4], ub[4];
|
---|
1069 |
|
---|
1070 | lb[0]=lb[1]=lb[2]=lb[3]=0.0;
|
---|
1071 |
|
---|
1072 | ub[0]=ub[2]=ub[3]=DBL_MAX;
|
---|
1073 | ub[1]=0.8;
|
---|
1074 |
|
---|
1075 | ret=dlevmar_bc_der(hatfldb, jachatfldb, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1076 | }
|
---|
1077 | break;
|
---|
1078 |
|
---|
1079 | case 14:
|
---|
1080 | /* hatfldc problem */
|
---|
1081 | m=4; n=4;
|
---|
1082 | p[0]=p[1]=p[2]=p[3]=0.9;
|
---|
1083 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1084 | //ret=dlevmar_der(hatfldc, jachatfldc, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1085 | {
|
---|
1086 | double lb[4], ub[4];
|
---|
1087 |
|
---|
1088 | lb[0]=lb[1]=lb[2]=lb[3]=0.0;
|
---|
1089 |
|
---|
1090 | ub[0]=ub[1]=ub[2]=ub[3]=10.0;
|
---|
1091 |
|
---|
1092 | ret=dlevmar_bc_der(hatfldc, jachatfldc, p, x, m, n, lb, ub, NULL, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1093 | }
|
---|
1094 | break;
|
---|
1095 |
|
---|
1096 | case 15:
|
---|
1097 | /* equilibrium combustion problem */
|
---|
1098 | m=5; n=5;
|
---|
1099 | p[0]=p[1]=p[2]=p[3]=p[4]=0.0001;
|
---|
1100 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1101 | //ret=dlevmar_der(combust, jaccombust, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1102 | {
|
---|
1103 | double lb[5], ub[5];
|
---|
1104 |
|
---|
1105 | lb[0]=lb[1]=lb[2]=lb[3]=lb[4]=0.0001;
|
---|
1106 |
|
---|
1107 | ub[0]=ub[1]=ub[2]=ub[3]=ub[4]=100.0;
|
---|
1108 |
|
---|
1109 | ret=dlevmar_bc_der(combust, jaccombust, p, x, m, n, lb, ub, NULL, 5000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
|
---|
1110 | }
|
---|
1111 | break;
|
---|
1112 |
|
---|
1113 | #ifdef HAVE_LAPACK
|
---|
1114 | case 16:
|
---|
1115 | /* Hock - Schittkowski modified #1 problem 52 */
|
---|
1116 | m=5; n=4;
|
---|
1117 | p[0]=2.0; p[1]=2.0; p[2]=2.0;
|
---|
1118 | p[3]=2.0; p[4]=2.0;
|
---|
1119 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1120 |
|
---|
1121 | {
|
---|
1122 | double A[3*5]={1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -2.0, 0.0, 1.0, 0.0, 0.0, -1.0},
|
---|
1123 | b[3]={0.0, 0.0, 0.0};
|
---|
1124 |
|
---|
1125 | double lb[5], ub[5];
|
---|
1126 |
|
---|
1127 | double weights[5]={2000.0, 2000.0, 2000.0, 2000.0, 2000.0}; // penalty terms weights
|
---|
1128 |
|
---|
1129 | lb[0]=-0.09; lb[1]=0.0; lb[2]=-DBL_MAX; lb[3]=-0.2; lb[4]=0.0;
|
---|
1130 | ub[0]=DBL_MAX; ub[1]=0.3; ub[2]=0.25; ub[3]=0.3; ub[4]=0.3;
|
---|
1131 |
|
---|
1132 | ret=dlevmar_blec_der(mod1hs52, jacmod1hs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, analytic Jacobian
|
---|
1133 | //ret=dlevmar_blec_dif(mod1hs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
|
---|
1134 | }
|
---|
1135 | break;
|
---|
1136 |
|
---|
1137 | case 17:
|
---|
1138 | /* Schittkowski modified problem 235 */
|
---|
1139 | m=3; n=2;
|
---|
1140 | p[0]=-2.0; p[1]=3.0; p[2]=1.0;
|
---|
1141 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1142 |
|
---|
1143 | {
|
---|
1144 | double A[2*3]={1.0, 0.0, 1.0, 0.0, 1.0, -4.0},
|
---|
1145 | b[2]={-1.0, 0.0};
|
---|
1146 |
|
---|
1147 | double lb[3], ub[3];
|
---|
1148 |
|
---|
1149 | lb[0]=-DBL_MAX; lb[1]=0.1; lb[2]=0.7;
|
---|
1150 | ub[0]=DBL_MAX; ub[1]=2.9; ub[2]=DBL_MAX;
|
---|
1151 |
|
---|
1152 | ret=dlevmar_blec_der(mods235, jacmods235, p, x, m, n, lb, ub, A, b, 2, NULL, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, analytic Jacobian
|
---|
1153 | //ret=dlevmar_blec_dif(mods235, p, x, m, n, lb, ub, A, b, 2, NULL, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
|
---|
1154 | }
|
---|
1155 | break;
|
---|
1156 |
|
---|
1157 | case 18:
|
---|
1158 | /* Boggs & Tolle modified problem 7 */
|
---|
1159 | m=5; n=5;
|
---|
1160 | p[0]=-2.0; p[1]=1.0; p[2]=1.0; p[3]=1.0; p[4]=1.0;
|
---|
1161 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1162 |
|
---|
1163 | {
|
---|
1164 | double A[3*5]={1.0, 1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0},
|
---|
1165 | b[3]={1.0, 0.0, 0.5};
|
---|
1166 |
|
---|
1167 | double lb[5], ub[5];
|
---|
1168 |
|
---|
1169 | lb[0]=-DBL_MAX; lb[1]=-DBL_MAX; lb[2]=-DBL_MAX; lb[3]=-DBL_MAX; lb[4]=-0.3;
|
---|
1170 | ub[0]=0.7; ub[1]= DBL_MAX; ub[2]= DBL_MAX; ub[3]= DBL_MAX; ub[4]=DBL_MAX;
|
---|
1171 |
|
---|
1172 | ret=dlevmar_blec_der(modbt7, jacmodbt7, p, x, m, n, lb, ub, A, b, 3, NULL, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, analytic Jacobian
|
---|
1173 | //ret=dlevmar_blec_dif(modbt7, p, x, m, n, lb, ub, A, b, 3, NULL, 10000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
|
---|
1174 | }
|
---|
1175 | break;
|
---|
1176 |
|
---|
1177 | case 19:
|
---|
1178 | /* Hock - Schittkowski modified #2 problem 52 */
|
---|
1179 | m=5; n=5;
|
---|
1180 | p[0]=2.0; p[1]=2.0; p[2]=2.0;
|
---|
1181 | p[3]=2.0; p[4]=2.0;
|
---|
1182 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1183 |
|
---|
1184 | {
|
---|
1185 | double C[3*5]={1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -2.0, 0.0, -1.0, 0.0, 0.0, 1.0},
|
---|
1186 | d[3]={-1.0, -2.0, -7.0};
|
---|
1187 |
|
---|
1188 | ret=dlevmar_bleic_der(mod2hs52, jacmod2hs52, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, 3, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, analytic Jacobian
|
---|
1189 | //ret=dlevmar_bleic_dif(mod2hs52, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, 3, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, no Jacobian
|
---|
1190 | }
|
---|
1191 | break;
|
---|
1192 |
|
---|
1193 | case 20:
|
---|
1194 | /* Hock - Schittkowski modified problem 76 */
|
---|
1195 | m=4; n=4;
|
---|
1196 | p[0]=0.5; p[1]=0.5; p[2]=0.5; p[3]=0.5;
|
---|
1197 | for(i=0; i<n; i++) x[i]=0.0;
|
---|
1198 |
|
---|
1199 | {
|
---|
1200 | double A[1*4]={0.0, 1.0, 4.0, 0.0},
|
---|
1201 | b[1]={1.5};
|
---|
1202 |
|
---|
1203 | double C[2*4]={-1.0, -2.0, -1.0, -1.0, -3.0, -1.0, -2.0, 1.0},
|
---|
1204 | d[2]={-5.0, -0.4};
|
---|
1205 |
|
---|
1206 | double lb[4]={0.0, 0.0, 0.0, 0.0};
|
---|
1207 |
|
---|
1208 | ret=dlevmar_bleic_der(modhs76, jacmodhs76, p, x, m, n, lb, NULL, A, b, 1, C, d, 2, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, analytic Jacobian
|
---|
1209 | //ret=dlevmar_bleic_dif(modhs76, p, x, m, n, lb, NULL, A, b, 1, C, d, 2, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, no Jacobian
|
---|
1210 | /* variations:
|
---|
1211 | * if no lb is used, the minimizer is (-0.1135922 0.1330097 0.3417476 0.07572816)
|
---|
1212 | * if the rhs of constr2 is 4.0, the minimizer is (0.0, 0.166667, 0.333333, 0.0)
|
---|
1213 | */
|
---|
1214 | }
|
---|
1215 | break;
|
---|
1216 |
|
---|
1217 | #endif /* HAVE_LAPACK */
|
---|
1218 | } /* switch */
|
---|
1219 |
|
---|
1220 | printf("Results for %s:\n", probname[problem]);
|
---|
1221 | printf("Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);
|
---|
1222 | for(i=0; i<m; ++i)
|
---|
1223 | printf("%.7g ", p[i]);
|
---|
1224 | printf("\n\nMinimization info:\n");
|
---|
1225 | for(i=0; i<LM_INFO_SZ; ++i)
|
---|
1226 | printf("%g ", info[i]);
|
---|
1227 | printf("\n");
|
---|
1228 |
|
---|
1229 | return 0;
|
---|
1230 | }
|
---|