初学，如有问题欢迎交流~
问题详述
某公司经销汽车配件A，下设2个分别位于沈阳和郑州的加工厂。该公司每月需要把各产地生产的产品分别运往北京、上海、广州3个销售点。运输过程中，允许经天津、武汉2个中转站进行转运。
该公司在租用运输车辆时，租赁公司给出了如下的单位运价表。当两地之间的货运量低于120t时，采用第一行的单位运价；当货运量高于120t时，高于120t的部分采用第二行的单位运价。（单位：元/t）
问在考虑到产销地之间直接运输和非直接运输的各种可能方案的情况下，如何将2个加工厂每月生产的产品运往销售地，使总的运费最低？此外，如果对天津和武汉的两个中转站分别增加容量限制：100t、80t，试分析问题的最优解会发生什么样的变化。
模型假设
在对问题进行了深入分析后，为了合理地简化我们的模型，我们提出了以下假设：
1.运往中转站的产品一定会运出，在最终的运输方案中，中转站的运入量等于运出量。
2.产销地都不具有中转作用。也即，产地只能运出产品，销地只能运入产品。事实上，产销地是否可以中转在我们的两个模型中都是一个无关紧要的因素，它不会给模型带来实质性的改变。因此，为使模型的求解过程更加简单，我们假设产销地无中转作用。
线性规划
线性规划模型建立
使用变量xi表示两地间的运量，变量yi表示两地间货物运输的单位运价。并使用0-1变量bi来实现对两地间单位运价的分段。当两地间的运量低于120 t时，bi=1（i为奇数），xi不为0；当两地间的运量高于120 t时，bi+1=1（i为奇数），xi+1不为0。
根据上述变量约定，可得到如下约束：
相应的详细变量说明如表3.1.1、表3.1.2所示：
由于该问题要求求解在考虑到产销地之间直接运输和非直接运输的各种可能方案的情况下，使得总运费最低的运输方案，则可得到目标函数：
值得注意的是，在该模型中我们并没有约定变量xi必须为整数。也就是说，两地间运量的取值范围为正实数。
线性规划模型求解
在3.1中，我们建立了一个线性规划模型，考虑到Matlab软件有专门用于求解混合整数线性规划模型的函数intlinprog，因此我们选择了直接调用该函数来进行编程求解。
调用intlinprog时，应将待求解的模型调整为如下形式：
中转站增加容量限制
考虑到天津中转站至多容纳100 t货物，而武汉中转站至多容纳80 t货物，在原模型的基础上我们增加以下约束条件：
同样利用Matlab的intlinprog函数进行求解，得到在中转站有容量限制的情况下，最佳运输方案对应的最低运费为216900元，与无容量限制时相比，高出6300元。这与直观判断相符合。此时的最佳运输方案如图3.3.1所示：
网络流
最小费用最大流模型建立
分析问题特点可以发现，该问题能够转化为最小运费最大流模型来进行求解。首先，构造有向图D=(V,A)，其中V，A分别表示D中所有点的集合和所有弧的集合。vi∈V（i=1,2,…,7)为有向图D中的第i个点，这些点分别表示产地、销地及中转站。
由于单位运价是分段的，根据实际运量的大小，两地之间在运输时有两个不同的运价，因此任意两地之间也应该由两条不同的弧连通。在连通两地的两条弧中，其中一条弧的容量为120（单位：吨），它对应实际运量低于120吨时的单位运价；另一条弧容量为+∞，它对应实际运量高于120吨时，高于120吨部分的单位运价。
以沈阳和天津为例，在有向图中连接两地的弧应如图4.1.1所示：
在有向图D的基础上构造辅助网络[1]D’=（V’, A’)。在D’中，我们定义一个总的发点s和总的收点t，分别用v0和v8表示，使得所有表示产地的点（即点v1，v2）都发自v0，所有表示销地的点（即点v5，v6，v7）都收于v8。此外，弧集A’与A相比，还增加了两类弧。一类是以各产地为终点的弧，它们的容量为相应产地的产量，单位运价取0；另一类是以各销地为起点的弧，它们的容量为相应销地的销量，单位运价也取0。
基于上述构造原则，辅助网络图D’应如图4.1.2所示：
以辅助网络图D’为基础，建立数学模型。变量Cijk表示从i地到j地的第k条弧的容量（k=1,2），变量Mijk表示从i地到j地的第k条弧的单位运价，变量表示从i地到j地的第k条弧的流量。数学模型如下：
最小费用最大流模型求解
整体算法设计
在4.1中，我们建立了一个典型的网络流模型，它的求解算法本质上是一个不断调整的过程。在这里简单介绍程序设计的基本思想[2]：
1.从最简单的可行流——零流（流量为0）开始考虑。
2.寻找从总发点到总收点的增广路径，这条路径上的每条弧的流量都必须严格小于其容量。
3.对增广路径进行增广，计算每条弧的（容量-流量）中的最小值，记为delta。调整当前流的流量：增广路径上的每条弧的流量加上delta，其他弧的流量不变，这样可保证调整后的新流依然是可行流。
4.不断地从总发点出发寻找增广路径，每次都对其进行增广，直到总发点和总收点不再连通，也就是找不到增广路径为止。当找不到增广路径时，当前的流量就是最大流。
具体落实到程序的编写时，我们只使用一个数组来记录每条弧的容量，而不记录流量。当流量加上delta时，通过容量减去delta来等价表示。这样的转换可以简化程序的编写。
此外，注意到以上的算法在做增广路径时可能会阻塞后来的增广路径。也就是说，如果按照上述方法进行迭代求解，必须要在做增广路径时遵循一定的顺序，只有按照这一指定顺序，才能最终求得最大流。这样便产生了矛盾，因为我们在寻找增广路径时是完全任意的。为了修正这一问题，我们为每条弧构造了一条相应的反向弧，它的费用为对应正相弧费用的相反数，容量为0。在每次流量增加时，都将流量加入到相应的反向弧中，这样就允许后面的流进行自我调整，而不会阻塞以后的增广路径。
贝尔曼-福特算法（Bellman-Ford）及其队列优化算法SPFA
贝尔曼-福特算法（Bellman–Ford），是求解单源最短路径问题的一种基本算法。它的基本原理是对图进行松弛操作，得到所有可能的最短路径。依4.2.1中所述，待求解的有向图中存在着权值为负数的弧，这样在寻找增广路径时，Dijkstra算法就不适用了。因此，我们改用贝尔曼-福特算法来进行求解。下面对这种算法的步骤进行简单的介绍。
假定有向图中存在着n个点，m条边。
1.定义一个距离数组d[n]，除起点对应的d[i]初始化为0外，其余均初始化为正无穷。该数组用于记录当前从总发点到图中各点的距离。
2.考虑最坏的情况，也就是从总发点出发到达总收点所要经过的边数最多的情况，在这种情况下需要对d[n]进行n-1次的更新才能找到最优解。因此设定循环次数为n-1，保证无论是何种有向图，都能迭代得到最优解。
3.在每一次循环中，遍历有向图中所有的弧。在访问每条弧时，假设该弧的起点为u，终点为v，权值为w，判断d[u]+w是否小于d[v]，若小于，则更新d[v]为d[u]+w。
4.n-1次循环结束时，d数组中的值即为从总发点到图中各点的最短距离（在最小费用最大流问题中为最低费用）。
值得注意的是，贝尔曼-福特算法也有不可忽视的缺点。首先，由于以最坏的情况为基础进行考虑，因此对于较好的有向图也不能做出适应性的改变，只能盲目地循环n-1次。其次，每次循环都需要遍历有向图中的每一条弧，但对于在上一次循环中相应的距离数组没有改变的弧来说，本次循环中再做判断是没有任何意义的，判断的结果也是一定的。因此，我们选择维护一个备选节点队列，当且仅当有节点被松弛之后才会加入到队列之中。这就是贝尔曼-福特算法的队列优化算法——SPFA算法。
网络流方法的求解结果与线性规划方法的求解结果相同，因此这里不再重复说明。同时，这两种方法也彼此证明了二者的正确性。
中转站增加容量限制
中转站增加了容量限制后，我们在辅助网络图D’的基础上构造了辅助网络图D’’。
在D’’中，为实现对中转站容量的限制，我们将每一个中转站vi拆分为vi与vi’，并在这两点之间构造了一条全新的弧，这条弧的容量即为中转站的容量，而费用取0。此外，原来在D’中以vi为起始点的弧，改为以vi’为起始点；以vi为终点的弧则保持不变。
值得一提的是，由于我们设计的算法是针对弧进行操作的，所以其实有向图中点的标号并不重要。这些标号是否连续、是否有一定的顺序，都不会对结果产生任何影响。因此，在图中增添新点的操作是方便而简单的。
基于上述构造原则，在辅助网络图D’’中，天津与武汉的两个中转站的局部图如下：
拓展延伸
运输领域的许多其他问题，也可以使用该模型来进行求解。甚至对于某些问题而言，我们完全不需要修改算法设计及具体程序，或者只需要进行简单的修改即可。下面举几个方面的典型例子来进行说明。
产销不平衡问题
在问题详述中我们知道，这是一个产销平衡的问题。但事实上，模型本身及其算法设计并没有对数据做出这样的要求。也就是说，如果该问题由产销平衡问题变化为产销不平衡问题，本文所建立的模型依然是适用的。只要给出相应的产量与销量，我们仍然可以在不作出任何改变的情况下求得最优解。
由于产销不平衡时，问题的求解只不过是改变了输入矩阵，模型与算法都没有任何改变。因此，这里不再进行示例求解。
与传统的最小运费最大流问题的融合
如果对两地间运量的最大值加以限定，也就是考虑到铁路或公路等的承载能力时，我们的模型依然是适用的。在原题中，虽然两地间弧的容量被用以实现单位运价的分段，但是这并不影响我们对这个量进行“复用”。若单位运价仍是以120吨为界进行分段，在两地间运量约束高于120吨时，只需要将容量为无穷的弧的容量改为（约束量-120吨）即可；而在两地间运量约束低于120吨时，则只需要删除容量为无穷的弧，并将容量为120吨的弧的容量改为当前的约束量即可。
下面给出一个简单的求解示例。该示例在上文考虑中转站容量的基础上，考虑了铁路或公路等的承载能力。路径承载能力如表4.4.1.1所示：
从结果中可以发现，由于增加了承载能力这一限制条件，导致位于沈阳的产地有40吨货物积压，最大流也由800减小为760。基于这一点，我们进一步分析可以得出这样一个结论：当且仅当以某产地为起点的各弧的“容量和”大于等于该产地的产量或以某销地为终点的各弧的“容量和”大于等于该销地的销量时，才能保证所有产地的全部产品都能够按照计划运往各销地。否则，该公司需要考虑拓展其他的运输方式来满足运输需求。
线性规划与网络流的联系与对比
线性规划是运筹学中研究较早、发展较快、应用广泛、方法较成熟的一个重要分支。迄今为止，已有许多行之有效的求解算法，例如单纯形算法、Karmarkar算法及其各种形式的改进算法等等。但这些算法在求解大规模问题时却显得不够理想，计算量太大。
据统计，专门的、有效的网络流算法可比一般的线性规划算法快10~20倍[3]，而很多大规模的线性规划问题往往是稀疏的、易转化为网络流问题的。本文所研究的最小费用最大流问题就是一种可以转化为网络流问题的线性规划问题。
在线性规划的求解部分，我们使用了Matlab中的intlinprog函数。这是一个通用的函数，它适用于所有的混合整数规划问题。针对本题，它的运行时间约为10ms。然而，当我们将这个问题转化为网络流问题进行求解时，同样去计算程序运行的时间可以发现，每次得到的运行时间并不相同。这是因为网络流中增广路径的寻找是任意的。因此，为尽量精确地估算其运行时间，我们将程序运行了10次求其平均值[4]，得到结果约为3077.63us，为intlinprog函数求解时间的三分之一。这个结果进一步证明：专门的、有效的网络流算法的时间复杂度比一般的线性规划算法更低。因此，在求解大规模的线性规划问题时，考虑其是否能够转化为网络流问题就显得十分必要。
在上文的求解过程中，线性规划和网络流是几乎完全不相干的，我们是单纯从网络流的视角来建立的网络流模型。那么，如果立足于3中已经建立好的线性规划模型，又应该如何将它转化为网络流模型以降低算法的时间复杂度，并提高运行速度呢？这是一个有普适意义的问题，因此，下面对转化方法进行简单的阐述。
首先必须明确，并不是所有的线性规划问题都可以使用网络流的方法来进行求解，要将线性规划问题转化为网络流问题必须要满足这种转化的前提条件。也就是说，只有当线性规划能够转化成每个变量都出现两次，且系数分别为+1和-1的形式时，才可以使用网络流的方法巧妙地解决[5]。
进行转化时，我们首先需要将约束条件中所有的不等式都加上一个变量，化为等式。根据网络流的性质：除了源点和汇点外，所有点的流入量总等于流出量。因此，我们可以把每一个方程都当做一个点，并把有范围限制的变量看做弧。建图时，对于每个变量，把+1看成流出-1看成流入，找到对应方程的位置连弧即可。
在建图过程中，如果某变量的出现次数不为2，则可以进行适当的差分，变成每个变量只出现两次并且一次系数为+1，一次系数为-1的形式。此外，使用这种方法所建立的图可能是不连通的，而这时就需要增加一个超级源点来保证有向图的连通性。
综上所述，线性规划与网络流的关系十分密切：所有的网络流问题都可以转化为线性规划问题，而满足一定条件的线性规划问题才可能转化为网络流问题。鉴于专门的、有效的网络流算法的时间复杂度比一般的线性规划算法更低，与直接求解线性规划问题相比，将能够转化为网络流的线性规划问题先进行转化再求解是更高效的。
模型评价
对于本文所研究的主要问题来说，线性规划模型是比较容易建立的，且它对中转站容量的限制也是比较容易实现的。但在使用该模型时，由于变量较多，输入矩阵较大，数据的增、删、改都相对比较麻烦，降低了其实用性。
而对于网络流方法而言，辅助网络图的构造需要一定的经验和技巧，模型的建立较为困难。此外，算法的设计与程序的编写也较为复杂，需要时间打磨和调试。但其数据量较线性规划方法少，且增、删、改都比较容易，将程序封装后有较好的实用性。
但从另一个角度来讲，线性规划方法的普适性更好，在面对不同问题时的灵活度更大。虽然网络流也具有一定的灵活性，但较之线性规划却略显逊色，这也是它不可忽视的缺点。
参考文献
[1]董鹏,杨超,陈新.一类带容量限制的运输问题[J].海军工程大学学报,2004(05):96-99.
[2][WhiteJunior]的原创文章：最小费用最大流问题与算法实现(Bellman-Ford、SPFA、Dijkstra)，CC 4.0 BY-SA版权协议，原文链接：https://blog.csdn.net/lym940928/article/details/90209172.
[3]陈开周, 竺菊梅. 线性规划问题转化为网络流问题的条件和算法[J]. 西安电子科技大学学报, 1993(A12):63-69.
[4][Bamboo_Shui]的原创文章：C语言 计算程序运行时间（精确到毫秒/微秒），CC 4.0 BY-SA版权协议，原文链接：https://blog.csdn.net/Bamboo_shui/article/details/78757281
[5][fcwww]的原创文章：网络流解线性规划问题 BZOJ1061: [Noi2008]志愿者招募，原文链接：https://www.cnblogs.com/suika/p/9079315.html
借鉴的代码链接: link.
#include "stdio.h"
#include <iostream>
#include <algorithm>
#include <map>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <vector>
#include <queue>
#include <stack>
#define max 5000
#define MAXN 5050
#define INF 0x3f3f3f3f
using namespace std;
int n, m, s, t;
int u, v, c, w;
int maxFlow, minCost;
struct Edge
{
int from, to, flow, cap, cost;
};
bool vis[MAXN];
int p[MAXN], a[MAXN], d[MAXN],route[MAXN];
vector<int> g[MAXN];
vector<Edge> edges;
void init(int n)
{
for (int i = 0; i <= n; i++)
g[i].clear();
edges.clear();
}
void addedge(int from, int to, int cap, int cost)
{
Edge temp1 = { from, to, 0, cap, cost };
Edge temp2 = { to, from, 0, 0, -cost };//允许反向增广
edges.push_back(temp1);
edges.push_back(temp2); //加进来两条边
int len = edges.size();
g[from].push_back(len - 2);
g[to].push_back(len - 1);
}
//贝尔曼-福特算法实现
bool bellmanford(int s, int t)
{
for (int i = 0; i < MAXN; i++)
{
d[i] = INF; route[i] = 0;
}
d[s] = 0;
memset(vis, false, sizeof(vis));
memset(p, -1, sizeof(p));
p[s] = -1;
a[s] = INF;
queue<int> que;
que.push(s);
vis[s] = true;
while (!que.empty())
{
int u = que.front();
que.pop();
vis[u] = false;
for (int i = 0; i < g[u].size(); i++)
{
Edge& e = edges[g[u][i]];
if (e.cap > e.flow && d[e.to] > d[u] + e.cost )//进行松弛，寻找最短路径也就是最小费用
{
route[e.to] = e.from;
d[e.to] = d[u] + e.cost;
p[e.to] = g[u][i];
a[e.to] = min(a[u], e.cap - e.flow);
if (!vis[e.to])
{
que.push(e.to);
vis[e.to] = true;
}
}
}
}
if (d[t] == INF)
return false;
maxFlow += a[t];
minCost += d[t] * a[t];
for (int i = t; i != s; i = edges[p[i]].from)
{
edges[p[i]].flow += a[t];
edges[p[i] ^ 1].flow -= a[t];
}
return true;
}
void MCMF()
{
while (bellmanford(s, t))
continue;
cout << endl;
cout << "我要开始输出结果啦！" << endl;
for (int i = 0; i < 20; i = i+2)
cout << edges[i].flow << " ";
cout << endl;
for (int i = 20; i < 40; i=i+2)
cout << edges[i].flow << " ";
cout << endl;
for (int i = 40; i < 74; i=i+2)
cout << edges[i].flow << " ";
cout << endl;
cout << endl;
return;
}
int main()
{
cout << "节点数为："; cin >> n;
cout << "边数为："; cin >> m;
cout << "源点编号为："; cin >> s;
cout << "汇点编号为："; cin >> t;
/*int shuru[] = { 1,3,max,210, 1,4,120,185,
1,8,max,210, 1,7,120,200,
1,10,max,500, 1,9,120,410,
1,11,120,620,
1,6,max,470, 1,5,120,400,
0,1,500,0, 0,2,300,0,
2,3,max,230, 2,4,120,205,
2,8,max,210, 2,7,120,200,
2,10,max,280, 2,9,120,240,
2,11,120,350,
2,5,120,150, 2,6,max,160,
3,8,max,100, 4,7,120,80,
3,10,max,170, 4,9,120,150,
4,11,120,250,
4,5,120,135, 3,6,max,160,
5,4,120,140, 6,3,max,160,
5,7,120,270, 6,8,max,315,
5,9,120,110, 6,10,max,130,
5,11,120,130,
7,12,max,0, 8,12,max,0,
9,13,max,0, 10,13,max,0,
12,14,300,0, 13,14,400,0, 11,14,100,0,
3,4,max,0,4,3,max,0,
5,6,max,0,6,5,max,0
};*/
/*int shuru[] = { 1,3,max,210, 1,4,120,185,
1,8,max,210, 1,7,120,200,
1,10,max,500, 1,9,120,410,
1,11,120,620,
1,6,max,470, 1,5,120,400,
0,1,500,0, 0,2,300,0,
2,3,max,230, 2,4,120,205,
2,8,max,210, 2,7,120,200,
2,10,max,280, 2,9,120,240,
2,11,120,350,
2,5,120,150, 2,6,max,160,
17,8,max,100, 18,7,120,80,
17,10,max,170, 18,9,120,150,
18,11,120,250,
18,5,120,135, 17,6,max,160,
21,4,120,140, 22,3,max,160,
21,7,120,270, 22,8,max,315,
21,9,120,110, 22,10,max,130,
21,11,120,130,
7,12,max,0, 8,12,max,0,
9,13,max,0, 10,13,max,0,
12,14,300,0, 13,14,400,0, 11,14,100,0,
3,15,max,0, 4,15,max,0,
15,16,100,0,
16,17,max,0, 16,18,max,0,
5,19,max,0, 6,19,max,0,
19,20,80,0,
20,21,max,0, 20,22,max,0
};*/
//简化版 不带容量限制
/*int shuru[] = { 1,3,max,210, 1,3,120,185,
1,7,max,210, 1,7,120,200,
1,9,max,500, 1,9,120,410,
1,11,120,620,
1,5,max,470, 1,5,120,400,
0,1,500,0, 0,2,300,0,
2,3,max,230, 2,3,120,205,
2,7,max,210, 2,7,120,200,
2,9,max,280, 2,9,120,240,
2,11,120,350,
2,5,120,150, 2,5,max,160,
3,7,max,100, 3,7,120,80,
3,9,max,170, 3,9,120,150,
3,11,120,250,
3,5,120,135, 3,5,max,160,
5,3,120,140, 5,3,max,160,
5,7,120,270, 5,7,max,315,
5,9,120,110, 5,9,max,130,
5,11,120,130,
7,14,300,0, 9,14,400,0, 11,14,100,0,
};*/
int shuru[] = { 1,3,max,210, 1,3,120,185,
1,7,max,210, 1,7,120,200,
1,9,max,500, 1,9,120,410,
1,11,120,620,
1,5,max,470, 1,5,120,400,
0,1,500,0, 0,2,300,0,
2,3,max,230, 2,3,120,205,
2,7,max,210, 2,7,120,200,
2,9,max,280, 2,9,120,240,
2,11,120,350,
2,5,120,150, 2,5,max,160,
4,7,max,100, 4,7,120,80,
4,9,max,170, 4,9,120,150,
4,11,120,250,
4,5,120,135, 4,5,max,160,
6,3,120,140, 6,3,max,160,
6,7,120,270, 6,7,max,315,
6,9,120,110, 6,9,max,130,
6,11,120,130,
7,14,300,0, 9,14,400,0, 11,14,100,0,
3,4,100,0, 5,6,80,0
};
int i = 0;
while (m--)
{
u = shuru[i];
i = i + 1;
v = shuru[i];
i = i + 1;
c = shuru[i];
i = i + 1;
w = shuru[i];
i = i + 1;
addedge(u, v, c, w);
}
MCMF();
cout << endl;
cout << "最大流为：" << maxFlow << endl;
cout << "最小费用为" << minCost << endl;
cout << endl;
system("pause");
return 0;
}
扫一扫
1820
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
