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