運輸問題中表上做業法實現

Step1：初始基可行解肯定

方法1：西北角法

西北角法函數文件

function [B,X]=XiBeiJiao(demand,supply,cost)
    [m,n] = size(cost)
    B = zeros(m,n);X=zeros(m,n)
    B(1,1) = 1;
    i = 1,j = 1;
    while((i<=m)&&(j<=n))
        if  supply(i) < demand(j)
            X(i,j)=supply(i)
            demand(j) = demand(j)-supply(i)
            supply(i)=0 
            B(i,j) = 1 
            
            i = i + 1   
        else supply(i) > demand(j) 
            X(i,j)=demand(j)
            supply(i) = supply(i)-demand(j) 
            demand(j)=0 
            B(i,j) = 1
           
            j = j + 1     
        end
    end
    disp("初始基所在位置")
    disp(B)
    disp("初始調運方案")
    disp(X)
end

　　

實例1運行代碼

demand=[3;6;5;6]
supply=[7;4;9]
cost=[3,11,3,10;1,9,2,8;7,4,10,5]
XiBeiJiao(demand,supply,cost)

實例1運行結果

初始基所在位置
     1     1     0     0
     0     1     1     0
     0     0     1     1
​
初始調運方案
     3     4     0     0
     0     2     2     0
     0     0     3     6

實例2運行代碼

clc;clear all
demand=[3;2;3;2]
supply=[5;2;3]
cost= [3 7 6 4;2 4 3 2;4 3 8 5]
XiBeiJiao(demand,supply,cost)

實例2運行結果

初始基所在位置
     1     1     1     0
     0     0     1     0
     0     0     1     1
​
初始調運方案
     3     2     0     0
     0     0     2     0
     0     0     1     2

方法2：最小元素法

最小元素法函數文件

function [B,X]=ZuiXiaoYuanSu(demand,supply,cost)
%本程序經過最小元素法求解運輸問題的初始解
%B爲初始基可行解所在位置，X爲初始調運方案
   [m,n]=size(cost);I=[1:m];X=zeros(m,n);B=zeros(m,n);flag=0;
   while flag==0
        [row,col] = find(cost ==min(min(cost)))%尋找最小元素
        a =row(1,:);b=col(1,:) %防止最小元素有兩個以上
        if  supply(a) < demand(b)
            demand(b) = demand(b) - supply(a);
            B(a,b) = 1;%記爲基
            X(a,b) = supply(a);%進行調運
            supply(a) = 0;
            cost(a,:) = inf*ones(1,n)%刪除此行
        elseif supply(a) > demand(b)
            supply(a) = supply(a) - demand(b);
            B(a,b) = 1;
            X(a,b) = demand(b);
            demand(b) = 0;cost(:,b) = inf*ones(m,1)%刪除此列
        else supply(a) = demand(b)%退化狀況
             B(a,b) = 1;
             X(a,b) = supply(a);
            supply(a)=0;demand(b)=0;
            cost(a,:)=inf*ones(1,n)
             cost(:,b)=inf*ones(m,1)
             %補0
             if length(find(B==1)) < m+n-2
                I1 =find(I~=a)
                B(I1(end),b)=1
             end         
        end  
        if length(find(B==1))==m+n-1
            flag=1
        else
            flag =0
        end  
    end
disp("初始基所在位置")
disp(B)
disp("初始調運方案")
disp(X)
end    

實例1運行代碼

demand=[3;6;5;6]
supply=[7;4;9]
cost=[3,11,3,10;1,9,2,8;7,4,10,5]
ZuiXiaoYuanSu(demand,supply,cost)

實例1運行結果

初始基所在位置
     0     0     1     1
     1     0     1     0
     0     1     0     1
​
初始調運方案
     0     0     4     3
     3     0     1     0
     0     6     0     3

實例2運行代碼

demand=[3;2;3;2]
supply=[5;2;3]
cost= [3 7 6 4;2 4 3 2;4 3 8 5]
ZuiXiaoYuanSu(demand,supply,cost)

 

實例2運行結果

初始基所在位置
     1     0     1     1
     1     0     0     0
     0     1     1     0
初始調運方案
     1     0     2     2
     2     0     0     0
     0     2     1     0

方法3：元素差額法

元素差額法函數文件

function [C,X]=vogel(demand,supply,cost)
%本程序經過元素差額法求解運輸問題的初始解
%C爲初始基可行解所在位置，X爲初始調運方案
[m,n]=size(cost)
X=zeros(m,n);C=zeros(m,n);I=[1:m];flag=0;
while flag==0
    A=sort(cost,2,'ascend') 
    zd1=A(:,1)%行最小值
    cd1 = A(:,2)%行次小值
    h = cd1-zd1%行差額
    B=sort(cost,'ascend')
    zd2=B(1,:)';cd2=B(2,:)'
    l = cd2-zd2%列差額
    [z1,row]=max(h);[z2,col]=max(l);z=max([z1,z2])
    if z==z1
        [~,col]= min(cost(row,:))%此行中最小值，記錄列數
    else z==z2
        [~,row] = min(cost(:,col))%此列中最小值,記錄行數
    end        
    if demand(col)>supply(row)
       X(row,col)=supply(row);C(row,col)=1; 
       demand(col) = demand(col)-supply(row);supply(row) = 0;
       cost(row,:) = inf*ones(1,n)
    elseif demand(col) <supply(row)
         X(row,col)=demand(col);C(row,col)=1
         supply(row) = supply(row)-demand(col);demand(col) = 0;
         cost(:,col) = inf*ones(m,1)
    else demand(col) = supply(row)%退化
           C(row,col) = 1;X(row,col) = supply(row);
           supply(row)=0;demand(col)=0;
           cost(row,:)=inf*ones(1,n);cost(:,col)=inf*ones(m,1)
           %下面對退化狀況補0
           if length(find(C==1)) < m+n-2
                I1 =find(I~=row)
                B(I1(end),col)=1
           end
    end    
    if length(find(C==1))==m+n-1
       flag=1
    else
        flag =0
    end  
disp("初始基所在位置")
disp(C)
disp("初始調運方案")
disp(X)    
end

實例1運行代碼

clc;clear all
demand=[3;6;5;6]
supply=[7;4;9]
cost=[3,11,3,10;1,9,2,8;7,4,10,5]
volger(demand,supply,cost)

實例1運行結果

初始基所在位置
     0     0     1     1
     1     0     0     1
     0     1     0     1
初始調運方案
     0     0     5     2
     3     0     0     1
     0     6     0     3

實例2運行代碼

clc;clear all
demand=[3;2;3;2]
supply=[5;2;3]
cost= [3 7 6 4;2 4 3 2;4 3 8 5]
volger(demand,supply,cost)

實例2運行結果

初始基所在位置
     1     0     1     1
     0     0     1     0
     0     1     0     1
​
初始調運方案
     3     0     1     1
     0     0     2     0
     0     2     0     1

　　

Step2:解的檢驗（位勢法）

%利用位勢法求解非基變量檢驗數
%此程序須要輸入初始調運方案B
clc;clear all
demand=[3;6;5;6];supply=[7;4;9];cost=[3,11,3,10;1,9,2,8;7,4,10,5]
[m,n]=size(cost)
B=[0,0,3,10;1,0,2,0;0,4,0,5]
[BI,BJ]=find(B~=0)%記錄基變量下標
b=zeros(m+n-1,1)
A=zeros(m+n-1,m+n)%初始化係數矩陣
%給初值
for i=1:m+n-1
    A(i,BI(i))=1
    A(i,BJ(i)+m)=1%生成係數矩陣 
end
A1=A(:,1:3);A2=A(:,4:end)
A=[1,0,0,0,0,0,0;A]
b = [0; cost(sub2ind([m,n],BI,BJ))];
x=A \ b
u = x(1:m);v = x(m+1:end);
[nI,nJ] = find(B==0); 
Sigma = zeros(m,n);               
for i=1:length(nI)
    Sigma(nI(i),nJ(i)) = cost(nI(i),nJ(i)) - (u(nI(i)) + v(nJ(i)));
end
disp("非基變量檢驗數")
disp(Sigma)

　輸出結果

非基變量檢驗數
    1     2     0     0
    0     1     0    -1
   10     0    12     0

非基變量檢驗數存在負數狀況，須要對調運方案進一步調整