function [x,fval,exitflag,output,lambda]=linprog(f,A,B,Aeq,Beq,lb,ub,x0,options)
%LINPROG Linear programming.
%   X=LINPROG(f,A,b) attempts to solve the linear programming problem:
%
%            min f'*x    subject to:   A*x <= b
%             x
%
%   X=LINPROG(f,A,b,Aeq,beq) solves the problem above while additionally
%   satisfying the equality constraints Aeq*x = beq.
%
%   X=LINPROG(f,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=LINPROG(f,A,b,Aeq,beq,LB,UB,X0) sets the starting point to X0.  This
%   option is only available with the active-set algorithm.  The default
%   interior point algorithm will ignore any non-empty starting point.
%
%   X=LINPROG(f,A,b,Aeq,beq,LB,UB,X0,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.
%   Options are Display, Diagnostics, TolFun, LargeScale, MaxIter.
%   Currently, only 'final' and 'off' are valid values for the parameter
%   Display when LargeScale is 'off' ('iter' is valid when LargeScale is 'on').
%
%   [X,FVAL]=LINPROG(f,A,b) returns the value of the objective function at X:
%   FVAL = f'*X.
%
%   [X,FVAL,EXITFLAG] = LINPROG(f,A,b) returns an EXITFLAG that describes the
%   exit condition of LINPROG. Possible values of EXITFLAG and the corresponding
%   exit conditions are
%
%     1  LINPROG converged to a solution X.
%     0  Maximum number of iterations reached.
%    -2  No feasible point found.
%    -3  Problem is unbounded.
%    -4  NaN value encountered during execution of algorithm.
%    -5  Both primal and dual problems are infeasible.
%    -7  Search direction became too small; no further progress can be made.
%
%   [X,FVAL,EXITFLAG,OUTPUT] = LINPROG(f,A,b) returns a structure
%   OUTPUT with the number of iterations taken in OUTPUT.iterations, the type
%   of algorithm used in OUTPUT.algorithm, the number of conjugate gradient
%   iterations (if used) in OUTPUT.cgiterations, and the exit message in
%   OUTPUT.message.
%
%   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(f,A,b) returns the set of
%   Lagrangian multipliers LAMBDA, at the solution: LAMBDA.ineqlin for the
%   linear inequalities A, LAMBDA.eqlin for the linear equalities Aeq,
%   LAMBDA.lower for LB, and LAMBDA.upper for UB.
%
%   NOTE: the LargeScale (the default) version of LINPROG uses a primal-dual
%         method. Both the primal problem and the dual problem must be feasible
%         for convergence. Infeasibility messages of either the primal or dual,
%         or both, are given as appropriate.  The primal problem in standard
%         form is
%              min f'*x such that A*x = b, x >= 0.
%         The dual problem is
%              max b'*y such that A'*y + s = f, s >= 0.


%   $Revision.2 $  $Date: 2004/12/06 16:36:12 $

% If just 'defaults' passed in, return the default options in X

% Default MaxIter is set to [] because its value depends on the algorithm.
defaultopt = struct('Display','final',...
   'TolFun',[],'Diagnostics','off',...
   'LargeScale','on','MaxIter',[], ...
   'Simplex','off');

if nargin==1 && nargout <= 1 && isequal(f,'defaults')
   x = defaultopt;
   return
end

% Handle missing arguments
if nargin < 9, options = [];
   if nargin < 8, x0 = [];
      if nargin < 7, ub = [];
         if nargin < 6, lb = [];
            if nargin < 5, Beq = [];
               if nargin < 4, Aeq = [];
               end, end, end, end, end, end

% Check for non-double inputs
% SUPERIORFLOAT errors when superior input is neither single nor double;
% We use try-catch to override SUPERIORFLOAT's error message when input
% data type is integer.
try
    dataType = superiorfloat(f,A,B,Aeq,Beq,lb,ub,x0);
    if ~isequal('double', dataType)
        error('optim:linprog:NonDoubleInput', ...
            'LINPROG only accepts inputs of data type double.')
    end
catch
    error('optim:linprog:NonDoubleInput', ...
        'LINPROG only accepts inputs of data type double.')
end

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

% Options setup
if isfield(options,'Simplex')
    useSimplex = isequal(optimget(options,'Simplex',defaultopt,'fast'), 'on');
else
    useSimplex = isequal(defaultopt.Simplex, 'on');
end

largescale = isequal(optimget(options,'LargeScale',defaultopt,'fast'),'on');
diagnostics = isequal(optimget(options,'Diagnostics',defaultopt,'fast'),'on');
switch optimget(options,'Display',defaultopt,'fast')
case {'off','none'}
   verbosity = 0;
case 'iter'
   verbosity = 2;
case 'final'
   verbosity = 1;
otherwise
   verbosity = 1;
end

% Set the constraints up: defaults and check size
[nineqcstr,nvarsineq]=size(A);
[neqcstr, nvarseq]=size(Aeq);
nvars = max([length(f),nvarsineq,nvarseq]); % In case A is empty
ncstr = nineqcstr + neqcstr;

if isempty(f), f=zeros(nvars,1); end
if isempty(A), A=zeros(0,nvars); end
if isempty(B), B=zeros(0,1); end
if isempty(Aeq), Aeq=zeros(0,nvars); end
if isempty(Beq), Beq=zeros(0,1); end

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

if ~isequal(length(B),nineqcstr)
    error('optim:linprog:SizeMismatchRowsOfA', ...
          'The number of rows in A must be the same as the length of b.');
elseif ~isequal(length(Beq),neqcstr)
    error('optim:linprog:SizeMismatchRowsOfAeq', ...
          'The number of rows in Aeq must be the same as the length of beq.');
elseif ~isequal(length(f),nvarsineq) && ~isempty(A)
    error('optim:linprog:SizeMismatchColsOfA', ...
          'The number of columns in A must be the same as the length of f.');
elseif ~isequal(length(f),nvarseq) && ~isempty(Aeq)
    error('optim:linprog:SizeMismatchColsOfAeq', ...
        'The number of columns in Aeq must be the same as the length of f.');
end

[x0,lb,ub,msg] = checkbounds(x0,lb,ub,nvars);
if ~isempty(msg)
   exitflag = -2;
   x=x0; fval = []; lambda = [];
   output.iterations = 0;
   output.algorithm = ''; % not known at this stage
   output.cgiterations = [];
   output.message = msg;
   if verbosity > 0
      disp(msg)
   end
   return
end

caller = 'linprog';
ncstr = nineqcstr + neqcstr;

if largescale
   OUTPUT.algorithm = 'large-scale: interior point';
elseif useSimplex
   OUTPUT.algorithm = 'medium-scale: simplex';
else
   OUTPUT.algorithm  = 'medium-scale: active-set';
end

if diagnostics
   % Do diagnostics on information so far
   gradflag = []; hessflag = []; line_search=[];
   constflag = 0; gradconstflag = 0; non_eq=0;non_ineq=0;
   lin_eq=size(Aeq,1); lin_ineq=size(A,1); XOUT=ones(nvars,1);
   funfcn{1} = [];ff=[]; GRAD=[];HESS=[];
   confcn{1}=[];c=[];ceq=[];cGRAD=[];ceqGRAD=[];
   msg = diagnose('linprog',OUTPUT,gradflag,hessflag,constflag,gradconstflag,...
      line_search,options,defaultopt,XOUT,non_eq,...
      non_ineq,lin_eq,lin_ineq,lb,ub,funfcn,confcn,ff,GRAD,HESS,c,ceq,cGRAD,ceqGRAD);
end

if (largescale)
    if (useSimplex)
        warning('optim:linprog:IgnoreSimplexOn', ...
                ['Simplex method does not currently solve large-scale problems;\n' ...
                 'calling large-scale interior point method. (To run simplex, set\n' ...
                 'OPTIONS.LargeScale to ''off''.)'])
    end
    if ~isempty(x0) && verbosity > 0
        warning('optim:linprog:IgnoreStartPoint', ...
                ['Large scale (interior point) method uses a built-in starting point;\n' ...
                 'ignoring user-supplied X0.'])
    end
    % Set the default values of TolFun and MaxIter for this algorithm
    defaultopt.TolFun = 1e-8;
    defaultopt.MaxIter = 85;
    [x,fval,lambda,exitflag,output] = lipsol(f,A,B,Aeq,Beq,lb,ub,options,defaultopt,computeLambda);
elseif (useSimplex)
    if ~isempty(x0) && verbosity > 0
        warning('optim:linprog:IgnoreStartPoint', ...
               ['Simplex method uses a built-in starting point; ignoring user-supplied X0.'])
    end
    % Set the default values of TolFun and MaxIter for this algorithm
    defaultopt.TolFun = 1e-6;
    defaultopt.MaxIter = '100*NumberOfVariables';
    [x,fval,lambda,exitflag,output] = simplex(f,A,B,Aeq,Beq,lb,ub,options,defaultopt,computeLambda);
    % Remap exitflags if necessary
    if exitflag == -1
      exitflag = -2;
    elseif exitflag == -2
      exitflag = -3;
    end
else
    if ~largescale && (issparse(A) || issparse(Aeq) )% asked for medium-scale but sparse
        if verbosity > 0
            disp('The medium-scale (active-set) algorithm does not currently handle sparse matrices.')
            disp('Converting to full matrices to solve.')
        end
    end
    if isempty(x0), x0=zeros(nvars,1); end
    % Set the default value of MaxIter for this algorithm
    defaultopt.MaxIter = '10*max(NumberOfVariables,NumberOfInequalities+NumberOfBounds)';
    [x,lambdaqp,exitflag,output,dum1,dum2,msg]= ...
        qpsub([],full(f),full([Aeq;A]),full([Beq;B]),lb,ub,x0,neqcstr,verbosity,caller,ncstr, ...
              nvars,options,defaultopt);
    output.algorithm = 'medium-scale: active-set';
end

if isequal(OUTPUT.algorithm , 'medium-scale: active-set')
   fval = f'*x;
   if computeLambda
       llb = length(lb);
       lub = length(ub);
       lambda.lower = zeros(llb,1);
       lambda.upper = zeros(lub,1);
       arglb = ~isinf(lb); lenarglb = nnz(arglb);
       argub = ~isinf(ub); lenargub = nnz(argub);
       lambda.eqlin = lambdaqp(1:neqcstr,1);
       lambda.ineqlin = lambdaqp(neqcstr+1:neqcstr+nineqcstr,1);
       lambda.lower(arglb) = lambdaqp(neqcstr+nineqcstr+1:neqcstr+nineqcstr+lenarglb);
       lambda.upper(argub) = lambdaqp(neqcstr+nineqcstr+lenarglb+1:neqcstr+nineqcstr+lenarglb+lenargub);
   end
   output.cgiterations =[];

   if exitflag == 1
     normalTerminationMsg = sprintf('Optimization terminated.');
     if verbosity > 0
       disp(normalTerminationMsg)
     end
     if isempty(msg)
       output.message = normalTerminationMsg;
     else
       % append normal termination msg to current output msg
       output.message = sprintf('%s\n%s',msg,normalTerminationMsg);
     end
   else
     output.message = msg;
   end
end

Error using linprog (line 99)
LINPROG only accepts inputs of data type double.

Copyright 1990-2004 The MathWorks, Inc.
Published with MATLAB® 7.14