Mercurial > octave
annotate scripts/optimization/sqp.m @ 9051:1bf0ce0930be
Grammar check TexInfo in all .m files
Cleanup documentation sources to follow a few consistent rules.
Spellcheck was NOT done. (but will be in another changeset)
author | Rik <rdrider0-list@yahoo.com> |
---|---|
date | Fri, 27 Mar 2009 22:31:03 -0700 |
parents | eb63fbe60fab |
children | 923c7cb7f13f |
rev | line source |
---|---|
8920 | 1 ## Copyright (C) 2005, 2006, 2007, 2008 John W. Eaton |
5289 | 2 ## |
3 ## This file is part of Octave. | |
4 ## | |
5 ## Octave is free software; you can redistribute it and/or modify it | |
6 ## under the terms of the GNU General Public License as published by | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
5289 | 9 ## |
10 ## Octave is distributed in the hope that it will be useful, but | |
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 ## General Public License for more details. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
5289 | 18 |
19 ## -*- texinfo -*- | |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
20 ## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{iter}, @var{nf}, @var{lambda}] =} sqp (@var{x}, @var{phi}, @var{g}, @var{h}, @var{lb}, @var{ub}, @var{maxiter}, @var{tolerance}) |
5289 | 21 ## Solve the nonlinear program |
6741 | 22 ## @iftex |
23 ## @tex | |
24 ## $$ | |
25 ## \min_x \phi (x) | |
26 ## $$ | |
27 ## @end tex | |
28 ## @end iftex | |
29 ## @ifnottex | |
5289 | 30 ## |
31 ## @example | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
32 ## @group |
5289 | 33 ## min phi (x) |
34 ## x | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
35 ## @end group |
5289 | 36 ## @end example |
37 ## | |
6741 | 38 ## @end ifnottex |
39 ## subject to | |
5289 | 40 ## @iftex |
41 ## @tex | |
6741 | 42 ## $$ |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
43 ## g(x) = 0 \qquad h(x) \geq 0 \qquad lb \leq x \leq ub |
6741 | 44 ## $$ |
5289 | 45 ## @end tex |
46 ## @end iftex | |
6741 | 47 ## @ifnottex |
5289 | 48 ## |
49 ## @example | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
50 ## @group |
5289 | 51 ## g(x) = 0 |
52 ## h(x) >= 0 | |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
53 ## lb <= x <= ub |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
54 ## @end group |
5289 | 55 ## @end example |
6741 | 56 ## @end ifnottex |
5289 | 57 ## |
58 ## @noindent | |
59 ## using a successive quadratic programming method. | |
60 ## | |
61 ## The first argument is the initial guess for the vector @var{x}. | |
62 ## | |
7001 | 63 ## The second argument is a function handle pointing to the objective |
5289 | 64 ## function. The objective function must be of the form |
65 ## | |
66 ## @example | |
67 ## y = phi (x) | |
68 ## @end example | |
69 ## | |
70 ## @noindent | |
71 ## in which @var{x} is a vector and @var{y} is a scalar. | |
72 ## | |
73 ## The second argument may also be a 2- or 3-element cell array of | |
74 ## function handles. The first element should point to the objective | |
75 ## function, the second should point to a function that computes the | |
76 ## gradient of the objective function, and the third should point to a | |
77 ## function to compute the hessian of the objective function. If the | |
78 ## gradient function is not supplied, the gradient is computed by finite | |
79 ## differences. If the hessian function is not supplied, a BFGS update | |
80 ## formula is used to approximate the hessian. | |
81 ## | |
82 ## If supplied, the gradient function must be of the form | |
83 ## | |
84 ## @example | |
7031 | 85 ## g = gradient (x) |
5289 | 86 ## @end example |
87 ## | |
88 ## @noindent | |
89 ## in which @var{x} is a vector and @var{g} is a vector. | |
90 ## | |
91 ## If supplied, the hessian function must be of the form | |
92 ## | |
93 ## @example | |
7031 | 94 ## h = hessian (x) |
5289 | 95 ## @end example |
96 ## | |
97 ## @noindent | |
98 ## in which @var{x} is a vector and @var{h} is a matrix. | |
99 ## | |
100 ## The third and fourth arguments are function handles pointing to | |
101 ## functions that compute the equality constraints and the inequality | |
102 ## constraints, respectively. | |
103 ## | |
104 ## If your problem does not have equality (or inequality) constraints, | |
105 ## you may pass an empty matrix for @var{cef} (or @var{cif}). | |
106 ## | |
107 ## If supplied, the equality and inequality constraint functions must be | |
108 ## of the form | |
109 ## | |
110 ## @example | |
7031 | 111 ## r = f (x) |
5289 | 112 ## @end example |
113 ## | |
114 ## @noindent | |
115 ## in which @var{x} is a vector and @var{r} is a vector. | |
116 ## | |
117 ## The third and fourth arguments may also be 2-element cell arrays of | |
118 ## function handles. The first element should point to the constraint | |
119 ## function and the second should point to a function that computes the | |
120 ## gradient of the constraint function: | |
121 ## | |
6741 | 122 ## @iftex |
123 ## @tex | |
124 ## $$ | |
125 ## \Bigg( {\partial f(x) \over \partial x_1}, | |
126 ## {\partial f(x) \over \partial x_2}, \ldots, | |
127 ## {\partial f(x) \over \partial x_N} \Bigg)^T | |
128 ## $$ | |
129 ## @end tex | |
130 ## @end iftex | |
131 ## @ifnottex | |
5289 | 132 ## @example |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
133 ## @group |
5289 | 134 ## [ d f(x) d f(x) d f(x) ] |
135 ## transpose ( [ ------ ----- ... ------ ] ) | |
136 ## [ dx_1 dx_2 dx_N ] | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
137 ## @end group |
5289 | 138 ## @end example |
6741 | 139 ## @end ifnottex |
5289 | 140 ## |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
141 ## The fifth and sixth arguments are vectors containing lower and upper bounds |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
142 ## on @var{x}. These must be consistent with equality and inequality |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
143 ## constraints @var{g} and @var{h}. If the bounds are not specified, or are |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
144 ## empty, they are set to -@var{realmax} and @var{realmax} by default. |
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
145 ## |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
146 ## The seventh argument is max. number of iterations. If not specified, |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
147 ## the default value is 100. |
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
148 ## |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
149 ## The eighth argument is tolerance for stopping criteria. If not specified, |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
150 ## the default value is @var{eps}. |
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
151 ## |
5289 | 152 ## Here is an example of calling @code{sqp}: |
153 ## | |
154 ## @example | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
155 ## @group |
7031 | 156 ## function r = g (x) |
157 ## r = [ sumsq(x)-10; | |
158 ## x(2)*x(3)-5*x(4)*x(5); | |
159 ## x(1)^3+x(2)^3+1 ]; | |
160 ## endfunction | |
161 ## | |
162 ## function obj = phi (x) | |
163 ## obj = exp(prod(x)) - 0.5*(x(1)^3+x(2)^3+1)^2; | |
164 ## endfunction | |
5289 | 165 ## |
7031 | 166 ## x0 = [-1.8; 1.7; 1.9; -0.8; -0.8]; |
167 ## | |
168 ## [x, obj, info, iter, nf, lambda] = sqp (x0, @@phi, @@g, []) | |
5289 | 169 ## |
7031 | 170 ## x = |
171 ## | |
172 ## -1.71714 | |
173 ## 1.59571 | |
174 ## 1.82725 | |
175 ## -0.76364 | |
176 ## -0.76364 | |
5289 | 177 ## |
7031 | 178 ## obj = 0.053950 |
179 ## info = 101 | |
180 ## iter = 8 | |
181 ## nf = 10 | |
182 ## lambda = | |
183 ## | |
184 ## -0.0401627 | |
185 ## 0.0379578 | |
186 ## -0.0052227 | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
187 ## @end group |
5289 | 188 ## @end example |
189 ## | |
190 ## The value returned in @var{info} may be one of the following: | |
191 ## @table @asis | |
192 ## @item 101 | |
193 ## The algorithm terminated because the norm of the last step was less | |
194 ## than @code{tol * norm (x))} (the value of tol is currently fixed at | |
195 ## @code{sqrt (eps)}---edit @file{sqp.m} to modify this value. | |
196 ## @item 102 | |
197 ## The BFGS update failed. | |
198 ## @item 103 | |
199 ## The maximum number of iterations was reached (the maximum number of | |
200 ## allowed iterations is currently fixed at 100---edit @file{sqp.m} to | |
201 ## increase this value). | |
202 ## @end table | |
5642 | 203 ## @seealso{qp} |
5289 | 204 ## @end deftypefn |
205 | |
6768 | 206 function [x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif, lb, ub, maxiter, tolerance) |
5289 | 207 |
6768 | 208 global __sqp_nfun__; |
5289 | 209 global __sqp_obj_fun__; |
210 global __sqp_ce_fun__; | |
211 global __sqp_ci_fun__; | |
6768 | 212 global __sqp_cif__; |
213 global __sqp_cifcn__; | |
5289 | 214 |
6768 | 215 if (nargin >= 2 && nargin <= 8 && nargin != 5) |
5289 | 216 |
217 ## Choose an initial NxN symmetric positive definite Hessan | |
218 ## approximation B. | |
219 | |
220 n = length (x); | |
221 | |
222 ## Evaluate objective function, constraints, and gradients at initial | |
223 ## value of x. | |
224 ## | |
225 ## obj_fun | |
226 ## obj_grad | |
227 ## ce_fun -- equality constraint functions | |
228 ## ci_fun -- inequality constraint functions | |
229 ## A == [grad_{x_1} cx_fun, grad_{x_2} cx_fun, ..., grad_{x_n} cx_fun]^T | |
230 | |
231 obj_grd = @fd_obj_grd; | |
232 have_hess = 0; | |
233 if (iscell (objf)) | |
234 if (length (objf) > 0) | |
235 __sqp_obj_fun__ = obj_fun = objf{1}; | |
236 if (length (objf) > 1) | |
237 obj_grd = objf{2}; | |
238 if (length (objf) > 2) | |
239 obj_hess = objf{3}; | |
240 have_hess = 1; | |
241 endif | |
242 endif | |
243 else | |
244 error ("sqp: invalid objective function"); | |
245 endif | |
246 else | |
247 __sqp_obj_fun__ = obj_fun = objf; | |
248 endif | |
249 | |
250 ce_fun = @empty_cf; | |
251 ce_grd = @empty_jac; | |
252 if (nargin > 2) | |
253 ce_grd = @fd_ce_jac; | |
254 if (iscell (cef)) | |
255 if (length (cef) > 0) | |
256 __sqp_ce_fun__ = ce_fun = cef{1}; | |
257 if (length (cef) > 1) | |
258 ce_grd = cef{2}; | |
259 endif | |
260 else | |
261 error ("sqp: invalid equality constraint function"); | |
262 endif | |
263 elseif (! isempty (cef)) | |
264 ce_fun = cef; | |
265 endif | |
266 endif | |
267 __sqp_ce_fun__ = ce_fun; | |
268 | |
269 ci_fun = @empty_cf; | |
270 ci_grd = @empty_jac; | |
6768 | 271 |
5289 | 272 if (nargin > 3) |
6768 | 273 ## constraint function given by user with possibly gradient |
274 __sqp_cif__ = cif; | |
275 ## constraint function given by user without gradient | |
276 __sqp_cifcn__ = @empty_cf; | |
277 if (iscell (__sqp_cif__)) | |
278 if (length (__sqp_cif__) > 0) | |
279 __sqp_cifcn__ = __sqp_cif__{1}; | |
280 endif | |
281 elseif (! isempty (__sqp_cif__)) | |
282 __sqp_cifcn__ = __sqp_cif__; | |
283 endif | |
284 | |
285 if (nargin < 5) | |
286 ci_grd = @fd_ci_jac; | |
287 if (iscell (cif)) | |
288 if (length (cif) > 0) | |
289 __sqp_ci_fun__ = ci_fun = cif{1}; | |
290 if (length (cif) > 1) | |
291 ci_grd = cif{2}; | |
292 endif | |
293 else | |
294 error ("sqp: invalid equality constraint function"); | |
5289 | 295 endif |
6768 | 296 elseif (! isempty (cif)) |
297 ci_fun = cif; | |
298 endif | |
299 else | |
300 global __sqp_lb__; | |
301 if (isvector (lb)) | |
302 __sqp_lb__ = lb; | |
303 elseif (isempty (lb)) | |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
304 if (isa (x, "single")) |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
305 __sqp_lb__ = -realmax ("single"); |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
306 else |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
307 __sqp_lb__ = -realmax; |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
308 endif |
5289 | 309 else |
6768 | 310 error ("sqp: invalid lower bound"); |
5289 | 311 endif |
6768 | 312 |
313 global __sqp_ub__; | |
314 if (isvector (ub)) | |
315 __sqp_ub__ = ub; | |
316 elseif (isempty (lb)) | |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
317 if (isa (x, "single")) |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
318 __sqp_ub__ = realmax ("single"); |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
319 else |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
320 __sqp_ub__ = realmax; |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7404
diff
changeset
|
321 endif |
6768 | 322 else |
323 error ("sqp: invalid upper bound"); | |
324 endif | |
325 | |
326 if (lb > ub) | |
327 error ("sqp: upper bound smaller than lower bound"); | |
328 endif | |
329 __sqp_ci_fun__ = ci_fun = @cf_ub_lb; | |
330 ci_grd = @cigrad_ub_lb; | |
331 endif | |
332 __sqp_ci_fun__ = ci_fun; | |
333 endif | |
334 | |
335 iter_max = 100; | |
336 if (nargin > 6 && ! isempty (maxiter)) | |
6921 | 337 if (isscalar (maxiter) && maxiter > 0 && round (maxiter) == maxiter) |
6768 | 338 iter_max = maxiter; |
339 else | |
340 error ("sqp: invalid number of maximum iterations"); | |
5289 | 341 endif |
342 endif | |
343 | |
6768 | 344 tol = sqrt (eps); |
345 if (nargin > 7 && ! isempty (tolerance)) | |
346 if (isscalar (tolerance) && tolerance > 0) | |
347 tol = tolerance; | |
348 else | |
349 error ("sqp: invalid value for tolerance"); | |
350 endif | |
351 endif | |
5289 | 352 |
353 iter = 0; | |
354 | |
355 obj = feval (obj_fun, x); | |
6768 | 356 __sqp_nfun__ = 1; |
5289 | 357 |
358 c = feval (obj_grd, x); | |
359 | |
6382 | 360 if (have_hess) |
361 B = feval (obj_hess, x); | |
362 else | |
363 B = eye (n, n); | |
364 endif | |
365 | |
5289 | 366 ce = feval (ce_fun, x); |
367 F = feval (ce_grd, x); | |
368 | |
369 ci = feval (ci_fun, x); | |
370 C = feval (ci_grd, x); | |
371 | |
372 A = [F; C]; | |
373 | |
374 ## Choose an initial lambda (x is provided by the caller). | |
375 | |
376 lambda = 100 * ones (rows (A), 1); | |
377 | |
378 qp_iter = 1; | |
379 alpha = 1; | |
380 | |
381 ## report (); | |
382 | |
6768 | 383 ## report (iter, qp_iter, alpha, __sqp_nfun__, obj); |
5289 | 384 |
6527 | 385 info = 0; |
386 | |
5289 | 387 while (++iter < iter_max) |
388 | |
389 ## Check convergence. This is just a simple check on the first | |
390 ## order necessary conditions. | |
391 | |
392 ## IDX is the indices of the active inequality constraints. | |
393 | |
394 nr_f = rows (F); | |
395 | |
396 lambda_e = lambda((1:nr_f)'); | |
397 lambda_i = lambda((nr_f+1:end)'); | |
398 | |
399 con = [ce; ci]; | |
400 | |
401 t0 = norm (c - A' * lambda); | |
402 t1 = norm (ce); | |
403 t2 = all (ci >= 0); | |
404 t3 = all (lambda_i >= 0); | |
405 t4 = norm (lambda .* con); | |
406 | |
407 if (t2 && t3 && max ([t0; t1; t4]) < tol) | |
408 break; | |
409 endif | |
410 | |
411 ## Compute search direction p by solving QP. | |
412 | |
413 g = -ce; | |
414 d = -ci; | |
415 | |
416 ## Discard inequality constraints that have -Inf bounds since those | |
417 ## will never be active. | |
418 idx = isinf (d) & d < 0; | |
419 d(idx) = []; | |
420 C(idx,:) = []; | |
421 | |
422 [p, obj_qp, INFO, lambda] = qp (x, B, c, F, g, [], [], d, C, | |
423 Inf * ones (size (d))); | |
424 | |
425 info = INFO.info; | |
426 | |
427 ## Check QP solution and attempt to recover if it has failed. | |
428 | |
429 ## Choose mu such that p is a descent direction for the chosen | |
430 ## merit function phi. | |
431 | |
432 [x_new, alpha, obj_new] = linesearch_L1 (x, p, obj_fun, obj_grd, | |
433 ce_fun, ci_fun, lambda, obj); | |
434 | |
435 ## Evaluate objective function, constraints, and gradients at | |
436 ## x_new. | |
437 | |
438 c_new = feval (obj_grd, x_new); | |
439 | |
440 ce_new = feval (ce_fun, x_new); | |
441 F_new = feval (ce_grd, x_new); | |
442 | |
443 ci_new = feval (ci_fun, x_new); | |
444 C_new = feval (ci_grd, x_new); | |
445 | |
446 A_new = [F_new; C_new]; | |
447 | |
448 ## Set | |
449 ## | |
450 ## s = alpha * p | |
451 ## y = grad_x L (x_new, lambda) - grad_x L (x, lambda}) | |
452 | |
453 y = c_new - c; | |
454 | |
455 if (! isempty (A)) | |
456 t = ((A_new - A)'*lambda); | |
457 y -= t; | |
458 endif | |
459 | |
460 delx = x_new - x; | |
461 | |
462 if (norm (delx) < tol * norm (x)) | |
463 info = 101; | |
464 break; | |
465 endif | |
466 | |
467 if (have_hess) | |
468 | |
469 B = feval (obj_hess, x); | |
470 | |
471 else | |
472 | |
473 ## Update B using a quasi-Newton formula. | |
474 | |
475 delxt = delx'; | |
476 | |
477 ## Damped BFGS. Or maybe we would actually want to use the Hessian | |
478 ## of the Lagrangian, computed directly. | |
479 | |
480 d1 = delxt*B*delx; | |
481 | |
482 t1 = 0.2 * d1; | |
483 t2 = delxt*y; | |
484 | |
485 if (t2 < t1) | |
486 theta = 0.8*d1/(d1 - t2); | |
487 else | |
488 theta = 1; | |
489 endif | |
490 | |
491 r = theta*y + (1-theta)*B*delx; | |
492 | |
493 d2 = delxt*r; | |
494 | |
495 if (d1 == 0 || d2 == 0) | |
496 info = 102; | |
497 break; | |
498 endif | |
499 | |
500 B = B - B*delx*delxt*B/d1 + r*r'/d2; | |
501 | |
502 endif | |
503 | |
504 x = x_new; | |
505 | |
506 obj = obj_new; | |
507 | |
508 c = c_new; | |
509 | |
510 ce = ce_new; | |
511 F = F_new; | |
512 | |
513 ci = ci_new; | |
514 C = C_new; | |
515 | |
516 A = A_new; | |
517 | |
6768 | 518 ## report (iter, qp_iter, alpha, __sqp_nfun__, obj); |
5289 | 519 |
520 endwhile | |
521 | |
522 if (iter >= iter_max) | |
523 info = 103; | |
524 endif | |
525 | |
6768 | 526 nf = __sqp_nfun__; |
5289 | 527 |
528 else | |
529 | |
6046 | 530 print_usage (); |
5289 | 531 |
532 endif | |
533 | |
7399 | 534 endfunction |
5289 | 535 |
536 | |
537 function [merit, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu) | |
538 | |
6768 | 539 global __sqp_nfun__; |
5289 | 540 |
541 ce = feval (ce_fun, x); | |
542 ci = feval (ci_fun, x); | |
543 | |
544 idx = ci < 0; | |
545 | |
546 con = [ce; ci(idx)]; | |
547 | |
548 if (isempty (obj)) | |
549 obj = feval (obj_fun, x); | |
6768 | 550 __sqp_nfun__++; |
5289 | 551 endif |
552 | |
553 merit = obj; | |
554 t = norm (con, 1) / mu; | |
555 | |
556 if (! isempty (t)) | |
557 merit += t; | |
558 endif | |
559 | |
7399 | 560 endfunction |
5289 | 561 |
562 | |
563 function [x_new, alpha, obj] = linesearch_L1 (x, p, obj_fun, obj_grd, | |
564 ce_fun, ci_fun, lambda, obj) | |
565 | |
566 ## Choose parameters | |
567 ## | |
568 ## eta in the range (0, 0.5) | |
569 ## tau in the range (0, 1) | |
570 | |
571 eta = 0.25; | |
572 tau = 0.5; | |
573 | |
574 delta_bar = sqrt (eps); | |
575 | |
576 if (isempty (lambda)) | |
577 mu = 1 / delta_bar; | |
578 else | |
579 mu = 1 / (norm (lambda, Inf) + delta_bar); | |
580 endif | |
581 | |
582 alpha = 1; | |
583 | |
584 c = feval (obj_grd, x); | |
585 ce = feval (ce_fun, x); | |
586 | |
587 [phi_x_mu, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu); | |
588 | |
589 D_phi_x_mu = c' * p; | |
590 d = feval (ci_fun, x); | |
591 ## only those elements of d corresponding | |
592 ## to violated constraints should be included. | |
593 idx = d < 0; | |
594 t = - norm ([ce; d(idx)], 1) / mu; | |
595 if (! isempty (t)) | |
596 D_phi_x_mu += t; | |
597 endif | |
598 | |
599 while (1) | |
600 [p1, obj] = phi_L1 ([], obj_fun, ce_fun, ci_fun, x+alpha*p, mu); | |
601 p2 = phi_x_mu+eta*alpha*D_phi_x_mu; | |
602 if (p1 > p2) | |
603 ## Reset alpha = tau_alpha * alpha for some tau_alpha in the | |
604 ## range (0, tau). | |
605 tau_alpha = 0.9 * tau; ## ?? | |
606 alpha = tau_alpha * alpha; | |
607 else | |
608 break; | |
609 endif | |
610 endwhile | |
611 | |
612 ## Set x_new = x + alpha * p; | |
613 | |
614 x_new = x + alpha * p; | |
615 | |
7399 | 616 endfunction |
5289 | 617 |
618 | |
619 function report (iter, qp_iter, alpha, nfun, obj) | |
620 | |
621 if (nargin == 0) | |
622 printf (" Itn ItQP Step Nfun Objective\n"); | |
623 else | |
624 printf ("%5d %4d %8.1g %5d %13.6e\n", iter, qp_iter, alpha, nfun, obj); | |
625 endif | |
626 | |
7399 | 627 endfunction |
5289 | 628 |
629 | |
630 function grd = fdgrd (f, x) | |
631 | |
632 if (! isempty (f)) | |
633 y0 = feval (f, x); | |
634 nx = length (x); | |
635 grd = zeros (nx, 1); | |
636 deltax = sqrt (eps); | |
637 for i = 1:nx | |
638 t = x(i); | |
639 x(i) += deltax; | |
640 grd(i) = (feval (f, x) - y0) / deltax; | |
641 x(i) = t; | |
642 endfor | |
643 else | |
644 grd = zeros (0, 1); | |
645 endif | |
646 | |
7399 | 647 endfunction |
5289 | 648 |
649 | |
650 function jac = fdjac (f, x) | |
6768 | 651 nx = length (x); |
5289 | 652 if (! isempty (f)) |
653 y0 = feval (f, x); | |
654 nf = length (y0); | |
655 nx = length (x); | |
656 jac = zeros (nf, nx); | |
657 deltax = sqrt (eps); | |
658 for i = 1:nx | |
659 t = x(i); | |
660 x(i) += deltax; | |
661 jac(:,i) = (feval (f, x) - y0) / deltax; | |
662 x(i) = t; | |
663 endfor | |
664 else | |
665 jac = zeros (0, nx); | |
666 endif | |
667 | |
7399 | 668 endfunction |
5289 | 669 |
670 | |
671 function grd = fd_obj_grd (x) | |
672 | |
673 global __sqp_obj_fun__; | |
674 | |
675 grd = fdgrd (__sqp_obj_fun__, x); | |
676 | |
7399 | 677 endfunction |
5289 | 678 |
679 | |
680 function res = empty_cf (x) | |
681 | |
682 res = zeros (0, 1); | |
683 | |
7399 | 684 endfunction |
5289 | 685 |
686 | |
687 function res = empty_jac (x) | |
688 | |
689 res = zeros (0, length (x)); | |
690 | |
7399 | 691 endfunction |
5289 | 692 |
693 | |
694 function jac = fd_ce_jac (x) | |
695 | |
696 global __sqp_ce_fun__; | |
697 | |
698 jac = fdjac (__sqp_ce_fun__, x); | |
699 | |
7399 | 700 endfunction |
5289 | 701 |
702 | |
703 function jac = fd_ci_jac (x) | |
704 | |
6768 | 705 global __sqp_cifcn__; |
706 ## __sqp_cifcn__ = constraint function without gradients and lb or ub | |
707 jac = fdjac (__sqp_cifcn__, x); | |
708 | |
7399 | 709 endfunction |
6768 | 710 |
7017 | 711 |
6768 | 712 function res = cf_ub_lb (x) |
5289 | 713 |
6768 | 714 ## combine constraint function with ub and lb |
715 global __sqp_cifcn__ __sqp_lb__ __sqp_ub__ | |
716 | |
717 res = [x-__sqp_lb__; __sqp_ub__-x]; | |
718 | |
719 if (! isempty (__sqp_cifcn__)) | |
720 res = [feval(__sqp_cifcn__,x); x-__sqp_lb__; __sqp_ub__-x]; | |
721 endif | |
5289 | 722 |
7399 | 723 endfunction |
6768 | 724 |
7017 | 725 |
6768 | 726 function res = cigrad_ub_lb (x) |
727 | |
728 global __sqp_cif__ | |
729 | |
730 res = [eye(numel(x)); -eye(numel(x))]; | |
731 | |
732 cigradfcn = @fd_ci_jac; | |
733 | |
734 if (iscell (__sqp_cif__) && length (__sqp_cif__) > 1) | |
735 cigradfcn = __sqp_cif__{2}; | |
736 endif | |
737 | |
738 if (! isempty (cigradfcn)) | |
739 res = [feval(cigradfcn,x); eye(numel(x)); -eye(numel(x))]; | |
740 endif | |
741 | |
7399 | 742 endfunction |
7361 | 743 |
7371 | 744 %!function r = g (x) |
745 %! r = [sumsq(x)-10; | |
746 %! x(2)*x(3)-5*x(4)*x(5); | |
747 %! x(1)^3+x(2)^3+1 ]; | |
7361 | 748 %! |
7371 | 749 %!function obj = phi (x) |
750 %! obj = exp(prod(x)) - 0.5*(x(1)^3+x(2)^3+1)^2; | |
7361 | 751 %! |
752 %!test | |
753 %! x0 = [-1.8; 1.7; 1.9; -0.8; -0.8]; | |
754 %! | |
7381 | 755 %! [x, obj, info, iter, nf, lambda] = sqp (x0, @phi, @g, []); |
7361 | 756 %! |
757 %! x_opt = [-1.717143501952599; | |
758 %! 1.595709610928535; | |
759 %! 1.827245880097156; | |
760 %! -0.763643103133572; | |
761 %! -0.763643068453300]; | |
762 %! | |
7371 | 763 %! obj_opt = 0.0539498477702739; |
7361 | 764 %! |
8047
d54f113aa983
Increase test script tolerances.
Thomas Treichl <Thomas.Treichl@gmx.net>
parents:
7795
diff
changeset
|
765 %! assert (all (abs (x-x_opt) < 5*sqrt (eps)) && abs (obj-obj_opt) < sqrt (eps)); |