Mercurial > hg > octave-nkf
annotate scripts/optimization/sqp.m @ 11371:c767bb1afa03
sqp.m: Fix indexing error in sqp bounds selection
author | Olaf Till <olaf.till@uni-jena.de> |
---|---|
date | Thu, 16 Dec 2010 07:39:04 -0800 |
parents | 2ae0ca4ee36b |
children | 1740012184f9 |
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 -*- | |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
20 ## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{iter}, @var{nf}, @var{lambda}] =} sqp (@var{x0}, @var{phi}) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
21 ## @deftypefnx {Function File} {[@dots{}] =} sqp (@var{x0}, @var{phi}, @var{g}) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
22 ## @deftypefnx {Function File} {[@dots{}] =} sqp (@var{x0}, @var{phi}, @var{g}, @var{h}) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
23 ## @deftypefnx {Function File} {[@dots{}] =} sqp (@var{x0}, @var{phi}, @var{g}, @var{h}, @var{lb}, @var{ub}) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
24 ## @deftypefnx {Function File} {[@dots{}] =} sqp (@var{x0}, @var{phi}, @var{g}, @var{h}, @var{lb}, @var{ub}, @var{maxiter}) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
25 ## @deftypefnx {Function File} {[@dots{}] =} sqp (@var{x0}, @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 |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10678
diff
changeset
|
57 ## |
6741 | 58 ## @end ifnottex |
5289 | 59 ## @noindent |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
60 ## using a sequential quadratic programming method. |
5289 | 61 ## |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
62 ## The first argument is the initial guess for the vector @var{x0}. |
5289 | 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: | |
6741 | 122 ## @tex |
123 ## $$ | |
124 ## \Bigg( {\partial f(x) \over \partial x_1}, | |
125 ## {\partial f(x) \over \partial x_2}, \ldots, | |
126 ## {\partial f(x) \over \partial x_N} \Bigg)^T | |
127 ## $$ | |
128 ## @end tex | |
129 ## @ifnottex | |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10678
diff
changeset
|
130 ## |
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 |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10678
diff
changeset
|
138 ## |
6741 | 139 ## @end ifnottex |
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. |
11290
646063a22a35
sqp.m: Use correct stopping tolerance in documentation.
Rik <octave@nomad.inbox5.com>
parents:
11120
diff
changeset
|
152 ## The default value is @code{sqrt(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: |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10678
diff
changeset
|
155 ## |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
156 ## @table @asis |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
157 ## @item 101 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
158 ## The algorithm terminated normally. |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
159 ## 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
|
160 ## @tex |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
161 ## $\Delta x,$ |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
162 ## @end tex |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
163 ## @ifnottex |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
164 ## delta @var{x}, |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
165 ## @end ifnottex |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
166 ## is less than @code{tol * norm (x)}. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10678
diff
changeset
|
167 ## |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
168 ## @item 102 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
169 ## The BFGS update failed. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10678
diff
changeset
|
170 ## |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
171 ## @item 103 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
172 ## 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
|
173 ## @end table |
8167
17352ccd860e
describe additional arguments in sqp() documentation string
Ivan Sutoris <ivan.sutoris@gmail.com>
parents:
8047
diff
changeset
|
174 ## |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
175 ## An example of calling @code{sqp}: |
5289 | 176 ## |
177 ## @example | |
7031 | 178 ## function r = g (x) |
179 ## r = [ sumsq(x)-10; | |
180 ## x(2)*x(3)-5*x(4)*x(5); | |
181 ## x(1)^3+x(2)^3+1 ]; | |
182 ## endfunction | |
183 ## | |
184 ## function obj = phi (x) | |
185 ## obj = exp(prod(x)) - 0.5*(x(1)^3+x(2)^3+1)^2; | |
186 ## endfunction | |
5289 | 187 ## |
7031 | 188 ## x0 = [-1.8; 1.7; 1.9; -0.8; -0.8]; |
189 ## | |
190 ## [x, obj, info, iter, nf, lambda] = sqp (x0, @@phi, @@g, []) | |
5289 | 191 ## |
7031 | 192 ## x = |
193 ## | |
194 ## -1.71714 | |
195 ## 1.59571 | |
196 ## 1.82725 | |
197 ## -0.76364 | |
198 ## -0.76364 | |
5289 | 199 ## |
7031 | 200 ## obj = 0.053950 |
201 ## info = 101 | |
202 ## iter = 8 | |
203 ## nf = 10 | |
204 ## lambda = | |
205 ## | |
206 ## -0.0401627 | |
207 ## 0.0379578 | |
208 ## -0.0052227 | |
5289 | 209 ## @end example |
210 ## | |
5642 | 211 ## @seealso{qp} |
5289 | 212 ## @end deftypefn |
213 | |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
214 function [x, obj, info, iter, nf, lambda] = sqp (x0, objf, cef, cif, lb, ub, maxiter, tolerance) |
5289 | 215 |
6768 | 216 global __sqp_nfun__; |
5289 | 217 global __sqp_obj_fun__; |
218 global __sqp_ce_fun__; | |
219 global __sqp_ci_fun__; | |
6768 | 220 global __sqp_cif__; |
221 global __sqp_cifcn__; | |
5289 | 222 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
223 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
|
224 print_usage (); |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
225 endif |
5289 | 226 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
227 if (!isvector (x0)) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
228 error ("sqp: X0 must be a vector"); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
229 endif |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
230 if (rows (x0) == 1) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
231 x0 = x0'; |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
232 endif |
5289 | 233 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
234 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
|
235 have_hess = 0; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
236 if (iscell (objf)) |
11120
a44f979a35ce
style fixes for some .m files
John W. Eaton <jwe@octave.org>
parents:
10821
diff
changeset
|
237 switch (numel (objf)) |
a44f979a35ce
style fixes for some .m files
John W. Eaton <jwe@octave.org>
parents:
10821
diff
changeset
|
238 case 1 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
239 obj_fun = objf{1}; |
11120
a44f979a35ce
style fixes for some .m files
John W. Eaton <jwe@octave.org>
parents:
10821
diff
changeset
|
240 case 2 |
10678
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}; |
11120
a44f979a35ce
style fixes for some .m files
John W. Eaton <jwe@octave.org>
parents:
10821
diff
changeset
|
243 case 3 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
244 obj_fun = objf{1}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
245 obj_grd = objf{2}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
246 obj_hess = objf{3}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
247 have_hess = 1; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
248 otherwise |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
249 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
|
250 endswitch |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
251 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
252 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
|
253 endif |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
254 __sqp_obj_fun__ = obj_fun; |
5289 | 255 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
256 ce_fun = @empty_cf; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
257 ce_grd = @empty_jac; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
258 if (nargin > 2) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
259 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
|
260 if (iscell (cef)) |
11120
a44f979a35ce
style fixes for some .m files
John W. Eaton <jwe@octave.org>
parents:
10821
diff
changeset
|
261 switch (numel (cef)) |
a44f979a35ce
style fixes for some .m files
John W. Eaton <jwe@octave.org>
parents:
10821
diff
changeset
|
262 case 1 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
263 ce_fun = cef{1}; |
11120
a44f979a35ce
style fixes for some .m files
John W. Eaton <jwe@octave.org>
parents:
10821
diff
changeset
|
264 case 2 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
265 ce_fun = cef{1}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
266 ce_grd = cef{2}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
267 otherwise |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
268 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
|
269 endswitch |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
270 elseif (! isempty (cef)) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
271 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
|
272 endif |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
273 endif |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
274 __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
|
275 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
276 ci_fun = @empty_cf; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
277 ci_grd = @empty_jac; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
278 if (nargin > 3) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
279 ## 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
|
280 __sqp_cif__ = cif; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
281 ## 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
|
282 __sqp_cifcn__ = @empty_cf; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
283 if (iscell (cif)) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
284 if (length (cif) > 0) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
285 __sqp_cifcn__ = cif{1}; |
5289 | 286 endif |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
287 elseif (! isempty (cif)) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
288 __sqp_cifcn__ = cif; |
5289 | 289 endif |
290 | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
291 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
|
292 ## 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
|
293 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
|
294 if (iscell (cif)) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
295 switch length (cif) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
296 case {1} |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
297 ci_fun = cif{1}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
298 case {2} |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
299 ci_fun = cif{1}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
300 ci_grd = cif{2}; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
301 otherwise |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
302 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
|
303 endswitch |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
304 elseif (! isempty (cif)) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
305 ci_fun = cif; # No cell array, only constraint inequality function set |
5289 | 306 endif |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
307 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
308 ## 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
|
309 global __sqp_lb__; |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
310 lb_idx = ub_idx = true (size (x0)); |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
311 ub_grad = - (lb_grad = eye (rows (x0))); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
312 if (isvector (lb)) |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
313 __sqp_lb__ = tmp_lb = lb(:); |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
314 lb_idx(:) = tmp_idx = (lb != -Inf); |
11371
c767bb1afa03
sqp.m: Fix indexing error in sqp bounds selection
Olaf Till <olaf.till@uni-jena.de>
parents:
11348
diff
changeset
|
315 __sqp_lb__ = __sqp_lb__(tmp_idx, 1); |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
316 lb_grad = lb_grad(lb_idx, :); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
317 elseif (isempty (lb)) |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
318 if (isa (x0, "single")) |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
319 __sqp_lb__ = tmp_lb = -realmax ("single"); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
320 else |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
321 __sqp_lb__ = tmp_lb = -realmax; |
10549 | 322 endif |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
323 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
324 error ("sqp: invalid lower bound"); |
6768 | 325 endif |
326 | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
327 global __sqp_ub__; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
328 if (isvector (ub)) |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
329 __sqp_ub__ = tmp_ub = ub(:); |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
330 ub_idx(:) = tmp_idx = (ub != Inf); |
11371
c767bb1afa03
sqp.m: Fix indexing error in sqp bounds selection
Olaf Till <olaf.till@uni-jena.de>
parents:
11348
diff
changeset
|
331 __sqp_ub__ = __sqp_ub__(tmp_idx, 1); |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
332 ub_grad = ub_grad(ub_idx, :); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
333 elseif (isempty (ub)) |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
334 if (isa (x0, "single")) |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
335 __sqp_ub__ = tmp_ub = realmax ("single"); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
336 else |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
337 __sqp_ub__ = tmp_ub = realmax; |
10549 | 338 endif |
6768 | 339 else |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
340 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
|
341 endif |
6768 | 342 |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
343 if (any (tmp_lb > tmp_ub)) |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
344 error ("sqp: upper bound smaller than lower bound"); |
6768 | 345 endif |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
346 bounds_grad = [lb_grad; ub_grad]; |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
347 ci_fun = @ (x) cf_ub_lb (x, lb_idx, ub_idx); |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
348 ci_grd = @ (x) cigrad_ub_lb (x, bounds_grad); |
6768 | 349 endif |
5289 | 350 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
351 __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
|
352 endif # if (nargin > 3) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
353 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
354 iter_max = 100; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
355 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
|
356 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
|
357 iter_max = maxiter; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
358 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
359 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
|
360 endif |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
361 endif |
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 tol = sqrt (eps); |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
364 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
|
365 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
|
366 tol = tolerance; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
367 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
368 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
|
369 endif |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
370 endif |
5289 | 371 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
372 ## Initialize variables for search loop |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
373 ## Seed x with initial guess and evaluate objective function, constraints, |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
374 ## and gradients at initial value x0. |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
375 ## |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
376 ## obj_fun -- objective function |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
377 ## obj_grad -- objective gradient |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
378 ## 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
|
379 ## 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
|
380 ## A == [grad_{x_1} cx_fun, grad_{x_2} cx_fun, ..., grad_{x_n} cx_fun]^T |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
381 x = x0; |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
382 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
383 obj = feval (obj_fun, x0); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
384 __sqp_nfun__ = 1; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
385 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
386 c = feval (obj_grd, x0); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
387 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
388 ## Choose an initial NxN symmetric positive definite Hessian approximation B. |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
389 n = length (x0); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
390 if (have_hess) |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
391 B = feval (obj_hess, x0); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
392 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
393 B = eye (n, n); |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
394 endif |
5289 | 395 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
396 ce = feval (ce_fun, x0); |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
397 F = feval (ce_grd, x0); |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
398 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
399 ci = feval (ci_fun, x0); |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
400 C = feval (ci_grd, x0); |
10678
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 A = [F; C]; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
403 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
404 ## 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
|
405 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
|
406 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
407 qp_iter = 1; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
408 alpha = 1; |
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 info = 0; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
411 iter = 0; |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
412 # report (); # Called with no arguments to initialize reporting |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
413 # report (iter, qp_iter, alpha, __sqp_nfun__, obj); |
5289 | 414 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
415 while (++iter < iter_max) |
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 ## 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
|
418 ## order necessary conditions. |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
419 nr_f = rows (F); |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
420 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
421 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
|
422 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
|
423 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
424 con = [ce; ci]; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
425 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
426 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
|
427 t1 = norm (ce); |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
428 t2 = all (ci >= 0); |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
429 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
|
430 t4 = norm (lambda .* con); |
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 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
|
433 info = 101; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
434 break; |
6382 | 435 endif |
436 | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
437 ## 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
|
438 g = -ce; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
439 d = -ci; |
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. |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
450 [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
|
451 ce_fun, ci_fun, lambda, obj); |
5289 | 452 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
453 ## Evaluate objective function, constraints, and gradients at x_new. |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
454 c_new = feval (obj_grd, x_new); |
5289 | 455 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
456 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
|
457 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
|
458 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
459 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
|
460 C_new = feval (ci_grd, x_new); |
5289 | 461 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
462 A_new = [F_new; C_new]; |
5289 | 463 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
464 ## Set |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
465 ## |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
466 ## s = alpha * p |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
467 ## y = grad_x L (x_new, lambda) - grad_x L (x, lambda}) |
6527 | 468 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
469 y = c_new - c; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
470 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
471 if (! isempty (A)) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
472 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
|
473 y -= t; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
474 endif |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
475 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
476 delx = x_new - x; |
5289 | 477 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
478 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
|
479 info = 101; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
480 break; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
481 endif |
5289 | 482 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
483 if (have_hess) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
484 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
485 B = feval (obj_hess, x); |
5289 | 486 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
487 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
488 ## Update B using a quasi-Newton formula. |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
489 delxt = delx'; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
490 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
491 ## Damped BFGS. Or maybe we would actually want to use the Hessian |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
492 ## of the Lagrangian, computed directly? |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
493 d1 = delxt*B*delx; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
494 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
495 t1 = 0.2 * d1; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
496 t2 = delxt*y; |
5289 | 497 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
498 if (t2 < t1) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
499 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
|
500 else |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
501 theta = 1; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
502 endif |
5289 | 503 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
504 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
|
505 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
506 d2 = delxt*r; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
507 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
508 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
|
509 info = 102; |
10549 | 510 break; |
5289 | 511 endif |
512 | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
513 B = B - B*delx*delxt*B/d1 + r*r'/d2; |
5289 | 514 |
515 endif | |
516 | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
517 x = x_new; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
518 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
519 obj = obj_new; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
520 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
521 c = c_new; |
5289 | 522 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
523 ce = ce_new; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
524 F = F_new; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
525 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
526 ci = ci_new; |
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 A = A_new; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
530 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
531 # 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
|
532 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
533 endwhile |
5289 | 534 |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
535 if (iter >= iter_max) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
536 info = 103; |
5289 | 537 endif |
538 | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
539 nf = __sqp_nfun__; |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
540 |
7399 | 541 endfunction |
5289 | 542 |
543 | |
544 function [merit, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu) | |
545 | |
6768 | 546 global __sqp_nfun__; |
5289 | 547 |
548 ce = feval (ce_fun, x); | |
549 ci = feval (ci_fun, x); | |
550 | |
551 idx = ci < 0; | |
552 | |
553 con = [ce; ci(idx)]; | |
554 | |
555 if (isempty (obj)) | |
556 obj = feval (obj_fun, x); | |
6768 | 557 __sqp_nfun__++; |
5289 | 558 endif |
559 | |
560 merit = obj; | |
561 t = norm (con, 1) / mu; | |
562 | |
563 if (! isempty (t)) | |
564 merit += t; | |
565 endif | |
566 | |
7399 | 567 endfunction |
5289 | 568 |
569 | |
570 function [x_new, alpha, obj] = linesearch_L1 (x, p, obj_fun, obj_grd, | |
10549 | 571 ce_fun, ci_fun, lambda, obj) |
5289 | 572 |
573 ## Choose parameters | |
574 ## | |
575 ## eta in the range (0, 0.5) | |
576 ## tau in the range (0, 1) | |
577 | |
578 eta = 0.25; | |
579 tau = 0.5; | |
580 | |
581 delta_bar = sqrt (eps); | |
582 | |
583 if (isempty (lambda)) | |
584 mu = 1 / delta_bar; | |
585 else | |
586 mu = 1 / (norm (lambda, Inf) + delta_bar); | |
587 endif | |
588 | |
589 alpha = 1; | |
590 | |
591 c = feval (obj_grd, x); | |
592 ce = feval (ce_fun, x); | |
593 | |
594 [phi_x_mu, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu); | |
595 | |
596 D_phi_x_mu = c' * p; | |
597 d = feval (ci_fun, x); | |
598 ## only those elements of d corresponding | |
599 ## to violated constraints should be included. | |
600 idx = d < 0; | |
601 t = - norm ([ce; d(idx)], 1) / mu; | |
602 if (! isempty (t)) | |
603 D_phi_x_mu += t; | |
604 endif | |
605 | |
606 while (1) | |
607 [p1, obj] = phi_L1 ([], obj_fun, ce_fun, ci_fun, x+alpha*p, mu); | |
608 p2 = phi_x_mu+eta*alpha*D_phi_x_mu; | |
609 if (p1 > p2) | |
610 ## Reset alpha = tau_alpha * alpha for some tau_alpha in the | |
611 ## range (0, tau). | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
612 tau_alpha = 0.9 * tau; # ?? |
5289 | 613 alpha = tau_alpha * alpha; |
614 else | |
615 break; | |
616 endif | |
617 endwhile | |
618 | |
619 x_new = x + alpha * p; | |
620 | |
7399 | 621 endfunction |
5289 | 622 |
623 | |
624 function grd = fdgrd (f, x) | |
625 | |
626 if (! isempty (f)) | |
627 y0 = feval (f, x); | |
628 nx = length (x); | |
629 grd = zeros (nx, 1); | |
630 deltax = sqrt (eps); | |
631 for i = 1:nx | |
632 t = x(i); | |
633 x(i) += deltax; | |
634 grd(i) = (feval (f, x) - y0) / deltax; | |
635 x(i) = t; | |
636 endfor | |
637 else | |
638 grd = zeros (0, 1); | |
639 endif | |
640 | |
7399 | 641 endfunction |
5289 | 642 |
643 | |
644 function jac = fdjac (f, x) | |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
645 |
6768 | 646 nx = length (x); |
5289 | 647 if (! isempty (f)) |
648 y0 = feval (f, x); | |
649 nf = length (y0); | |
650 nx = length (x); | |
651 jac = zeros (nf, nx); | |
652 deltax = sqrt (eps); | |
653 for i = 1:nx | |
654 t = x(i); | |
655 x(i) += deltax; | |
656 jac(:,i) = (feval (f, x) - y0) / deltax; | |
657 x(i) = t; | |
658 endfor | |
659 else | |
660 jac = zeros (0, nx); | |
661 endif | |
662 | |
7399 | 663 endfunction |
5289 | 664 |
665 | |
666 function grd = fd_obj_grd (x) | |
667 | |
668 global __sqp_obj_fun__; | |
669 | |
670 grd = fdgrd (__sqp_obj_fun__, x); | |
671 | |
7399 | 672 endfunction |
5289 | 673 |
674 | |
675 function res = empty_cf (x) | |
676 | |
677 res = zeros (0, 1); | |
678 | |
7399 | 679 endfunction |
5289 | 680 |
681 | |
682 function res = empty_jac (x) | |
683 | |
684 res = zeros (0, length (x)); | |
685 | |
7399 | 686 endfunction |
5289 | 687 |
688 | |
689 function jac = fd_ce_jac (x) | |
690 | |
691 global __sqp_ce_fun__; | |
692 | |
693 jac = fdjac (__sqp_ce_fun__, x); | |
694 | |
7399 | 695 endfunction |
5289 | 696 |
697 | |
698 function jac = fd_ci_jac (x) | |
699 | |
6768 | 700 global __sqp_cifcn__; |
701 ## __sqp_cifcn__ = constraint function without gradients and lb or ub | |
702 jac = fdjac (__sqp_cifcn__, x); | |
703 | |
7399 | 704 endfunction |
6768 | 705 |
7017 | 706 |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
707 function res = cf_ub_lb (x, lbidx, ubidx) |
5289 | 708 |
6768 | 709 ## combine constraint function with ub and lb |
710 global __sqp_cifcn__ __sqp_lb__ __sqp_ub__ | |
711 | |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
712 if (isempty (__sqp_cifcn__)) |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
713 res = [x(lbidx,1)-__sqp_lb__; __sqp_ub__-x(ubidx,1)]; |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
714 else |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
715 res = [feval(__sqp_cifcn__,x); \ |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
716 x(lbidx,1)-__sqp_lb__; __sqp_ub__-x(ubidx,1)]; |
6768 | 717 endif |
5289 | 718 |
7399 | 719 endfunction |
6768 | 720 |
7017 | 721 |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
722 function res = cigrad_ub_lb (x, bgrad) |
6768 | 723 |
724 global __sqp_cif__ | |
725 | |
726 cigradfcn = @fd_ci_jac; | |
727 | |
728 if (iscell (__sqp_cif__) && length (__sqp_cif__) > 1) | |
729 cigradfcn = __sqp_cif__{2}; | |
730 endif | |
10549 | 731 |
11347
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
732 if (isempty (cigradfcn)) |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
733 res = bgrad; |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
734 else |
2726132f77f6
sqp.m: Remove never violated Inf bounds from computation (bug #31742)
Rik <octave@nomad.inbox5.com>
parents:
11290
diff
changeset
|
735 res = [feval(cigradfcn,x); bgrad]; |
6768 | 736 endif |
737 | |
7399 | 738 endfunction |
7361 | 739 |
11348
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
740 # Utility function used to debug sqp |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
741 function report (iter, qp_iter, alpha, nfun, obj) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
742 |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
743 if (nargin == 0) |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
744 printf (" Itn ItQP Step Nfun Objective\n"); |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
745 else |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
746 printf ("%5d %4d %8.1g %5d %13.6e\n", iter, qp_iter, alpha, nfun, obj); |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
747 endif |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
748 |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
749 endfunction |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
750 |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
751 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
752 %% Test Code |
2ae0ca4ee36b
sqp.m: Change docstring to refer to x0 as the initial seed vector
Rik <octave@nomad.inbox5.com>
parents:
11347
diff
changeset
|
753 |
7371 | 754 %!function r = g (x) |
755 %! r = [sumsq(x)-10; | |
756 %! x(2)*x(3)-5*x(4)*x(5); | |
757 %! x(1)^3+x(2)^3+1 ]; | |
7361 | 758 %! |
7371 | 759 %!function obj = phi (x) |
760 %! obj = exp(prod(x)) - 0.5*(x(1)^3+x(2)^3+1)^2; | |
7361 | 761 %! |
762 %!test | |
10678
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
763 %! |
7361 | 764 %! x0 = [-1.8; 1.7; 1.9; -0.8; -0.8]; |
765 %! | |
7381 | 766 %! [x, obj, info, iter, nf, lambda] = sqp (x0, @phi, @g, []); |
7361 | 767 %! |
768 %! x_opt = [-1.717143501952599; | |
769 %! 1.595709610928535; | |
770 %! 1.827245880097156; | |
771 %! -0.763643103133572; | |
772 %! -0.763643068453300]; | |
773 %! | |
7371 | 774 %! obj_opt = 0.0539498477702739; |
7361 | 775 %! |
8047
d54f113aa983
Increase test script tolerances.
Thomas Treichl <Thomas.Treichl@gmx.net>
parents:
7795
diff
changeset
|
776 %! 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
|
777 |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
778 %% Test input validation |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
779 %!error sqp () |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
780 %!error sqp (1) |
35338deff753
Guarantee equivalent results if sqp called with or wihout bounds
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
781 %!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
|
782 %!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
|
783 %!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
|
784 %!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
|
785 %!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
|
786 %!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
|
787 %!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
|
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),1,-1) |
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),[],[],ones(2,2)) |
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),[],[],-1) |
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.5) |
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),[],[],[],ones(2,2)) |
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),[],[],[],-1) |