【数学建模】物资调度的线性约束条件优化模型

1. 问题描述

某地区各县市的分布情况如图所示，图中的边表示连接各县市的公路，边权值表示车辆通过该路段所需的时间。

现在 D1, D2, D3 三地突发自然灾害，每天需要应急物资分别为 100 吨、80 吨和 60 吨；而能提供应急物资的县市为 A1, A2, …, A12，它们每天能提供 (30, 15, 15, 20, 35, 40, 30, 20, 10, 25, 25, 30) 吨的应急物资。
如果 A1, A2, …, A12 各地所提供的应急物资，其成本分别为 (2, 4, 5, 3, 1, 6, 2, 3, 6, 1, 5, 5)
万元/ 吨，请设计一份经济快速的应急物资运送方案，使得每天都能保证 D1, D2, D3 三地对应急物资的需求。

2. 基本假设与符号约定

为了简化问题和方便讨论，除问题中给出的假设外，我们进一步做如下的假设和说明：
(1)假设所有物资对于灾区人民的价值是相同的。
在此，我们也约定文中所用符号如下：

运输物资的吨数：

运输吨数	28($A_1$)	34($A_2$)	21($A_3$)	30($A_4$)	17($A_5$)	25($A_6$)	39($A_7$)	40($A_8$)	7($A_9$)	14($A_{10}$)	10($A_{11}$)	3($A_{12}$)
$D_1$	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$	$x_8$	$x_9$	$x_{10}$	$x_{11}$	$x_{12}$
$D_2$	$x_{13}$	$x_{14}$	$x_{15}$	$x_{16}$	$x_{17}$	$x_{18}$	$x_{19}$	$x_{20}$	$x_{21}$	$x_{22}$	$x_{23}$	$x_{24}$
$D_3$	$x_{25}$	$x_{26}$	$x_{27}$	$x_{28}$	$x_{29}$	$x_{30}$	$x_{31}$	$x_{32}$	$x_{33}$	$x_{34}$	$x_{35}$	$x_{36}$

3. 问题的分析与模型的建立

• 问题的分析

显然，给定的问题是一个优化问题，目标是经济快速的应急物资运送方案。分析题意可知，运送的总代价是时间代价T和成本代价W的和，而时间的代价是由每个物品运达的时间总和决定的，也就是县市到受灾地t运输的吨数x与所花费的时间的乘积x*t。成本代价就是运输吨数x与每吨的成本w的乘积x*w。这两部分相互制约：时间代价越低，成本代价可能就会越高，反之依然。因此，原问题要求在时间代价和成本代价之间找到一种平衡，使得总代价最小。
我们通过函数F(x)=aT(x)+bW(x)（a+b=1,a0,b0）来对总代价也就是经济快速进行量化，结合实际情况对a,b参数的值进行调整，a的值越大越侧重于节约时间，b值越大越侧重于节约成本。
需要注意的是，因A12提供物资不能于当日到达,所以第一天需要单独考虑。

（1）节点重编码
对题目所给路径模型各节点进行编码，编码后网络图如下：

（2）求最短时间
通过floyd算法(源程序见附录1)求解提供应急物资的县市对应结点到达突发自然灾害的地区对应结点的最短时间，所得结果如下表格所示：

最短时间	28($A_1$)	34($A_2$)	21($A_3$)	30($A_4$)	17($A_5$)	25($A_6$)	39($A_7$)	40($A_8$)	7($A_9$)	14($A_{10}$)	10($A_{11}$)	3($A_{12}$)
$D_1$	7	8	4	6	11	12.5	17.5	11.5	16	22	22.5	28
$D_2$	11	12	4	6	3	7	16	13	13.5	19.5	20	25.5
$D_3$	7	4	8	2	11	7	11	5	20	21	16	32

（3）数学模型建立

$$令：x=(x_1,x_2...x_{36})$$
$$总代价：F(x)=aT(x)+bW(x)$$
$$时间代价：T(x)=\sum_{i=1}^{36}x_i*t(i)$$
$$成本代价：W(x)=\sum_{i=1}^{12}(x_1+x_{i+12}+x_{i+24})*w(i)$$
等式约束条件：
$$\begin{cases} x_1+x_2...+x_{12}=100\\ x_{13}+x_{13}...+x_{24}=100\\ x_{25}+x_{24}...+x_{36}=100 \end{cases}$$
不等式约束条件：
$$\begin{cases} x_{1}+x_{13}+x_{25}\le 30\\ x_{2}+x_{14}+x_{26}\le 15\\ x_{3}+x_{15}+x_{27}\le 15\\ x_{4}+x_{16}+x_{28}\le 20\\ x_{5}+x_{17}+x_{29}\le 35\\ x_{6}+x_{18}+x_{30}\le 40\\ x_{7}+x_{19}+x_{31}\le 30\\ x_{8}+x_{20}+x_{32}\le 20\\ x_{9}+x_{21}+x_{33}\le 10\\ x_{10}+x_{22}+x_{34}\le 25\\ x_{11}+x_{23}+x_{35}\le 25\\ x_{12}+x_{24}+x_{36}\le 30\\ \end{cases}$$

（4）编程求解结果
a=1,b=0时：

运输吨数	28($A_1$)	34($A_2$)	21($A_3$)	30($A_4$)	17($A_5$)	25($A_6$)	39($A_7$)	40($A_8$)	7($A_9$)	14($A_{10}$)	10($A_{11}$)	3($A_{12}$)
$D_1$	30	15	15	20	0	0	0	0	7.9986	12.0014	0	0
$D_2$	0	0	0	0	35	40	0	0	2.0014	2.9986	0	0
$D_3$	0	0	0	0	0	0	30	20	0	0	10	0

F(x)=1962.5

a=0,b=1时：

运输吨数	28($A_1$)	34($A_2$)	21($A_3$)	30($A_4$)	17($A_5$)	25($A_6$)	39($A_7$)	40($A_8$)	7($A_9$)	14($A_{10}$)	10($A_{11}$)	3($A_{12}$)
$D_1$	12.1789	6.4604	5.7443	8.3709	14.1984	0	12.1789	8.3709	0	10.4323	9.7661	12.2990
$D_2$	10.0588	5.0131	4.4851	6.6978	11.7359	0	10.0588	6.6978	0	8.3688	7.6832	9.2007
$D_3$	7.7623	3.5265	3.2069	4.9314	9.0657	0	7.7623	4.9314	0	6.1989	5.7452	6.8694

F(x)=685

a=0.5,b=0.5时：

运输吨数	28($A_1$)	34($A_2$)	21($A_3$)	30($A_4$)	17($A_5$)	25($A_6$)	39($A_7$)	40($A_8$)	7($A_9$)	14($A_{10}$)	10($A_{11}$)	3($A_{12}$)
$D_1$	30	10.3667	15	14.6333	0	0	0	0	7.8151	22.1849	0	0
$D_2$	0	0	0	0	35	40	0	0	2.1849	2.8151	0	0
$D_3$	0	4.6333	0	5.3667	0	0	30	20	0	0	0	0

F(x)=1358.8

a	b	总代价F(x)
1	0	1962.5
0	1	710
0.5	0.5	1358.8

（5）总结分析
上图表格分别展示了侧重于考虑时间，侧重于考虑成本，以及综合考虑时间与成本的分配方案。其中a=0.5,b=0.5仅为综合考虑时间与成本时的一种情况，要获得具体的分配方案，则需依据实际情况对a，b参数进行赋值。
通过表格1可以观察到A12提供应急物资到三个受灾地区所需时间均大于24小时，其余地区都可在当天内运达。因A12提供物资不能于当日到达，A12所提供物资将计入次日提供物资量。由于灾害发生后第一天内不能接收到A12的物资，所以第一天需单独考虑。
因此计算第一天分配方案时，则需要除去A12所对应的三个变量，得到部分总代价值如下表:

a	b	F(x)
1	0	1962.5
0	1	710
0.5	0.5	1358.8

floyd:

#include<cstdio>
#include<iostream>
#include<cstdlib>
using namespace std;
const int INF=99999999;
const int maxn=100;
float dis[100][100];
int n,m;
void Init()
{
      for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
            {
                 if(i==j)dis[i][j]=0;
                 else dis[i][j]=INF;
            }
}
void Floyd()
{
      for(int k=1;k<=n;k++)
           for(int i=1;i<=n;i++)
              for(int j=1;j<=n;j++)
                   if(dis[i][j]>dis[i][k]+dis[k][j])
                      dis[i][j]=dis[i][k]+dis[k][j];
}
int main()
{
      scanf("%d%d",&n,&m);
      Init();
      int a,b;
      float c;
      for(int i=1;i<=m;i++)
      {
            scanf("%d%d%f",&a,&b,&c);
            dis[a][b]=c;
            dis[b][a]=c;
      }
      Floyd();
      cout<<"check:"<<endl;
      while(1)
      {
           scanf("%d%d",&n,&m);
           if(n==0 && m==0)break;
           cout<<dis[n][m]<<endl;
      }
}

MATLAB求解线性规划模型:

   a=1;
   b=0;
   f=[7*a+2*b;   8*a+4*b;   4*a+5*b;   6*a+3*b;   11*a+1*b;   12.5*a+6*b;   17.5*a+2*b;      11.5*a+3*b;   16*a+6*b;   22*a+1*b;   22.5*a+5*b;   
     11*a+2*b;  12*a+4*b;   4*a+5*b;   6*a+3*b;    3*a+1*b;      7*a+6*b;     16*a+2*b;     13*a+3*b; 13.5*a+6*b; 19.5*a+1*b;     20*a+5*b;
      7*a+2*b;   4*a+4*b;   8*a+5*b;   2*a+3*b;   11*a+1*b;      7*a+6*b;     11*a+2*b;      5*a+3*b;   20*a+6*b;   21*a+1*b;     16*a+5*b; ];
   A=[1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 
      0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 
      0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 
      0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 
      0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 
      0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 
      0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 
      0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 
      0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0  
      0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 
      0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 
      ];
   b=[30 15 15 20 35 40 30 20 10 25 25 ];
   Aeq=[1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 ];
   beq=[100 80 60];
   lb=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
   [x,fval]=linprog(f,A,b,Aeq,beq,lb,[])