Python小白的数学建模课20。网络流优化案例
在实际工作和数模竞赛中,网络最大流问题、最小费用流问题和最小费用最大流问题都有很多延伸和推广的应用。本文介绍了常见的网络最大流问题、最小费用流问题的应用与推广,及其解题方法。本文选择多源多汇物流转运问题、多商品流问题案例,详细介绍网络流问题的分析方法和解决方案,并使用线性规划方法建模和编程。1. 网络流问题的应用与推广
流在生活中十分常见,例如交通系统中的人流、车流、物流,供水管网中的水流,金融系统中的现金流,网络中的信息流。网络流优化问题是基本的网络优化问题,应用非常广泛,遍及通讯、运输、电力、工程规划、任务分派、设备更新以及计算机辅助设计等领域。
网络流优化问题最重要的指标是每条边的成本和容量限制,既要考虑成本最低(最短路径问题),又要满足容量限制(最大流问题),由此产生了网络最大流问题、最小费用流问题和最小费用最大流问题。在《Python小白的数学建模课-19.网络流优化》一文中,我们详细介绍了网络最大流问题、最小费用流问题和最小费用最大流问题的算法、编程和案例。
不论在实际工作中,还是在数模竞赛中,网络最大流问题、最小费用流问题和最小费用最大流问题都有很多延伸和推广的应用。小白同学对于网络图的概念和思维逻辑比较生疏,面对变换和深化的题型往往对题目能有基本判断和思路、但找不到具体的解决方法,觉得应该能做出来但就是差一点。
本文选择多源多汇物流转运问题、多商品流问题案例,详细介绍网络流问题的分析方法和解决方案,并使用线性规划方法建模和编程。
对于小白同学来说,这两个案例的建模和编程都有些复杂。别急,也别怕,问题分析和程序说明都很详细,只要静下来慢慢读,将程序说明、源程序和运行结果对照着读,相信小白同学也能看懂。2. 网络最大流问题的应用与推广2.1 指定流量的可行流
我们在《19.网络流优化》一文讨论了网络最大流的算法。如果不是计算最大流,而是要计算指定流量 v 的可行流,可以用如下图所示的增广路(相当于辅助线)方法处理。
(1)在容量网络 G 上增加一个顶点 t’ 和一条边 (t,t’),构建一个新的容量网络 G’,令边 (t,t’) 的容量为 v;
(2)对于容量网络 G‘,计算从源点 s 到汇点 t’ 的最大流 f m a x 显然 f m a x v ;
(3)若 v < f m a x ,则容量网络 G 不存在满足流量 v 的可行流;若 v = f m a x ,则从最大流 f m a x的路径中去掉边 (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() 函数求出最小费用流。
该问题也是一个典型的线性规划问题,可以用线性规划方法求解。
式中: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 # Demo21 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(" 最小运输费用 = ", 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。所有线路的商品总流量不超过边的容量限制。
该问题也是一个典型的线性规划问题,可以用如下模型描述:
式中: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) 的运输费用之和 m i , j E w i j f i j 。多商品流的模型主要包含流守恒约束和容量限制。
本案例的问题中,不同商品 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)定义约束条件:
(4)定义约束条件:
多商品流的模型主要包含边的容量约束和顶点的流守恒约束和两个方面。
1.边的容量约束,即每条边上各种商品流量之和不大于边的容量。
2.顶点的流守恒约束,分为源点(供应点)、汇点(需求点)和中间节点三种情况:通过 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 # Demo21 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:{} Supply 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(" 最小运输费用 = ", 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
【本节完】
————————————————版权声明:本文为CSDN博主「youcans」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
原文链接:https://blog.csdn.net/youcans/article/details/118882174