Contents

function [x,Out] = yzsimplex(A,b,c,bas,nbs,opts)
% simplex method for standard LP: min{c'*x: A*x=b, x >= 0}
% INPUT:
%               A = constraint coefficient matrix
%               b = constraint right-hand side vector
%               c = objective coefficient vector
%             bas = initial basic indices
%             nbs = initial nonbasic indices
%            opts = parameter options structure
%               opts.maxit = maximum iterations allowed
%               opts.prule = pivot rule for choosing entering variables
%                         prule = 1: min-value (most negative) rule
%                         prule = 2: smallest subscript rule
%               opts.frep = prequency of printout
%               opts.x = initial basic feasible solution
%
% OUTPUT:
%                 x = final basic feasible solution if it exists
%             Out = output options structure
%                   Out.info = character string containing status info:
%                       info(1) = 'm' (e.g., info = 'maxiter reached')
%                       info(1) = 'u' (e.g., info = 'unbounded LP')
%                       info(1) = 'o' (e.g., info = 'optimum obtained')
%                   1st letter must be the above 3 (lower or upper case)
%                   Out.iter  = iterations taken
%                   Out.bas = final basic indices
%                   Out.nbs = final nonbasic indices

[m,n] = size(A);
if isfield(opts,'maxit'); maxit = opts.maxit; else maxit = m*n; end
if isfield(opts,'prule'); prule = opts.prule; else prule = 1; end
if isfield(opts,'freq');  freq  = opts.freq;  else freq = inf; end
if isfield(opts,'x'); x = opts.x; else x = []; end
info = 'maxiter reached';

for iter = 1:maxit

Extract and factorize basic submatrix

A_B = A(:,bas);
[L, U] = lu(A_B);
if isempty(x);
    x = zeros(n,1);
    x(bas) = U \ (L \ b);
    x(nbs) = 0;
end

computer y and sN

%y = A_B' \ c_B and sN = c_N - A_N'*y
y = L' \ (U' \ c(bas));
sN = c(nbs) - A(:,nbs)'*y;
if ~mod(iter,freq)
    fprintf('Iter %4i:  Obj = %.8e',iter,c'*x);
    fprintf('   basis: '); fprintf('%d,', bas(1:3));
    fprintf('\n')
end

check optimality

if all(sN >= -eps*100)
   if mod(iter,freq)
        fprintf('Iter %4i:  Obj = %.8e\n',iter,c'*x);
   end
   info = 'optimum obtained';
   break;
end

find entering variable local/global indices

[ie_local,ie_global] = get_enter(sN,nbs,prule);

compute direction and check unboundedness

%d = A_B \ A(:,ie_global);
d = U \ (L \ A(:,ie_global));
if all(d <= eps)
    info = 'unbounded'; break;
end

find leaving variable local/global indices and step

[il_local,il_global,step] = get_leave(x(bas),d,bas);

update x

x(bas) = x(bas) - step*d;
x(ie_global) = step;
x(il_global) = 0;

update partition

bas(il_local) = ie_global;
nbs(ie_local) = il_global;

end

if freq; fprintf('%s\n',info); end
Out.info   = info;
Out.iter   = iter;
Out.basis  = bas;
Out.nonbas = nbs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ie_local,ie_global] = get_enter(sN, nbs, prule)
% Find entering variable's local and global indices
% INPUT:
%           sN  = sN vector in lecture notes
%           nbs = indices for the nonbasics
%         prule = pivot rule
% OUTPUT:
%            ie_local  = local index for entering variable in bas
%            ie_global = global index for entering variable in x

switch prule
    case 1; % min-value rule
        [~,ie_local] = min(sN);
    case 2  % smallest-subscript rule
        I_neg = find(sN < 0);
        [~, imin] = min(nbs(I_neg));
        ie_local = I_neg(imin);
    otherwise
        error('pivot rule must be 1 or 2');
end

% find global index
ie_global = nbs(ie_local);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [il_local,il_global,step] = get_leave(xB,d,bas)
% performing ratio test to determine leaving variable
% local and global indices, and the step length
%
% INPUT:
%            xB  = basic variables
%            d   = direction vector
%            bas = basic indices
% OUTPUT:
%            il_local  = local index for leaving variable in bas
%            il_global = global index for leaving variable in x
%            step = step length (how large entering variable can be)

tol = eps;                          % tolerance
ipos = find(d > tol);               % find indices for which d > 0
ratio = xB(ipos)./d(ipos);          % compute ratios for those d > 0
[step, imin] = min(ratio);          % find the min ratio and index
il_local = ipos(imin);              % leaving variable local index
il_global = bas(il_local);          % leaving variable blobal index

Published with MATLAB® R2013b