引入0-1规划的运输学问题

122 阅读 0 评论 81 点赞

我是靠谱客的博主追寻书包，这篇文章主要介绍引入0-1规划的运输学问题，现在分享给大家，希望可以做个参考。

产销平衡的运输学问题

产销平衡的运输学问题是运筹学里面的一个重要部分。常规解法是通过表上作业法求解：首先通过最小元素法或者伏格尔法找到初始基可行解，之后用闭回路法或者位势法进行解的最优性判断，最后通过闭回路法对结果进行调整。这一方法虽然在手动计算时简化了计算量，但是在通过计算机进行计算时难以写出代码。考虑到表上作业法本质上是一种单纯形法，因此通过对matlab内置的linprog函数进行求解，代码如下：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
function x = ys(a,beq)%先输入右侧列向量，再输入下面向量
[m,n]=size(a);
aeq=zeros(m+n,m*n);
for i=1:m
    for j=1:n
        aeq(i,(i-1)*n+j)=1;
    end     
end
for i=m+1:m+n
    for j=1:m
        aeq(i,(j-1)*n+i-m)=1;
    end
end
c=zeros(1,m*n);
for i=1:m
    for j=1:n
        c(1,(i-1)*n+j)=a(i,j);
    end
end
lb=zeros(m*n,1);
ub=zeros(m*n,1);
for i=1:m*n
    ub(i,1)=Inf;
end
[b,fval]=linprog(c,[],[],aeq,beq,lb,ub);
x=zeros(m,n);
for i=1:m
    for j=1:n
       x(i,j)=b((i-1)*n+j,1);
    end
end

通过这一方法求解下面的问题：

	B1	B2	B3	B4	产量
A1	3	11	3	10	7
A2	1	9	2	8	4
A3	7	4	10	5	9
销量	3	6	5	6	20

这是我所使用的运筹学教材中的一道例题。输入情况如下：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
a=[3 11 3 10;1 9 2 8;7 4 10 5]
a =
     3    11     3    10
     1     9     2     8
     7     4    10     5
>> b=[7 4 9 3 6 5 6]
b =
     7     4     9     3     6     5     6
>> ys(a,b)

Optimal solution found.

ans =
     0     0     5     2
     3     0     0     1
     0     6     0     3

和书中给出的正确答案一致。

引入了0-1规划的运输学问题

在实际生活中，有时候会遇到这样的问题，即一个工厂最好能只给一个销地供货。在这种情况下理应允许一部分销地运量略大于实际销量，另一部分销地运量略小于实际销量。在这一问题中不具备普遍性，但是在考虑病人的就医问题等情况时，为了使病人得到及时的救治，可以令一部分医院多收治一些病人。为了方便对接，常常令一个街道办事处对接一家医院，此时这一研究便有了实际意义。为对这一问题进行控制，我规定了一个e，代表销地可以比平常多出售或少出售的商品占实际能出售商品的比例。即在产销平衡运输问题中，销量为b，那么在这一问题中，其销量为[(1-e)*b,(1+e)*b]这么一个区间，下面附上代码：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
function x = ya(a,beq,e)%先输入右侧列向量，再输入下面向量
%考虑0-1规划模型运输问题算法
[m,n]=size(a);
s=m*n;
A=zeros(n,s);
aeq=zeros(m,s);
a0=beq(1:m,1);
b0=beq(m+1:m+n,1);
b=[(1+e)*b0;-(1-e)*b0];
beq=ones(m,1);
lb=zeros(s,1);
ub=ones(s,1);
for i=1:s
    intcon(i)=i;
end
for i =1:m
    for j=1:n
        a(i,j)=a0(i,1)*a(i,j);
        aeq(i,(i-1)*n+j)=1;
    end
end
for i=1:n
    for j=1:m
        A(i,(j-1)*n+i)=a0(j,1);
        A(i+n,(j-1)*n+i)=-a0(j,1);
        f((j-1)*n+i,1)=a(j,i);
    end
end 
[x0,f1,f2]=intlinprog(f,intcon,A,b,aeq,beq,lb,ub); 
x=zeros(m,n);
for i=1:m
    for j=1:n
       x(i,j)=x0((i-1)*n+j,1);
    end
end

需要注意的是，在该代码中，运价矩阵的行数应大于等于列数。如果不满足这一条件，将矩阵转置即可。下面给出该问题的实例：
首先转置矩阵：

	A1	A2	A3	销量
B1	3	1	7	3
B2	11	9	4	6
B3	3	2	10	5
B4	10	8	5	6
产量	7	4	9	20
矩阵输入的方法与上一个代码相同。首先看到在e取0.1时会报错，说明无可行解：

复制代码

1
2
3
4
5
6
7
8
9
10
 ya(a,b,0.1)

No feasible solution found.

Intlinprog stopped because no point satisfies the constraints.

位置 1 处的索引超出数组边界。
出错 ya (line 33)
       x(i,j)=x0((i-1)*n+j,1);

接下来取e=0.5，得到结果：

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ya(a,b,0.5)
LP:                Optimal objective value is 70.500000.                                            


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0 (the
default value). The intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05 (the default value).

ans =
         0    1.0000         0
         0         0    1.0000
    1.0000         0         0
         0         0    1.0000