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