標準形の線形計画問題
minimize z = c'x
subject to Ax = b, x >= 0
を解くためのOctave関数を紹介します．
octave:1> A = [2 5 3 -1 0 0;
> 3 2.5 8 0 -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,basis] = 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
basis =
3 1 5
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
endfunction
simplex.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;
endfunction
pivot.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 0 -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