Mercurial > octave

comparison scripts/optimization/sqp.m @ 31653:972dcc46bb41

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
sqp.m: Ignore missing lower/upper bounds (bug #32008) * sqp.m: Avoid setting lower and upper bounds to +/-realmax when unspecified in the function call. Set bounds to +/-Inf in main function. Update constraint subfunction to correctly return empty distance to constraint when bounds are unspecified.
author Julien Bect <julien.bect@supelec.fr>
date Thu, 08 Dec 2022 15:23:29 -0500
parents c05ef94a2bbc
children 5f11de0e7440
comparison
equal deleted inserted replaced
31652:a18897e4c7b5 31653:972dcc46bb41
119 ## The fifth and sixth arguments, @var{lb} and @var{ub}, contain lower and 119 ## The fifth and sixth arguments, @var{lb} and @var{ub}, contain lower and
120 ## upper bounds on @var{x}. These must be consistent with the equality and 120 ## upper bounds on @var{x}. These must be consistent with the equality and
121 ## inequality constraints @var{g} and @var{h}. If the arguments are vectors 121 ## inequality constraints @var{g} and @var{h}. If the arguments are vectors
122 ## then @var{x}(i) is bound by @var{lb}(i) and @var{ub}(i). A bound can also 122 ## then @var{x}(i) is bound by @var{lb}(i) and @var{ub}(i). A bound can also
123 ## be a scalar in which case all elements of @var{x} will share the same 123 ## be a scalar in which case all elements of @var{x} will share the same
124 ## bound. If only one bound (lb, ub) is specified then the other will 124 ## bound.
125 ## default to (-@var{realmax}, +@var{realmax}).
126 ## 125 ##
127 ## The seventh argument @var{maxiter} specifies the maximum number of 126 ## The seventh argument @var{maxiter} specifies the maximum number of
128 ## iterations. The default value is 100. 127 ## iterations. The default value is 100.
129 ## 128 ##
130 ## The eighth argument @var{tolerance} specifies the tolerance for the stopping 129 ## The eighth argument @var{tolerance} specifies the tolerance for the stopping
286 endif 285 endif
287 else 286 else
288 ## constraint inequality function with bounds present 287 ## constraint inequality function with bounds present
289 lb_idx = ub_idx = true (size (x0)); 288 lb_idx = ub_idx = true (size (x0));
290 ub_grad = - (lb_grad = eye (rows (x0))); 289 ub_grad = - (lb_grad = eye (rows (x0)));
290
291 ## if unspecified set ub and lb to +/-Inf, preserving single type
292 if (isempty (lb))
293 if (isa (x0, "single"))
294 lb = - inf (size (x0), "single");
295 else
296 lb = - inf (size (x0));
297 endif
298 endif
299
291 if (isvector (lb)) 300 if (isvector (lb))
292 globals.lb = tmp_lb = lb(:); 301 lb = lb(:);
293 lb_idx(:) = tmp_idx = (lb != -Inf); 302 lb_idx(:) = (lb != -Inf);
294 globals.lb = globals.lb(tmp_idx, 1); 303 globals.lb = lb(lb_idx, 1);
295 lb_grad = lb_grad(lb_idx, :); 304 lb_grad = lb_grad(lb_idx, :);
296 elseif (isempty (lb))
297 if (isa (x0, "single"))
298 globals.lb = tmp_lb = -realmax ("single");
299 else
300 globals.lb = tmp_lb = -realmax ();
301 endif
302 else 305 else
303 error ("sqp: invalid lower bound"); 306 error ("sqp: invalid lower bound");
304 endif 307 endif
305 308
309 if (isempty (ub))
310 if (isa (x0, "single"))
311 ub = inf (size (x0), "single");
312 else
313 ub = inf (size (x0));
314 endif
315 endif
316
306 if (isvector (ub)) 317 if (isvector (ub))
307 globals.ub = tmp_ub = ub(:); 318 ub = ub(:);
308 ub_idx(:) = tmp_idx = (ub != Inf); 319 ub_idx(:) = (ub != Inf);
309 globals.ub = globals.ub(tmp_idx, 1); 320 globals.ub = ub(ub_idx, 1);
310 ub_grad = ub_grad(ub_idx, :); 321 ub_grad = ub_grad(ub_idx, :);
311 elseif (isempty (ub))
312 if (isa (x0, "single"))
313 globals.ub = tmp_ub = realmax ("single");
314 else
315 globals.ub = tmp_ub = realmax ();
316 endif
317 else 322 else
318 error ("sqp: invalid upper bound"); 323 error ("sqp: invalid upper bound");
319 endif 324 endif
320 325
321 if (any (tmp_lb > tmp_ub)) 326 if (any (globals.lb > globals.ub))
322 error ("sqp: upper bound smaller than lower bound"); 327 error ("sqp: upper bound smaller than lower bound");
323 endif 328 endif
324 bounds_grad = [lb_grad; ub_grad]; 329 bounds_grad = [lb_grad; ub_grad];
325 ci_fcn = @(x) cf_ub_lb (x, lb_idx, ub_idx, globals); 330 ci_fcn = @(x) cf_ub_lb (x, lb_idx, ub_idx, globals);
326 ci_grd = @(x) cigrad_ub_lb (x, bounds_grad, globals); 331 ci_grd = @(x) cigrad_ub_lb (x, bounds_grad, globals);
703 708
704 endfunction 709 endfunction
705 710
706 711
707 function res = cf_ub_lb (x, lbidx, ubidx, globals) 712 function res = cf_ub_lb (x, lbidx, ubidx, globals)
708 713 ## function returning constraint evaluated at x, and distance between ub/lb
709 ## combine constraint function with ub and lb 714 ## and x, when they are specified (otherwise return empty)
710 if (isempty (globals.cifcn)) 715
711 res = [x(lbidx,1)-globals.lb; globals.ub-x(ubidx,1)]; 716 ## inequality constraints
712 else 717 if (! isempty (globals.cifcn))
713 res = [feval(globals.cifcn,x); x(lbidx,1)-globals.lb; 718 ci = feval (globals.cifcn, x);
714 globals.ub-x(ubidx,1)]; 719 else
715 endif 720 ci = [];
721 endif
722
723 ## lower bounds
724 if (! isempty (globals.lb))
725 lb = x(lbidx, 1) - globals.lb;
726 else
727 lb = [];
728 endif
729
730 ## upper bounds
731 if (! isempty (globals.ub))
732 ub = globals.ub - x(ubidx, 1);
733 else
734 ub = [];
735 end
736
737 res = [ci; lb; ub];
716 738
717 endfunction 739 endfunction
718 740
719 741
720 function res = cigrad_ub_lb (x, bgrad, globals) 742 function res = cigrad_ub_lb (x, bgrad, globals)