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

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

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

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

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

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时会报错，说明无可行解：

 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，得到结果：

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

取值为1说明建立运输关系。