📄 solvopt.m
字号:
function [x,f,options]=solvopt(x,fun,grad,options,func,gradc)
% Usage:
% [x,f,options]=solvopt(x,fun,grad,options,func,gradc)
% The function SOLVOPT performs a modified version of Shor's r-algorithm in
% order to find a local minimum resp. maximum of a nonlinear function
% defined on the n-dimensional Euclidean space
% or
% a solution of a nonlinear constrained problem:
% min { f(x): g(x) (<)= 0, g(x) in R(m), x in R(n) }
% Arguments:
% x is the n-vector (row or column) of the coordinates of the starting
% point,
% fun is the name of an M-file (M-function) which computes the value
% of the objective function <fun> at a point x,
% synopsis: f=fun(x)
% grad is the name of an M-file (M-function) which computes the gradient
% vector of the function <fun> at a point x,
% synopsis: g=grad(x)
% func is the name of an M-file (M-function) which computes the MAXIMAL
% RESIDUAL(!) for a set of constraints at a point x,
% synopsis: fc=func(x)
% gradc is the name of an M-file (M-function) which computes the gradient
% vector for the maximal residual consyraint at a point x,
% synopsis: gc=gradc(x)
% options is a row vector of optional parameters:
% options(1)= H, where sign(H)=-1 resp. sign(H)=+1 means minimize
% resp. maximize <fun> (valid only for unconstrained problem)
% and H itself is a factor for the initial trial step size
% (options(1)=-1 by default),
% options(2)= relative error for the argument in terms of the
% infinity-norm (1.e-4 by default),
% options(3)= relative error for the function value (1.e-6 by default),
% options(4)= limit for the number of iterations (15000 by default),
% options(5)= control of the display of intermediate results and error
% resp. warning messages (default value is 0, i.e., no intermediate
% output but error and warning messages, see more in the manual),
% options(6)= admissible maximal residual for a set of constraints
% (options(6)=1e-8 by default, see more in the manual),
% *options(7)= the coefficient of space dilation (2.5 by default),
% *options(8)= lower bound for the stepsize used for the difference
% approximation of gradients (1e-12 by default, see more in the manual).
% (* ... changes should be done with care)
% Returned values:
% x is the optimizer (row resp. column),
% f is the optimum function value,
% options returns the values of the counters
% options(9), the number of iterations, if positive,
% or an abnormal stop code, if negative (see more in the manual),
% options(10), the number of objective
% options(11), the number of gradient evaluations,
% options(12), the number of constraint function evaluations,
% options(13), the number of constraint gradient evaluations.
% ____________________________________________________________________________
% strings: ----{
errmes='SolvOpt error:';
wrnmes='SolvOpt warning:';
error1='No function name and/or starting point passed to the function.';
error2='Argument X has to be a row or column vector of dimension > 1.';
error30='<fun> returns an empty string.';
error31='Function value does not exist (NaN is returned).';
error32='Function equals infinity at the point.';
error40='<grad> returns an improper matrix. Check the dimension.';
error41='Gradient does not exist (NaN is returned by <grad>).';
error42='Gradient equals infinity at the starting point.';
error43='Gradient equals zero at the starting point.';
error50='<func> returns an empty string.';
error51='<func> returns NaN at the point.';
error52='<func> returns infinite value at the point.';
error60='<gradc> returns an improper vector. Check the dimension';
error61='<gradc> returns NaN at the point.';
error62='<gradc> returns infinite vector at the point.';
error63='<gradc> returns zero vector at an infeasible point.';
error5='Function is unbounded.';
error6='Choose another starting point.';
warn1= 'Gradient is zero at the point, but stopping criteria are not fulfilled.';
warn20='Normal re-setting of a transformation matrix.' ;
warn21='Re-setting due to the use of a new penalty coefficient.' ;
warn4= 'Iterations limit exceeded.';
warn31='The function is flat in certain directions.';
warn32='Trying to recover by shifting insensitive variables.';
warn09='Re-run from recorded point.';
warn08='Ravine with a flat bottom is detected.';
termwarn0='SolvOpt: Normal termination.';
termwarn1='SolvOpt: Termination warning:';
appwarn='The above warning may be reasoned by inaccurate gradient approximation';
endwarn=[...
'Premature stop is possible. Try to re-run the routine from the obtained point. ';...
'Result may not provide the optimum. The function apparently has many extremum points. ';...
'Result may be inaccurate in the coordinates. The function is flat at the optimum. ';...
'Result may be inaccurate in a function value. The function is extremely steep at the optimum.'];
% ----}
% ARGUMENTS PASSED ----{
if nargin<2 % Function and/or starting point are not specified
options(9)=-1; disp(errmes); disp(error1); return
end
if nargin<3, app=1; % No user-supplied gradients
elseif isempty(grad), app=1;
else, app=0; % Exact gradients are supplied
end
% OPTIONS ----{
doptions=[-1,1.e-4,1.e-6,15000,0,1.e-8,2.5,1e-11];
if nargin<4, options=doptions;
elseif isempty(options), options=doptions;
else,
% Replace default options by user specified options:
ii=find(options~=0);doptions(ii)=options(ii);
options=doptions;
end
% Check the values:
options([2:4,6:8])=abs(options([2:4,6:8]));
options(2:3)=max(options(2:3),[1.e-12,1.e-12]);
options(2)=max(options(8)*1.e2,options(2));
options(2:3)=min(options(2:3),[1,1]);
options(6)=max(options(6),1e-12);
options(7)=max([options(7),1.5]);
options(8)=max(options(8),1e-11);
% ----}
if nargin<5, constr=0; % Unconstrained problem
elseif isempty(func), constr=0;
else, constr=1; % Constrained problem
if nargin<6, appconstr=1; t=3; % No user-supplied gradients for constraints
elseif isempty(gradc),
appconstr=1;
else, appconstr=0; % Exact gradients of constraints are supplied
end
end
% ----}
% STARTING POINT ----{
if max(size(x))<=1, disp(errmes); disp(error2);
options(9)=-2; return
elseif size(x,2)==1, n=size(x,1); x=x'; trx=1;
elseif size(x,1)==1, n=size(x,2); trx=0;
else, disp(errmes); disp(error2);
options(9)=-2; return
end
% ----}
% WORKING CONSTANTS AND COUNTERS ----{
options(10)=0; options(11)=0; % function and gradient calculations
if constr
options(12)=0; options(13)=0; % same for constraints
end
epsnorm=1.e-15;epsnorm2=1.e-30; % epsilon & epsilon^2
if constr, h1=-1; % NLP: restricted to minimization
cnteps=options(6); % Max. admissible residual
else, h1=sign(options(1)); % Minimize resp. maximize a function
end
k=0; % Iteration counter
wdef=1/options(7)-1; % Default space transf. coeff.
%Gamma control ---{
ajb=1+.1/n^2; % Base I
ajp=20;
ajpp=ajp; % Start value for the power
ajs=1.15; % Base II
knorms=0; gnorms=zeros(1,10); % Gradient norms stored
%---}
%Display control ---{
if options(5)<=0, dispdata=0;
if options(5)==-1, dispwarn=0; else, dispwarn=1; end
else, dispdata=round(options(5)); dispwarn=1;
end, ld=dispdata;
%---}
%Stepsize control ---{
dq=5.1; % Step divider (at f_{i+1}>gamma*f_{i})
du20=2;du10=1.5;du03=1.05; % Step multipliers (at certain steps made)
kstore=3;nsteps=zeros(1,kstore); % Steps made at the last 'kstore' iterations
if app, des=6.3; % Desired number of steps per 1-D search
else, des=3.3; end
mxtc=3; % Number of trial cycles (steep wall detect)
%---}
termx=0; limxterm=50; % Counter and limit for x-criterion
ddx =max(1e-11,options(8)); % stepsize for gradient approximation
low_bound=-1+1e-4; % Lower bound cosine used to detect a ravine
ZeroGrad=n*1.e-16; % Lower bound for a gradient norm
nzero=0; % Zero-gradient events counter
% Lower bound for values of variables taking into account
lowxbound=max([options(2),1e-3]);
% Lower bound for function values to be considered as making difference
lowfbound=options(3)^2;
krerun=0; % Re-run events counter
detfr=options(3)*100; % relative error for f/f_{record}
detxr=options(2)*10; % relative error for norm(x)/norm(x_{record})
warnno=0; % the number of warn.mess. to end with
kflat=0; % counter for points of flatness
stepvanish=0; % counter for vanished steps
stopf=0;
% ----} End of setting constants
% ----} End of the preamble
% COMPUTE THE FUNCTION ( FIRST TIME ) ----{
if trx, f=feval(fun,x');
else, f=feval(fun,x); end
options(10)=options(10)+1;
if isempty(f), if dispwarn,disp(errmes);disp(error30);end
options(9)=-3; if trx, x=x';end, return
elseif isnan(f), if dispwarn,disp(errmes);disp(error31);disp(error6);end
options(9)=-3; if trx, x=x';end, return
elseif abs(f)==Inf, if dispwarn,disp(errmes);disp(error32);disp(error6);end
options(9)=-3; if trx, x=x';end, return
end
xrec=x; frec=f; % record point and function value
% Constrained problem
if constr, fp=f; kless=0;
if trx, fc=feval(func,x');
else, fc=feval(func,x); end
if isempty(fc),
if dispwarn,disp(errmes);disp(error50);end
options(9)=-5; if trx, x=x';end, return
elseif isnan(fc),
if dispwarn,disp(errmes);disp(error51);disp(error6);end
options(9)=-5; if trx, x=x';end, return
elseif abs(fc)==Inf,
if dispwarn,disp(errmes);disp(error52);disp(error6);end
options(9)=-5; if trx, x=x';end, return
end
options(12)=options(12)+1;
PenCoef=1; % first rough approximation
if fc<=cnteps, FP=1; fc=0; % feasible point
else, FP=0; % infeasible point
end
f=f+PenCoef*fc;
end
% ----}
% COMPUTE THE GRADIENT ( FIRST TIME ) ----{
if app, deltax=h1*ddx*ones(size(x));
if constr, if trx, g=apprgrdn(x',fp,fun,deltax',1);
else, g=apprgrdn(x ,fp,fun,deltax,1); end
else, if trx, g=apprgrdn(x',f,fun,deltax',1);
else, g=apprgrdn(x ,f,fun,deltax,1); end
end, options(10)=options(10)+n;
else, if trx, g=feval(grad,x');
else, g=feval(grad,x); end
options(11)=options(11)+1;
end
if size(g,2)==1, g=g'; end, ng=norm(g);
if size(g,2)~=n, if dispwarn,disp(errmes);disp(error40);end
options(9)=-4; if trx, x=x';end, return
elseif isnan(ng), if dispwarn,disp(errmes);disp(error41);disp(error6);end
options(9)=-4; if trx, x=x';end, return
elseif ng==Inf, if dispwarn,disp(errmes);disp(error42);disp(error6);end
options(9)=-4; if trx, x=x';end, return
elseif ng<ZeroGrad, if dispwarn,disp(errmes);disp(error43);disp(error6);end
options(9)=-4; if trx, x=x';end, return
end
if constr, if ~FP
if appconstr,
deltax=sign(x); idx=find(deltax==0);
deltax(idx)=ones(size(idx)); deltax=ddx*deltax;
if trx, gc=apprgrdn(x',fc,func,deltax',0);
else, gc=apprgrdn(x ,fc,func,deltax ,0); end
options(12)=options(12)+n;
else, if trx, gc=feval(gradc,x');
else, gc=feval(gradc,x); end
options(13)=options(13)+1;
end
if size(gc,2)==1, gc=gc'; end, ngc=norm(gc);
if size(gc,2)~=n,
if dispwarn,disp(errmes);disp(error60);end
options(9)=-6; if trx, x=x';end, return
elseif isnan(ngc),
if dispwarn,disp(errmes);disp(error61);disp(error6);end
options(9)=-6; if trx, x=x';end, return
elseif ngc==Inf,
if dispwarn,disp(errmes);disp(error62);disp(error6);end
options(9)=-6; if trx, x=x';end, return
elseif ngc<ZeroGrad,
if dispwarn,disp(errmes);disp(error63);end
options(9)=-6; if trx, x=x';end, return
end
g=g+PenCoef*gc; ng=norm(g);
end, end
grec=g; nng=ng;
% ----}
% INITIAL STEPSIZE
h=h1*sqrt(options(2))*max(abs(x)); % smallest possible stepsize
if abs(options(1))~=1,
h=h1*max(abs([options(1),h])); % user-supplied stepsize
else,
h=h1*max(1/log(ng+1.1),abs(h)); % calculated stepsize
end
% RESETTING LOOP ----{
while 1,
kcheck=0; % Set checkpoint counter.
kg=0; % stepsizes stored
kj=0; % ravine jump counter
B=eye(n); % re-set transf. matrix to identity
fst=f; g1=g; dx=0;
% ----}
% MAIN ITERATIONS ----{
while 1,
k=k+1;kcheck=kcheck+1;
laststep=dx;
% ADJUST GAMMA --{
gamma=1+max([ajb^((ajp-kcheck)*n),2*options(3)]);
gamma=min([gamma,ajs^max([1,log10(nng+1)])]);
% --}
gt=g*B; w=wdef;
% JUMPING OVER A RAVINE ----{
if (gt/norm(gt))*(g1'/norm(g1))<low_bound
if kj==2, xx=x; end, if kj==0, kd=4; end,
kj=kj+1; w=-.9; h=h*2; % use large coef. of space dilation
if kj>2*kd, kd=kd+1; warnno=1;
if any(abs(x-xx)<epsnorm*abs(x)), % flat bottom is detected
if dispwarn,disp(wrnmes);disp(warn08); end
end
end
else, kj=0;
end
% ----}
% DILATION ----{
z=gt-g1;
nrmz=norm(z);
if(nrmz>epsnorm*norm(gt))
z=z/nrmz;
g1=gt+w*(z*gt')*z; B=B+w*(B*z')*z;
else
z=zeros(1,n); nrmz=0; g1=gt;
end
d1=norm(g1); g0=(g1/d1)*B';
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -