We show two Octave functions for solving the standard linear programming problem:
minimize z = c'x
subject to Ax = b, x >= 0
octave:1> A = [2 5 3 -1 0 0; > 3 2.5 8 -1 0; > 8 10 4 0 0 -1] A = 2.00000 5.00000 3.00000 -1.00000 0.00000 0.00000 3.00000 2.50000 8.00000 0.00000 -1.00000 0.00000 8.00000 10.00000 4.00000 0.00000 0.00000 -1.00000 octave:2> b = [185;155;600] b = 185 155 600 octave:3> c = [4 8 3 0 0 0]' c = 4 8 3 0 0 0 octave:4> [optx,zmin,is_bounded,sol] = lp(c,A,b) optx = 66.25000 0.00000 17.50000 0.00000 183.75000 0.00000 zmin = 317.50 is_bounded = 1 sol = 1
function [optx,zmin,is_bounded,sol,basis] = lp(c,A,b) % [optx,zmin,is_bounded,sol,basis] = lp(c,A,b) % This program finds a solution of the standard linear programming problem: % minimize z = c'x % subject to Ax = b, x >= 0 % using the two phase method, where the simplex method is used at each stage. % optx: an optimal solution. % zmin: the optimal value. % is_bounded: 1 if the solution is bounded; 0 if unbounded. % sol: 1 if the problem is solvable; 0 if unsolvable. % basis: indices of the basis of the solution. [m,n] = size(A); % Phase one for i = 1:m if b(i) < 0 A(i,:) = -A(i,:); b(i) = -b(i); endif endfor if m == 1 d = -A; else d = -sum(A); endif w0 = sum(b); H = [A b;c' 0;d -w0]; % The initial simplex table of phase one index = 1:n; basis = n+1:n+m; [H,basis,is_bounded] = simplex(H,basis,index,1); if H(m+2,n+1) < -1e-10 sol = 0; disp('unsolvable') optx = []; zmin = []; is_bounded = []; else sol = 1; j = 0; for i = 1:n j = j+1; if H(m+2,j) > 1e-10 H(:,j) = []; index(j) = []; j = j-1; endif endfor H(m+2,:) = []; if length(index) > 0 % Phase two [H,basis,is_bounded] = simplex(H,basis,index,2); if is_bounded == 1 optx = zeros(n+m,1); [n1,n2] = size(H); for i = 1:m optx(basis(i)) = H(i,n2); endfor optx(n+1:n+m,1) = []; zmin = -H(n1,n2); else optx = []; zmin = -Inf; endif else optx = zeros(n+m,1); zmin = 0; endif endif endfunctionsimplex.m
function [H1,basis,is_bounded] = simplex(H,basis,index,s) % [H1,basis,is_bounded] = simplex(H,basis,index,s) % H: simplex table. % basis: the indices of basis. % index: the indices of x. % s: 1 for phase one; 2 for phase two. % H1: new simplex table. % is_bounded: 1 if the solution is bounded; 0 if unbounded. switch s case 1 s0 = 2; case 2 s0 = 1; endswitch [n1,n2] = size(H); sol = 0; while sol == 0 [fm,jp] = min(H(n1,1:n2-1)); if fm >= 0 is_bounded = 1; % bounded solution sol = 1; else [hm,ip] = max(H(1:n1-s0,jp)); if hm <= 0 is_bounded = 0; % unbounded solution sol = 1; else for i = 1:n1-s0 if H(i,jp) > 0 h1(i) = H(i,n2)/H(i,jp); else h1(i) = Inf; endif endfor [minh1,ip] = min(h1); basis(ip) = index(jp); H = pivot(H,ip,jp); endif endif endwhile H1 = H; endfunctionpivot.m
function H1 = pivot(H,ip,jp) % H1 = pivot(H,ip,jp) [n,m] = size(H); piv = H(ip,jp); if piv == 0 disp('singular') else H(ip,:) = H(ip,:)/piv; for i = 1:n if i != ip H(i,:) = H(i,:) - H(i,jp)*H(ip,:); endif endfor endif H1 = H; endfunction
octave:1> A = [2 5 3 -1 0 0; > 3 2.5 8 -1 0; > 8 10 4 0 0 -1] A = 2.00000 5.00000 3.00000 -1.00000 0.00000 0.00000 3.00000 2.50000 8.00000 0.00000 -1.00000 0.00000 8.00000 10.00000 4.00000 0.00000 0.00000 -1.00000 octave:2> b = [185;155;600] b = 185 155 600 octave:3> c = [4 8 3 0 0 0]' c = 4 8 3 0 0 0 octave:4> [xn,sn,ln,mn,kn,nxs,rp,rd] = pclp(A,b,c); octave:5> xn xn = 6.6250e+01 5.9031e-11 1.7500e+01 2.0661e-10 1.8375e+02 2.7548e-10 octave:6> zmin = c'*xn zmin = 317.50