Mercurial > octave

comparison scripts/optimization/sqp.m @ 6768:40e1255eda0e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-06-29 20:11:58 by jwe]
author jwe
date Fri, 29 Jun 2007 20:12:04 +0000
parents 00116015904d
children 2ae5d4353d0b
comparison
equal deleted inserted replaced
6767:a6c8000f113e 6768:40e1255eda0e
180 ## increase this value). 180 ## increase this value).
181 ## @end table 181 ## @end table
182 ## @seealso{qp} 182 ## @seealso{qp}
183 ## @end deftypefn 183 ## @end deftypefn
184 184
185 function [x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif) 185 function [x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif, lb, ub, maxiter, tolerance)
186 186
187 global nfun; 187 global __sqp_nfun__;
188 global __sqp_obj_fun__; 188 global __sqp_obj_fun__;
189 global __sqp_ce_fun__; 189 global __sqp_ce_fun__;
190 global __sqp_ci_fun__; 190 global __sqp_ci_fun__;
191 191 global __sqp_cif__;
192 if (nargin >= 2 && nargin <= 4) 192 global __sqp_cifcn__;
193
194 if (nargin >= 2 && nargin <= 8 && nargin != 5)
193 195
194 ## Choose an initial NxN symmetric positive definite Hessan 196 ## Choose an initial NxN symmetric positive definite Hessan
195 ## approximation B. 197 ## approximation B.
196 198
197 n = length (x); 199 n = length (x);
243 endif 245 endif
244 __sqp_ce_fun__ = ce_fun; 246 __sqp_ce_fun__ = ce_fun;
245 247
246 ci_fun = @empty_cf; 248 ci_fun = @empty_cf;
247 ci_grd = @empty_jac; 249 ci_grd = @empty_jac;
250
248 if (nargin > 3) 251 if (nargin > 3)
249 ci_grd = @fd_ci_jac; 252 ## constraint function given by user with possibly gradient
250 if (iscell (cif)) 253 __sqp_cif__ = cif;
251 if (length (cif) > 0) 254 ## constraint function given by user without gradient
252 __sqp_ci_fun__ = ci_fun = cif{1}; 255 __sqp_cifcn__ = @empty_cf;
253 if (length (cif) > 1) 256 if (iscell (__sqp_cif__))
254 ci_grd = cif{2}; 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");
255 endif 274 endif
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;
256 else 284 else
257 error ("sqp: invalid equality constraint function"); 285 error ("sqp: invalid lower bound");
258 endif 286 endif
259 elseif (! isempty (cif)) 287
260 ci_fun = cif; 288 global __sqp_ub__;
261 endif 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;
262 endif 304 endif
263 __sqp_ci_fun__ = ci_fun;
264 305
265 iter_max = 100; 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");
312 endif
313 endif
314
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
266 323
267 iter = 0; 324 iter = 0;
268 325
269 obj = feval (obj_fun, x); 326 obj = feval (obj_fun, x);
270 nfun = 1; 327 __sqp_nfun__ = 1;
271 328
272 c = feval (obj_grd, x); 329 c = feval (obj_grd, x);
273 330
274 if (have_hess) 331 if (have_hess)
275 B = feval (obj_hess, x); 332 B = feval (obj_hess, x);
292 qp_iter = 1; 349 qp_iter = 1;
293 alpha = 1; 350 alpha = 1;
294 351
295 ## report (); 352 ## report ();
296 353
297 ## report (iter, qp_iter, alpha, nfun, obj); 354 ## report (iter, qp_iter, alpha, __sqp_nfun__, obj);
298 355
299 info = 0; 356 info = 0;
300 357
301 while (++iter < iter_max) 358 while (++iter < iter_max)
302 359
316 t1 = norm (ce); 373 t1 = norm (ce);
317 t2 = all (ci >= 0); 374 t2 = all (ci >= 0);
318 t3 = all (lambda_i >= 0); 375 t3 = all (lambda_i >= 0);
319 t4 = norm (lambda .* con); 376 t4 = norm (lambda .* con);
320 377
321 tol = sqrt (eps);
322
323 if (t2 && t3 && max ([t0; t1; t4]) < tol) 378 if (t2 && t3 && max ([t0; t1; t4]) < tol)
324 break; 379 break;
325 endif 380 endif
326 381
327 ## Compute search direction p by solving QP. 382 ## Compute search direction p by solving QP.
429 ci = ci_new; 484 ci = ci_new;
430 C = C_new; 485 C = C_new;
431 486
432 A = A_new; 487 A = A_new;
433 488
434 ## report (iter, qp_iter, alpha, nfun, obj); 489 ## report (iter, qp_iter, alpha, __sqp_nfun__, obj);
435 490
436 endwhile 491 endwhile
437 492
438 if (iter >= iter_max) 493 if (iter >= iter_max)
439 info = 103; 494 info = 103;
440 endif 495 endif
441 496
442 nf = nfun; 497 nf = __sqp_nfun__;
443 498
444 else 499 else
445 500
446 print_usage (); 501 print_usage ();
447 502
450 ### endfunction 505 ### endfunction
451 506
452 507
453 function [merit, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu) 508 function [merit, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu)
454 509
455 global nfun; 510 global __sqp_nfun__;
456 511
457 ce = feval (ce_fun, x); 512 ce = feval (ce_fun, x);
458 ci = feval (ci_fun, x); 513 ci = feval (ci_fun, x);
459 514
460 idx = ci < 0; 515 idx = ci < 0;
461 516
462 con = [ce; ci(idx)]; 517 con = [ce; ci(idx)];
463 518
464 if (isempty (obj)) 519 if (isempty (obj))
465 obj = feval (obj_fun, x); 520 obj = feval (obj_fun, x);
466 nfun++; 521 __sqp_nfun__++;
467 endif 522 endif
468 523
469 merit = obj; 524 merit = obj;
470 t = norm (con, 1) / mu; 525 t = norm (con, 1) / mu;
471 526
562 617
563 ### endfunction 618 ### endfunction
564 619
565 620
566 function jac = fdjac (f, x) 621 function jac = fdjac (f, x)
567 622 nx = length (x);
568 if (! isempty (f)) 623 if (! isempty (f))
569 y0 = feval (f, x); 624 y0 = feval (f, x);
570 nf = length (y0); 625 nf = length (y0);
571 nx = length (x); 626 nx = length (x);
572 jac = zeros (nf, nx); 627 jac = zeros (nf, nx);
616 ### endfunction 671 ### endfunction
617 672
618 673
619 function jac = fd_ci_jac (x) 674 function jac = fd_ci_jac (x)
620 675
621 global __sqp_ci_fun__; 676 global __sqp_cifcn__;
622 677 ## __sqp_cifcn__ = constraint function without gradients and lb or ub
623 jac = fdjac (__sqp_ci_fun__, x); 678 jac = fdjac (__sqp_cifcn__, x);
624 679
625 ### endfunction 680 ### endfunction
681
682 function res = cf_ub_lb (x)
683
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
692
693 ### endfunction
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