运筹学 matlab 表上作业法 运输问题
文章目录
一，实验内容
二，建立运输问题的数学模型
根据题意是求预期盈利最大的采购方案，而我们所学的运输问题模型是已知运价求最小运价的运输方案，则我们需要将表中的利润转换为运价，需满足利润越高的则运价越低，可以定义单位运价表中的运价为表中最大的利润减去原来的利润，满足利润越高则运价越低。
单位运价表为：
| A | B | C | D | 产量 | |
|---|---|---|---|---|---|
| 城市1 | 0 | 5 | 4 | 3 | 2500 |
| 城市2 | 2 | 8 | 3 | 4 | 2500 |
| 城市3 | 1 | 7 | 6 | 2 | 5000 |
| 销量 | 1500 | 2000 | 3000 | 3500 |
m i n z = 0 x 11 + 5 x 12 + 4 x 13 + 3 x 14 + 2 x 21 + 8 x 22 + 3 x 23 + 4 x 24 + 1 x 31 + 7 x 32 + 6 x 33 + 2 x 34 minz=0x_{11}+5x_{12}+4x_{13}+3x_{14}+2x_{21}+8x_{22}+3x_{23}+4x_{24}+1x_{31}+7x_{32}+6x_{33}+2x_{34} minz=0x11+5x12+4x13+3x14+2x21+8x22+3x23+4x24+1x31+7x32+6x33+2x34
{ x 11 + x 12 + x 13 + x 14 = 2500 x 21 + x 22 + x 23 + x 24 = 2500 x 31 + x 32 + x 33 + x 34 = 5000 x 11 + x 21 + x 31 = 1500 x 12 + x 22 + x 32 = 2000 x 13 + x 23 + x 33 = 3000 x 14 + x 24 + x 34 = 3500 x 11 , x 12 , x 13 , x 14 , x 21 , x 22 , x 23 , x 24 , x 31 , x 32 , x 33 , x 34 ≥ 0 \left\{\begin{matrix} x_{11}+x_{12}+x_{13}+x_{14}&= 2500\\ x_{21}+x_{22}+x_{23}+x_{24}&= 2500\\ x_{31}+x_{32}+x_{33}+x_{34}&= 5000\\ x_{11}+x_{21}+x_{31}&= 1500\\ x_{12}+x_{22}+x_{32}&= 2000\\ x_{13}+x_{23}+x_{33}&= 3000\\ x_{14}+x_{24}+x_{34}&= 3500\\ x_{11},x_{12},x_{13},x_{14},x_{21},x_{22},x_{23},x_{24},x_{31},x_{32},x_{33},x_{34} &\ge 0 \end{matrix}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧x11+x12+x13+x14x21+x22+x23+x24x31+x32+x33+x34x11+x21+x31x12+x22+x32x13+x23+x33x14+x24+x34x11,x12,x13,x14,x21,x22,x23,x24,x31,x32,x33,x34=2500=2500=5000=1500=2000=3000=3500≥0
三，实验方法与步骤
1，定义输入，利用表上作业法求解
function [sigma]=Transport(N,out,in)
%N:运价表 out:每个产地产量(按行输入) in:每个销地销量(按行输入)
%x:最优运输方案
2，判断是何种运输问题，即是否产销平衡
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
3，求出初始基可行解
%求初始方案(西北角法)
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
4，用位势法求检验数
%求检验数(位势法)
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
5，判断是否是最优解（无负检验数）
if sigma>0
disp("有唯一最优方案,最优方案为(表中-1表示空格):")
6，闭回路调整法改进
%找闭回路
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
7，完整matlab实现
function [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表示空格):")
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表示空格):")
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
8，matlab函数输入
N=[ 0 5 4 3
2 8 3 4
1 7 6 2];
in=[1500 2000 3000 3500];
out=[2500 2500 5000];
[sigma]=Transport(N,out,in)
四，实验结果
结果为：
该问题为产销平衡问题
有唯一最优方案,最优方案为(表中-1表示空格):
x =
0 2000 500 -1
-1 -1 2500 -1
1500 -1 -1 3500
最小运价为：
sum_min = 28000
sigma =
Inf Inf Inf 2
3 4 Inf 4
Inf 1 1 Inf
以下是matlab截图：
![[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-188s5MLP-1648973333417)(C:\Users\86150\AppData\Roaming\Typora\typora-user-images\image-20220403155217361.png)]](https://i-blog.csdnimg.cn/blog_migrate/1821aa0bc9005e04a78b5b37594d8bab.png)
转换为原问题：
最优解为：(预期赢利最大的采购方案)
| A | B | C | D | |
|---|---|---|---|---|
| 城市1 | 2000 | 500 | ||
| 城市2 | 2500 | |||
| 城市3 | 1500 | 3500 |
最大赢利为：10*(2500+2500+5000)-28000=72000(元)
五，实验结果分析
经过计算，实验结果正确，并使用多个样例测试程序皆有正确结果。
六，附录：部分其余样例
%P94运筹学第四版（清华大学出版社）
N=[ 3 11 3 10
1 9 2 8
7 4 10 5];
in=[3 6 5 6];
out=[7 4 9];
[sigma]=Transport(N,out,in)
扫一扫
2012
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
