OD需求和连通可靠性不确定下的连续航运网络设计方法

OD需求和连通可靠性不确定下的连续航运网络设计方法

本文针对航运网络中由于供需的不确定性会对航运网络规划带来影响的问题，在剖析已有研究基础上提出了一种考虑OD需求不确定和连通可靠性不确定的航运网络设计方法，该问题本质上属于连续交通网络设计范畴。本文首先给出不确定网络连通可靠性的量化方法，构建考虑两类不确定因素的连续航运网络设计双层规划模型，并基于鲁棒优化的对偶理论对双层规划模型进行了转化处理使其容易求解。然后，针对该双层规划模型，设计了一种基于松弛算法（Relaxation Algorithm，RA）和带精英保留策略的快速非支配排序遗传算法（A Fast and Elitist multiobjective Genetic Algorithm，NSGA-II）进行求解。最后，本文通过以中国原油进口航运网络为案例分析对象，给出两种优化结果下的航运网络设计方案，并进行算法性能对比实验，验证本文提出模型的科学性和算法的优越性。结果表明了在航运网络设计问题中考虑需求和连通可靠性不确定的求解结果更加贴近实际，本文研究成果可以为中国原油进口航运网络设计提供政策性指导意义。

交通网络设计；连续航运网络设计；供需不确定；连通可靠性；双层规划模型；松弛算法和带精英保留策略的快速非支配排序遗传算法（R-NSGA-II）

当前的全球供应链网络体系正在面临诸多风险不确定性，这将直接导致未来规划和构建运输网络、优化运输通道时，必须要有新的视角和新的方法。2020年以来，受新冠肺炎疫情的影响，国际局势风云突变，出现百年未有之大变局，中美贸易摩擦加剧，国际运输通道需求受限，20201年，长荣海运集装箱“长赐号”因机舱突然失电，横向堵塞苏伊士运河，造成了40年以来该运河最严重的中断事故，该事故令全球贸易每周损失60亿美元至100亿美元，从中国到欧洲的集装箱价格比一年前高出近3倍[1]，国际海上运输通道供给能力失衡。因此，从应对非常规突发事件角度思考，现有的运输网络规划理论与方法对运营风险不确定性考虑不足，传统以网络运输效率最高、运输能力最大、运输成本最低的规划思想需要拓展创新。不可忽视的是，需要在新的综合运输网络规划中考虑应急管理和风险防控因素，以推动全球供应链体系发展行稳致远。本文研究的重点聚焦航运网络设计问题，通常的航运网络设计问题普遍都是考虑运输成本最小、运输时间最短、污染物排放最小、服务可靠性最高为单目标或多目标进行研究[2-4]，忽略了航运网络设计中的不确定性，比如节点间OD需求、连通可靠性、航运网络建设成本等因素的不确定，并且鲜有研究将网络配流均衡考虑在模型构建中，本文就是将不确定因素和网络配流均衡同时考虑在航运网络设计问题中。总结来看，本文的创新点和贡献主要集中在三个方面：（1）提出了一种量化航运网络连通可靠性不确定性的表征方法，构建了考虑OD需求和连通可靠性不确定的连续航运网络设计双层规划模型，采用鲁棒优化中的椭球集合来代表不确定需求，并基于鲁棒优化对偶理论对模型进行转化处理；（2）创新了一种求解连续航运网络设计双层规划模型的现代启发式算法，该算法结合了松弛算法和快速非支配排序遗传算法，对航运网络设计双层规划模型求解具有很好的适用性，比一般现代启发式算法和求解器的求解性能更优越；（3）将本文模型和算法应用在了中国原油进口航运网络设计实际案例分析中，通过该案例分析进一步论证了模型和算法的科学合理性，并为中国原油进口航运网络规划提供政策性指导建议。

本文所研究的航运网络设计问题其本质上是属于交通网络设计问题，交通网络设计问题（Network Design Problem, NDP）是交通学科研究的热点问题之一[5]，研究通常将NDP问题分为三种类型：连续交通网络设计问题（Continuous Network Design Problem, CNDP），离散交通网络设计问题（Discrete Network Design Problem, DNDP）和混合交通网络设计问题（Mixed Network Design Problem, MNDP）。它们三者之间的区别在于，连续交通网络设计问题主要研究改善现有路段通行能力，决策变量连续，比如，在进行道路扩建时，由于道路扩建所新增的路段通行能力即为连续变量；离散交通网络设计问题主要研究在已有网络中新建或改扩建路段，决策变量离散，决策变量为0-1变量；混合交通网络设计主要同时采用在已有网络中新建路段和改善已有路段通行能力两种手段，决策变量部分连续部分离散[6]。1973年Morlok首次提出交通网络设计问题[7]，四十多年以来，学者们在交通网络设计问题研究中取得较为丰富和成熟的研究成果。需要说明的是，本文研究的航运网络设计其实是连续交通网络设计问题，不考虑在网络中新增航道，在航运网络规划中一般都是考虑对原有航道进行扩容加宽处理以提高通航能力，如若新增一条航道不仅投资成本不可估计，更主要的是，不同于道路交通网络设计中新建道路后即可投入使用，新增航道能否开通运营需要开展很长周期的试验，因此为了符合大多实际航运网络规划，本研究选择以连续航运网络设计为研究对象。

连续交通网络设计适用范围包括道路改建、信号控制、匝道合并、车速限制等问题，国内外关于此类问题研究较多。Fitzpatrick等[8]属于较早一批研究连续交通网络设计的代表学者，他们研究了道路设计中的设计速度和限制速度的关系，提出了相关车速管理政策；Yang等[9]构建了一种用于缓解交通拥堵问题的双层规划模型，并提出用启发式算法求解模型，后来研究共识性地认为交通网络设计问题核心之一是构建双层规划模型；Meng等[10]针对连续交通网络设计双层规划模型中求解算法较少的问题，提出了一种基于局部收敛的增广拉格朗日算法，并验证了该算法求解的优越性；Tsao和Linh[11]将环境影响考虑进连续交通网络设计中，研究设计了一个考虑碳排放最小和系统总成本最小的海陆港网络；Wang和Szeto[12]将出行者出行的风险规避和路径选择偏好行为引入到网络设计建模中，构建了一个基于网络可靠性的连续交通网络设计双层规划模型。国内学者高自有等[13]详细综述城市道路中的连续交通网络设计建模方法，并进一步根据双层规划模型中下层用户均衡分配模型特征，设计了一种基于相继平均算法的全局收敛算法，避免模型结果陷入局部最优，并简化了启发式算法求解流程[14]；孙华等[15]研究OD需求不确定下的连续交通网络设计鲁棒优化问题，构建了一个带均衡约束的双层规划数学模型；孙飞和龙建成[16]提出了一个连续交通网络设计双目标双层规划模型，以网络总投资成本和尾气排放最小为目标，并设计非支配排序遗传算法求解模型；Han等[17]提出了一种将路径引导和连续交通网络设计相结合的流量优化模型，并设计改进模拟退火算法求解；占莉莎[18]提出了一个考虑道路交通拥挤收费策略的连续交通网络双层规划模型，结合免疫算法和网络配流理论进行求解。总的来讲，关于连续交通网络设计问题的研究国内外大多聚焦优化道路车辆行驶速度、控制道路交通流量、制定交通管理措施、保护交通生态环境等角度，通过改进模型目标、增加约束条件以及优化求解算法来拓展此类问题的研究。上述研究可以为本文航运网络设计提供理论和方法借鉴，但是这些研究大多忽略了交通网络设计中的不确定问题，一旦网络中存在不确定因素，则会对网络整体规划和设计带来影响，这是本文需要重点考虑的。

事实上，随着人们对交通网络设计问题研究的不断深入，交通系统的不确定性以及交通网络的可靠性早已引起了运输研究学者的极大关注。在实际网络中确实存在大量不确定因素，如需求层面的不确定性、供给层面的不确定性以及人们在出行行为方面的不确定性，忽视交通系统的不确定性将会给道路网络规划建设带来风险，比如可能出现一些次优的投资决策或导致交通拥堵问题的进一步恶化。本研究重点关注的是交通网络设计中的需求不确定情形和连通可靠性不确定情形。其中，关于交通网络不确定问题研究，最具有代表性的学者是Ukkusuri等[19]在最早的研究中将始发地和目的地出行OD矩阵作为概率分布已知的随机变量，提出了一种适应各类需求的鲁棒网络设计模型，引入稳健性来衡量网络稳定程度，但是假定需求量呈一定概率分布，很难贴近实际应用，因此他们将需求变化分成多个时间阶段，研究单阶段的网络设计决策，将问题进行灵活处理[20]。然而，这样研究划分时间阶段的主观性较强，且假定单阶段需求固定不变，在最新的研究中，他们直接将网络设计问题认为是柔性网络设计问题，网络中的需求、运力和成本都是不确定的，提出了一种考虑需求随机性和需求弹性的多级柔性网络设计方法，其本质是带平衡约束的随机数学规划问题[21]；国内学者陆化普等[22,23]研究以OD需求不确性为基本前提，以随机双层规划和均值方差理论为基础，建立OD需求不确定的离散交通网络设计模型；进一步的，他们考虑到交通规划中供给需求及决策准则的不确定性，提出了一种基于不确定理论的离散交通网络设计模型。研究假定供给和需求为随机变量，使用半标准差和风险价值来度量不确定性，建立了考虑决策者风险态度的多目标双层优化模型[24]。Cao等[25]研究了运输成本和需求不确定的连续交通网络设计问题，提出了一个基于鲁棒优化的双层规划模型；Wang和Xie[26]在研究中既考虑了需求动态日变化性，又考虑了需求概率分布和参数估计的不确定性，提出了一种不确定需求变异性的连续交通网络设计模型；Yang[27]提出了一个两阶段随机规划模型来衡量航空需求呈季节性变化的情形，第一阶段是枢纽选址决策，第二阶段是根据第一阶段的枢纽位置和不确定需求确定航线和流量分配。上述研究在刻画不确定因素中都用到了鲁棒优化和不确定理论的思想，可以为本文刻画不确定连通可靠性和不确定OD需求提供理论支撑，但航运网络设计还是有别于交通网络设计，交通网络设计中重点关注的是道路通行能力提升和道路网络中节点连通可靠性，航运网络具有比道路网络更加复杂的拓扑结构，在航运网络中并不是每个节点连通的航线都能提升通航能力，航运网络中的通航能力提升主要指的是航道扩容，但还要考虑连接航道的其他航线，加之，当考虑网络配流均衡后使得此类问题研究将会变得更加复杂。

遗憾的是，几乎没有研究在航运网络设计中考虑网络配流问题和多目标不确定问题，关于航运网络设计，赵宇哲等[28]考虑OD需求不确定构建基于轴辐理论的集装箱海运网络设计模型，借助机会约束规划将模型转化为了等价确定模型，在此基础上他们又以路径连通为决策变量，构建考虑竞争环境下的航运网络设计模型[29]；张祺等[30]以集装箱海运网络总成本最小为目标，考虑港口通过能力约束，应用混沌优化方法结合混合优化算法进行求解；朱小林等[31]针对建立区域整体贸易运输网络，构建多枢纽单分配多级轴辐式航运网络设计模型，研究枢纽选址对航运网络可靠性的影响。上述研究都是从货运网络规划和设计角度进行研究[32]，几乎没有从交通视角进行航运网络设计，因此，这是本文需要重点聚焦的研究视角。

综上所述，对已有文献总结来看，目前已有研究对支撑本文连续航运网络设计问题还存在三个方面的局限性：（1）虽然连续交通网络设计问题研究已经非常成熟和应用广泛，但是连续航运网络设计问题不完全等同于连续交通网络设计问题，其中的建模方法、约束条件和求解算法还需要进一步改进和创新，以符合本研究背景；（2）在交通网络设计问题中考虑运输需求、建设成本、服务可靠性等的不确定性研究也比较广泛，但是鲜有同时考虑多个不确定因素的研究，而在我们本研究航运网络设计建模问题中需要考虑多目标多因素不确定，因此复杂性相对较高；（3）求解算法实现难度较大，对于双层规划模型的求解算法本身难以实现，尤其在双层规划模型中既考虑基于鲁棒优化的不确定因素量化和网络配流问题，算法设计比较复杂，当前在交通网络设计中求解双层规划模型的算法研究存在瓶颈。

本文研究重点考虑OD对间需求不确定和港口（节点）连通可靠性不确定的连续航运网络设计方法，鲁棒优化方法能够较好地刻画数学建模问题中的不确定因素，借鉴交通网络设计理论中的连续交通网络设计方法可以构建连续航运网络设计模型，因此本节围绕鲁棒优化和连续交通网络设计进行基本理论分析。

连通交通网络设计（CNDP）顾名思义就是CNDP模型中的决策变量是连续的，这些连续变量主要包括了道路通行能力改善[33]、道路收费政策[34]、信号配时控制[35]、匝道控制限流[9]等。基于上述文献综述的讨论，简单描述连续交通网络设计问题的数学模型，连续交通网络设计模型是一个双层规划数学模型，上层模型是从政策制定决策者角度考虑改善路段通行能力条件，以提高路段通过能力；下层模型是从路网使用者角度，即用户（出行者）层面考虑出行行为下的网络配流特征，可以对上层模型起到反馈作用。假定存在一个交通网络 ，其中 为网络节点集合； 为网络边的集合； 为网络中的OD对，通过不同的上层优化目标选择和下层网络配流模型组合，可以得到不同的连续交通网络设计模型，这里选择连续交通网络设计一般建模方法，上层规划模型以网络总阻抗值最小，下层规划模型采用基于用户均衡的网络配流模型，则CDNP问题的数学模型可以写成：

                                  (1)

                                             (2)

                                             (3)

其中， 是关于 的隐函数，由下层模型决定：

                                   (4)

                                              (5)

                                         (6)

                                            (7)

上述连续交通网络设计双层规划模型中， 表示路段 上新增加的通行能力， 表示路段 上的走行时间， 为OD对 之间路径 上的流量， 表示OD对 之间的路径集合， 表示路段 上的流量， 表示如果OD对 之间的路段 在路径 上，那么 ，否则为0。 为OD对 之间的总运输需求量， 为路段 的投资函数， 是总投资预算。上层模型目标函数式（1）表示的是最小化网络中总阻抗值，约束（2）为投资预算约束，约束（3）为新增通行能力非负约束。下层模型目标函数式（4）表示用户均衡的目标函数，约束（5）为流量守恒约束，表示当网络达到用户均衡状态时，OD对 之间的总需求量等于该OD对之间所有路径的流量总和，约束（6）表示路段流量和路径流量之间的关系，当网络达到均衡状态时，路段 上的流量等于所有经过该路段的路径 流量总和，约束（7）为路径流量的非负约束。

在鲁棒优化问题中，对应于不确定优化问题的鲁棒模型被称为鲁棒对等模型，它通常是一个半无限规划，鲁棒对等模型的可行解和最优解被称为不确定优化问题的鲁棒可行解和鲁棒最优解。鲁棒优化按照对不确定因素的描述方法大致可以分为：基于情景的鲁棒优化和基于集合的鲁棒优化。基于情景的鲁棒优化就是通过有限多个以一定概率实现的离散片段来描述优化模型中的不确定参数，找到一个近似最优解，通常用期望-方差模型作为目标函数。基于集合的鲁棒优化方法就是假定不确定参数在一个有界闭集内变化，找到一个解，使得它对不确定参数在这个有界闭集内的所有取值可行。需要说明的是，本文研究考虑不确定因素为需求和连通可靠性，且不确定需求下不考虑需求服从随机概率分布，与基于情景的鲁棒优化无关，因此本节主要介绍基于集合的鲁棒优化方法。为了便于介绍基于集合的鲁棒优化方法，给出如下不确定优化问题。

                                               (8)

                                               (9)

为目标函数 的约束函数，参数 为不确定参数，基于集合的鲁棒优化方法不同于基于情景的鲁棒优化方法，不需要知道不确定参数 的概率分布，而是假定不确定参数 属于一个有界闭集 ，基于集合的鲁棒优化目的是求得这样一个解，对于不确定参数 的任何取值，约束条件均满足，并且使得最坏情况下的目标函数值达到最优，参数 是决策变量。其核心思想是将原始问题以一定近似程度转化成一个具有多项式计算复杂度的凸优化问题，关键是建立鲁棒对等模型。通常鲁棒对等模型（Robust Counterparts, RC）可以写成：

                                     (10)

                                          (11)

它等价于如下模型：

                                        (12)

                                                (13)                                         

                                          (14)

从上述模型可以看出，鲁棒对等问题是一个半无限规划，一般很难求。为了求解鲁棒对等模型，通常是利用对偶理论将其等价转化为可求解的鲁棒对等模型。基于集合的鲁棒优化模型存在一个默许的假设，就是所有决策变量在不确定参数被观测到之前必须做出决策[36,37]。然而，现实中有些优化问题可能不满足这个假设，一部分决策变量能够在不确定参数被观测到之后再进行调整。针对这种情况，Ben-Tal等[38]研究了线性规划模型的可调整鲁棒等价模型（Adjustable Robust Counterparts, ARC），Takeda等[39]进一步将ARC推广到非线性规划上去：

                                        (15)

 

这里的决策变量 在不确定参数 被观测到之前必须做出决策，称为不可调整的决策变量，而变量 可以在不确定参数 被观测到之后再进行调整，称为可调整决策变量。显然的，可调整的鲁棒对等模型可行域要大于原鲁棒对等模型可行域，因此可调整的鲁棒对等模型更加灵活，但是可调整的鲁棒对等模型求解很困难，从下列ARC等价模型可以看出。

                                         (16)

目标函数式（16）描述的是一个极小极大极小的问题，但是几乎没有研究方法能够求解此类问题，Ben-Tal等[40,41]通过限制可调整变量 和不确定参数 之间存在映射关系，即 ，这里的参数 和 是不可调整变量。于是，可调整鲁棒对等模型（Affinely Adjustable Robust Counterparts, AARC）可以写为：

                                          (17)

                                     (18)

                                       (19)

通过这样的模型变化目的是为了便于求解，Ben-Tal已经证明了上述AARC模型可以通过等价转化成比较容易处理的优化问题，如线性规划、线性半定规划问题等。

本文研究不考虑新建或新增航道，仅考虑对现有航道进行改扩建从而使得航运网络系统性能最优，借鉴连续交通网络设计方法构建连续航运网络设计双层规划模型，对于连续交通网络设计问题通常研究构建上层规划模型以总阻抗值最小为目标[12,42]，本研究为了更好描述航运网络中的不确定风险因素，将港口节点间连通可靠性引入上层规划模型。构建的上层规划模型是以系统总走行时间（总阻抗值）最小和连通可靠性最大的双目标，下层规划模型则是基于用户均衡的交通配流模型。针对连续航运网络设计双层规划建模问题，本文首先研究基于不确定理论的连通可靠性量化，然后构建考虑不确定连通可靠性的双目标上层规划模型和考虑不确定需求的下层交通配流模型，并基于鲁棒优化对偶理论进行转化处理，最后推导出连续航运网络设计可求解的双层规划模型。

为了便于本文建模，将模型参数符号定义和说明梳理如下表1所示。

表1 模型参数符号定义

参数名称	参数说明
	在一个航运网络 中， 为航运网络节点集合， 为航道（边）集合
	航道 上通行的交通量
	航道 上新增交通量
	航道 上的最大通航能力
	航道改扩建总投资建设成本预算规模
	不确定需求集合
	航道 上密度为0时的班轮自由流行驶时间
	BPR函数中的固定参数
	航道 扩建成本的常数值
	航运网络的连通可靠性
	航运网络中路径的连通可靠性
	航运网络中的航线路径
	航运网络中的航线路径数量
	航运网络中OD对 之间的路径 上的走行时间函数
	OD对 间的最小走行时间
	OD对 之间的总需求量
	OD对 间路径 上的流量
	OD对 之间的路径集合
	若航道 在路径 上，则取值为1，否则为0

航运网络的连通可靠性是指OD对之间至少有一条路径连通的概率，如果所有路径都不连通，则该航运网络被认为是不连通的，为了描述网络的连通状态，引入下面的函数：

                                 (20)

借鉴吕靖等[43]研究基于网络最小不交回路的连通可靠性度量方法，航运网络的连通可靠性 可以通过加总所有路径可靠性来计算：

                                     (21)

每条路径连通可靠性可以根据路径构成的节点连通可靠性进行计算，则路径连通可靠性可以表示为：

                                         (22)

式中， 为第 个节点的连通可靠性，是不确定变量， 为构成路径 的节点个数。期望值是用来衡量不确定变量的平均值，表示不确定变量的大小，本文采用 的期望值来表示其可靠性的大小。根据不确定理论中的定理：设 是不确定变量并具有不确定分布 ，如果 的期望值存在，则计算得到第 个节点的连通可靠性 ，进而计算路径的连通可靠性：

                                         (23)

式中， 称为 的逆不确定分布， 为置信水平。如果不确定变量 有不确定分布：

                            (24)

则称不确定变量 为Zigzag不确定变量，记为 ，其中 为实数，且 。

Zigzag分布式正则的，其逆不确定分布为：

                               (25)

将式（22）~式（25）代入式（21）即可得到航运网络的连通可靠性。

上层模型是一个考虑不确定连通可靠性的双目标优化问题。上层规划模型如下所示：

                                                 (26)

                                                   (27)

                                                      (28)

                                          (29)

                                                      (30)

                                                    (31)

                               (32)

                                        (33)

                                         (34)

                                                 (35)

                                              (36)

                                                 (37)

上层规划模型中，式（26）和式（27）表示以系统总走行时间（系统总阻抗值）最小和连通可靠性最大为双目标函数；式（28）是投资建设成本约束；式（29）假设航道 上的走行时间服从BPR函数，在该BPR函数中考虑了航道新增通行能力 ；式（30）表示BPR函数中的参数修订，借鉴已有文献[44,45]，认为用参数 修订BPR函数更符合中国国情；式（31）是投资建设成本的函数表达式；式（32）是基于不确定理论的连通可靠性函数表达式；式（33）是对航道 原始通航能力和新增通航能力的非负约束；式（34）表示OD对 之间关于路径流量和航道通航能力增加量之间的路径走行时间函数满足航道总走行时间之和；式（35）和式（36）表示路径走行时间和OD对 之间的最小走行时间之间关系；式（37）表示最小化航运网络的走行时间约束。

下层规划模型是一个考虑不确定需求的交通配流问题。下层规划模型如下所示：

                                            (38)

                                                   (39)

                                             (40)

                                                  (41)

下层规划模型中，式（38）是基于用户均衡的交通配流目标函数表达式；式（39）表示当航运网络达到均衡状态时，网络中的总需求量与路径流量之间关系；式（40）表示当航运网络达到均衡状态时，航道上的流量与路径流量之间的关系；式（41）是对路径流量的非负约束条件。

从上文可以看出，连续航运网络设计双层规划模型实际上是一个半无限规划，一般来讲该类问题是一个非凸的并且在任何可行点均不满足约束要求，所以求解起来非常困难，在此类问题求解中最具有代表性的就是Yin和Lawphonganich[46]提出了一种类似割平面的需求生成算法来求解半无限规划模型，他们证明了在合适条件下，算法能够收敛得到全局最优解，但是该算法要去条件局限性太强，且求得的解过于保守。本节研究基于对偶理论转化鲁棒对等模型。根据Ben-Tal和Nemirovski等[38,47]的研究，将4.3中的双层规划模型写成RC（上层模型）和RC（下层模型）：

RC（上层规划模型）：

                         (42)

                                  (43)

                                             (44)

                                      (45)

                                 (46)

                                     (47)

                                (48)

RC（下层规划模型）：

                                      (49)

                                          (50)

                                   (51)

                                         (52)

本文借鉴Ben-Tal的经典研究“鲁棒优化不确定转化”采用如下形式的椭球集合 来代表不确定需求集合：

                            (53)

式（53）中的 为OD对 的“名义”需求， 为OD对 的需求可能偏差， 为不确定参数向量， 为不确定度，本研究中的不确定度的取值假设为 。则根据Houska和Diehl[48]提出的基于拉格朗日算子的鲁棒对等模型偏导转化为近似ARC模型：

ARC（上层规划模型）：

                                (54)

                                     (55)

                                                 (56)

                        (57)

                                         (58)

                      (59)

                                        (60)

                           (61)

                                          (62)

                                           (63)

                                                (64)

ARC（下层规划模型）：

                                         (65)

                                  (66)

                                        (67)

                             (68)

                                        (69)

                                    (70)

                                          (71)

                                               (72)

上面模型中的 是在约束（45）~（47）和约束（50）~（52）下，式（53）所对应的拉格朗日乘子。研究需要进一步对上面模型进行可解性转化，因此对上述偏导函数式还可以进行如下转化为AARC模型。

                         (73)

                 (74)

          (75)

                            (76)

                      (77)

                                (78)

                               (79)

                                 (80)

从式（73）~式（80）可以看出，上层模型和下层模型的关键是航道流量 ，路径流量 和OD对 之间的最小走行时间 对网络总需求 的偏导数，即 ，由于总需求 是不确定的，将椭球集合 中的不确定需求 表达式代入上述三个偏导数中，则有：

                    (81)

根据式（81）可以看出，只需要确定不确定参数 即可，本文研究假定不确定参数取值为 ，且“名义”需求 和偏差需求 之间的关系满足 。

综合上文模型推导研究，给出可用于求解的连续航运网络设计双层规划模型。

上层规划模型：

                                       (82)

                                            (83)

  

下层规划模型：

                                             (84)

其中，上层规划模型中的约束条件式（56）~式（64）和下层规划模型中的约束条件式（66）~式（72）在求解中等效转化为式（73）~式（81）中的可求解式进行求解。

本文构建的是一个多目标非线性规划模型，显然，各目标函数各自的最优解并不能满足多目标的整体最优，也不存在可以同时使所有目标达到最优的解决方案，这是一个NP难问题。目前解决多目标规划问题的方法主要有两类：经典的解决方法是将多目标统一量纲后赋予权重转化为单个目标问题进行求解，此外，并非每次求解都能得到唯一解的结果[49,50]；另一种是采用支配关系对目标向量部分排序，多目标规划一般会产生由许多次优解组成的决策组合，针对所有的目标来说，决策组合所包含的解释不能比较优劣的，其特点是不能进一步满足任一目标优化程度而不减小其他任一目标的优化程度，一般称此类解为帕累托最优解[51,52]。一些多目标优化算法被相继提出来，其能够在单次迭代后找到多个非支配解。带精英保留策略的快速非支配排序遗传算法（NSGA-II）是目前在求解多目标优化问题中应用最为广泛的多目标遗传算法之一，具有算法计算复杂程度相对较低，帕累托最优解收敛效果较好，易于编程实现，具有潜在的并行性和分布式等优点[53]。在本文研究的连续航运网络设计问题本质上属于交通网络设计问题，该算法也被学者广泛采用[54-56]。但是仅用NSGA-II算法是无法完成本文构建考虑不确定需求和连通可靠性的双层规划模型，最主要的问题在于需求量无法准确输入导致不能通过下层交通配流模型得到航道流量，从而优化上层目标函数。因此，本文提出一种松弛算法（Relaxation Algorithm，RA）和带精英保留策略的快速非支配排序遗传算法（A Fast and Elitist multiobjective Genetic Algorithm，NSGA-II）的求解方法，下文称之为R- NSGA-II。

首先应用松弛算法求解不确定需求下的RC模型，步骤如下：

步骤1：初始化参数，初始迭代次数 ，当 ，进行基于用户均衡的交通配流得到初始配流参数 和 ；

步骤2：基于 的交通配流结果，根据鲁棒优化的对偶转化方法可以得到 和 ；

步骤3：根据步骤2中计算得到的 和 ，引入松弛算法求解RC模型。具体如下：

步骤3.1：设置初始的辅助参数 ，设置迭代次数限制 和更新因子 ，迭代次数从 开始；

步骤3.2：基于当前的辅助参数 ，求解如下模型：

                                    (85)

                                         (86)

                                         (87)

                                         (88)

                                         (89)

                                         (90)

                                         (91)

模型中的其他约束见式（57）、式（59）、式（63）~式（64）、式（66）、式（68）、式（70）和式（72）。

步骤3.3：如果 ，那么 ，按照 更新辅助参数。返回步骤2，否则转到步骤4。

步骤4：通过上面求解RC模型可以得到 和 。

将松弛算法求得的不确定需求 和通航能力增量 作为输入变量，接下来执行NSGA-II的算法步骤。

步骤5：初始化NSGA-II算法参数。定义参数变量包括染色体编码方式，种群规模 ，适应度值和交叉、变异概率，以及算法中的最大迭代次数 ，给定OD需求 ，随机产生父代种群，令进化迭代开始为 。

步骤6：航道流量分配和适应度值计算，具体如下：

步骤6.1：对给定需求 采用Frank-wolf算法求解下层用户均衡的交通配流模型，得到每个OD对 之间的航道流量 。Frank-wolf算法是专门用来求解固定需求下的用户均衡模型[57]，即当需求为固定时的模型（49）~（52），该算法求解配流的思路就是反复在最短路径上使用“全有全无”流量加载法，由于该算法应用已经很经典[58]，本文不涉及对此算法的改进，因此不再赘述Frank-wolf算法执行步骤。

步骤6.2：将各航道流量 代入上层规划模型目标函数中，求解适应值。

步骤7：采用NSGA-II算法更新种群，具体如下：

步骤7.1：将所有个体按非支配规则排序，令其适应值等于非支配序 。

步骤7.2：通过轮盘赌的方法选择交配池中的个体，进行交叉、变异产生初始父代种群 ，交叉概率为 ，变异概率为 ；

步骤8：基于NSGA-II算法流程执行操作，具体如下：

步骤8.1：快速非支配排序操作。快速非支配是根据个体的非劣解水平对种群分层，其作用是指搜索向帕累托最优解集方向进行。它是一个循环的适应值分层评估过程。首先对非支配解得定义进行说明：假设任意两个解S1和S2对所有目标而言，S1优于S2，则称S1支配S2，若S1的解没有被其他解所支配，则称S1为非支配解，即帕累托最优。在算法中，需要找出群体中非支配解集，记为第一非支配层Level1，将其所有个体赋予非支配序值 （其中 是个体 的非支配排序值），并从整个种群中除去，然后继续在余下群体中寻找非支配解集，记为第二层支配层Level2，个体被赋予非支配序值 ；照此下去，直到整个种群被分层，同一分层个体具有相同的非支配序值 。

步骤8.2：个体拥挤距离计算。通过快速非支配排序后，种群会被分层，每一层中会存在多个个体，然而，这些解得分布情况是无法判断的，若解的分布集中在某一空间，则会导致帕累托前沿面分布不均匀，难以保证种群的多样性。为了解决这一问题，个体拥挤距离 的排序十分重要。个体 的拥挤距离是指：在同一非支配层内，将个体按照任意一个目标函数值进行升序排列，此时个体 与相邻的个体 与 之间的距离被称为拥挤距离。通过选择拥挤距离最大的个体，以保证种群分布的均匀性，显然同层内部分在边缘的个体选择优势最大，因此给定一个足够大的数确保边缘的个体能够被选中 ，再利用下式（92）对中间个体求拥挤距离。根据上述描述即可求得所有非支配层中的个体拥挤距离。

                                    (92)

步骤8.3：精英保留策略。精英保留策略即保留父代中的优良个体直接进入子代，以防止获得的帕累托最优解丢失。精英保留策略选择算子由父代和子代合成的种群 进行优选，以组成新的父代种群 。按照非支配序值 从低到高进行排序，将整层种群依次放入 ，直到出现 大小超过种群规模限制 的情况。

步骤9：判断是否达到最大迭代次数。若 ，则转入步骤8；否则转入步骤10。

步骤10：输出最优帕累托前沿面，算法终止。

根据上述步骤1~步骤10，我们给出算法简单流程示意图如下图1所示。

图1 R-NSGA-II算法流程图

6 案例分析

6.1 案例背景及数据整理

本文以中国进口原油航运网络为例进行案例分析。本文首先根据中国进口原油来源国分布情况，获取实际的进口原油航运网络，案例数据来源中国学者王爽、吕靖等[43,59]的研究。利用本文构建的航运网络设计双层规划模型确定连通可靠性最高的航线路径和扩建航道方案。

根据联合国国际贸易中心的统计数据显示，2020年中国进口原油3.35亿吨，中东、非洲、俄罗斯以及拉丁美洲是中国的四大原油主要进口地区，分别占有51.1%、18.8%、12.6%和12.5%，共计约占全部原油进口95.0%的市场份额[60]。除了从俄罗斯、哈萨克斯坦和蒙古进口的原油经由陆地运输，其余全部依靠海上运输，结合中国近年来主要原油进口来源国的分布情况，给出中国原油进口海上运输网络拓扑示意图，如下图2所示。

图2 中国进口原油航运网络示意图

从图2可以看出，中国进口原油航运网络中的主要进口来源国包括中东地区的沙特阿拉伯、阿曼、伊拉克、伊朗、阿联酋和科威特，非洲地区的安哥拉、刚果、苏丹和利比亚，以及拉丁美洲的委内瑞拉、哥伦比亚和巴西。原油进口主要海峡运河包括：霍尔木兹海峡、马六甲海峡、巽他海峡、龙目海峡、直布罗陀海峡、台湾海峡和巴拿马运河。本文以上述13个原油进口国和中国宁波港为代表，根据图2的航运网络示意图构建21个节点的原油进口航运网络拓扑结构图，如图3所示，21个节点包括了原油进口来源国港口、宁波港以及海峡运河构成，共有13个OD对，图3中的各节点名称如表2所示。（说明：由于资料查找范围受限，节点1~13用进口来源国的名称表示这些进口来源国的实际港口节点）

本案例分析中的算法参数设置为种群规模P=200 ，最大迭代次数imax ，交叉概率pic=0.7 ，变异概率pim=0.05 ，连通可靠性参数置信水平取值υ=0.2 ，BPR函数中的参数取值，航道建设总投资成本预算为B=100 （单位为亿美元），案例研究假设航道每扩宽2km，投资为1000万美元。

图3 中国进口原油航运网络拓扑示意图

表2 图3中的各节点名称

节点编号	节点名称	节点编号	节点名称
1	沙特阿拉伯	12	哥伦比亚
2	伊拉克	13	巴西
3	伊朗	14	霍尔木兹海峡
4	阿联酋	15	马六甲海峡
5	科威特	16	巽他海峡
6	阿曼	17	龙目海峡
7	安哥拉	18	直布罗陀海峡
8	刚果	19	台湾海峡
9	苏丹	20	巴拿马运河
10	利比亚	21	宁波港
11	委内瑞拉

需要说明的是，图3不同于道路交通网络设计，航运网络设计中主要考虑的是对航道的通航能力提升，在道路交通网络设计中任意两个节点之间连接的边都会有容量限制，但是在图3中两个节点连接的边为航运线路，本文假设在海上运输线路上没有容量限制，真正带来通航能力限制的是航运网络中的航道（海峡、运河），因此本文重点研究的是航道通航能力提升和连通可靠性优化问题，即表2中的节点14~节点20。根据图3还可以看出，所有OD对中的非海峡、运河节点，达到其他任意节点都必须要通过上述7个航道节点。根据已有资料显示，7个航道的现状通航能力特征如下表3所示，本研究航运网络中的航道通航能力与一般意义上的通航能力有所不同，一般意义上的航道通航能力是指可容纳邮轮的最大运载能力，而本文研究航运网络设计问题，关注航道中的邮轮量，因此本研究中的航道通航能力是指航道每天可以最大容纳邮轮的数量。

表3 七个航道通航能力现状特征

序号	名称	航道长度（千米）	航道平均宽度（千米）	自由通航时间（小时）	现状通行能力（艘/天）	航道容量拓展上线（艘/天）
1	霍尔木兹海峡	150	70	5.5	144	250
2	马六甲海峡	1080	140	40	275	380
3	巽他海峡	120	60	4.4	80	200
4	龙目海峡	80	50	3.0	120	300
5	直布罗陀海峡	90	43	3.3	220	350
6	台湾海峡	500	300	18.5	180	400
7	巴拿马运河	65	55	2.4	120	180

注：表3中假设各航道非拥挤状态下的自由通航时间为15节（1节=1.85公里/小时）

各海峡和运河的不确定分布参数借鉴王爽、吕靖等[59]的研究中关于航运网络关键海峡和运河的安全评价结果以及世界银行发布的数据确定，整理如下表4所示。（ 的单位为百万吨）

表4 不确定参数分布设置

节点
沙特阿拉伯	0.70	0.75	0.85	50
伊拉克	0.65	0.73	0.78	40
伊朗	0.68	0.78	0.80	20
阿联酋	0.80	0.88	0.95	20
科威特	0.70	0.73	0.82	20
阿曼	0.72	0.75	0.78	10
安哥拉	0.65	0.71	0.82	15
刚果	0.65	0.68	0.82	15
苏丹	0.60	0.70	0.75	10
利比亚	0.60	0.68	0.75	15
委内瑞拉	0.62	0.75	0.82	15
哥伦比亚	0.80	0.82	0.90	20
巴西	0.68	0.75	0.80	20
霍尔木兹海峡	0.60	0.78	0.88	0
马六甲海峡	0.60	0.70	0.75	0
巽他海峡	0.65	0.80	0.85	0
龙目海峡	0.70	0.75	0.85	0
直布罗陀海峡	0.70	0.80	0.83	0
台湾海峡	0.78	0.80	0.85	0
巴拿马运河	0.75	0.82	0.85	0
宁波港	0.95	0.95	0.95	0

6.2 案例结果分析

图3所示的航运网络拓扑结构图中，以宁波港为终点的13个OD对之间一共有30条运输路径，由于论文篇幅限制，此处不再全部罗列连通路径。基于R-NSGA-II算法求解关于30条路径的运输时间和连通可靠性的帕累托前沿面，得到图4和图5所示的帕累托解分布（运输时间单位：天）。

 根据图4和图5，从宏观层面可以看出不同于单目标规划只能得到一个最优解，在权衡系统总走行时间和连通可靠性两个目标后，得到了不同的优化方案，找到最佳决策组合方案。图4和图5还可以看出，系统总走行时间和网络连通可靠性之间存在一定相关关系。从图5可以观察到，随着系统总走行时间增加，网络连通可靠性呈现上下波动的趋势，说明在航运网络中若总走行时间改变会影响网络连通可靠性；根据图4可以观察到，随着连通可靠性的增加，网络总走行时间也是呈现上下波动趋势。由此证明了：航运网络中连通可靠性越高，系统总走行时间不一定最少，系统总走行时间最少，不一定连通可靠性最高。这与实际航运网络中的运营环境相符合，比如有的航道距离短，运输时间较少，但是运营突发风险较大，有的航道连通可靠性较好，但是运输时间较长。则实际决策中根据连通可靠性和运输时间偏好来选择运输路径，这里分别给出基于连通可靠性结果对应走行时间和基于鲁棒优化的最优走行时间总和对应的连通可靠性方案。

基于上文路径连通可靠性计算方法来计算中国进口原油航运网络的最高连通可靠性 。首先识别网络中的所有连通路径关系，并将这些路径连接起来组成联合路径集。然后按照式（22）~式（26）及表4中的不确定参数取值计算每个节点的连通可靠性，最后计算得到航运网络的连通可靠性最高为0.875。将连通可靠性代入上层规划模型计算，通过R-NSGA-II可求得对应的一个系统总走行时间为21天。

计算不同投资预算规模下的系统总走行时间，投资预算规模上限为 （单位为亿美元），求解连续航运网络设计双层规划RC模型，可以得到最优的投资组合方案，如下图6所示，图6对比了求解双层规划RC模型和只求解上层规划模型的求解结果（即是否考虑网络配流问题）。从下图6可以看出，仅求解上层规划模型对应的投资预算不仅较大，且系统总走行时间也较长，系统总走行时间随着投资预算增加呈现递减趋势，但是求解双层规划RC模型可以找到最佳的投资预算方案和最优系统总走行时间，即当投资预算为41亿美元时，求得的系统总走行时间最小为15天，在投资预算规模为50亿美元范围内时，仅求解上层规划模型无法得到最佳的投资预算规模和系统总走行时间方案。则当系统总走行时间为15天时，求解上层规划模型，通过R-NSGA-II可求得对应的网络连通可靠性是0.8。

根据以上两种组合方案，可以得到连续航运网络航道通航能力扩容结果，如下表5所示。

根据表5可以看出，两种最优方案对应的航道通航能力扩容结果不同，方案2对于航道通航能力增加量更接近于容量上限，而方案1对应航道通航能力增加量相对较小，这是由于方案2对于系统总走行时间优化较大，因此更倾向于扩大航道容量以满足更多航船通行，提升航道运输效率。

本文设计的R-NSGA-II是用来求解双目标连续航运网络设计双层规划鲁棒优化模型，由于目前已有的单一现代启发式算法和求解器无法完成对本文提出的目标函数进行求解，为了体现本文设计R-NSGA-II的优越性，只能采用部分模型求解结果对比，即以上层规划模型中的系统总走行时间为目标函数，将R-NSGA-II与模拟退火算法、标准遗传算法和CPLEX求解器进行对比，对于算法参数设置规定：遗传算法参数设置与本案例设置相同，模拟退火算法初始温度设置为500，终止温度设置为0.5，链长设置为10，降温系数为0.9，规定若温度低于0.5，即可停止迭代，四种求解方法的结果对比如下表6和图7所示。

根据表6和图7中四种求解方法运行过程和求解结果可以看出，标准遗传算法虽然迭代次数少，算法运行时间短，但是求解质量是最差的，这是由于遗传算法具有本身容易过早收敛从而陷入局部最优的特征。模拟退火算法在启发式算法中，迭代次数最大，算法运行时间最长，算法稳定性也较差，因此不适合此类问题求解。CPLEX求解器具有可以精确求解的能力，虽然求解结果较好，但是算法运行时间太长，计算耗时约为R-NSGA-II的2倍。综合四种算法可以看出，R-NSGA-II具备遗传算法求解速度快的优势，且求解质量最优，因此可以证明R-NSGA-II在求解此类问题中具有适用性和优越性。

本文提出了一种考虑OD需求不确定和连通可靠性不确定下的连续航运网络设计方法，通过构建双层规划模型、设计启发式算法、案例分析验证了本文模型和算法的合理性，具体而言，通过本文研究得到以下结论：

（1）在航运网络设计问题中，考虑网络配流因素构建的双层规划模型比传统的以系统总阻抗值最小的网络设计模型更具有实际指导意义，主要体现在对目标函数优化而言，考虑网络配流因素构建的双层规划模型能够得到投资规模预算下的最优目标函数值，而不考虑网络配流因素则不能得到最优目标函数值，且从系统总走行时间来看，考虑网络配流因素构建的双层规划模型能够使目标函数优化效率最大提升25%；

（2）在航运网络设计问题中，若以网络连通可靠性最大和系统总走行时间最小为双目标，可以通过本文算法得到帕累托最优前沿面，两者具有相关关系，但是不易量化，因此当分别考虑两者之一最优得到另一个值时，航运网络设计方案具有差异，差异结果表明了当倾向于系统总走行时间最优时，航道通航能力增加量更贴近通航能力上限，通航能力最高可以达到上限的96%；当倾向于连通可靠性最大时，航道通航能力增加量相对较少，通航能力最高达到上限的88%；

（3）本文提出的基于鲁棒优化的航运网络设计双层规划模型可以较好地解决不确定条件下的网络设计问题，R-NSGA-II对求解此类问题具有较好的适用性和优越性，R-NSGA-II相较于标准遗传算法、模拟退火算法和CPLEX求解器算法性能分别提升60%、30%和50%；

（4）通过本文研究可以为中国原油航运网络鲁棒性提升提供政策性指导建议：霍尔木兹海峡、马六甲海峡、巽他海峡、龙目海峡、直布罗陀海峡、台湾海峡和巴拿马运河七条航道为中国原油进口主通道，由七条航道构成的航运网络连通可靠性最大为0.9，在航运网络设计中建议重点提升这几条航道通航能力，以适应未来非常规突发事件下由于不确定环境因素带来的挑战。

还需要进一步研究的是关于航运网络设计的下层规划模型是基于用户均衡的网络配流模型，没有充分考虑决策者对航线决策行为偏好，实际上，未来应对全球供应链不确定风险中，不确定因素将会很大程度上影响人们的决策行为。此外，探索更多求解航运网络设计双层规划模型的算法，从而更加全面论证本文研究提出的算法适用性。