FMINSEARCHCON: Extension of FMINSEARCHBND with general inequality constraints
usage: x=FMINSEARCHCON(fun,x0)
usage: x=FMINSEARCHCON(fun,x0,LB)
usage: x=FMINSEARCHCON(fun,x0,LB,UB)
usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b)
usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b,nonlcon)
usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b,nonlcon,options)
usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b,nonlcon,options,p1,p2,...)
usage: [x,fval,exitflag,output]=FMINSEARCHCON(fun,x0,...)
arguments:
fun, x0, options - see the help for FMINSEARCH
x0 MUST be a feasible point for the linear and nonlinear
inequality constraints. If it is not inside the bounds
then it will be moved to the nearest bound. If x0 is
infeasible for the general constraints, then an error will
be returned.
LB - lower bound vector or array, must be the same size as x0
If no lower bounds exist for one of the variables, then
supply -inf for that variable.
If no lower bounds at all, then LB may be left empty.
Variables may be fixed in value by setting the corresponding
lower and upper bounds to exactly the same value.
UB - upper bound vector or array, must be the same size as x0
If no upper bounds exist for one of the variables, then
supply +inf for that variable.
If no upper bounds at all, then UB may be left empty.
Variables may be fixed in value by setting the corresponding
lower and upper bounds to exactly the same value.
A,b - (OPTIONAL) Linear inequality constraint array and right
hand side vector. (Note: these constraints were chosen to
be consistent with those of fmincon.)
A*x <= b
nonlcon - (OPTIONAL) general nonlinear inequality constraints
NONLCON must return a set of general inequality constraints.
These will be enforced such that nonlcon is always <= 0.
nonlcon(x) <= 0
Notes:
If options is supplied, then TolX will apply to the transformed
variables. All other FMINSEARCH parameters should be unaffected.
Variables which are constrained by both a lower and an upper
bound will use a sin transformation. Those constrained by
only a lower or an upper bound will use a quadratic
transformation, and unconstrained variables will be left alone.
Variables may be fixed by setting their respective bounds equal.
In this case, the problem will be reduced in size for FMINSEARCH.
The bounds are inclusive inequalities, which admit the
boundary values themselves, but will not permit ANY function
evaluations outside the bounds. These constraints are strictly
followed.
If your problem has an EXCLUSIVE (strict) constraint which will
not admit evaluation at the bound itself, then you must provide
a slightly offset bound. An example of this is a function which
contains the log of one of its parameters. If you constrain the
variable to have a lower bound of zero, then FMINSEARCHCON may
try to evaluate the function exactly at zero.
Inequality constraints are enforced with an implicit penalty
function approach. But the constraints are tested before
any function evaluations are ever done, so the actual objective
function is NEVER evaluated outside of the feasible region.
Example usage:
rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
Fully unconstrained problem
fminsearchcon(rosen,[3 3])
ans =
1.0000 1.0000
lower bound constrained
fminsearchcon(rosen,[3 3],[2 2],[])
ans =
2.0000 4.0000
x(2) fixed at 3
fminsearchcon(rosen,[3 3],[-inf 3],[inf,3])
ans =
1.7314 3.0000
simple linear inequality: x(1) + x(2) <= 1
fminsearchcon(rosen,[0 0],[],[],[1 1],.5)
ans =
0.6187 0.3813
general nonlinear inequality: sqrt(x(1)^2 + x(2)^2) <= 1
fminsearchcon(rosen,[0 0],[],[],[],[],@(x) norm(x)-1)
ans =
0.78633 0.61778
Of course, any combination of the above constraints is
also possible.
See test_main.m for other examples of use.
See also: fminsearch, fminspleas, fminsearchbnd
Author: John D'Errico
E-mail: woodchips@rochester.rr.com
Release: 1.0
Release date: 12/16/060001 function [x,fval,exitflag,output]=fminsearchcon(fun,x0,LB,UB,A,b,nonlcon,options,varargin) 0002 % FMINSEARCHCON: Extension of FMINSEARCHBND with general inequality constraints 0003 % usage: x=FMINSEARCHCON(fun,x0) 0004 % usage: x=FMINSEARCHCON(fun,x0,LB) 0005 % usage: x=FMINSEARCHCON(fun,x0,LB,UB) 0006 % usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b) 0007 % usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b,nonlcon) 0008 % usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b,nonlcon,options) 0009 % usage: x=FMINSEARCHCON(fun,x0,LB,UB,A,b,nonlcon,options,p1,p2,...) 0010 % usage: [x,fval,exitflag,output]=FMINSEARCHCON(fun,x0,...) 0011 % 0012 % arguments: 0013 % fun, x0, options - see the help for FMINSEARCH 0014 % 0015 % x0 MUST be a feasible point for the linear and nonlinear 0016 % inequality constraints. If it is not inside the bounds 0017 % then it will be moved to the nearest bound. If x0 is 0018 % infeasible for the general constraints, then an error will 0019 % be returned. 0020 % 0021 % LB - lower bound vector or array, must be the same size as x0 0022 % 0023 % If no lower bounds exist for one of the variables, then 0024 % supply -inf for that variable. 0025 % 0026 % If no lower bounds at all, then LB may be left empty. 0027 % 0028 % Variables may be fixed in value by setting the corresponding 0029 % lower and upper bounds to exactly the same value. 0030 % 0031 % UB - upper bound vector or array, must be the same size as x0 0032 % 0033 % If no upper bounds exist for one of the variables, then 0034 % supply +inf for that variable. 0035 % 0036 % If no upper bounds at all, then UB may be left empty. 0037 % 0038 % Variables may be fixed in value by setting the corresponding 0039 % lower and upper bounds to exactly the same value. 0040 % 0041 % A,b - (OPTIONAL) Linear inequality constraint array and right 0042 % hand side vector. (Note: these constraints were chosen to 0043 % be consistent with those of fmincon.) 0044 % 0045 % A*x <= b 0046 % 0047 % nonlcon - (OPTIONAL) general nonlinear inequality constraints 0048 % NONLCON must return a set of general inequality constraints. 0049 % These will be enforced such that nonlcon is always <= 0. 0050 % 0051 % nonlcon(x) <= 0 0052 % 0053 % 0054 % Notes: 0055 % 0056 % If options is supplied, then TolX will apply to the transformed 0057 % variables. All other FMINSEARCH parameters should be unaffected. 0058 % 0059 % Variables which are constrained by both a lower and an upper 0060 % bound will use a sin transformation. Those constrained by 0061 % only a lower or an upper bound will use a quadratic 0062 % transformation, and unconstrained variables will be left alone. 0063 % 0064 % Variables may be fixed by setting their respective bounds equal. 0065 % In this case, the problem will be reduced in size for FMINSEARCH. 0066 % 0067 % The bounds are inclusive inequalities, which admit the 0068 % boundary values themselves, but will not permit ANY function 0069 % evaluations outside the bounds. These constraints are strictly 0070 % followed. 0071 % 0072 % If your problem has an EXCLUSIVE (strict) constraint which will 0073 % not admit evaluation at the bound itself, then you must provide 0074 % a slightly offset bound. An example of this is a function which 0075 % contains the log of one of its parameters. If you constrain the 0076 % variable to have a lower bound of zero, then FMINSEARCHCON may 0077 % try to evaluate the function exactly at zero. 0078 % 0079 % Inequality constraints are enforced with an implicit penalty 0080 % function approach. But the constraints are tested before 0081 % any function evaluations are ever done, so the actual objective 0082 % function is NEVER evaluated outside of the feasible region. 0083 % 0084 % 0085 % Example usage: 0086 % rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2; 0087 % 0088 % Fully unconstrained problem 0089 % fminsearchcon(rosen,[3 3]) 0090 % ans = 0091 % 1.0000 1.0000 0092 % 0093 % lower bound constrained 0094 % fminsearchcon(rosen,[3 3],[2 2],[]) 0095 % ans = 0096 % 2.0000 4.0000 0097 % 0098 % x(2) fixed at 3 0099 % fminsearchcon(rosen,[3 3],[-inf 3],[inf,3]) 0100 % ans = 0101 % 1.7314 3.0000 0102 % 0103 % simple linear inequality: x(1) + x(2) <= 1 0104 % fminsearchcon(rosen,[0 0],[],[],[1 1],.5) 0105 % 0106 % ans = 0107 % 0.6187 0.3813 0108 % 0109 % general nonlinear inequality: sqrt(x(1)^2 + x(2)^2) <= 1 0110 % fminsearchcon(rosen,[0 0],[],[],[],[],@(x) norm(x)-1) 0111 % ans = 0112 % 0.78633 0.61778 0113 % 0114 % Of course, any combination of the above constraints is 0115 % also possible. 0116 % 0117 % See test_main.m for other examples of use. 0118 % 0119 % 0120 % See also: fminsearch, fminspleas, fminsearchbnd 0121 % 0122 % 0123 % Author: John D'Errico 0124 % E-mail: woodchips@rochester.rr.com 0125 % Release: 1.0 0126 % Release date: 12/16/06 0127 0128 % size checks 0129 xsize = size(x0); 0130 x0 = x0(:); 0131 n=length(x0); 0132 0133 if (nargin<3) || isempty(LB) 0134 LB = repmat(-inf,n,1); 0135 else 0136 LB = LB(:); 0137 end 0138 if (nargin<4) || isempty(UB) 0139 UB = repmat(inf,n,1); 0140 else 0141 UB = UB(:); 0142 end 0143 0144 if (n~=length(LB)) || (n~=length(UB)) 0145 error 'x0 is incompatible in size with either LB or UB.' 0146 end 0147 0148 % defaults for A,b 0149 if (nargin<5) || isempty(A) 0150 A = []; 0151 end 0152 if (nargin<6) || isempty(b) 0153 b = []; 0154 end 0155 nA = []; 0156 nb = []; 0157 if (isempty(A)&&~isempty(b)) || (isempty(b)&&~isempty(A)) 0158 error 'Sizes of A and b are incompatible' 0159 elseif ~isempty(A) 0160 nA = size(A); 0161 b = b(:); 0162 nb = size(b,1); 0163 if nA(1)~=nb 0164 error 'Sizes of A and b are incompatible' 0165 end 0166 if nA(2)~=n 0167 error 'A is incompatible in size with x0' 0168 end 0169 end 0170 0171 % defaults for nonlcon 0172 if (nargin<7) || isempty(nonlcon) 0173 nonlcon = []; 0174 end 0175 0176 % test for feasibility of the initial value 0177 % against any general inequality constraints 0178 if ~isempty(A) 0179 if any(A*x0>b) 0180 error 'Infeasible starting values (linear inequalities failed).' 0181 end 0182 end 0183 if ~isempty(nonlcon) 0184 if any(feval(nonlcon,(reshape(x0,xsize)),varargin{:})>0) 0185 error 'Infeasible starting values (nonlinear inequalities failed).' 0186 end 0187 end 0188 0189 % set default options if necessary 0190 if (nargin<8) || isempty(options) 0191 options = optimset('fminsearch'); 0192 end 0193 0194 % stuff into a struct to pass around 0195 params.args = varargin; 0196 params.LB = LB; 0197 params.UB = UB; 0198 params.fun = fun; 0199 params.n = n; 0200 params.xsize = xsize; 0201 0202 params.OutputFcn = []; 0203 0204 params.A = A; 0205 params.b = b; 0206 params.nonlcon = nonlcon; 0207 0208 % 0 --> unconstrained variable 0209 % 1 --> lower bound only 0210 % 2 --> upper bound only 0211 % 3 --> dual finite bounds 0212 % 4 --> fixed variable 0213 params.BoundClass = zeros(n,1); 0214 for i=1:n 0215 k = isfinite(LB(i)) + 2*isfinite(UB(i)); 0216 params.BoundClass(i) = k; 0217 if (k==3) && (LB(i)==UB(i)) 0218 params.BoundClass(i) = 4; 0219 end 0220 end 0221 0222 % transform starting values into their unconstrained 0223 % surrogates. Check for infeasible starting guesses. 0224 x0u = x0; 0225 k=1; 0226 for i = 1:n 0227 switch params.BoundClass(i) 0228 case 1 0229 % lower bound only 0230 if x0(i)<=LB(i) 0231 % infeasible starting value. Use bound. 0232 x0u(k) = 0; 0233 else 0234 x0u(k) = sqrt(x0(i) - LB(i)); 0235 end 0236 0237 % increment k 0238 k=k+1; 0239 case 2 0240 % upper bound only 0241 if x0(i)>=UB(i) 0242 % infeasible starting value. use bound. 0243 x0u(k) = 0; 0244 else 0245 x0u(k) = sqrt(UB(i) - x0(i)); 0246 end 0247 0248 % increment k 0249 k=k+1; 0250 case 3 0251 % lower and upper bounds 0252 if x0(i)<=LB(i) 0253 % infeasible starting value 0254 x0u(k) = -pi/2; 0255 elseif x0(i)>=UB(i) 0256 % infeasible starting value 0257 x0u(k) = pi/2; 0258 else 0259 x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1; 0260 % shift by 2*pi to avoid problems at zero in fminsearch 0261 % otherwise, the initial simplex is vanishingly small 0262 x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k)))); 0263 end 0264 0265 % increment k 0266 k=k+1; 0267 case 0 0268 % unconstrained variable. x0u(i) is set. 0269 x0u(k) = x0(i); 0270 0271 % increment k 0272 k=k+1; 0273 case 4 0274 % fixed variable. drop it before fminsearch sees it. 0275 % k is not incremented for this variable. 0276 end 0277 0278 end 0279 % if any of the unknowns were fixed, then we need to shorten 0280 % x0u now. 0281 if k<=n 0282 x0u(k:n) = []; 0283 end 0284 0285 % were all the variables fixed? 0286 if isempty(x0u) 0287 % All variables were fixed. quit immediately, setting the 0288 % appropriate parameters, then return. 0289 0290 % undo the variable transformations into the original space 0291 x = xtransform(x0u,params); 0292 0293 % final reshape 0294 x = reshape(x,xsize); 0295 0296 % stuff fval with the final value 0297 fval = feval(params.fun,x,params.args{:}); 0298 0299 % fminsearchbnd was not called 0300 exitflag = 0; 0301 0302 output.iterations = 0; 0303 output.funcCount = 1; 0304 output.algorithm = 'fminsearch'; 0305 output.message = 'All variables were held fixed by the applied bounds'; 0306 0307 % return with no call at all to fminsearch 0308 return 0309 end 0310 0311 % Check for an outputfcn. If there is any, then substitute my 0312 % own wrapper function. 0313 if ~isempty(options.OutputFcn) 0314 params.OutputFcn = options.OutputFcn; 0315 options.OutputFcn = @outfun_wrapper; 0316 end 0317 0318 % now we can call fminsearch, but with our own 0319 % intra-objective function. 0320 [xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params); 0321 0322 % undo the variable transformations into the original space 0323 x = xtransform(xu,params); 0324 0325 % final reshape 0326 x = reshape(x,xsize); 0327 0328 % Use a nested function as the OutputFcn wrapper 0329 function stop = outfun_wrapper(x,varargin); 0330 % we need to transform x first 0331 xtrans = xtransform(x,params); 0332 0333 % then call the user supplied OutputFcn 0334 stop = params.OutputFcn(xtrans,varargin{1:(end-1)}); 0335 0336 end 0337 0338 end % mainline end 0339 0340 % ====================================== 0341 % ========= begin subfunctions ========= 0342 % ====================================== 0343 function fval = intrafun(x,params) 0344 % transform variables, test constraints, then call original function 0345 0346 % transform 0347 xtrans = xtransform(x,params); 0348 0349 % test constraints before the function call 0350 0351 % First, do the linear inequality constraints, if any 0352 if ~isempty(params.A) 0353 % Required: A*xtrans <= b 0354 if any(params.A*xtrans(:) > params.b) 0355 % linear inequality constraints failed. Just return inf. 0356 fval = inf; 0357 return 0358 end 0359 end 0360 0361 % resize xtrans to be the correct size for the nonlcon 0362 % and objective function calls 0363 xtrans = reshape(xtrans,params.xsize); 0364 0365 % Next, do the nonlinear inequality constraints 0366 if ~isempty(params.nonlcon) 0367 % Required: nonlcon(xtrans) <= 0 0368 cons = feval(params.nonlcon,xtrans,params.args{:}); 0369 if any(cons(:) > 0) 0370 % nonlinear inequality constraints failed. Just return inf. 0371 fval = inf; 0372 return 0373 end 0374 end 0375 0376 % we survived the general inequality constraints. Only now 0377 % do we evaluate the objective function. 0378 0379 % append any additional parameters to the argument list 0380 fval = feval(params.fun,xtrans,params.args{:}); 0381 0382 end % sub function intrafun end 0383 0384 % ====================================== 0385 function xtrans = xtransform(x,params) 0386 % converts unconstrained variables into their original domains 0387 0388 xtrans = zeros(1,params.n); 0389 % k allows some variables to be fixed, thus dropped from the 0390 % optimization. 0391 k=1; 0392 for i = 1:params.n 0393 switch params.BoundClass(i) 0394 case 1 0395 % lower bound only 0396 xtrans(i) = params.LB(i) + x(k).^2; 0397 0398 k=k+1; 0399 case 2 0400 % upper bound only 0401 xtrans(i) = params.UB(i) - x(k).^2; 0402 0403 k=k+1; 0404 case 3 0405 % lower and upper bounds 0406 xtrans(i) = (sin(x(k))+1)/2; 0407 xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i); 0408 % just in case of any floating point problems 0409 xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i))); 0410 0411 k=k+1; 0412 case 4 0413 % fixed variable, bounds are equal, set it at either bound 0414 xtrans(i) = params.LB(i); 0415 case 0 0416 % unconstrained variable. 0417 xtrans(i) = x(k); 0418 0419 k=k+1; 0420 end 0421 end 0422 0423 end % sub function xtransform end 0424 0425 0426 0427 0428