Python小白的数学建模课-20.网络流优化案例-编程知识网

本站消息

出租广告位,需要合作请联系站长

今日名言-想象你自己对困难作出的反应，不是逃避或绕开它们，而是面对它们，同它们打交道，以一种进取的和明智的方式同它们奋斗。——马克斯威尔·马尔兹

今日名言-用谅解、宽恕的目光和心理看人、待人。人就会觉得葱笼的世界里，春意盎然，到处充满温暖。——蔡文甫

小酷狗

文章

24403

访问

+关注

分类

python(8)

java(2)

日期归档

暂无数据

Python小白的数学建模课-20.网络流优化案例

发布于2021-07-25 06:30 阅读(686) 评论(0) 点赞(23) 收藏(2)

在实际工作和数模竞赛中，网络最大流问题、最小费用流问题和最小费用最大流问题都有很多延伸和推广的应用。
本文介绍了常见的网络最大流问题、最小费用流问题的应用与推广，及其解题方法。
本文选择多源多汇物流转运问题、多商品流问题案例，详细介绍网络流问题的分析方法和解决方案，并使用线性规划方法建模和编程。
『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。

1. 网络流问题的应用与推广

流在生活中十分常见，例如交通系统中的人流、车流、物流，供水管网中的水流，金融系统中的现金流，网络中的信息流。网络流优化问题是基本的网络优化问题，应用非常广泛，遍及通讯、运输、电力、工程规划、任务分派、设备更新以及计算机辅助设计等领域。

网络流优化问题最重要的指标是每条边的成本和容量限制，既要考虑成本最低（最短路径问题），又要满足容量限制（最大流问题），由此产生了网络最大流问题、最小费用流问题和最小费用最大流问题。在《Python小白的数学建模课-19.网络流优化》一文中，我们详细介绍了网络最大流问题、最小费用流问题和最小费用最大流问题的算法、编程和案例。

不论在实际工作中，还是在数模竞赛中，网络最大流问题、最小费用流问题和最小费用最大流问题都有很多延伸和推广的应用。小白同学对于网络图的概念和思维逻辑比较生疏，面对变换和深化的题型往往对题目能有基本判断和思路、但找不到具体的解决方法，觉得应该能做出来但就是差一点。

本文选择多源多汇物流转运问题、多商品流问题案例，详细介绍网络流问题的分析方法和解决方案，并使用线性规划方法建模和编程。

对于小白同学来说，这两个案例的建模和编程都有些复杂。别急，也别怕，问题分析和程序说明都很详细，只要静下来慢慢读，将程序说明、源程序和运行结果对照着读，相信小白同学也能看懂。

2. 网络最大流问题的应用与推广

2.1 指定流量的可行流

我们在《19.网络流优化》一文讨论了网络最大流的算法。如果不是计算最大流，而是要计算指定流量 v 的可行流，可以用如下图所示的增广路（相当于辅助线）方法处理。

（1）在容量网络 G 上增加一个顶点 t’ 和一条边 (t,t’)，构建一个新的容量网络 G’，令边 (t,t’) 的容量为 v；

（2）对于容量网络 G‘，计算从源点 s 到汇点 t’ 的最大流 $f_{max}$ ，显然 $f_{max}\leq v$ ；

（3）若 $v<f_{max}$ ，则容量网络 G 不存在满足流量 v 的可行流；若 $v=f_{max}$ ，则从最大流 $f_{max}$ 的路径中去掉边 (t,t’) 即得到容量网络 G 满足流量 v 的可行流。

在这里插入图片描述

另一种思路是，设置源点 s /汇点 t 的流量为 v，再按最小费用流即可。这种方法不仅得到了流量 v 的可行流，而且是最小费用流，但该方法的计算量较大。主要程序如下，完整程序详见《Python小白的数学建模课-19.网络流优化》。

G2.add_node("s", demand=-v)  # nx.min_cost_flow() 的设置要求
G2.add_node("t", demand=v)  # 设置源点/汇点的流量
# 求最小费用流(demand=v)
minFlowCost = nx.min_cost_flow_cost(G2)  # 求最小费用流的费用
minFlowDict = nx.min_cost_flow(G2)  # 求最小费用流

2.2 多源点多汇点的最大流

《19.网络流优化》一文讨论的网络流优化问题，都是单源点、单汇点的网络。对于多源点单汇点、单源点多汇点、多源点多汇点的容量网络的最大流问题，可以通过如下图所示的增加一个超级源点和/或超级汇点的方法处理。

（1）在容量网络 G 上增加一个超级源点 s 和一个超级源点 t；从超级源点 s 向网络 G 中的每个源点连一条容量为 $+\infty$ 的边，从网络 G 中的每个汇点向超级汇点 t 连一条容量为 $+\infty$ 的边。于是得到了一个新的容量网络 G’。

（2）对于容量网络 G‘，计算从超级源点 s 到超级汇点 t 的最大流 $f_{max}$ ，即为多源多汇容量网络 G 的最大流。

在这里插入图片描述

使用 NetworkX 工具包，也有另一种特殊方法可以解决多源多汇最大流问题。

NetworkX 中求解费用最小流问题的函数 min_cost_flow()、min_cost_flow_cost() 是求解费用最小流问题的函数，并不是设定供应点、需求点，而是通过设置顶点属性 ‘demand’ 确定供应点、需求点及各顶点的净流量，因而允许网络中存在多个供应点、需求点。

min_cost_flow(G, demand=‘demand’, capacity=‘capacity’, weight=‘weight’)
min_cost_flow_cost(G, demand=‘demand’, capacity=‘capacity’, weight=‘weight’)

因此，对于指定流量的多源多汇网络可行流问题，只要在每个源点、汇点的顶点属性 ‘demand’ 设置对应的供应量、需求量，再按最小费用流即可。这种方法不仅得到了流量 v 的可行流，而且是最小费用流。类似地，对于多源多汇网络的最大流问题，可以参考《19.网络流优化》中“3.4 案例：运输费用”的方法，v 从 1 逐渐递增，计算所有流量下的最小费用流，直到达到网络容量的极限，即得到容量网络的最大费用流。该方法不需要构造新的容量网络，编程实现简单，而且可以获得最大流量的最小费用流，但计算量很大。

2.3 带顶点容量约束的最大流

标准形式的网络流优化问题，只有边的容量约束，没有顶点的容量约束。在实际问题中，顶点所表示的工厂、仓库、销售点，都存在最大储存量，因而带有顶点的容量约束。

对于顶点容量约束的网络流问题，可以通过如下图所示的将每个带约束的顶点 N 分拆为两个顶点（入顶点 Nin 和出顶点 Nout）的方法处理：

（1）将指向顶点 N 的边，改为指向入顶点 Nin；

（2）将顶点 N 所指向的边，改为出顶点 Nout 所指向的边；

（3）增加一条从入顶点 Nin 指向出顶点 Nout 的边，边的容量为顶点 N 的容量，单位费用为 0。

于是得到了一个新的容量网络 G’。由此，原有网络 G 的顶点 N 的容量限制，转化为网络 G’ 的边 (Nin,Nout) 的容量限制。带顶点约束的容量网络 G 的最大流问题，转化为新的标准形式的容量网络 G’ 的最大流问题。

在这里插入图片描述

3. 最小费用流问题的应用与推广

3.1 运输问题

有出发地（供应点、供应量）和目的地（需求点、需求量），没有转运点，没有路段（边）的容量限制，目标是总运输成本最小或总利润最大。

3.2 指派问题

出发地（供应点、供应量=1）是人，目的地（需求点、需求量=1）是任务，没有转运点，没有路段（边）的容量限制，目标是总指派成本最小或总利润最大。

3.3 转运问题

有出发地（供应点、供应量）和目的地（需求点、需求量），有转运点，没有路段（边）的容量限制（或带有容量限制），目标是总流量费用最小或总利润最大。

3.4 最大流问题

最大流问题是特殊的最小费用流问题：有供应点、需求点、转运点，有路段（边）的容量限制，没有供应量和需求量的限制，目标是通过网络到目的地的总流量最大。

3.5 最短路问题

有供应点（供应量=1）、需求点（需求量=1）、转运点，没有路段（边）的容量限制，目标是通过网络到目的地的总距离最短。

4. 案例：多源多汇的物流转运问题

4.1 问题描述

如图所示，某公司的工厂（供应点）位于F1、F2，仓库（中转点）位于W1、W2，零售商（需求点）位于R1、R2、R3、R4。从工厂生产的产品先送到仓库，再由仓库发到零售商。图中给出了各供应点和需求点的流量、每条发货线路的最大流量、单位流量的运输成本，要求在产销平衡的条件下，找出总运输费用最小的运输方案。

在这里插入图片描述

4.2 问题分析

这是一个运输路径规划问题，也是一个多源点、多汇点的最小费用流问题。

对于多源多汇的最小费用流问题，构造一个费用网络 G，网络 G 的各条边没有容量限制，单位流量的运输成本为边的权值 w。F1、F2 是网络 G 的源点（供应点），R1~R4 是汇点（需求点），W1、W2 是中间节点。F1、F2 的供应量与 R1~R4 的需求量是平衡的。因此可以用 NetworkX 的 maximum_flow() 函数求出最小费用流。

该问题也是一个典型的线性规划问题，可以用线性规划方法求解。
$min\;\sum_{i,j\in E} w_{ij}f_{ij}\\ s.t.:\;$

{\begin{cases} \sum_{j} f_{i j} \leq - v_{i}, i \in S \\ \sum_{j} f_{j i} = v_{i}, i \in T \\ \sum_{i} f_{i k} - \sum_{j} f_{k j} = 0, k \notin S, T \end{cases}

$\begin{cases} \sum_j f_{ij} \leq -v_i \;,i \in S\\ \sum_j f_{ji}=v_i \;,i \in T\\ \sum_i f_{ik} - \sum_j f_{kj} = 0 \;,k\notin S,T\\ \end{cases}$

m i n i, j \in E \sum w_{i j} f_{i j} s . t . : ⎩ ⎪ ⎨ ⎪ ⎧ \sum_{j} f_{i j} \leq - v_{i}, i \in S \sum_{j} f_{j i} = v_{i}, i \in T \sum_{i} f_{i k} - \sum_{j} f_{k j} = 0, k \in / S, T

式中：f(i,j) 表示路段 (i,j) 的流量，w(i,j) 表示路段 (i,j) 的单位流量运输成本，S 表示源点（供应点），T 表示汇点（需求点）。

对于上式描述的线性规划问题，可以用 PuLP工具包求解（参见《Python小白的数学建模课-03.线性规划》），编程步骤如下：

（1）导入 PuLP库函数；
（2）定义一个线性规划问题；
（3）定义决策变量，决策变量的上下限：

    f11：F1-W1 流量, 0<=f11<=600
    f12：F1-W2 流量, 0<=f12<=600
    f21：F2-W1 流量, 0<=f21<=400
    f22：F2-W2 流量, 0<=f22<=400
    w11：W1-R1 流量, 0<=w11<=200
    w12：W1-R2 流量, 0<=w12<=150
    w13：W1-R3 流量, 0<=w13<=350
    w14：W1-R4 流量, 0<=w14<=300
    w21：W2-R1 流量, 0<=w21<=200
    w22：W2-R2 流量, 0<=w22<=150
    w23：W2-R3 流量, 0<=w23<=350
    w24：W2-R4 流量, 0<=w24<=300

（4）定义目标函数：

    min cost = 2*f11 + 3*f12 + 3*f21 +   f22 + 2*w11 + 6*w12
             + 3*w13 + 6*w14 + 4*w21 + 4*w22 + 6*w23 + 5*w24

（5）定义约束条件：

    f11 + f12 <= 600
    f21 + f22 <= 400
    f11 + f21 = w11 + w12 + w13 + w14
    f12 + f22 = w21 + w22 + w23 + w24
    w11 + w21 = 200
    w12 + w22 = 150
    w13 + w23 = 350
    w14 + w24 = 300

注意，f11+f12<=600, f21+f22<=400，意味着允许存在产销不平衡的情况。

（6）求解规划问题，输出优化结果

本文给出了基于 NetworkX 工具包和基于 PuLP 工具包的两种方法的程序，以便小白理解图论中的最小费用流问题与线性规划问题的联系。

4.3 Python 例程

# mathmodel21_v1.py
# Demo20 of mathematical modeling algorithm
# Demo of network flow problem optimization with NetworkX
# Copyright 2021 YouCans, XUPT
# Crated：2021-07-18

import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
import networkx as nx  # 导入 NetworkX 工具包
import pulp      # 导入 pulp库

# # 4. 多源多汇最小费用流问题 (Capacity network with multi source and multi sink)
# # 4.1 费用网络
# 创建有向图
G1 = nx.DiGraph()  # 创建一个有向图 DiGraph
G1.add_edges_from([('F1','W1',{'capacity': 9e5, 'weight': 2}), # F1~F2：工厂
                  ('F1','W2',{'capacity': 9e5, 'weight': 3}),  # W1~W2：仓库
                  ('F2','W1',{'capacity': 9e5, 'weight': 3}),
                  ('F2','W2',{'capacity': 9e5, 'weight': 1}),
                  ('W1','R1',{'capacity': 9e5, 'weight': 2}),  # R1~R4：零售店
                  ('W1','R2',{'capacity': 9e5, 'weight': 6}),
                  ('W1','R3',{'capacity': 9e5, 'weight': 3}),
                  ('W1','R4',{'capacity': 9e5, 'weight': 6}),
                  ('W2','R1',{'capacity': 9e5, 'weight': 4}),
                  ('W2','R2',{'capacity': 9e5, 'weight': 4}),
                  ('W2','R3',{'capacity': 9e5, 'weight': 6}),
                  ('W2','R4',{'capacity': 9e5, 'weight': 5})]) # 添加边的属性 'capacity', 'weight'
G1.add_node("F1", demand=-600)  # nx.min_cost_flow() 的设置要求
G1.add_node("F2", demand=-400)  # 设置源点的流量，负值表示净流出
G1.add_node("R1", demand=200)  # 设置汇点的流量，正值表示净流入
G1.add_node("R2", demand=150)
G1.add_node("R3", demand=350)
G1.add_node("R4", demand=300)
pos={'F1':(2,6.5),'F2':(2,3.5),'W1':(5,6.5),'W2':(5,3.5),'R1':(9,8),'R2':(9,6),'R3':(9,4),'R4':(9,2)}  # 指定顶点绘图位置

# # 4.2 用 NetworkX 求最小费用流
# 求最小费用流(demand=v)
minFlowCost = nx.min_cost_flow_cost(G1)  # 求最小费用流的费用
minFlowDict = nx.min_cost_flow(G1)  # 求最小费用流
# minFlowCost, minFlowDict = nx.network_simplex(G2)  # 求最小费用流--与上行等效

# 整理边的标签，用于绘图显示
# 数据格式转换
edgeCapacity = nx.get_edge_attributes(G1, 'weight')
edgeLabel = {}  # 边的标签
for i in edgeCapacity.keys():  # 整理边的标签，用于绘图显示
    edgeLabel[i] = f'w={edgeCapacity[i]:}'  # 边的容量
edgeLists = []
for i in minFlowDict.keys():
    for j in minFlowDict[i].keys():
        edgeLabel[(i,j)] += ',f=' + str(minFlowDict[i][j])  # 取出每条边流量信息存入边显示值
        if minFlowDict[i][j] > 0:
            edgeLists.append((i, j))

print("1. NetworkX 网络与图（最小费用流优化结果）：")  # NetworkX 工具包
print("最小费用:{}".format(minFlowCost))  # 输出最小费用的值
print("最小费用流的路径及流量: ", minFlowDict)  # 输出最小费用流的途径和各路径上的流量
print("最小费用流的路径：", edgeLists)  # 输出最小费用流的途径

# 绘制有向网络图
fig, ax = plt.subplots(figsize=(8,6))
ax.text(1.2,6.4,"600",color='navy')
ax.text(1.2,3.4,"400",color='navy')
ax.text(9.3,7.9,"200",color='navy')
ax.text(9.3,5.9,"150",color='navy')
ax.text(9.3,3.9,"350",color='navy')
ax.text(9.3,1.9,"300",color='navy')
ax.set_title("Capacity network with multi source and multi sink")
nx.draw(G1,pos,with_labels=True,node_color='skyblue',node_size=300,font_size=10)   # 绘制有向图，显示顶点标签
edgeLabel1 = nx.get_edge_attributes(G1, 'weight')
nx.draw_networkx_nodes(G1, pos, nodelist=['F1','F2'], node_color='orange')  # 设置指定顶点的颜色、宽度
nx.draw_networkx_nodes(G1, pos, nodelist=['W1','W2'], node_color='c')  # 设置指定顶点的颜色、宽度
nx.draw_networkx_edge_labels(G1,pos,edgeLabel,font_size=10,label_pos=0.25)  # 显示边的标签：'capacity','weight' + minCostFlow
nx.draw_networkx_edges(G1,pos,edgelist=edgeLists,edge_color='m',width=2)  # 设置指定边的颜色、宽度
plt.xlim(0.5, 10.5)
plt.ylim(1, 9)
plt.axis('on')
# plt.show()

# # 4.3 用线性规划求最小费用流
# (1) 建立优化问题 TransportCost: 求最小值(LpMinimize)
transportLP = pulp.LpProblem("TransportCostProblem", sense=pulp.LpMinimize)  # 定义问题，求最小值
# (2) 定义决策变量 f11~f22, w11~w24
f11 = pulp.LpVariable('F1-W1', lowBound=0, upBound=600, cat='Continuous')  # 定义 f11
f12 = pulp.LpVariable('F1-W2', lowBound=0, upBound=600, cat='Continuous')  # 定义 f12
f21 = pulp.LpVariable('F2-W1', lowBound=0, upBound=400, cat='Continuous')  # 定义 f21
f22 = pulp.LpVariable('F2-W2', lowBound=0, upBound=400, cat='Continuous')  # 定义 f22
w11 = pulp.LpVariable('W1-R1', lowBound=0, upBound=200, cat='Continuous')  # 定义 w11
w21 = pulp.LpVariable('W2-R1', lowBound=0, upBound=200, cat='Continuous')  # 定义 w21
w12 = pulp.LpVariable('W1-R2', lowBound=0, upBound=150, cat='Continuous')  # 定义 w12
w22 = pulp.LpVariable('W2-R2', lowBound=0, upBound=150, cat='Continuous')  # 定义 w22
w13 = pulp.LpVariable('W1-R3', lowBound=0, upBound=350, cat='Continuous')  # 定义 w13
w23 = pulp.LpVariable('W2-R3', lowBound=0, upBound=350, cat='Continuous')  # 定义 w23
w14 = pulp.LpVariable('W1-R4', lowBound=0, upBound=300, cat='Continuous')  # 定义 w14
w24 = pulp.LpVariable('W2-R4', lowBound=0, upBound=300, cat='Continuous')  # 定义 w24
# (3) 定义目标函数 cost
transportLP += (2*f11 + 3*f12 + 3*f21 +   f22 + 2*w11 + 6*w12
             +  3*w13 + 6*w14 + 4*w21 + 4*w22 + 6*w23 + 5*w24)  # 运输费用
# (4) 设置约束条件
transportLP += (f11 + f12 <= 600)  # 等式约束
transportLP += (f21 + f22 <= 400)  # 等式约束
transportLP += (w11 + w21 == 200)  # 等式约束
transportLP += (w12 + w22 == 150)  # 等式约束
transportLP += (w13 + w23 == 350)  # 等式约束
transportLP += (w14 + w24 == 300)  # 等式约束
transportLP += (f11 +f21 == w11 + w12 + w13 + w14)  # 等式约束
transportLP += (f12 +f22 == w21 + w22 + w23 + w24)  # 等式约束
# (5) 求解线性规划问题
transportLP.solve()
# (6) 输出优化结果
print("2. PuLP 线性规划（最小费用的优化结果）：")  # PuLP 工具包
print(transportLP)  # 输出问题设定参数和条件
print("Optimization result")  # PuLP 工具包
for v in transportLP.variables():
    print(v.name, " = ", v.varValue)  # 输出每个变量的最优值
print("\n最小运输费用 = ", pulp.value(transportLP.objective))  # 输出最优解的目标函数值

plt.show()

4.4 程序运行结果

1. NetworkX 网络与图（最小费用流优化结果）：
最小费用:5200
最小费用流的路径及流量:  {'F1': {'W1': 600, 'W2': 0}, 'W1': {'R1': 200, 'R2': 0, 'R3': 350, 'R4': 50}, 'W2': {'R1': 0, 'R2': 150, 'R3': 0, 'R4': 250}, 'F2': {'W1': 0, 'W2': 400}, 'R1': {}, 'R2': {}, 'R3': {}, 'R4': {}}
最小费用流的路径： [('F1', 'W1'), ('W1', 'R1'), ('W1', 'R3'), ('W1', 'R4'), ('W2', 'R2'), ('W2', 'R4'), ('F2', 'W2')]

2. PuLP 线性规划（最小费用的优化结果）：
MINIMIZE
2*F1_W1 + 3*F1_W2 + 3*F2_W1 + 1*F2_W2 + 2*W1_R1 + 6*W1_R2 + 3*W1_R3 + 6*W1_R4 + 4*W2_R1 + 4*W2_R2 + 6*W2_R3 + 5*W2_R4 + 0
SUBJECT TO
_C1: F1_W1 + F1_W2 <= 600
_C2: F2_W1 + F2_W2 <= 400
_C3: W1_R1 + W2_R1 = 200
_C4: W1_R2 + W2_R2 = 150
_C5: W1_R3 + W2_R3 = 350
_C6: W1_R4 + W2_R4 = 300
_C7: F1_W1 + F2_W1 - W1_R1 - W1_R2 - W1_R3 - W1_R4 = 0
_C8: F1_W2 + F2_W2 - W2_R1 - W2_R2 - W2_R3 - W2_R4 = 0
VARIABLES
F1_W1 <= 600 Continuous
F1_W2 <= 600 Continuous
F2_W1 <= 400 Continuous
F2_W2 <= 400 Continuous
W1_R1 <= 200 Continuous
W1_R2 <= 150 Continuous
W1_R3 <= 350 Continuous
W1_R4 <= 300 Continuous
W2_R1 <= 200 Continuous
W2_R2 <= 150 Continuous
W2_R3 <= 350 Continuous
W2_R4 <= 300 Continuous
Optimization result
F1_W1  =  550.0
F1_W2  =  50.0
F2_W1  =  0.0
F2_W2  =  400.0
W1_R1  =  200.0
W1_R2  =  -0.0
W1_R3  =  350.0
W1_R4  =  -0.0
W2_R1  =  0.0
W2_R2  =  150.0
W2_R3  =  0.0
W2_R4  =  300.0
最小运输费用 =  5200.0

在这里插入图片描述

5. 案例：多商品流问题

多商品流问题（Multi-Commodity Flow Prolem）是多种商品（物流、人流、数据流等）在具有容量限制的网络中从供应点流到需求点的网络流问题，通常是 NP 问题。多商品流问题的优化目标是以最小的成本实现商品在网络中的流通，约束条件主要是每条边的容量。

多商品流路径优化问题具有广泛的实际应用，如交通运输、物流网络、电信网络以及产销系统。一些实际问题带有特殊要求，如在电信网络中要考虑时间延迟和可靠性问题，在零担货运中需要考虑配货拼车，有时需要综合考虑运输的时间和成本。

5.1 问题描述

将城市 v0 生产的商品 A、B 通过铁路网分别运送到城市 v6、v5，商品 A、B 的需求量分别为 6.0、5.0 单位。铁路网的结构如图所示，每条线路（边）上的数值表示该线路的容量和单位运费。求完成商品运输任务的最小费用的运输方案。

在这里插入图片描述

5.2 问题分析

这是一个多商品流问题，具有 1 个供应点 V0 和 2个需求点 V5、V6。

对于多源多汇的多商品流最小费用流问题，构造一个费用网络 G，网络 G 的各条边具有容量限制，单位流量的运输成本为边的权值 w。源点（供应点）为流量净输出，汇点（需求点）为流量净输入；中间节点不带存储，每种商品的净流量都为 0。所有线路的商品总流量不超过边的容量限制。

该问题也是一个典型的线性规划问题，可以用如下模型描述：
$min\;\sum_m\sum_{i,j\in E} w_{ij}f_{ij,m}\\ s.t.:\;$

{\begin{cases} \sum_{m} f_{i j, m} \leq= c_{i j}, m \in M \\ \sum_{j} f_{i j, m} = - v_{i, m}, i \in S \\ \sum_{j} f_{j i, m} = v_{i, m}, i \in T \\ \sum_{i} f_{i k, m} - \sum_{j} f_{k j, m} = 0, k \notin S, T \end{cases}

$\begin{cases} \sum_m f_{ij,m} \leq = c_{ij}\;,m \in M\\ \sum_j f_{ij,m}=-v_{i,m} \;,i \in S\\ \sum_j f_{ji,m}=v_{i,m} \;,i \in T\\ \sum_i f_{ik,m} - \sum_j f_{kj,m} = 0 \;,k\notin S,T\\ \end{cases}$

m i n m \sum i, j \in E \sum w_{i j} f_{i j, m} s . t . : ⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ \sum_{m} f_{i j, m} \leq = c_{i j}, m \in M \sum_{j} f_{i j, m} = - v_{i, m}, i \in S \sum_{j} f_{j i, m} = v_{i, m}, i \in T \sum_{i} f_{i k, m} - \sum_{j} f_{k j, m} = 0, k \in / S, T

式中：M 表示商品品种，S 表示源点（供应点），T 表示汇点（需求点），w(i,j) 表示路段 (i,j) 的单位流量运输成本，f(i,j,m) 表示商品 m 在路段 (i,j) 的流量，

NetworkX 工具包中没有多商品流问题的解决方法，在 Python 其它图与网络工具包也没有找到多商品流问题的函数。网络流问题本质上是线性规划问题，对于上式描述的线性规划问题，可以用 PuLP工具包求解（参见《Python小白的数学建模课-03.线性规划》）。

多商品流的目标是各种商品从供应点运送到需求点的总费用最小。定义各商品 m 在各路段 (i,j) 的流量 f(i,j,m) 为决策变量，则目标函数为商品 m 在各路段 (i,j) 的运输费用之和 $\sum_m\sum_{i,j\in E} w_{ij}f_{ij,m}$ 。多商品流的模型主要包含流守恒约束和容量限制。

本案例的问题中，不同商品 A、B 是从同一供应点发出，运送到不同需求点。但本案例的数学模型和例程，并不限定从不同商品从同一供应点发出，也就是说可以适用于多源多汇的多商品流问题。另外，如果如果不同商品在各条边的运费单价不同，只要相应修改例程的目标函数即可处理。

5.3 程序说明

本程序并不复杂，但有些地方的处理可能会令小白感到困惑，详细说明如下：

图的输入：稀疏有向图，使用 nx.DiGraph() 定义有向图，使用G.add_weighted_edges_from() 函数以列表向图中添加多条赋权边，每个赋权边以元组 (node1,node2,{‘capacity’:c1, ‘weight’:w1}) 定义属性 ‘capacity’ 和 ‘weight’。注意必须以关键字 ‘capacity’ 表示容量，以 ‘weight’ 表示单位流量的成本。
supplyA、supplyB 表示商品 A、B 的供应量，来自题目要求。注意网络的最大流量必须大于等于各种商品的总供应量或总需求量，否则将没有可行解。
绘制网络图，以边的标签显示边的容量 ‘capacity’ 和单位流量的成本 ‘weight’ 。
用 PuLP 求解多商品流线性规划问题：

（1）定义一个线性规划问题 MCFproblem。
（2）定义决策变量 fA、fB ，并设置决策变量类型和上下限：

fA、fB 都是一维数组，表示商品 A、B 在网络 G2 各条边的实际流量。根据题意将决策变量设为连续变量，如果问题要求整车配货等条件，可以相应地设为整数变量。
fA、fB 的下限为 0，上限应为各条边的容量。但由于 pulp.LpVariable.dicts 函数中定义上限 upBound 为实数类型，不能是数组，无法对每个决策变量设定不同的上限值。
因此，在定义决策变量 fA、fB 时，将上限设为所有边的容量的最大值 maxCapacity。每条边的具体容量约束，则作为约束条件处理。
G2.edges() 表示遍历网络 G2 的边，因此 fA、fB 是 A、B 在所有边 edges 的流量：

{('v0', 'v1'): FlowA_('v0',_'v1'), ('v0', 'v2'): FlowA_('v0',_'v2'), ('v1', 'v3'): FlowA_('v1',_'v3'), ('v2', 'v1'): FlowA_('v2',_'v1'), ('v2', 'v4'): FlowA_('v2',_'v4'), ('v3', 'v4'): FlowA_('v3',_'v4'), ('v3', 'v5'): FlowA_('v3',_'v5'), ('v3', 'v6'): FlowA_('v3',_'v6'), ('v4', 'v1'): FlowA_('v4',_'v1'), ('v4', 'v5'): FlowA_('v4',_'v5'), ('v4', 'v6'): FlowA_('v4',_'v6')}
{('v0', 'v1'): FlowB_('v0',_'v1'), ('v0', 'v2'): FlowB_('v0',_'v2'), ('v1', 'v3'): FlowB_('v1',_'v3'), ('v2', 'v1'): FlowB_('v2',_'v1'), ('v2', 'v4'): FlowB_('v2',_'v4'), ('v3', 'v4'): FlowB_('v3',_'v4'), ('v3', 'v5'): FlowB_('v3',_'v5'), ('v3', 'v6'): FlowB_('v3',_'v6'), ('v4', 'v1'): FlowB_('v4',_'v1'), ('v4', 'v5'): FlowB_('v4',_'v5'), ('v4', 'v6'): FlowB_('v4',_'v6')}

（3）定义目标函数：

目标函数是 A、B 在所有边 edges 的总运费，即流量 fA、fB 与单位流量成本 ‘weight’的乘积的总和。

需要说明的是：（1）通过 edgeWeight = nx.get_edge_attributes(G2, ‘weight’) 将所有边的属性 ‘weight’ 转换为字典类型，就可以用 for edge in G2.edges 遍历操作。（2）如果不同商品在各条边的运费单价不同，只要将目标函数中的 $KaTeX parse error: Undefined control sequence: \* at position 8: w(edge)\̲*̲ [fA(edge)+fB(e…$ 修改为 $KaTeX parse error: Undefined control sequence: \* at position 9: wA(edge)\̲*̲fA(edge)+wB(edg…$ 即可。

# (3). 设置目标函数
MCFproblem += pulp.lpSum([edgeWeight[edge] * (fA[edge]+fB[edge]) for edge in G2.edges])  # 总运输费用

（4）定义约束条件：

多商品流的模型主要包含边的容量约束和顶点的流守恒约束和两个方面。

边的容量约束，即每条边上各种商品流量之和不大于边的容量。
顶点的流守恒约束，分为源点（供应点）、汇点（需求点）和中间节点三种情况：
- 通过 for node in G2.nodes 可以遍历网络的所有顶点。
- in_degree(node) 计算顶点 node 的入度，入度为 0 的顶点是源点；out_degree(node) 计算顶点 node 的出度，出度为 0 的顶点是汇点。由此可以判断网络的源点、汇点和中间节点；也可以根据题目直接判断各中商品的供应点、需求点进行定义。
- in_edges(nbunch=node) 返回顶点 node 的所有入边，out_edges(nbunch=node) 返回顶点 node 的所有出边。
- 源点：对于每种商品，所有输出边的流量之和不大于该顶点该商品的供应量。对于供应量、需求量相等的问题，也可以将该约束条件设为输出流量等于供应量。
- 汇点：对于每种商品，所有输入边的流量之和等于该顶点该商品的供应量。案例问题有多个需求点，并对应了不同品种的商品，因此需要分别进行设置。
- 中间节点：对于每种商品，中间节点的总流入量与总流出量相等。
- 注意，如果供应点、需求点也是流通线路而不是源点、汇点，则不必分作三种情况，都按照中间节点处理，但对每个节点需要设置每种商品的净流量。

    源点：v0, 出边：[('v0', 'v1'), ('v0', 'v2')]
    中间点：v1, 入边：[('v0', 'v1'), ('v2', 'v1'), ('v4', 'v1')], 出边：[('v1', 'v3')]
    中间点：v2, 入边：[('v0', 'v2')], 出边：[('v2', 'v1'), ('v2', 'v4')]
    中间点：v3, 入边：[('v1', 'v3')], 出边：[('v3', 'v4'), ('v3', 'v5'), ('v3', 'v6')]
    中间点：v4, 入边：[('v2', 'v4'), ('v3', 'v4')], 出边：[('v4', 'v1'), ('v4', 'v5'), ('v4', 'v6')]
    汇点：v5, 入边：[('v3', 'v5'), ('v4', 'v5')]
    汇点：v6, 入边：[('v3', 'v6'), ('v4', 'v6')]

（6）求解规划问题，输出优化结果

5.4 Python 例程

# mathmodel21_v1.py
# Demo20 of mathematical modeling algorithm
# Demo of network flow problem optimization with NetworkX
# Copyright 2021 YouCans, XUPT
# Crated：2021-07-22

import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
import networkx as nx  # 导入 NetworkX 工具包
import pulp      # 导入 pulp库

# 5. 多商品流问题 (Multi-commodity Flow Problem)
# 5.1 创建有向图
G2 = nx.DiGraph()  # 创建一个有向图 DiGraph
G2.add_edges_from([('v0','v1',{'capacity': 7, 'weight': 4}),
                  ('v0','v2',{'capacity': 8, 'weight': 4}),
                  ('v1','v3',{'capacity': 9, 'weight': 1}),
                  ('v2','v1',{'capacity': 5, 'weight': 5}),
                  ('v2','v4',{'capacity': 9, 'weight': 4}),
                  ('v3','v4',{'capacity': 6, 'weight': 2}),
                  ('v3','v5',{'capacity': 10, 'weight': 6}),
                  ('v3','v6',{'capacity': 10, 'weight': 6}),
                  ('v4','v1',{'capacity': 2, 'weight': 1}),
                  ('v4','v5',{'capacity': 5, 'weight': 2}),
                  ('v4','v6',{'capacity': 5, 'weight': 2})]) # 添加边的属性 'capacity', 'weight'
pos={'v0':(0,5),'v1':(4,2),'v2':(4,8),'v3':(10,2),'v4':(10,8),'v5':(15,3),'v6':(15,7)}  # 指定顶点绘图位置

supplyA = 6.0  # 商品 A 供应量
supplyB = 5.0  # 商品 B 供应量
print("Supply A:{}\tSupply B:{}".format(supplyA,supplyB))

# 整理边的标签
edgeLabel1 = nx.get_edge_attributes(G2, 'capacity')
edgeLabel2 = nx.get_edge_attributes(G2, 'weight')
edgeLabel = {}
for i in edgeLabel1.keys():
    edgeLabel[i] = f'({edgeLabel1[i]:},{edgeLabel2[i]:})'  # 边的(容量，成本)

# 5.2 绘制有向网络图
fig, ax = plt.subplots(figsize=(8,6))
nx.draw(G2,pos,with_labels=True,node_color='skyblue',node_size=400,font_size=10)   # 绘制有向图，显示顶点标签
nx.draw_networkx_nodes(G2, pos, nodelist=['v0'], node_color='orange',node_size=400)  # 设置指定顶点的颜色、宽度
nx.draw_networkx_nodes(G2, pos, nodelist=['v5','v6'], node_color='c',node_size=400)  # 设置指定顶点的颜色、宽度
nx.draw_networkx_edge_labels(G2,pos,edgeLabel,font_size=10)  # 显示边的标签：'capacity','weight'
ax.set_title("Multi-commodity Flow Problem by youcans@xupt")
ax.text(-1.8,5.2,"A:{}".format(supplyA),color='m')
ax.text(-1.8,4.6,"B:{}".format(supplyB),color='navy')
ax.text(15.8,7.0,"A:{}".format(supplyA),color='m')
ax.text(15.8,2.8,"B:{}".format(supplyB),color='navy')
plt.xlim(-3, 18)
plt.ylim(1, 9)
plt.axis('on')
# plt.show(YouCans-XUPT)

# 5.3  用 PuLP 求解多商品流最小费用问题 (Multi-commodity Flow Problem YouCans-XUPT)
edgeWeight = nx.get_edge_attributes(G2, 'weight')  # 'weight', 单位流量的成本
edgeCapacity = nx.get_edge_attributes(G2, 'capacity')  # 'capacity', 边的容量
maxCapacity = max(edgeCapacity.values())  # 边的容量的最大值
print(edgeWeight)
print("max(Weight)",max(edgeWeight.values()))
print("max(Capacity)",max(edgeCapacity.values()))

# (1) 建立优化问题 MCFproblem: 求最小值(LpMinimize)
MCFproblem = pulp.LpProblem("MultiCommodityFlowProb", sense=pulp.LpMinimize)  # 定义问题，求最小值

# (2) 定义决策变量 fA(edges), fB(edges)
# itemsG2 = ["A","B"]  # 商品种类
fA = pulp.LpVariable.dicts("FlowA", G2.edges(), lowBound=0.0, upBound=maxCapacity, cat='Continuous')
fB = pulp.LpVariable.dicts("FlowB", G2.edges(), lowBound=0.0, upBound=maxCapacity, cat='Continuous')
print(fA)
print(fB)

# (3). 设置目标函数
MCFproblem += pulp.lpSum([edgeWeight[edge] * (fA[edge]+fB[edge]) for edge in G2.edges])  # 总运输费用

# (4) 设置约束条件
for edge in G2.edges:  # 边的最大流量约束
    MCFproblem += (fA[edge] + fB[edge] - edgeCapacity[edge] <= 0)  # edgeCapacity[edge], 边的容量

for node in G2.nodes:  # 顶点的净流量约束
    if G2.in_degree(node) == 0:  # 入度为 0，判断是否为源点
        print("源点：{}, 出边：{}".format(node,G2.out_edges(nbunch=node)))
        MCFproblem += (sum(fA[edge] for edge in G2.out_edges(nbunch=node)) <= supplyA)  # A 供应量约束
        MCFproblem += (sum(fB[edge] for edge in G2.out_edges(nbunch=node)) <= supplyB)  # B 供应量约束
    elif G2.out_degree(node) == 0:  # 出度为 0，判断是否为汇点
        print("汇点：{}, 入边：{}".format(node,G2.in_edges(nbunch=node)))
        if node=='v6':  # 题目条件, v6 需求为 B
            MCFproblem += (sum(fA[edge] for edge in G2.in_edges(nbunch=node)) == supplyA)  # A 需求量约束
        if node=='v5':  # 题目条件, v5 需求为 A
            MCFproblem += (sum(fB[edge] for edge in G2.in_edges(nbunch=node)) == supplyB)  # B 需求量约束
    else:  # 中间节点，每种商品都是流量平衡
        print("中间点：{}, 入边：{}, 出边：{}".format(node,G2.in_edges(nbunch=node),G2.out_edges(nbunch=node)))
        MCFproblem += (sum(fA[edge1] for edge1 in G2.out_edges(nbunch=node))
                       - sum(fA[edge2] for edge2 in G2.in_edges(nbunch=node)) == 0)  # 总流出 = 总流入
        MCFproblem += (sum(fB[edge1] for edge1 in G2.out_edges(nbunch=node))
                       - sum(fB[edge2] for edge2 in G2.in_edges(nbunch=node)) == 0)  # 总流出 = 总流入
# (5) 求解线性规划问题
MCFproblem.solve()

# (6) 输出优化结果
print("2. PuLP 线性规划（最小费用的优化结果）：")  # PuLP 工具包
print(MCFproblem)  # 输出问题设定参数和条件
print("Optimization result")  # PuLP 工具包
for v in MCFproblem.variables():
    print(v.name, " = ", v.varValue)  # 输出每个变量的最优值
print("\n最小运输费用 = ", pulp.value(MCFproblem.objective))  # 输出最优解的目标函数值

plt.show()

5.5 程序运行结果

2. PuLP 线性规划（最小费用的优化结果）：
MultiCommodityFlowProb:
MINIMIZE
4*FlowA_('v0',_'v1') + 4*FlowA_('v0',_'v2') + 1*FlowA_('v1',_'v3') + 5*FlowA_('v2',_'v1') + 4*FlowA_('v2',_'v4') + 2*FlowA_('v3',_'v4') + 6*FlowA_('v3',_'v5') + 6*FlowA_('v3',_'v6') + 1*FlowA_('v4',_'v1') + 2*FlowA_('v4',_'v5') + 2*FlowA_('v4',_'v6') + 4*FlowB_('v0',_'v1') + 4*FlowB_('v0',_'v2') + 1*FlowB_('v1',_'v3') + 5*FlowB_('v2',_'v1') + 4*FlowB_('v2',_'v4') + 2*FlowB_('v3',_'v4') + 6*FlowB_('v3',_'v5') + 6*FlowB_('v3',_'v6') + 1*FlowB_('v4',_'v1') + 2*FlowB_('v4',_'v5') + 2*FlowB_('v4',_'v6') + 0
SUBJECT TO
_C1: FlowA_('v0',_'v1') + FlowB_('v0',_'v1') <= 7
_C2: FlowA_('v0',_'v2') + FlowB_('v0',_'v2') <= 8
_C3: FlowA_('v1',_'v3') + FlowB_('v1',_'v3') <= 9
_C4: FlowA_('v2',_'v1') + FlowB_('v2',_'v1') <= 5
_C5: FlowA_('v2',_'v4') + FlowB_('v2',_'v4') <= 9
_C6: FlowA_('v3',_'v4') + FlowB_('v3',_'v4') <= 6
_C7: FlowA_('v3',_'v5') + FlowB_('v3',_'v5') <= 10
_C8: FlowA_('v3',_'v6') + FlowB_('v3',_'v6') <= 10
_C9: FlowA_('v4',_'v1') + FlowB_('v4',_'v1') <= 2
_C10: FlowA_('v4',_'v5') + FlowB_('v4',_'v5') <= 5
_C11: FlowA_('v4',_'v6') + FlowB_('v4',_'v6') <= 5
_C12: FlowA_('v0',_'v1') + FlowA_('v0',_'v2') <= 6
_C13: FlowB_('v0',_'v1') + FlowB_('v0',_'v2') <= 5
_C14: - FlowA_('v0',_'v1') + FlowA_('v1',_'v3') - FlowA_('v2',_'v1')
 - FlowA_('v4',_'v1') = 0
_C15: - FlowB_('v0',_'v1') + FlowB_('v1',_'v3') - FlowB_('v2',_'v1')
 - FlowB_('v4',_'v1') = 0
_C16: - FlowA_('v0',_'v2') + FlowA_('v2',_'v1') + FlowA_('v2',_'v4') = 0
_C17: - FlowB_('v0',_'v2') + FlowB_('v2',_'v1') + FlowB_('v2',_'v4') = 0
_C18: - FlowA_('v1',_'v3') + FlowA_('v3',_'v4') + FlowA_('v3',_'v5')
 + FlowA_('v3',_'v6') = 0
_C19: - FlowB_('v1',_'v3') + FlowB_('v3',_'v4') + FlowB_('v3',_'v5')
 + FlowB_('v3',_'v6') = 0
_C20: - FlowA_('v2',_'v4') - FlowA_('v3',_'v4') + FlowA_('v4',_'v1')
 + FlowA_('v4',_'v5') + FlowA_('v4',_'v6') = 0
_C21: - FlowB_('v2',_'v4') - FlowB_('v3',_'v4') + FlowB_('v4',_'v1')
 + FlowB_('v4',_'v5') + FlowB_('v4',_'v6') = 0
_C22: FlowB_('v3',_'v5') + FlowB_('v4',_'v5') = 5
_C23: FlowA_('v3',_'v6') + FlowA_('v4',_'v6') = 6
VARIABLES
FlowA_('v0',_'v1') <= 10 Continuous
FlowA_('v0',_'v2') <= 10 Continuous
FlowA_('v1',_'v3') <= 10 Continuous
FlowA_('v2',_'v1') <= 10 Continuous
FlowA_('v2',_'v4') <= 10 Continuous
FlowA_('v3',_'v4') <= 10 Continuous
FlowA_('v3',_'v5') <= 10 Continuous
FlowA_('v3',_'v6') <= 10 Continuous
FlowA_('v4',_'v1') <= 10 Continuous
FlowA_('v4',_'v5') <= 10 Continuous
FlowA_('v4',_'v6') <= 10 Continuous
FlowB_('v0',_'v1') <= 10 Continuous
FlowB_('v0',_'v2') <= 10 Continuous
FlowB_('v1',_'v3') <= 10 Continuous
FlowB_('v2',_'v1') <= 10 Continuous
FlowB_('v2',_'v4') <= 10 Continuous
FlowB_('v3',_'v4') <= 10 Continuous
FlowB_('v3',_'v5') <= 10 Continuous
FlowB_('v3',_'v6') <= 10 Continuous
FlowB_('v4',_'v1') <= 10 Continuous
FlowB_('v4',_'v5') <= 10 Continuous
FlowB_('v4',_'v6') <= 10 Continuous

Optimization result
FlowA_('v0',_'v1')  =  6.0
FlowA_('v0',_'v2')  =  -0.0
FlowA_('v1',_'v3')  =  6.0
FlowA_('v2',_'v1')  =  0.0
FlowA_('v2',_'v4')  =  0.0
FlowA_('v3',_'v4')  =  5.0
FlowA_('v3',_'v5')  =  0.0
FlowA_('v3',_'v6')  =  1.0
FlowA_('v4',_'v1')  =  0.0
FlowA_('v4',_'v5')  =  0.0
FlowA_('v4',_'v6')  =  5.0
FlowB_('v0',_'v1')  =  1.0
FlowB_('v0',_'v2')  =  4.0
FlowB_('v1',_'v3')  =  1.0
FlowB_('v2',_'v1')  =  0.0
FlowB_('v2',_'v4')  =  4.0
FlowB_('v3',_'v4')  =  1.0
FlowB_('v3',_'v5')  =  0.0
FlowB_('v3',_'v6')  =  0.0
FlowB_('v4',_'v1')  =  0.0
FlowB_('v4',_'v5')  =  5.0
FlowB_('v4',_'v6')  =  0.0

最小运输费用 =  105.0

【本节完】

欢迎关注『Python小白的数学建模课 @ Youcans』原创作品

原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/118882174)。

Crated：2021-07-22

欢迎关注『Python小白的数学建模课 @ Youcans』系列，持续更新
Python小白的数学建模课-01.新手必读
 Python小白的数学建模课-02.数据导入
 Python小白的数学建模课-03.线性规划
 Python小白的数学建模课-04.整数规划
 Python小白的数学建模课-05.0-1规划
 Python小白的数学建模课-06.固定费用问题
 Python小白的数学建模课-07.选址问题
 Python小白的数学建模课-09.微分方程模型
 Python小白的数学建模课-10.微分方程边值问题
 Python小白的数学建模课-12.非线性规划
 Python小白的数学建模课-15.图论的基本概念
 Python小白的数学建模课-16.最短路径算法
 Python小白的数学建模课-17.条件最短路径算法
 Python小白的数学建模课-18.最小生成树问题
 Python小白的数学建模课-19.网络流优化问题
 Python小白的数学建模课-20.网络流优化案例
 Python小白的数学建模课-A1.国赛赛题类型分析
 Python小白的数学建模课-A2.2021年数维杯C题探讨
 Python小白的数学建模课-A3.12个新冠疫情数模竞赛赛题及短评
 Python小白的数学建模课-B2. 新冠疫情 SI模型
 Python小白的数学建模课-B3. 新冠疫情 SIS模型
 Python小白的数学建模课-B4. 新冠疫情 SIR模型
 Python小白的数学建模课-B5. 新冠疫情 SEIR模型
 Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型
 Python数模笔记-PuLP库
 Python数模笔记-StatsModels统计回归
 Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法