货运服务网络设计：经典文献阅读笔记（3）复现Netplan

最新推荐文章于 2024-04-15 14:51:03 发布

vapour3223

最新推荐文章于 2024-04-15 14:51:03 发布

阅读量549

点赞数 1

分类专栏：货运服务网络设计 Gurobi学习文章标签： python

本文链接：https://blog.csdn.net/vapour3223/article/details/116396762

版权

货运服务网络设计同时被 2 个专栏收录

3 篇文章 0 订阅

订阅专栏

Gurobi学习

3 篇文章 0 订阅

订阅专栏

本文介绍了针对零担货运服务网络设计的Netplan模型复现过程，基于Jacques Roy & Louis Delorme的1989年论文。模拟了货运需求、服务网络和相关成本，利用Python和Gurobipy构建优化模型，旨在最小化运输成本。通过生成不超过两个中转节点的货运服务路径，考虑了服务次数和货运量的平衡，最终求解得出服务频次和货运分配方案，展示了成本降低的效果。

摘要由CSDN通过智能技术生成

**货运服务网络设计：经典文献阅读笔记（2）**提到说要把Crainic T G（1984年）文献使用的模型复现一下，但是文章给出的通用框架还是太笼统，在尝试后决定使用Jacques Roy & Louis Delorme在Crainic T G（1984年）基础上提出的面向零担货运的Netplan（1989年）作为基础学习的蓝本来实现复现。

Netplan: A Network Optimization Model For Tactical Planning In The
Less-Than-Truckload Motor-Carrier Industry, INFOR: Information Systems and Operational Research, 27:1, 22-35, Jacques Roy & Louis Delorme (1989)
DOI: 10.1080/03155986.1989.11732079

一论文内容概述

加拿大公路货运市场管制放松，市场竞争激烈，零担货运服务企业谋求更加精细化地管理，从而促使了Netplan这一决策工具的产生。Netplan源于货运服务网络设计模型和算法的通用框架，融入了零担货运服务中长期规划的特点。

1.1 零担货运业务流程示意图

零担货运业务流程示意图

1.2 零担货运服务网络设计中三大问题

运载计划（Load Plan）：
如何制定站点之间的服务，在网络上运输（分开和合并）货物的策略，使得运输策略保证最小化运输和处理成本，同时改善或至少保持现有的服务水平——Powell（1986)

运载计划（Load Plan）必须考虑以下三个相互关联的子问题：

1.服务网络设计（Service Network Design），即选择提供运营商服务的路线，并确定每项服务的特征，如方式、速度和频率（每个规划期内给定路线上的车辆发车次数）
2.货运量分配（The Traffic Distribution），即对于具有正货运需求的每个出发地-目的地对（或货运类），将使用哪个服务序列来移动货物以及货物将通过哪个终端
3.空车移动（The Empty Vehicle Movement），即确定空车移动所需的数量，以在网络中每个始发终端的出发车辆和到达车辆数量之间建立平衡

第三个问题在论文中并没有详细说明，因此在这里也暂时忽略。

1.3 货运网络决策的影响范围

二复现论文模型

Netplan模型的主要任务是帮助零担货运企业在货运网络层面做出中长期的运载计划（Load Plan)。使用的语言是python，求解器使用的是gurobipy，还有绘图的networkx

import sys
import math
import random
import gurobipy as gp
from gurobipy import GRB
from gurobipy import *
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

2.1 输入数据模拟

在复现论文时遇到的第一个问题是没有数据，因此所有的数据都需要按照一定逻辑去模拟得到。那么下面先来看一下，会需要用到的数据有哪些？

2.1.1 输入数据分析

2.1.1.1 成本数据

城际运输成本（intercity transportation costs）：随着所使用的服务类型及其频率而变化
货运中转成本（the freight consolidation costs）：随着在散杂货站（Breakbulk）处理的货运量而变化
容量惩罚成本（the capacity penalty costs）：当每个运输服务中拖车容量的过度使用时发生的罚款
服务罚款成本（the service penalty costs）：当服务绩效标准未达到时产生的罚款，同时考虑每个货运类的计划运输时间的平均值和方差

以上列出了相关的成本以及成本单位，循序渐进，本文先只考虑第一个成本。

2.1.1.2 基础性信息

在这里先把需要的信息都列在这里，由于暂时先考虑一个成本，因此部分已知信息暂且不会用到（也就不需要在这次进行模拟生成）
在这里插入图片描述
约束信息

2.1.1.3 待决策信息

在这里插入图片描述

2.1.2 输入数据生成的代码

# 算例生成 for the article NETPLAN(1989)
# 节点集合,编号+坐标

N = dict()
# 每个节点的中转成本
Cn = dict()
for n in range(6):
    N.update({n:np.array([random.randint(0,20),random.randint(0,20)],dtype = np.float32)}) #编号+坐标
    Cn.update({n:random.randint(5,10)}) # 每个节点的中转成本

# 随机生成OD对的货运需求,暂时不考虑货运需求惩罚金与运输时间要求
D = dict()
for i in range(30):
    o = random.randint(0,len(N)-1)
    d = random.randint(0,len(N)-1)
    while o == d:
        d = random.randint(0,len(N)-1)
    D.update({(o,d):random.randint(30,200)})

我在这里模拟时，选择了需求推动的思路。先随机生成零担货运需求集合，货运需求由OD对和对应的货运量构成，采用了字典类型来存储。

货运服务的生成非常麻烦，为了简化模拟，所以仅用路（path）来定义不同的服务，即不同的弧构成的序列即为不同的服务。这样一来，需要制定路径搜索的规则，而用程序实现这些路径搜索规则前，需要确定网络的表示方式。

货运需求字典也可以看作是一个货运需求网络，是一个用弧来表示的有向网络。【ps，边一般指代无向的连接，弧指代有向的连接）

但是我觉得在搜索路径时用弧来表示的网络并不好用，因此我在代码中重新生成了snet 这一字典类型变量，它的表示方式是以节点为键，节点指向节点列表为值。

def find_degree(n, edges,out):
    '''
    在用边集表示的图中寻找从点n出发后到达的点
    或到达点n的出发点
    The Function is used to find the arriving node or 
    the starting node for the node n in the net which is 
    described by the edges set.
    thus parameter n indicates the node;
    edges is the net which is described by edges set;
    out indicates whether it is to find the arriving node (out =1) 
    or the starting node (out =0).
    '''
    degree_n = []
    if out == 0:
        for i in edges:
            if n == i[0]:
                degree_n.append(i[1])
    else:
        for i in edges:
            if n == i[1]:
                degree_n.append(i[0])
    return degree_n

# 生成需求网络
snet = dict()
for i in N.keys():
    deout = find_degree(i, D.keys(),0)
    snet.update({i:deout})

基于snet来生成单纯地由路径定义的服务

# 生成服务,服务经过的弧
L = dict()
#服务的容量
alpha = dict()
# 服务的运输单位成本
Cl = dict()
# 基于需求网络生成经过不超过2个中转节点的路径，然后随机生成服务容量以及服务成本
k = 0
for o in range(len(N)):
    for d in range(o+1,len(N)):
        for l in searchGraph(snet, o, d): # 使用的是不经过2个以上中转节点的路径
            # 满足坐标的三角形关系
            if len(l) == 3:
                l01 = ((N[l[0]][0] - N[l[1]][0])**2 + (N[l[0]][1] - N[l[1]][1])**2 )**0.5
                l02 = ((N[l[0]][0] - N[l[2]][0])**2 + (N[l[0]][1] - N[l[2]][1])**2 )**0.5
                l12 = ((N[l[1]][0] - N[l[2]][0])**2 + (N[l[1]][1] - N[l[2]][1])**2 )**0.5
                indicator = l01 + l12 > l02
                if indicator == True:
                    L.update({k:l}) 
                    Cl.update({k:random.randint(10,20)})
                    alpha.update({k:random.randint(2,4)*100})
                    k = k+1
            elif len(l) == 4:
                l01 = ((N[l[0]][0] - N[l[1]][0])**2 + (N[l[0]][1] - N[l[1]][1])**2 )**0.5
                l02 = ((N[l[0]][0] - N[l[2]][0])**2 + (N[l[0]][1] - N[l[2]][1])**2 )**0.5
                l03 = ((N[l[0]][0] - N[l[3]][0])**2 + (N[l[0]][1] - N[l[3]][1])**2 )**0.5
                l12 = ((N[l[1]][0] - N[l[2]][0])**2 + (N[l[1]][1] - N[l[2]][1])**2 )**0.5
                l13 = ((N[l[1]][0] - N[l[3]][0])**2 + (N[l[1]][1] - N[l[3]][1])**2 )**0.5
                l23 = ((N[l[2]][0] - N[l[3]][0])**2 + (N[l[2]][1] - N[l[3]][1])**2 )**0.5 
                l03 = ((N[l[0]][0] - N[l[3]][0])**2 + (N[l[0]][1] - N[l[3]][1])**2 )**0.5 
                
                ind1 = (l01 + l12) > l02
                ind2 = (l12 + l23) > l13
                ind3 = l03 > (l01+l12+l23)
                ind4
                ind5 = (l23 >= (0.3*(l01+l12+l23)))
                indicator1 = (ind1 & ind2 & ind3 & ind4)
                if indicator1 == True:
                    L.update({k:l}) 
                    Cl.update({k:random.randint(10,20)})
                    alpha.update({k:random.randint(2,4)*100})
                    k = k+1
                
            else:
                L.update({k:l}) 
                Cl.update({k:random.randint(10,20)})
                alpha.update({k:random.randint(2,4)*100})
                k = k+1

上面服务生成的代码虽然计算了路径的长度，也使得路径保持三角平衡，但是还是会生成看上去就很不合理的路径，但是因为只是模拟，所以就暂时放弃一下进一步细化路径生成的规则。
看起来就不合理的货运服务
看起来合理的货运服务
所以在路径生成上需要更多的规则来规范生成的服务。当然实际情况下，服务类型是货运企业给定，或者也是一个小的寻优问题（比如，不同步长下的最短路径）

至此基本的数据可以算是模拟好了，下面就是在这些数据上做进一步的计算。

2.2 模型数据分析

Netplan数学模型
学习要一步步来，本文就只考虑了使用每个服务的运输成本。一个约束是指，要在每个货运类m可以选择的服务集合中去分配货运量；第二个约束是指货运量不能为负，第三个约束是服务次数是有限的，可以不使用某一个服务，而且服务次数是整数。可以发现除了第一个约束外，其他的约束都是对决策变量取值范围的约束，这在编程中直接在定义决策变量的函数中就可以实现。
下面是Gurobipy构建模型的代码。

# 模型构建
Model = gp.Model('fsnd') 

# 添加变量

# 服务l的频率
F_l = Model.addVars(L,lb=0, vtype =GRB.INTEGER, name='Fre')

#每个货运流m的可选择路线k上的运输量
X_mk = dict()
for m in m_path:
    for k in m_path[m]:
        X_mk[m,k] = Model.addVar(lb=0,vtype =GRB.INTEGER,name=f'x_{m}{k}')
        
# 最小化运输成本    
objective = Model.setObjective(gp.quicksum(F_l[l] * Cl[l] for l in L),GRB.MINIMIZE)

# 添加约束
# 对任一货运类m在k条路径上的运输总量等于运输需求
i = 0
for m in m_path:
    Model.addConstr(gp.quicksum(X_mk[m,k] for k in m_path[m].keys())== D[m],name=f'C{i}')
    i = i + 1

# 对于任一服务的次数要大于所需的运输量
for l in L:
    Model.addConstr(gp.quicksum(X_mk[m,l] for m in l_mset[l])<=F_l[l]*alpha[l], name=f'C{i}')
    i = i + 1

Model.optimize()

# 模型错误检查
if Model.status == GRB.Status.INFEASIBLE:
        print('Optimization was stopped with status %d' % Model.status)
        # do IIS, find infeasible constraints
        Model.computeIIS()
        for c in Model.getConstrs():
            if c.IISConstr:
                print('%s' % c.constrName)

决策变量生成中，值得一提的是，其实每个货运类的备选服务数量是不同，在构建模型前需要将每个货运类m的备选服务集给筛选出来。

# 对每个货运需求，可选服务l集合
m_path = dict()
for m in D:
    path = dict()
    for l in L:
        o_ind = findindex(m[0],L[l])
        d_ind = findindex(m[1],L[l])
        if 0 <= o_ind < d_ind:
            path.update({l:L[l]})
    m_path.update({m:path})

这个模型的约束在实现的时候其实在目标函数的超载罚款成本上也去约束了服务次数 $\frac{X_l}{\alpha_l F_l}$ 。这里的 $X_l$ 是一个中间变量，因为我在这里暂时没有考虑超载的成本，因此我就在代码实现中增加了这个约束.

写成公式的话是
$\sum_{m \in M(l)}{X^m_k } <=F_l\alpha_l \quad \forall l$

$X_l = \sum_{m \in M(l)}{X^m_k } \quad \forall l$
$M (l)$ 是对于每个货运服务l适用的货运类m的集合。

下面是生成 $M (l)$ 的代码

# 对每个服务l,可能服务的货运类m集合
l_mset = dict()
for l in L:
    mdict = dict()
    for m in m_path:
        if l in m_path[m].keys():
            mdict.update({m:0})
    l_mset.update({l:mdict})

约束对应的代码

# 对于任一服务的次数要大于所需的运输量
for l in L:
    Model.addConstr(gp.quicksum(X_mk[m,l] for m in l_mset[l])<=F_l[l]*alpha[l], name=f'C{i}')
    i = i + 1

2.3 数据输出

数据输出的部分主要是对gurobi求解后的结果进行可视化。

这里放一个一次算例的log

Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 76 rows, 259 columns and 458 nonzeros
Model fingerprint: 0x6b2ec18d
Variable types: 0 continuous, 259 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e+00, 4e+02]
  Objective range  [1e+01, 2e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [8e+01, 2e+02]
Found heuristic solution: objective 218.0000000
Presolve time: 0.00s
Presolved: 76 rows, 259 columns, 458 nonzeros
Variable types: 0 continuous, 259 integer (17 binary)

Root relaxation: objective 6.456032e+01, 131 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   64.56032    0    9  218.00000   64.56032  70.4%     -    0s
H    0     0                     121.0000000   64.56032  46.6%     -    0s
H    0     0                     111.0000000   64.56032  41.8%     -    0s
H    0     0                      98.0000000   64.56032  34.1%     -    0s
H    0     0                      95.0000000   64.56032  32.0%     -    0s
H    0     0                      92.0000000   64.56032  29.8%     -    0s
     0     0   68.23180    0   21   92.00000   68.23180  25.8%     -    0s
H    0     0                      77.0000000   68.23180  11.4%     -    0s
     0     0   68.37794    0   19   77.00000   68.37794  11.2%     -    0s
     0     0   68.38825    0   25   77.00000   68.38825  11.2%     -    0s
     0     0   70.02606    0   18   77.00000   70.02606  9.06%     -    0s
     0     0   70.02606    0    8   77.00000   70.02606  9.06%     -    0s
     0     0   70.02606    0   11   77.00000   70.02606  9.06%     -    0s
     0     0   70.02606    0   21   77.00000   70.02606  9.06%     -    0s
     0     0   70.09949    0   18   77.00000   70.09949  8.96%     -    0s
     0     0   71.47247    0   31   77.00000   71.47247  7.18%     -    0s
     0     0   71.47286    0   31   77.00000   71.47286  7.18%     -    0s
     0     0   71.62734    0   18   77.00000   71.62734  6.98%     -    0s
     0     0   71.66762    0   38   77.00000   71.66762  6.93%     -    0s
     0     0   71.68363    0   38   77.00000   71.68363  6.90%     -    0s
H    0     0                      76.0000000   71.68363  5.68%     -    0s
     0     0   71.71158    0   15   76.00000   71.71158  5.64%     -    0s
     0     0   71.71158    0   19   76.00000   71.71158  5.64%     -    0s
     0     0   71.74907    0   45   76.00000   71.74907  5.59%     -    0s
     0     0   71.81435    0   28   76.00000   71.81435  5.51%     -    0s
H    0     0                      74.0000000   71.81435  2.95%     -    0s
     0     0   72.17717    0   24   74.00000   72.17717  2.46%     -    0s
     0     0   72.17717    0   19   74.00000   72.17717  2.46%     -    0s
     0     0   72.17717    0    8   74.00000   72.17717  2.46%     -    0s
     0     0   72.17717    0   17   74.00000   72.17717  2.46%     -    0s
     0     0   72.63333    0   22   74.00000   72.63333  1.85%     -    0s
     0     0   73.18750    0   26   74.00000   73.18750  1.10%     -    0s

Cutting planes:
  Gomory: 3
  MIR: 13
  Relax-and-lift: 2

Explored 1 nodes (554 simplex iterations) in 0.13 seconds
Thread count was 8 (of 8 available processors)

Solution count 9: 74 76 77 ... 218

Optimal solution found (tolerance 1.00e-04)
Best objective 7.400000000000e+01, best bound 7.400000000000e+01, gap 0.0000%

输出的结果

货运类(2, 0)_58：81.0:[4, 2, 0, 5]
货运类(1, 4)_15：41.0:[1, 0, 4, 2]
货运类(1, 4)_19：39.0:[1, 4, 2, 3]
货运类(1, 4)_30：109.0:[1, 3, 4, 5]
货运类(4, 2)_15：189.0:[1, 0, 4, 2]
货运类(1, 0)_15：149.0:[1, 0, 4, 2]
货运类(1, 3)_19：94.0:[1, 4, 2, 3]
货运类(5, 2)_4：100.0:[0, 4, 5, 2]
货运类(2, 5)_54：78.0:[3, 2, 1, 5]
货运类(0, 4)_4：188.0:[0, 4, 5, 2]
货运类(2, 3)_19：167.0:[1, 4, 2, 3]
货运类(4, 5)_4：64.0:[0, 4, 5, 2]
货运类(4, 5)_30：23.0:[1, 3, 4, 5]
货运类(0, 5)_4：48.0:[0, 4, 5, 2]
货运类(0, 5)_58：58.0:[4, 2, 0, 5]
货运类(4, 0)_58：161.0:[4, 2, 0, 5]
货运类(1, 5)_30：73.0:[1, 3, 4, 5]
货运类(1, 5)_54：54.0:[3, 2, 1, 5]
货运类(3, 4)_30：195.0:[1, 3, 4, 5]
货运类(2, 1)_54：88.0:[3, 2, 1, 5]
货运类(3, 2)_54：180.0:[3, 2, 1, 5]
服务编号：4, 服务路径：[0, 4, 5, 2], 服务容量：400,服务次数：1.0, 服务货运量：49
服务编号：15, 服务路径：[1, 0, 4, 2], 服务容量：400,服务次数：1.0, 服务货运量：53
服务编号：19, 服务路径：[1, 4, 2, 3], 服务容量：300,服务次数：1.0, 服务货运量：45
服务编号：30, 服务路径：[1, 3, 4, 5], 服务容量：400,服务次数：1.0, 服务货运量：61
服务编号：54, 服务路径：[3, 2, 1, 5], 服务容量：400,服务次数：1.0, 服务货运量：48
服务编号：58, 服务路径：[4, 2, 0, 5], 服务容量：300,服务次数：1.0, 服务货运量：35