fmincon.m

数学建模的源代码

function [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = fmincon(FUN,X,A,B,Aeq,Beq,LB,UB,NONLCON,options,varargin)
%求解下列形式非线性规划问题：
%     min f(x)  s.t.  Ax<= b,  Aeqx=beq, 
%                   c(x)<=0, ceq(x) = 0, lb<=x<=ub
%调用格式为：
% x=fmincon(fun, 初值,A,b,Aeq,beq,lb,ub,nonlcon)
%      此时当约束条件中缺A和b、Aeq和beq或lb和ub时，
%       相关项可用[ ]代替以表示省略。
%      fun写成如下的M-函数形式 (fun.m) ：
%                function  f = fun (x)
%                        f = f(x); 
%        非线性约束条件写成如下的M-函数形式 (nonlcon.m) ：
%         function [c,ceq]=nonlcon(x)
%        c = c(x);ceq=ceq(x);
% 注意：方程变量必须拼成一个向量变量，即用x(1),x(2),...
%
%例题
%   max x*y*z
%   -x+x*y+2*z>=0
%   x+2*y+2*z<=72
%   10<=y<=20
%    x-y=10
%  先化为
%   min -x(1)*x(2)*x(3)
%   x(1)+2*x(2)+2*x(3)<=72
%   x(1)-x(2)=10
%   x(1)-x(1)*x(2)-2*x(3)<=0
%   10<=x(2)<=20
%  写M函数 fconfun.m
%          function  f=fconfun(x)
%          f=-x(1)*x(2)*x(3)
%  再写M函数 fconfun2.m
%          function  [g,geq]=fconfun2(x)
%          g=x(1)-x(1)*x(2)-2*x(3);
%          geq=0;
%  求解
%   x0=[10,10,10];
%   A=[1 2 2];b=72;
%   Aeq=[1 -1 0];beq=10;
%   [x,f]=fmincon('fconfun',x0,A,b,Aeq,beq,[-inf,10,-inf]',[inf,20,inf],'fconfun2')
%
%FMINCON Finds the constrained minimum of a function of several variables.
%   FMINCON solves problems of the form:
%       min F(X)  subject to:  A*X  <= B, Aeq*X  = Beq (linear constraints)
%        X                       C(X) <= 0, Ceq(X) = 0   (nonlinear constraints)
%                                LB <= X <= UB            
%                                                           
%   X=FMINCON(FUN,X0,A,B) starts at X0 and finds a minimum X to the function
%   described in FUN, subject to the linear inequalities A*X <= B. X0 may be a 
%   scalar, vector or matrix. The function FUN (usually an M-file or inline object) 
%   should return a scalar function value F evaluated at X when called with 
%   feval: F=feval(FUN,X).
%
%   X=FMINCON(FUN,X0,A,B,Aeq,Beq) minimizes FUN subject to the linear equalities
%   Aeq*X = Beq as well as A*X <= B. (Set A=[] and B=[] if no inequalities exist.)
%
%   X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB) defines a set of lower and upper
%   bounds on the design variables, X, so that the solution is in 
%   the range LB <= X <= UB. Use empty matrices for LB and UB
%   if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; 
%   set UB(i) = Inf if X(i) is unbounded above.
%
%   X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON) subjects the minimization to the 
%   constraints defined in NONLCON. The function NONLCON should return the vectors
%   C and Ceq, representing the nonlinear inequalities and equalities respectively, 
%   when called with feval: [C, Ceq] = feval(NONLCON,X). FMINCON minimizes
%   FUN such that C(X)<=0 and Ceq(X)=0. (Set LB=[] and/or UB=[] if no bounds exist.)
%
%   X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS) minimizes with the 
%   default optimization parameters replaced by values in the structure OPTIONS, 
%   an argument created with the OPTIMSET function.  See OPTIMSET for details.  Used
%   options are Display, TolX, TolFun, TolCon, DerivativeCheck, Diagnostics, GradObj, 
%   GradConstr, Hessian, MaxFunEvals, MaxIter, DiffMinChange and DiffMaxChange, 
%   LargeScale, MaxPCGIter, PrecondBandWidth, TolPCG, TypicalX, HessPattern. 
%   Use the GradObj option to specify that FUN may be called with two output 
%   arguments where the second, G, is the partial derivatives of the
%   function df/dX, at the point X: [F,G] = feval(FUN,X). Use the GradConstr
%   option to specify that NONLCON may be called with four output arguments:
%   [C,Ceq,GC,GCeq] = feval(NONLCON,X) where GC is the partial derivatives of the 
%   constraint vector of inequalities C an GCeq is the partial derivatives of the 
%   constraint vector of equalities Ceq. Use OPTIONS = [] as a place holder if 
%   no options are set.
%  
%   X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS,P1,P2,...) passes the 
%   problem-dependent parameters P1,P2,... directly to the functions FUN 
%   and NONLCON: feval(FUN,X,P1,P2,...) and feval(NONLCON,X,P1,P2,...).  Pass
%   empty matrices for A, B, Aeq, Beq, OPTIONS, LB, UB, and NONLCON to use the 
%   default values.
%
%   [X,FVAL]=FMINCON(FUN,X0,...) returns the value of the objective 
%   function FUN at the solution X.
%
%   [X,FVAL,EXITFLAG]=FMINCON(FUN,X0,...) returns a string EXITFLAG that 
%   describes the exit condition of FMINCON.  
%   If EXITFLAG is:
%      > 0 then FMINCON converged to a solution X.
%      0   then the maximum number of function evaluations was reached.
%      < 0 then FMINCON did not converge to a solution.
%   
%   [X,FVAL,EXITFLAG,OUTPUT]=FMINCON(FUN,X0,...) returns a structure
%   OUTPUT with the number of iterations taken in OUTPUT.iterations, the number
%   of function evaluations in OUTPUT.funcCount, the algorithm used in 
%   OUTPUT.algorithm, the number of CG iterations (if used) in OUTPUT.cgiterations, 
%   and the first-order optimality (if used) in OUTPUT.firstorderopt.
%
%   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=FMINCON(FUN,X0,...) returns the Lagrange multipliers
%   at the solution X: LAMBDA.lower for LB, LAMBDA.upper for UB, LAMBDA.ineqlin is
%   for the linear inequalities, LAMBDA.eqlin is for the linear equalities,
%   LAMBDA.ineqnonlin is for the nonlinear inequalities, and LAMBDA.eqnonlin
%   is for the nonlinear equalities.
%
%   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD]=FMINCON(FUN,X0,...) returns the value of 
%   the gradient of FUN at the solution X.
%
%   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]=FMINCON(FUN,X0,...) returns the 
%   value of the HESSIAN of FUN at the solution X.
%
%   See also OPTIMSET, FMINUNC, FMINBND, FMINSEARCH.

%   Copyright (c) 1990-98 by The MathWorks, Inc.
%   $Revision: 1.19 $  $Date: 1998/10/22 20:11:07 $

defaultopt = optimset('display','final','LargeScale','on', ...
   'TolX',1e-6,'TolFun',1e-6,'TolCon',1e-6,'DerivativeCheck','off',...
   'Diagnostics','off',...
   'GradObj','off','GradConstr','off','Hessian','off','MaxFunEvals','100*numberOfVariables',...
   'DiffMaxChange',1e-1,'DiffMinChange',1e-8,...
   'PrecondBandWidth',0,'TypicalX','ones(numberOfVariables,1)','MaxPCGIter','max(1,floor(numberOfVariables/2))', ...
   'TolPCG',0.1,'MaxIter',400,'HessPattern',[]);
% If just 'defaults' passed in, return the default options in X
if nargin==1 & nargout <= 1 & isequal(FUN,'defaults')
   X = defaultopt;
   return
end

large = 'large-scale';
medium = 'medium-scale';

if nargin < 4, error('FMINCON requires at least four input arguments'); end
if nargin < 10, options=[];
   if nargin < 9, NONLCON=[];
      if nargin < 8, UB = [];
         if nargin < 7, LB = [];
            if nargin < 6, Beq=[];
               if nargin < 5, Aeq =[];
               end, end, end, end, end, end
if isempty(NONLCON) & isempty(A) & isempty(Aeq) & isempty(UB) & isempty(LB)
   error('FMINCON is for constrained problems. Use FMINUNC for unconstrained problems.')
end

if nargout > 4
   computeLambda = 1;
else 
   computeLambda = 0;
end

caller='constr';
lenVarIn = length(varargin);
XOUT=X(:);
numberOfVariables=length(XOUT);

options = optimset(defaultopt,options);
switch optimget(options,'display')
case {'off','none'}
   verbosity = 0;
case 'iter'
   verbosity = 2;
case 'final'
   verbosity = 1;
otherwise
   verbosity = 1;
end

% Set to column vectors
B = B(:);
Beq = Beq(:);

[XOUT,l,u,msg] = checkbounds(XOUT,LB,UB,numberOfVariables);
if ~isempty(msg)
   EXITFLAG = -1;
   [FVAL,OUTPUT,LAMBDA,GRAD,HESSIAN] = deal([]);
   X(:)=XOUT;
   if verbosity > 0
      disp(msg)
   end
   return
end
lFinite = l(~isinf(l));
uFinite = u(~isinf(u));


meritFunctionType = 0;

diagnostics = isequal(optimget(options,'diagnostics','off'),'on');
gradflag =  strcmp(optimget(options,'GradObj'),'on');
hessflag = strcmp(optimget(options,'Hessian'),'on');
if isempty(NONLCON)
   constflag = 0;
else
   constflag = 1;
end
gradconstflag =  strcmp(optimget(options,'GradConstr'),'on');
line_search = strcmp(optimget(options,'largescale','off'),'off'); % 0 means trust-region, 1 means line-search

% Convert to inline function as needed
if ~isempty(FUN)  % will detect empty string, empty matrix, empty cell array
   [funfcn, msg] = fprefcnchk(FUN,'fmincon',length(varargin),gradflag,hessflag);
else
   errmsg = sprintf('%s\n%s', ...
      'FUN must be a function name, valid string expression, or inline object;', ...
      ' or, FUN may be a cell array that contains these type of objects.');
   error(errmsg)
end

if constflag % NONLCON is non-empty
   [confcn, msg] = fprefcnchk(NONLCON,'fmincon',length(varargin),gradconstflag,[],1);
else
   confcn{1} = '';
end

[rowAeq,colAeq]=size(Aeq);
% if only l and u then call sfminbx
if ~line_search & isempty(NONLCON) & isempty(A) & isempty(Aeq) & gradflag
   OUTPUT.algorithm = large;
   % if only Aeq beq and Aeq has as many columns as rows, then call sfminle
elseif ~line_search & isempty(NONLCON) & isempty(A) & isempty(lFinite) & isempty(uFinite) & gradflag ...
      & colAeq >= rowAeq
   OUTPUT.algorithm = large;
elseif ~line_search
   warning(['Trust region method does not currently solve this type of problem,',...
         sprintf('\n'), 'switching to line search.'])
   if isequal(funfcn{1},'fungradhess')
      funfcn{1}='fungrad';
      warning('Hessian provided by user will be ignored in line search algorithm')
      
   elseif  isequal(funfcn{1},'fun_then_grad_then_hess')
      funfcn{1}='fun_then_grad';
      warning('Hessian provided by user will be ignored in line search algorithm')
   end    
   hessflag = 0;
   OUTPUT.algorithm = medium;
elseif line_search
   OUTPUT.algorithm = medium;
   if issparse(Aeq) | issparse(A)
      warning('can not do sparse with line_search, converting to full')
   end
   
   % else call nlconst
else
   error('Unrecognized combination of OPTIONS flags and calling sequence.')
end


lenvlb=length(l);
lenvub=length(u);

if isequal(OUTPUT.algorithm,medium)
   CHG = 1e-7*abs(XOUT)+1e-7*ones(numberOfVariables,1);
   i=1:lenvlb;
   lindex = XOUT(i)<l(i);
   if any(lindex),
      XOUT(lindex)=l(lindex)+1e-4; 
   end
   i=1:lenvub;
   uindex = XOUT(i)>u(i);
   if any(uindex)
      XOUT(uindex)=u(uindex);
      CHG(uindex)=-CHG(uindex);
   end
   X(:) = XOUT;
else
   arg = (u >= 1e10); arg2 = (l <= -1e10);
   u(arg) = inf*ones(length(arg(arg>0)),1);
   l(arg2) = -inf*ones(length(arg2(arg2>0)),1);
   if min(min(u-XOUT),min(XOUT-l)) < 0, 
      XOUT = startx(u,l);
      X(:) = XOUT;
   end
end

% Evaluate function
GRAD=zeros(numberOfVariables,1);
HESS = [];

switch funfcn{1}
case 'fun'
   try
      f = feval(funfcn{3},X,varargin{:});
   catch
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
case 'fungrad'
   try
      [f,GRAD(:)] = feval(funfcn{3},X,varargin{:});
   catch
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
case 'fungradhess'
   try
      [f,GRAD(:),HESS] = feval(funfcn{3},X,varargin{:});
   catch
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
case 'fun_then_grad'
   try
      f = feval(funfcn{3},X,varargin{:});
   catch
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
   try
      GRAD(:) = feval(funfcn{4},X,varargin{:});
   catch
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective gradient function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
case 'fun_then_grad_then_hess'
   try
      f = feval(funfcn{3},X,varargin{:});
   catch
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
   try
      GRAD(:) = feval(funfcn{4},X,varargin{:});
   catch
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective gradient function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
   try
      HESS = feval(funfcn{5},X,varargin{:});
   catch 
      errmsg = sprintf('%s\n%s\n\n%s',...
         'FMINCON cannot continue because user supplied objective Hessian function', ...
         ' failed with the following error:', lasterr);
      error(errmsg);
   end
otherwise
   error('Undefined calltype in FMINCON');
end

% Evaluate constraints
switch confcn{1}
case 'fun'
   try