プログラミングで学ぶ線形計画法（LP)

% --------------------------
% LP12.m 最適基底解を求める（２段階単体法の２段階目、MATLAB_FreeMat）
global A Binv b c xi xn c0 sw
A=csvread('A.csv'); % Excelのcsvファイルから読み込む
b=csvread('B.csv'); c=csvread('C.csv')
a1=diag(ones(1,size(A,1))); % 単位行列を準備する
A=[A,a1] % 標準型に書き直した場合の係数行列
c=[c,zeros(1,size(A,1))]; c=c'; b=b'
x=zeros(size(A,2),1); sw=0;
xn=[1,2] % xnは非基底の番号を入れる場所
for rep=1:999 % 改善されなくなるとreturnで終了する
c0=c; B=A; % Bを以下で基底行列にする
% 非基底変数に対応する係数を削除（降順に削除する）
xi=1:size(A,2); % xiは基底の番号を入れる場所
for j=size(xn,2):-1:1 % 降順、xn(j)は非基底の番号
B(:,xn(j))=[]; xi(xn(j))=[]; c0(xn(j))=[];
end
Binv=B^(-1); % 基底行列の逆行列
Kaisu=rep % 回数を表示する
BasicEx; % 基底変数の入れ替えfunctionへ（BasicEx.m）
if sw==1; return; end
end

% --------------------------------------------
% function BasicEx
function BasicEx
global A Binv b c xi xn c0 sw
RC=[];
for k=1:size(xn,2)
RC(k)=c0'*Binv*A(:,xn(k))-c(xn(k)); % リデューストコスト
end
ReducedCost=RC % リデューストコストの表示
[r,i]=min(RC); KiteiNiIreru=xn(i) % 基底に入れる変数の番号
if r>=0 % 最適な基底変数が選ばれている時
sw=1; SaitekiKitei=999999
return;
end
xi=[xi xn(i)]; c0=[c0; c(xn(i))]; xn(i)=[];
beta=Binv*b; alpha=Binv*A(:,KiteiNiIreru);
for j=size(alpha,1):-1:1
if alpha(j)<=0 alpha(j)=1; beta(j)=1e+10; end
end
[theta,ss]=min(beta./alpha); % 基底変数の中でss番目を非基底変数にする
% 非基底変数x(Ireru)を基底に入れてΘまで増やせる
c0(ss)=[]; xn=[xn xi(ss)]; KiteiKaraDasu=xi(ss)
xi(ss)=[];
xx=Binv*(b-theta*A(:,KiteiNiIreru)); % 新基底解β-Θ*α
x=zeros(size(A,2),1);
for j=1:size(xx,1)
x(xi(j))=xx(j); % 基底も非基底もxに入れる
x(KiteiNiIreru)=theta;
end
JikkoKanoKai=x' % xを表示
MokutekiKansuTi=c'*x % 目的関数値
xn=sort(xn); xi=sort(xi); % 番号順に並べなおす

% 逆行列のピボット演算、縮退の処理はやらない