lipsol.m

数学建模的源代码

function [xsol,fval,lambda,exitflag,output] = lipsol(f,Aineq,bineq,Aeq,beq,lb,ub,options,computeLambda)
%LIPSOL  Linear programming Interior-Point SOLver.
%   X = LIPSOL(f,A,b) will solves the Linear Programming problem
%
%          min f'*x    subject to:   A*x <= b
%           x
%
%   X = LIPSOL(f,A,b,Aeq,beq) solves the problem above while additionally
%   satisfying the equality constraints Aeq*x = beq.
%
%   X = LIPSOL(f,A,b,Aeq,beq,LB,UB) defines a set of lower and upper bounds
%   on the design variables, X, so that the solution X is always in the
%   range LB <= X <= UB.  If LB is [] or if any of the entries in LB are -Inf,
%   that is taken to mean the corresponding variable is not bounded below.
%   If UB is [] or if any of the entries in UB are Inf, that is
%   taken to mean the corresponding variable is not bounded above.
%
%   X = LIPSOL(f,A,b,Aeq,beq,LB,UB,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, TolFun, LargeScale, MaxIter. Use OPTIONS = [] as a 
%   place holder if no options are set.
%
%   [X,FVAL] = LIPSOL(f,A,b) returns the value of the objective
%   function at X: FVAL = f' * X.
%
%   [X,FVAL,LAMBDA] = LIPSOL(f,A,b) returns the set of Lagrange multipliers,
%   LAMBDA, at the solution, X.
%   LAMBDA.ineqlin are the Lagrange multipliers for the inequalities.
%   LAMBDA.eqlin are the Lagrange multipliers for the equalities.
%   LAMBDA.lb the Lagrange multipliers for the lower bounds.
%   LAMBDA.ub are the Lagrange multipliers for the upper bounds.
%
%   [X,FVAL,LAMBDA,EXITFLAG] = LIPSOL(f,A,b) describes the exit conditions:
%   If EXITFLAG is:
%      > 0 then LIPSOL converged with a solution X.
%      = 0 then LIPSOL did not converge within the maximum number of iterations.
%      < 0 then the problem was infeasible.
%
%   [X,FVAL,LAMBDA,EXITFLAG,OUTPUT] = LIPSOL(f,A,b) returns a structure:
%     output.iterations: number of iterations.
%     OUTPUT.cgiterations: total number of PCG iterations.

%   Original author Yin Zhang 1995.
%   Modified by The MathWorks, Inc. since 1997.
%   Copyright (c) 1990-98 The University of Maryland Baltimore County and The MathWorks, Inc.
%   $Revision: 1.10 $ $Date: 1998/09/14 19:07:56 $

if (nargin < 3)
   error('Not enough input arguments.')
end

n = size(f,1);
n_orig = n;
if (size(f,2) ~= 1)
   error('First input to lipsol must be a column vector of weights.')
end
f_orig = f;

if (nargin < 8)
   options = [];
   if ( (nargin < 5) | isempty(Aeq) )
      Aeq = sparse(0,n);
      beq = zeros(0,1);
   end
end

if isempty(Aineq)
   Aineq = sparse(0,n);
   bineq = zeros(0,1);
end

tol = optimget(options,'TolFun');
maxiter = optimget(options,'maxiter');
switch optimget(options,'display')
case 'none'
   diagnostic_level = 0;
case 'iter'
   diagnostic_level = 2;
case 'final'
   diagnostic_level = 1;
case 'testing'
   diagnostic_level = 4;
otherwise
   diagnostic_level = 1;
end

showstat = optimget(options,'ShowStatusWindow','off');
switch showstat
case {'iter','iterplus','iterplusbounds'}
   showstat = 1;
case {'final','none','off'}
   showstat = 0;
otherwise
   showstat = 0;
end

if ((nargin < 7) | isempty(ub))
   ub = repmat(Inf,n,1);
   if (diagnostic_level >= 4)
      fprintf('  Assuming all variables are unbounded above.\n');
   end
end
if ((nargin < 6) | isempty(lb))
   lb = repmat(-Inf,n,1);
   if (diagnostic_level >= 4)
      fprintf('  Assuming all variables are unbounded below.\n');
   end
end
lbinf = ~isfinite(lb);
nlbinf = sum(lbinf);

% Collect some information for lagrange multiplier calculations
if computeLambda
   lbIn = lb;
   ubIn = ub;
   lbnotinf = ~lbinf;
   ubnotinf = isfinite(ub);
   numRowsLb = sum(lbnotinf);
   numRowsUb = sum(ubnotinf);
   numRowsAineq = size(Aineq,1);
   numRowsAeq = size(Aeq,1);
   % Indices for lagrange multipliers and initialize multipliers
   Aindex = 1:numRowsAineq+numRowsAeq;
   Aineqstart = 1; Aineqend = 0;  % end value is zero since not in A yet
   Aeqstart = 1; Aeqend = 0;      % end value is zero since not in A yet
   lbindex = 1:length(lb);
   ubindex = 1:length(ub);
   lbend = length(lb);  % we don't remove the inf ones, we just ignore them
   ubend = length(ub);  % we don't remove the inf ones, we just ignore them
   lambda.ineqlin = zeros(numRowsAineq,1);
   lambda.eqlin = zeros(numRowsAeq,1);
   lambda.upper = zeros(n_orig,1);
   lambda.lower = zeros(n_orig,1);
   extraAineqvars = 0;
   ilbinfubfin = [];
   Aineq_orig = Aineq; Aeq_orig = Aeq;
else
   lambda.ineqlin = [];
   lambda.eqlin = [];
   lambda.upper = [];
   lambda.lower = [];
end

% IF one of the variables is unbounded below: -Inf == lb(i) <= x(i) <= ub(i)
% THEN create a new variable splitting the old x(i) into its difference:
% x(i) -> x(i) - x(n+1) such that both are positive:
% 0 <= x(i) <= Inf and 0 <= x(n+1) <= Inf
% and a new inequality constraint reflecting the old upper bound:
% x(i) - x(n+1) <= ub(i), unless ub(i) == Inf.
% The inequality and equality constraints must be padded to reflect the new
% variable and the function f must also be adjusted.
if (nlbinf > 0)
   f = [f; -f(lbinf)];
   Aineq = [Aineq -Aineq(:,lbinf)];
   Aeq = [Aeq -Aeq(:,lbinf)];
   lbinfubfin = lbinf & isfinite(ub);
   ilbinfubfin = find(lbinfubfin);
   nlbinfubfin = sum(lbinfubfin);
   Aineq = [Aineq; ...
         sparse(1:nlbinfubfin,find(lbinfubfin),1,nlbinfubfin,n) ...
         sparse(1:nlbinfubfin,find(lbinfubfin(lbinf)),-1,nlbinfubfin,nlbinf)];
   bineq = [bineq; ub(lbinfubfin)];
   lb(lbinf) = 0;
   lb = [lb; zeros(nlbinf,1)];
   ub(lbinf) = Inf;
   ub = [ub; repmat(Inf,nlbinf,1)];
   n = n + nlbinf;
   if computeLambda
      extraAineqvars = nlbinfubfin;
      Aindex = 1:length(Aindex)+ extraAineqvars;
      ubindex = 1:length(ub);
      lbindex = 1:length(lb);
   end
   if (diagnostic_level >= 4)
      fprintf(['  Splitting %d (unbounded below) variables into the' ...
            ' difference of two strictly\n  positive variables:' ...
            ' x(i) = xpos-xneg, such that xpos >= 0, xneg >= 0.\n'], nlbinf);
      if (nlbinfubfin > 0)
         fprintf(['  Adding %d inequality constraints of the form:' ...
               ' xpos - xneg <= ub(i) for the finite upper bounds.\n'], ...
            nlbinfubfin);
      end
   end
end

nslack = n;
if isempty(Aineq) % Equality constraints only; no xpos/xneg variables
   A = Aeq;
   b = beq;
   m = size(A,1);
   if computeLambda
      Aeqstart = 1; Aeqend = m;
   end
   if (diagnostic_level >= 4)
      fprintf(['  Linear Programming problem has %d equality' ...
            ' constraints on %d variables.\n'],m,n);
   end
else
   mineq = size(Aineq,1);
   meq = size(Aeq,1);
   m = mineq + meq;
   f = [f; sparse(mineq,1)];
   A = [Aineq speye(mineq); Aeq sparse(meq,mineq)];
   b = [bineq; beq];
   lb = [lb; sparse(mineq,1)];
   ub = [ub; repmat(Inf,mineq,1)];
   if computeLambda
      Aineqstart = 1; Aineqend = numRowsAineq;
      Aeqstart = Aineqend + extraAineqvars + 1; Aeqend = Aeqstart + numRowsAeq - 1;
      lbindex = 1:length(lb);
      ubindex = 1:length(ub);
   end
   if (diagnostic_level >= 4)
      fprintf(['  Adding %d slack variables, one for each inequality' ...
            ' constraint, to the existing\n  set of %d variables,'],mineq,n);
   end
   n = n + mineq;
   if (diagnostic_level >= 4)
      fprintf([' resulting in %d equality constraints on %d variables.\n'], ...
         mineq,n);
      if  ~isempty(Aeq)
         fprintf(['  Combining the %d (formerly inequality) constraints' ...
               ' with %d equality constraints,\n  resulting in %d equality' ...
               ' constraints on %d variables.\n'],mineq,meq,m,n);
      end
   end
end

% cholinc(A,'inf') only accepts sparse matrices.
if (~issparse(A))
   A = sparse(A);
end
f = sparse(f);
data_changed = 0;

% Perform preliminary solving, rearranging and checking for infeasibility

if computeLambda
   % We might delete variables from ub/lb (and so w/z) so we need
   % a separate index into w/z 
   windex = zeros(length(ub),1);
   windex(ubindex) = 1;
   zindex = zeros(length(lb),1);
   zindex(lbindex) = 1;
   yindex = ones(size(A,1),1);
   % Save f at this point -- may be needed to compute multipliers
   f_before_deletes = f;
end

% delete fixed variables
fixed = (lb == ub);
Fixed_exist = any(fixed);
if (Fixed_exist)
   ifix = find(fixed); 
   notfixed = ~fixed;
   infx = find(notfixed);
   if computeLambda
      % Decide which bound multiplier to be nonzero according to -f(i) in case
      % we later find no other constraints are active involving this x(i);
      %  leave the other zero.
      % If we find that one of the other constraints on x(i) is active, then
      %   it is ok for both these multipliers to remain zero.
      fixedlower = ( ( f >= 0 ) & fixed );
      fixedupper = ( ( f < 0) &  fixed );
      lbindex = lbindex(infx);
      zindex = zindex(infx); 
      ubindex = ubindex(infx);
      windex = windex(infx);
   end
   xfix = lb(ifix);
   f = f(infx);
   b = b - A(:,ifix)*xfix;
   A = A(:,infx);
   lb = lb(infx);
   ub = ub(infx);
   data_changed = 1;
   if (diagnostic_level >= 4)
      fprintf(['  Assigning %i equal values of' ...
            ' lower and upper bounds to solution.\n'],length(ifix))
   end
end

% delete zero rows
rnnzct = sum(spones(A'));
if ~isempty(A) & (any(rnnzct == 0))
   zrows = (rnnzct == 0);
   izrows = find(zrows);
   if (any(b(izrows) ~= 0))
      if (diagnostic_level >= 1)
         fprintf(['Exiting due to infeasibility:   an all zero row in the constraint' ...
               ' matrix does not have a zero in corresponding right hand size entry.\n']);
      end
      xsol = [];
      fval = [];
      lambda.ineqlin = [];
      lambda.eqlin = [];
      lambda.upper = [];
      lambda.lower = [];
      exitflag = -1;
      output = [];
      return
   else   
      nzrows = ~zrows;
      inzrows = find(nzrows);
      A = A(inzrows,:);
      b = b(inzrows); 
      rnnzct = rnnzct(inzrows);
      if (diagnostic_level >= 4)
         fprintf(['  Deleting %i all zero rows from' ...
               ' constraint matrix.\n'],length(izrows));
      end
      if computeLambda
         % multipliers for zero rows are zero
         rnnzct == 0;
         Aindex = Aindex(inzrows);
         yindex = yindex(inzrows);
      end
      data_changed = 1;
      [m,n] = size(A);
   end
end

% make A structurally "full rank"
sprk = sprank(A');
Row_deleted = 0;
if (sprk < m)
   Row_deleted = 1;
   [dmp, tmp] = dmperm(A);
   irow = dmp(1:sprk);          % Note permutation of rows might occur 
   idelrow = dmp(sprk+1:m);
   Adel = A(idelrow,:);
   bdel = b(idelrow);
   A = A(irow,:);
   b = b(irow);
   if computeLambda
      Aindex = Aindex(irow);
      yindex = yindex(irow);
   end
   rnnzct = rnnzct(irow);
   if (diagnostic_level >= 4)
      fprintf(['  Deleting %i dependent rows from' ...
            ' constraint matrix.\n'],m-sprk);
   end
   data_changed = 1;
end

% delete zero columns
Zrcols_exist = 0;
if isempty(A)
   zrcol = 0;
else
   zrcol = (max(abs(A)) == 0)';  % max behaves differently if A only has one row:
                                 %    so zero cols won't be deleted in this case.
end
if (any(zrcol == 1))
   Zrcols_exist = 1;
   izrcol = find(zrcol);
   if ((any(f(izrcol) < 0) & (ub(izrcol) == 0)))
      if (diagnostic_level >= 1)
         fprintf('Exiting due to infeasibility:   objective f''*x is unbounded below.\n');
      end
      xsol = [];
      fval = [];
      lambda.ineqlin = [];
      lambda.eqlin = [];
      lambda.upper = [];