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