MATLAB表上作业法解决运输问题

算例：

产销平衡算例：

计算结果：

产大于销算例：

计算结果：

销大于产算例：

计算结果：

代码如下：

%表上作业法求解运输问题
%要求每个产地和每个销地之间都可达（即运价不可以存在inf）

function [x,sigma]=Transport(N,out,in)
%N:运价表   out:每个产地产量(按行输入)   in:每个销地销量(按行输入)
%x:最优运输方案

sum_out=sum(out);
sum_in=sum(in);

if sum_out>sum_in          %产大于销的情况，转换为产销平衡问题
    [old_row,old_col]=size(N);
    new=zeros(old_row,1);
    N=[N,new];
    in=[in,sum_out-sum_in];
    disp("该问题产大于销,方案最后一列为虚拟销地")
elseif sum_in>sum_out     %销大于产的情况，转换为产销平衡问题
    [old_row1,old_col1]=size(N);
    new1=zeros(1,old_col1);
    N=[N;new1];
    out=[out,sum_in-sum_out];
    disp("该问题销大于产,方案最后一行为虚拟产地")
else
    disp("该问题为产销平衡问题")
end

[row,col]=size(N);
sigma=zeros(row,col);%定义检验数表
x=-ones(row,col);%定义初始运输表
ui=zeros(row,1);%定义ui
vj=zeros(1,col);%定义vj
out1=out;    %求初始方案的运算过程的产量out1
in1=in;       %求初始方案的运算过程的销量in1

%求初始方案(西北角法)
for i=1:row
    for j=1:col
        if out1(i)==0
            out1(i)=-1;%若该产地产量为0，则置为-1，即划掉该行
        end
        if in1(j)==0
            in1(j)=-1; %若该销地销量为0，则置为-1，即划掉该列
        end
        if out1(i)>0&&in1(j)>0
          if out1(i)>in1(j)        
            x(i,j)=in1(j);           %若产量大于销量，则该销地销量填入方案
            out1(i)=out1(i)-x(i,j);  %该产地剩余的产量
            in1(j)=in1(j)-x(i,j);    %该销地剩余的销量
          else                       
            x(i,j)=out1(i);          %若销量大于产量，则直接将产地产量填入方案
            in1(j)=in1(j)-x(i,j);    %该产地剩余的产量
            out1(i)=out1(i)-x(i,j);   %该销地剩余的销量
            break;
          end
        end
    end
end

%迭代过程
while 1
%求检验数(位势法)
for i=2:row  %初始化ui，除了第一个元素为0，其他所有ui元素变为inf
    ui(i)=inf;
end
for i=1:col  %初始化vj,所有vj元素变为inf
   vj(i)=inf; 
end
while 1    %求ui,vj
    checku=find(ui==inf);
    checkv=find(vj==inf);
    if isempty(checku)&&isempty(checkv)  %当ui,vj都不为inf时,ui和vj计算完成,跳出循环
        break;
    end
  for i=1:row
    for j=1:col
       if x(i,j)~=-1
          if ui(i)==inf&&vj(j)==inf%若ui,vj全为inf,则先跳过计算
              continue;
          elseif vj(j)==inf
             vj(j)=N(i,j)-ui(i);%若vj为inf,ui不为inf,计算vj
          else
              ui(i)=N(i,j)-vj(j);%若ui为inf,vj不为inf,计算ui
          end
       end
    end    
  end
end
for i=1:row  %初始化检验数表，所有元素变为inf
   for j=1:col
      sigma(i,j)=inf; 
   end
end
for i=1:row  %计算检验数表
    for j=1:col
        if x(i,j)==-1
            sigma(i,j)=N(i,j)-ui(i)-vj(j);          
        end
    end
end


%判断是否得到最优方案
if sigma>0   
    disp("有唯一最优方案,最优方案为(表中-1表示空格):")
%  可将最优方案中的-1点变为0,更美观,若执行此操作,还应删去计算最小运价过程中的判断条件
%  当最优方案中存在运量为0时,此操作会影响对方案的阅读
%     for i=1:row   
%        for j=1:col    
%           if x(i,j)==-1
%               x(i,j)=0;
%           end
%        end
%     end
    x
    sum_min=0;
    for i=1:row   %计算最小运价
       for j=1:col
           if x(i,j)~=-1
              sum_min=sum_min+x(i,j)*N(i,j);
           end
       end
    end
    disp("最小运价为：")
    sum_min
     break;
elseif sigma>=0
    disp("最优方案不唯一,其中一个为(表中-1表示空格):")
%   可将最终结果中的-1点变为0,更美观,若执行此操作,还应删去计算最小运价过程中的判断条件
%   当最优方案中存在运量为0时,此操作会影响对方案的阅读
%      for i=1:row  
%        for j=1:col     
%           if x(i,j)==-1
%               x(i,j)=0;
%           end
%        end
%      end
     x
     sum_min=0;
    for i=1:row  %计算最小运价
       for j=1:col
          if x(i,j)~=-1
              sum_min=sum_min+x(i,j)*N(i,j);
          end
       end
    end 
    disp("最小运价为：")
    sum_min
    break;
end

 
%闭回路调整法
visit=x;    
for i=1:row        %初始化访问表,可以被访问的点标为0
    for j=1:col
        if visit(i,j)~=-1
            visit(i,j)=0;  
        end
    end
end

m=min(sigma(sigma<0));%找到小于零的最小检验数m
[r2,c2]=find(sigma==m);%找到m的位置

%记录m的行标和列标,由于可能出现检验数相同的情况,我们取其中第一个
r=r2(1);
c=c2(1);
r1=r2(1);          
c1=c2(1); 

visit(r,c)=2;  %标记m已被访问,记为2   
circle=-ones(row+col+1,3);  %定义闭回路路径表
%circle表的结构:[行标  列标  运量]

%将起点(m点)填入路径表
circle(1,1)=r1;
circle(1,2)=c1;
p=2;


%找闭回路
while 1  
   for i=1:row         %找第c列中有无未被访问的点
      if visit(i,c)==0  %若该点未被访问,则访问该点,进行标记并存入路径表
         visit(i,c)=1;
         r=i;           %记录该点行标
         circle(p,1)=i; 
         circle(p,2)=c;
         circle(p,3)=x(i,c);
         p=p+1;
          break;
      end
   end
   
   for j=1:col         %找第r行中有无未被访问的点
      if visit(r,j)==0 %若该点未被访问,则访问该点,进行标记并存入路径表
         visit(r,j)=1;
         c=j;          %记录该点列标
         circle(p,1)=r;
         circle(p,2)=j;
         circle(p,3)=x(r,j);
         p=p+1;
          break;
      end
   end
   
   a=find(visit(r,:)==0);
   b=find(visit(:,c)==0);
   a1=find(visit(r,:)==2);
   b1=find(visit(:,c)==2);
   if isempty(a)&&isempty(b)  %判断该点所在行和列中有无未被访问的点
       if ~isempty(a1)||~isempty(b1)  %判断最后访问点是否与起始点在同一行或同一列，若是则跳出循环，找到闭回路
           break;
       else           %若不是,则将该点置为-1(此路不通,下次循环不走此路),重置访问表和路径表,开始下一次循环
           visit(r,c)=-1;
           for i=1:row
              for j=1:col
                  if visit(i,j)==1
                     visit(i,j)=0;
                  end
              end
           end
       
      r=r1;
      c=c1;
    
      circle=-ones(row+col,3);
      circle(1,1)=r1;
      circle(1,2)=c1;
       p=2;
       end
   end   
end

[rows,cols]=size(circle);
%定义circle表时我们给了足够大的维度,而由于闭回路可能无法包括所有点,现在要将多余行删去
for i=1:rows   
    if circle(i,1)==-1
        break;
    end
end
circle(i:rows,:)=[];

%开始找闭回路中的顶点表
add=circle(1,:);
circle=[circle;add];  %将起始点填入circle表,形成完整的闭回路
[row1,col1]=size(circle);
i=1;

%若闭回路中同一行或同一列中有两个以上的点,则删去中间点,只留下顶点
while 1             
     if i==row1-1
       break; 
     end
    i=i+1;
    if (circle(i-1,1)==circle(i+1,1))
        circle(i,:)=[];
        row1=row1-1;
        i=i-1;
    end
    if (circle(i-1,2)==circle(i+1,2))
        circle(i,:)=[];
        row1=row1-1; 
        i=i-1;
    end 
end

%由于起始点在闭回路中出现了两次,因此删去第二次出现的起始点
k=find(circle(:,3)==-1);
circle(k(2),:)=[];

[row2,col2]=size(circle);

%将顶点中奇数编号和偶数编号分开
%顶点表的结构:[行标  列标  运量]
if mod(row2,2)==0    %若顶点总数是偶数
    single=zeros(row2/2,3);  %定义奇数编号顶点表
    double=zeros(row2/2,3);  %定义偶数编号顶点表
end
if mod(row2,2)==1    %若顶点总数是奇数
   single=zeros(row2/2+0.5,3);   %定义奇数编号顶点表
   double=zeros(row2/2-0.5,3);   %定义偶数编号顶点表
end

j1=1;
j2=1;
%将闭回路中的点按编号分别存入奇数顶点表和偶数顶点表
for i=1:row2
   if mod(i,2)==1
       single(j1,:)=circle(i,:);
       j1=j1+1;
   end
   if mod(i,2)==0
      double(j2,:)=circle(i,:);
      j2=j2+1;
   end
end


%更新运量表
value=double(:,3);
[min_x,index]=min(value(value>=0));%找到偶数顶点的运量最小值及其位置
x(r1,c1)=min_x;       %把该运量最小值填入闭回路起始点的运量
x(double(index,1),double(index,2))=-1;  %把该偶数顶点的运量标记为-1
double(index,:)=[];    %将该偶数顶点从偶数顶点表中删去,以免影响后续计算

[row3,col3]=size(single);
[row4,col4]=size(double);
%将奇数编号顶点的运量加上min_x
for i=2:row3
   x(single(i,1),single(i,2))=single(i,3)+min_x;
end
%将偶数编号顶点的运量减去min_x
for i=1:row4
    x(double(i,1),double(i,2))=double(i,3)-min_x;
end
end