phaseii.doc

线性规划工具箱-performing the standard two phase simplex method on linear programming problems.具体介绍里面有

function [z,xbasic,ibasic,ienter,iter,PCOL,OPTEST,CYCTEST]=phaseii(A,b,c,ibasic);
%PHASEII performs phase II of the simplex method starting with the basic
%columns specified by the vector ibasic.  The output values are z, the last 
%value of the objective function; xbasic, the vector of basic variable values;
%ibasic, the vector of basic variable indices; ienter, the index of the entering %variable when an unbounded ray is detected; iter, the number of phase II %iterations performed; PCOL, the pivot column vector corresponding to an %unbounded ray when an unbounded ray is detected; OPTEST, a flag set to 1 when %the current solution is the optimal solution for the original LP; and CYCTEST, %a flag set to 1 when the cycling test detects possible cycling.
%See also PHASEI and LINPROG.
%
%Written by Jeff Stuart, Department of Mathematics, University of Southern %Mississippi, Hattiesburg, MS 39406. December, 1993. Revised, October, 1997.
[e-mail]
%
[m,n]=size(A);
PCOL=[];              %INITIALIZE PCOL 
iter=0;               %INITIALIZE PHASE II ITERATION COUNTER
cycle=0;              %INITIALIZE CYCLING COUNTER   
CYCTEST=0;            %INITIALIZE CYCLE TEST FLAGX=zeros(1,n);         %INITIALIZE VARIABLES STORAGE VECTOR
J=X;                  %CREATE BASIC VARIABLES SELECTION VECTOR 
J(ibasic)=ones(1,m);  %INITIALIZE BASIC VARIABLES SELECTION VECTOR
K=[1:n];
inon=K(~J);           %DETERMINE NONBASIC INDICES 
B=A(:,ibasic);        %DETERMINE BASIS MATRIX
xbasic=B\b;           %SOLVE FOR CURRENT BASIC VARIABLE VALUES
z=c(ibasic)*xbasic;   %CURRENT OBJECTIVE VALUE
if m<n;               %WHEN n=n, THERE IS A UNIQUE FEASIBLE POINT

  X(ibasic)=xbasic;     %STORE BASIC VARIABLE VALUES
  Cred=c(inon) - (c(ibasic)/B)*A(:,inon); %DETERMINE NONBASIC REDUCED COST
  OPTEST=1;             %TEST FOR OPTIMALITY/UNBOUNDEDNESS
  loop =1;
  while loop ==1;       %LOOP WHILE NONBASIC REDUCED COST HAS POSITIVE ENTRY
     if max(Cred) > 0;    
        iter=iter + 1;  %INCREMENT PHASE II ITERATION COUNTER
        [Maxcost,j]=max(Cred); %FIND VALUE & INDEX OF MAX NONBASIC REDUCED COST
        ienter=inon(j); %ENTERING BASIC INDEX
        PCOL=B\A(:,ienter);         %ENTERING COLUMN FROM CURRENT TABLEAU
        if PCOL <= 0 , OPTEST = 0;  %TEST FOR UNBOUNDEDNESS
           loop = 0;                    
        else                        %CONTINUE IF UNBOUNDEDNESS NOT DETECTED      
           J(ienter)=1;             %UPDATE BASIC VARIABLE SELECTION VECTOR
           TESTROWS=find(PCOL > 0); %FIND POSITIVE ENTRIES IN ENTERING COL
           TESTCOL=PCOL(TESTROWS);  %SUBCOL WITH POSITIVE ENTRIES
           [minrat,j]=min(xbasic(TESTROWS)./TESTCOL);%MINIMUM RATIO TEST
           if minrat <=0, cycle = cycle+1;%IF MINRATIO=0 INCREMENT CYCLE COUNTER
              if cycle > m; %IF CYCLE COUNTER EXCEEDS M, PROBABLE INFINITE CYCLE
                 disp('Algorithm terminated due to excessive cycling.')
                 disp('Restart algorithm from phase II using a perturbed')
                 disp(' RHS vector b and the current basis.')
                 disp(ibasic)
                 CYCTEST = 1;        %SET CYCLE TEST FLAG
                 break
              end
           else
               cycle = 0;  %IF MINRAT > 0 , RESET CYCLING COUNTER
           end
           iexit=ibasic(TESTROWS(j));    %EXITING BASIC INDEX
           J(iexit)=0;                   %UPDATE BASIC SELECTION VECTOR
           xbasic=xbasic - minrat*PCOL;  %UPDATE BASIC VARIABLE VALUES
           X(ibasic)=xbasic;             %STORE CURRENT BASIC VALUES  
           X(ienter)=minrat;             %STORE ENTERING BASIC VALUE
           X(iexit)=0;                   %SET EXITING BASIC VALUE TO ZERO 
           z=z + Maxcost*minrat;         %COMPUTE CURRENT OBJECTIVE VALUE
           J=logical(J);
           ibasic=K(J);                  %UPDATE BASIC INDICES
           inon=K(~J);                   %UPDATE NONBASIC INDICES 
           B=A(:,ibasic);                %UPDATE BASIS MATRIX    
           xbasic=X(ibasic)';            %UPDATE BASIC VARIABLES VECTOR
           Cred=c(inon) - (c(ibasic)/B)*A(:,inon);%COMPUTE NONBASIC RED.COST
        end
     else
      loop = 0;
     end
  end
end