lipsol.m

线性规划算法

function [xsol,fval,lambda,exitflag,output] = lipsol(f,Aineq,bineq,Aeq,beq,lb,ub,options,defaultopt,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 1990-2002 The University of Maryland Baltimore County and The MathWorks, Inc. 
%   $Revision: 1.21 $ $Date: 2002/03/12 20:36:18 $ 
 
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',defaultopt,'fast'); 
maxiter = optimget(options,'MaxIter',defaultopt,'fast'); 
switch optimget(options,'Display',defaultopt,'fast') 
case {'none','off'} 
   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');  % no default, so use old optimget 
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+i) such that both are positive: 
% 0 <= x(i) <= Inf and 0 <= x(n+i) <= Inf 
% and a new inequality constraint reflecting the old upper bound: 
% x(i) - x(n+i) <= 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; 
   % update n 
   n = size(f,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'), 1);   
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 = size(A,1); 
   end 
end 
 
% make A structurally "full rank" 
sprk = sprank(A'); 
Row_deleted = 0; 
if (sprk < m) & ~isempty(A) 
   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; 
   m = size(A,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 = []; 
      lambda.lower = []; 
      exitflag = -1; 
      output = [];