Mercurial > octave-nkf

comparison scripts/optimization/sqp.m @ 10549:95c3e38098bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Untabify .m scripts
author Rik <code@nomad.inbox5.com>
date Fri, 23 Apr 2010 11:28:50 -0700
parents 83de7b060e91
children 35338deff753
comparison
equal deleted inserted replaced
10548:479536c5bb10 10549:95c3e38098bf
222 222
223 obj_grd = @fd_obj_grd; 223 obj_grd = @fd_obj_grd;
224 have_hess = 0; 224 have_hess = 0;
225 if (iscell (objf)) 225 if (iscell (objf))
226 if (length (objf) > 0) 226 if (length (objf) > 0)
227 __sqp_obj_fun__ = obj_fun = objf{1}; 227 __sqp_obj_fun__ = obj_fun = objf{1};
228 if (length (objf) > 1) 228 if (length (objf) > 1)
229 obj_grd = objf{2}; 229 obj_grd = objf{2};
230 if (length (objf) > 2) 230 if (length (objf) > 2)
231 obj_hess = objf{3}; 231 obj_hess = objf{3};
232 have_hess = 1; 232 have_hess = 1;
233 endif 233 endif
234 endif 234 endif
235 else 235 else
236 error ("sqp: invalid objective function"); 236 error ("sqp: invalid objective function");
237 endif 237 endif
238 else 238 else
239 __sqp_obj_fun__ = obj_fun = objf; 239 __sqp_obj_fun__ = obj_fun = objf;
240 endif 240 endif
241 241
242 ce_fun = @empty_cf; 242 ce_fun = @empty_cf;
243 ce_grd = @empty_jac; 243 ce_grd = @empty_jac;
244 if (nargin > 2) 244 if (nargin > 2)
245 ce_grd = @fd_ce_jac; 245 ce_grd = @fd_ce_jac;
246 if (iscell (cef)) 246 if (iscell (cef))
247 if (length (cef) > 0) 247 if (length (cef) > 0)
248 __sqp_ce_fun__ = ce_fun = cef{1}; 248 __sqp_ce_fun__ = ce_fun = cef{1};
249 if (length (cef) > 1) 249 if (length (cef) > 1)
250 ce_grd = cef{2}; 250 ce_grd = cef{2};
251 endif 251 endif
252 else 252 else
253 error ("sqp: invalid equality constraint function"); 253 error ("sqp: invalid equality constraint function");
254 endif 254 endif
255 elseif (! isempty (cef)) 255 elseif (! isempty (cef))
256 ce_fun = cef; 256 ce_fun = cef;
257 endif 257 endif
258 endif 258 endif
259 __sqp_ce_fun__ = ce_fun; 259 __sqp_ce_fun__ = ce_fun;
260 260
261 ci_fun = @empty_cf; 261 ci_fun = @empty_cf;
262 ci_grd = @empty_jac; 262 ci_grd = @empty_jac;
263 263
264 if (nargin > 3) 264 if (nargin > 3)
265 ## constraint function given by user with possibly gradient 265 ## constraint function given by user with possibly gradient
266 __sqp_cif__ = cif; 266 __sqp_cif__ = cif;
267 ## constraint function given by user without gradient 267 ## constraint function given by user without gradient
268 __sqp_cifcn__ = @empty_cf; 268 __sqp_cifcn__ = @empty_cf;
269 if (iscell (__sqp_cif__)) 269 if (iscell (__sqp_cif__))
270 if (length (__sqp_cif__) > 0) 270 if (length (__sqp_cif__) > 0)
271 __sqp_cifcn__ = __sqp_cif__{1}; 271 __sqp_cifcn__ = __sqp_cif__{1};
272 endif 272 endif
273 elseif (! isempty (__sqp_cif__)) 273 elseif (! isempty (__sqp_cif__))
274 __sqp_cifcn__ = __sqp_cif__; 274 __sqp_cifcn__ = __sqp_cif__;
275 endif 275 endif
276 276
277 if (nargin < 5) 277 if (nargin < 5)
278 ci_grd = @fd_ci_jac; 278 ci_grd = @fd_ci_jac;
279 if (iscell (cif)) 279 if (iscell (cif))
280 if (length (cif) > 0) 280 if (length (cif) > 0)
281 __sqp_ci_fun__ = ci_fun = cif{1}; 281 __sqp_ci_fun__ = ci_fun = cif{1};
282 if (length (cif) > 1) 282 if (length (cif) > 1)
283 ci_grd = cif{2}; 283 ci_grd = cif{2};
284 endif 284 endif
285 else 285 else
286 error ("sqp: invalid equality constraint function"); 286 error ("sqp: invalid equality constraint function");
287 endif 287 endif
288 elseif (! isempty (cif)) 288 elseif (! isempty (cif))
289 ci_fun = cif; 289 ci_fun = cif;
290 endif 290 endif
291 else 291 else
292 global __sqp_lb__; 292 global __sqp_lb__;
293 if (isvector (lb)) 293 if (isvector (lb))
294 __sqp_lb__ = lb; 294 __sqp_lb__ = lb;
295 elseif (isempty (lb)) 295 elseif (isempty (lb))
296 if (isa (x, "single")) 296 if (isa (x, "single"))
297 __sqp_lb__ = -realmax ("single"); 297 __sqp_lb__ = -realmax ("single");
298 else 298 else
299 __sqp_lb__ = -realmax; 299 __sqp_lb__ = -realmax;
300 endif 300 endif
301 else 301 else
302 error ("sqp: invalid lower bound"); 302 error ("sqp: invalid lower bound");
303 endif 303 endif
304 304
305 global __sqp_ub__; 305 global __sqp_ub__;
306 if (isvector (ub)) 306 if (isvector (ub))
307 __sqp_ub__ = ub; 307 __sqp_ub__ = ub;
308 elseif (isempty (lb)) 308 elseif (isempty (lb))
309 if (isa (x, "single")) 309 if (isa (x, "single"))
310 __sqp_ub__ = realmax ("single"); 310 __sqp_ub__ = realmax ("single");
311 else 311 else
312 __sqp_ub__ = realmax; 312 __sqp_ub__ = realmax;
313 endif 313 endif
314 else 314 else
315 error ("sqp: invalid upper bound"); 315 error ("sqp: invalid upper bound");
316 endif 316 endif
317 317
318 if (lb > ub) 318 if (lb > ub)
319 error ("sqp: upper bound smaller than lower bound"); 319 error ("sqp: upper bound smaller than lower bound");
320 endif 320 endif
321 __sqp_ci_fun__ = ci_fun = @cf_ub_lb; 321 __sqp_ci_fun__ = ci_fun = @cf_ub_lb;
322 ci_grd = @cigrad_ub_lb; 322 ci_grd = @cigrad_ub_lb;
323 endif 323 endif
324 __sqp_ci_fun__ = ci_fun; 324 __sqp_ci_fun__ = ci_fun;
325 endif 325 endif
326 326
327 iter_max = 100; 327 iter_max = 100;
328 if (nargin > 6 && ! isempty (maxiter)) 328 if (nargin > 6 && ! isempty (maxiter))
329 if (isscalar (maxiter) && maxiter > 0 && round (maxiter) == maxiter) 329 if (isscalar (maxiter) && maxiter > 0 && round (maxiter) == maxiter)
330 iter_max = maxiter; 330 iter_max = maxiter;
331 else 331 else
332 error ("sqp: invalid number of maximum iterations"); 332 error ("sqp: invalid number of maximum iterations");
333 endif 333 endif
334 endif 334 endif
335 335
336 tol = sqrt (eps); 336 tol = sqrt (eps);
337 if (nargin > 7 && ! isempty (tolerance)) 337 if (nargin > 7 && ! isempty (tolerance))
338 if (isscalar (tolerance) && tolerance > 0) 338 if (isscalar (tolerance) && tolerance > 0)
339 tol = tolerance; 339 tol = tolerance;
340 else 340 else
341 error ("sqp: invalid value for tolerance"); 341 error ("sqp: invalid value for tolerance");
342 endif 342 endif
343 endif 343 endif
344 344
345 iter = 0; 345 iter = 0;
346 346
396 t3 = all (lambda_i >= 0); 396 t3 = all (lambda_i >= 0);
397 t4 = norm (lambda .* con); 397 t4 = norm (lambda .* con);
398 398
399 if (t2 && t3 && max ([t0; t1; t4]) < tol) 399 if (t2 && t3 && max ([t0; t1; t4]) < tol)
400 info = 101; 400 info = 101;
401 break; 401 break;
402 endif 402 endif
403 403
404 ## Compute search direction p by solving QP. 404 ## Compute search direction p by solving QP.
405 405
406 g = -ce; 406 g = -ce;
411 idx = isinf (d) & d < 0; 411 idx = isinf (d) & d < 0;
412 d(idx) = []; 412 d(idx) = [];
413 C(idx,:) = []; 413 C(idx,:) = [];
414 414
415 [p, obj_qp, INFO, lambda] = qp (x, B, c, F, g, [], [], d, C, 415 [p, obj_qp, INFO, lambda] = qp (x, B, c, F, g, [], [], d, C,
416 Inf (size (d))); 416 Inf (size (d)));
417 417
418 info = INFO.info; 418 info = INFO.info;
419 419
420 ## Check QP solution and attempt to recover if it has failed. 420 ## Check QP solution and attempt to recover if it has failed.
421 421
422 ## Choose mu such that p is a descent direction for the chosen 422 ## Choose mu such that p is a descent direction for the chosen
423 ## merit function phi. 423 ## merit function phi.
424 424
425 [x_new, alpha, obj_new] = linesearch_L1 (x, p, obj_fun, obj_grd, 425 [x_new, alpha, obj_new] = linesearch_L1 (x, p, obj_fun, obj_grd,
426 ce_fun, ci_fun, lambda, obj); 426 ce_fun, ci_fun, lambda, obj);
427 427
428 ## Evaluate objective function, constraints, and gradients at 428 ## Evaluate objective function, constraints, and gradients at
429 ## x_new. 429 ## x_new.
430 430
431 c_new = feval (obj_grd, x_new); 431 c_new = feval (obj_grd, x_new);
444 ## y = grad_x L (x_new, lambda) - grad_x L (x, lambda}) 444 ## y = grad_x L (x_new, lambda) - grad_x L (x, lambda})
445 445
446 y = c_new - c; 446 y = c_new - c;
447 447
448 if (! isempty (A)) 448 if (! isempty (A))
449 t = ((A_new - A)'*lambda); 449 t = ((A_new - A)'*lambda);
450 y -= t; 450 y -= t;
451 endif 451 endif
452 452
453 delx = x_new - x; 453 delx = x_new - x;
454 454
455 if (norm (delx) < tol * norm (x)) 455 if (norm (delx) < tol * norm (x))
456 info = 101; 456 info = 101;
457 break; 457 break;
458 endif 458 endif
459 459
460 if (have_hess) 460 if (have_hess)
461 461
462 B = feval (obj_hess, x); 462 B = feval (obj_hess, x);
463 463
464 else 464 else
465 465
466 ## Update B using a quasi-Newton formula. 466 ## Update B using a quasi-Newton formula.
467 467
468 delxt = delx'; 468 delxt = delx';
469 469
470 ## Damped BFGS. Or maybe we would actually want to use the Hessian 470 ## Damped BFGS. Or maybe we would actually want to use the Hessian
471 ## of the Lagrangian, computed directly. 471 ## of the Lagrangian, computed directly.
472 472
473 d1 = delxt*B*delx; 473 d1 = delxt*B*delx;
474 474
475 t1 = 0.2 * d1; 475 t1 = 0.2 * d1;
476 t2 = delxt*y; 476 t2 = delxt*y;
477 477
478 if (t2 < t1) 478 if (t2 < t1)
479 theta = 0.8*d1/(d1 - t2); 479 theta = 0.8*d1/(d1 - t2);
480 else 480 else
481 theta = 1; 481 theta = 1;
482 endif 482 endif
483 483
484 r = theta*y + (1-theta)*B*delx; 484 r = theta*y + (1-theta)*B*delx;
485 485
486 d2 = delxt*r; 486 d2 = delxt*r;
487 487
488 if (d1 == 0 || d2 == 0) 488 if (d1 == 0 || d2 == 0)
489 info = 102; 489 info = 102;
490 break; 490 break;
491 endif 491 endif
492 492
493 B = B - B*delx*delxt*B/d1 + r*r'/d2; 493 B = B - B*delx*delxt*B/d1 + r*r'/d2;
494 494
495 endif 495 endif
496 496
497 x = x_new; 497 x = x_new;
498 498
552 552
553 endfunction 553 endfunction
554 554
555 555
556 function [x_new, alpha, obj] = linesearch_L1 (x, p, obj_fun, obj_grd, 556 function [x_new, alpha, obj] = linesearch_L1 (x, p, obj_fun, obj_grd,
557 ce_fun, ci_fun, lambda, obj) 557 ce_fun, ci_fun, lambda, obj)
558 558
559 ## Choose parameters 559 ## Choose parameters
560 ## 560 ##
561 ## eta in the range (0, 0.5) 561 ## eta in the range (0, 0.5)
562 ## tau in the range (0, 1) 562 ## tau in the range (0, 1)
725 cigradfcn = @fd_ci_jac; 725 cigradfcn = @fd_ci_jac;
726 726
727 if (iscell (__sqp_cif__) && length (__sqp_cif__) > 1) 727 if (iscell (__sqp_cif__) && length (__sqp_cif__) > 1)
728 cigradfcn = __sqp_cif__{2}; 728 cigradfcn = __sqp_cif__{2};
729 endif 729 endif
730 730
731 if (! isempty (cigradfcn)) 731 if (! isempty (cigradfcn))
732 res = [feval(cigradfcn,x); eye(numel(x)); -eye(numel(x))]; 732 res = [feval(cigradfcn,x); eye(numel(x)); -eye(numel(x))];
733 endif 733 endif
734 734
735 endfunction 735 endfunction