论文复现20260914
论文题目:Integrated vessel traffic scheduling and berth allocation with restricted channel widths
作者:Duan Hongda et al.
问题归属:航道调度+内锚地分配+泊位分配
决策过程简述
船舶到港后首先到锚地临时停泊,VTS需要决策什么时候把船舶放进航道,亦或是是否拒绝在当前规划期为船舶提供服务,船舶进入港内水域后,港口管理人员则需要根据港内水域泊位的使用情况决定让船舶直接前往某个泊位还是先前往某个内锚地等待。作业结束后VTS需要根据航道的使用情况决定让船舶在前往某个内锚地等待还是直接进入航道。最后船舶离港。
有啥贡献?
- 问题层面:先前研究航道调度的文章对航道宽度和船舶宽度的刻画都是固定的,进而使得船舶在航道的航行规则是固定的,要么只能单向通行(单向航道)、要么只能双向通行(双向航道)、要么某些段可以单向某些段可以双向(复式航道),他这篇文章考虑船舶的在航道的航行模式并不是固定的,而是同时取决于变化的航道宽度和异质船舶宽度(单双向切换航道)。
- 模型层面:建模比较漂亮,频繁使用“把复杂因素隐藏到集合的构造”当中的技巧。
- 算法层面:使用列生成算法的过程中利用机器学习模型去缩小子问题的规模,思想比较有意思。
单双向切换航道有人做过,那他这个还有啥贡献?
刘老师的Short-term berth planning and ship scheduling for a busy seaport with channel restrictions以及Corry19年那篇TS也研究的单双向切换航道,但某种程度上是假设航道宽度固定,从而船舶是否可以双向通行只取决于船舶自身的宽度。如果一艘船舶禁止双向航行,在整个航道他都只能单向航行模式。但是Duan他们这篇文章,船舶在航道的航行模式是同时依赖于航道宽度和船舶宽度的。可能船舶在某个航段只能单向航行,但是在另一个航段又可以双向航行了。
模型推导
- 决策:船舶进入航道的时刻、内锚地选择、某个时刻是否占用内锚地、某个时刻是否占用泊位占用哪个泊位
- 目标:最小化船舶服务拒绝惩罚+靠泊延误惩罚+离港延误惩罚
- 考虑的因素(这里只列举耦合因素,就是把船舶耦合到一起的因素,对应于
后面没有 的那些约束) - 船舶安全裕度
- 船舶宽度和航段宽度对船舶航行模式的影响(这篇文章的特色因素)
- 泊位使用不重叠
- 内锚地使用不重叠
| Notation | Description |
|---|---|
| set of navigation channels. | |
| Set of channels accessible for vessel |
|
| set of time steps in the planning horizon. | |
| Set of time points when vessel |
|
| set of segments of channel |
|
| set of inner anchorages. | |
| set of dummy anchorages, where |
|
| set of berths. | |
| Set of berths eligible for serving vessel |
|
| set of vessels to be served during the planning horizon. | |
| Travel duration between inner end of channel |
|
| Travel duration between inner end of channel |
|
| Travel duration between inner anchorage |
|
| Travel duration to navigate through channel |
|
| Vessel safety clearance in channel |
|
| Width of the |
|
| Arrival time of vessel |
|
| Desired berthing time of vessel |
|
| Desired port departure time of vessel |
|
| Deadline on the port departure time of vessel |
|
| Service duration of vessel |
|
| Width of vessel |
|
| The cost incurred for the non-completion of vessel |
|
| The cost incurred for the berthing delay of vessel |
|
| The cost incurred for the departure delay of vessel |
|
这篇文章比较有技巧的地方是他把一些复杂的因素放到集合的构造中了,这样子模型看上去会特别的简洁(目前还不太清楚能不能带来性能上面的改善,因为这种把复杂因素放入集合构造中的技巧,本身会引入大量 0-1 变量,从集合里面选择合法元素,太多的01变量会不会让求解比较慢呢?)。比如要求船舶在航道的航行过程必须处于某个潮汐窗口内部,常见的做法是引入潮汐窗口选择约束变量,然后用一个不等式约束保证船舶在航道的整个航行过程一定位于潮汐窗口内部。他是直接把潮汐窗对船舶进航道时间的影响隐藏到集合
BIP模型咋推?
目标函数比较好理解,
约束部分,船舶进出港时间、内锚点和泊位使用都是通过01变量刻画的。
Part1:因此,对于一条船舶他的这些选择应该是唯一确定的,所以存在一批选择唯一性约束:
(1)船舶进港(出港)对航道和时刻的选择唯一
- 这篇文章船舶进港时刻指的是船舶处于航道的outer end的那个时刻
- 这篇文章船舶出港时刻指的是船舶处于航道的inner end的那个时刻
(2)船舶进港(出港)对内锚地的选择是唯一的
(3)船舶开始占用哪个泊位、开始占用时刻是唯一的
Part2: 泊位和内锚地是离散资源,可看做并行机,对于这些资源的占用要求使用时段不发生重叠。常见套路是利用资源使用次序变量结合大M约束把时间错开,可参考往期博客。但这篇文章船舶靠泊哪个泊位以及开始靠泊的时刻是通过01选择变量刻画的,所以写法不太一样了,套路是每个资源每个时刻最多被一个资源所占用:
但是泊位这块儿
这种写法有什么弊端?资源的占用时长必须得是已知的,如果说资源占用时长是变量,这种写法是不是就不行了呢?
Part3:然后是集成类问题里面常见的资源使用时序关系约束,我觉得这篇文章的BIP模型写的不对:
- 进港船开始使用内锚地的时刻在进港船到达航道inner end的那个时刻之后,缺失相关约束
- 出港船开始使用内锚地的时刻在船舶完成装卸作业之后,缺失相关约束
- 船舶开始靠泊的时刻在船舶结束使用内锚地的那个时刻+船舶在内锚地和泊位之间航行时间之后,缺失相关约束
文中有的相关约束仅仅是限制“船舶对锚地占用的时长等于离开锚地的时刻减去进入锚地的时刻”,“离开锚地的时刻在进入锚地的时刻之前”,“虚拟锚地的占用时长为0”,这些约束并不能保证1,2,3得到满足吧?
Part4:接着是航行安全的刻画,船舶在航道的航行必须保证最低安全裕度、以及船舶宽度和航段宽度对船舶航行模式的影响。
他这个考虑起来实际上比较简单,因为船舶是同航速的,所以不用考虑追击问题。所以同航向船舶保持最低安全裕度只需要保证进入航道的时间间隔保持最低安全裕度就ok了。
这句话等价于
作者文中定义了
作者应该是在创建变量的时候就没有创建非法的变量,所以写约束的时候也得用复杂的集合运算把合法的变量给取出来。
实际上我觉得可以创建变量的时候把所有可能变量创建出来,然后非法取值变量通过约束强制为0就好了,更好理解,比如:
非法变量特别多的时候用他这种写法可能性能还是会更加好一些
船舶宽度和航段宽度对船舶航行模式的影响这个约束比较难想,要求任意一个航道的任意一个航段上,如果同时存在进港和出港船舶,那么进港船宽度加上出港船的宽度要小于航段的宽度。难点在于时刻t我怎么知道船舶是不是位于航段s,没有变量直接表明这一点,所以需要转换:(1)进港船时刻t位于航段s等价于进港船什么时候进港?
论文里面把
Part 5:最后是船舶服务时间的要求,如果船舶服务完成的时间会超过预计离港时间,那这样的船直接拒绝在当前规划期服务
我自己第一遍推,写的是下面这个约束,他是不对的:
想表达,如果
MP模型咋推?
MP模型写的套路就是定义方案集合,以及配套的描述方案的参数(Indicator parameter),定义方案选择变量。MP模型的约束就是方案选择唯一性约束(如果一个对象只能选唯一一个方案的话)+BIP当中那些把对象耦合到一起的因素,对于这篇文章就是:(1)船舶安全裕度(2)船舶宽度和航段宽度对船舶航行模式的影响(这篇文章的特色因素)(3)泊位使用不重叠(4)内锚地使用不重叠
方案定义
关于
因为Colossus框架目前不支持方案集合分组的操作,只能引入Indicator parameter表征方案的服务对象。因此,下面得到的formulation和原文在形式上略有不同。此外,有了
定义
方案的定义包括——(1)定义方案集合(2)定义描述方案的参数。
关键点在于应该定义哪些描述方案的参数,我目前认为(2)的定义不是唯一的,只要包含原问题模型所需决策的信息,那方案的定义就是有效的。我在设计方案之前心里应该有一个预案,就是我子问题想要按照什么对象分解,比如这篇文章,我想按照船舶分解,那方案定义很显然就应该包含单一船舶作业方案的全部信息。
MP模型
定义方案选择变量
目标函数简单,就是最小化拒绝当前规划期对船舶进行服务的惩罚+方案成本。
其中方案成本
所以目标函数应该是:
约束就是方案选择或选择推迟唯一性约束、以及BIP中的那些把船舶i耦合到一起的约束:
推到这里算是知道为什么作者要对方案按照船舶分组了,且参数下标显含i!因为方案成本是跟船舶i相关的,如果不分组就得先用
把当前方案 对应的是哪个i给计算出来。
投降… 我选择给Colosus框架新增方案集合分组的功能,方案定义和MP模型还是按原文的方式写,这种MP中包含大量参数需要检索方案对应服务对象的情况真的做分组处理,不然表达上太麻烦了!
定价子问题
模型
先推一下变量
在这个问题中,把i耦合起来的约束全部放到主问题了,所以可以按照i对子问题分解,得到对应于船舶i的子问题:
注意,子问题本身就是生成方案的,所以
在子问题中不再充当索引。i由于咱们做了分解,所以在子问题中也固定下来了,不作为索引
关于定价子问题目标函数的一点看法
定价子问题目标函数不必直接是reduced cost,只要最优解和min reduced cost问题对应的最优解一致就可以了。比如它这篇文章子问题就包含决策是否推迟船舶服务的变量
time-expanded上重构成ssp
定价子问题转换成SSP这部分是这篇文章最精彩的内容。有必要搞清楚定价子问题到底是一个什么问题,这样以后我们在遇到类似的问题时也可以使用类似的模型转换思路!
问题界定:
- 归属:单船在港作业方案设计问题。即到港船舶只有一条船的情况下的船舶作业方案设计问题。
- 决策:是否推迟船舶的服务+MP中用到的indicator parameters(船舶进港的时刻、船舶选择哪个内锚地、某个时刻是否占用内锚地、开始靠泊的时刻、船舶出港的时刻)
- 目标:最小化推迟服务惩罚成本+靠泊延误成本+离港延误成本+和indicator parameter关联的对偶成本
- 考虑的因素:船舶使用哪些资源,这些资源的使用时序关系,是否推迟船舶服务
类比一下,以后遇到问题考虑的因素是单一对象使用哪些资源以及资源的使用时序关系,那这样的问题就可以用这篇文章的套路转换到time-expanded network上转换成ssp问题
时间拓展网络重构的关键就是恰当的定义Node和Arc,其中:
- Node代表某个时刻处于某个位置
- Arc是从某个Node指向某个Node的弧,它表达对象从某个时刻离开某个位置于另一个时刻到达新位置或仍处于当前位置,因此对弧的选择就对应了一个调度决策,进而要保证调度的合法性实际上就等价于保证弧的合法性,如果对某个弧的选择对应的调度会违反原问题的约束,那我们就应该在网络构造的时候就禁止相应的弧或节点的生成。比如是推迟服务和不推辞服务只能二选一,那从
出去的弧就只能有两类,一类对应推迟服务一类对应不推辞服务,这样推迟服务和不推辞服务只能二选一的约束自然得到满足。
那这俩玩意儿咋定义出来呢?
定义源点
然后咱们回归问题涉及的变量,思考怎么定义剩下的Node和Arc:
- 子问题需要决策船舶(进港)进入航道外侧的时刻以及所选航道
。我需要弧来表达,如果选中这个弧,就等价于我选择了某个航道以及某个进港时刻。我们把这样的弧集合定义为 。显然“进入航道外侧的时刻以及所选航道”是船舶调度所需做的第一个决策,因此,这个弧的起点就应该是 ,而终点应该是船舶可选航道以及对应的可进入航道的笛卡尔积中的某个元素。因此,我们定义Node集合 表达船舶i可用航道以及对应的可进入航道时刻的笛卡尔积构成的集合。 - 子问题还需要决策是否推迟船舶的服务
,我们把这样的弧集合定义为 ,这也是船舶调度一开始就要决策的,所以 中弧的起点的集合自然是 ,如果推迟服务,那船舶的调度实际上就终止了,所以终点自然是 。 - 子问题需要决策船舶(进港)选择哪个锚地
,这属于是分配决策,而弧的选择对应一个调度决策,但调度决策=分配+时间,所以分配决策的信息也可以包含在调度决策里面。注意到这样的分配决策应该发生在船舶(进港)选哪个航道什么时候进港这个调度之后,所以对应的弧的起点自然应该是 中的某个元素,终点可以定义成船舶可用锚地以及可使用时刻的笛卡尔积中的某个元素,不妨记作 。我们把对应的弧归结到集合 中。 - 子问题还需要决策某个时刻船舶(进港)是否占用锚地
。这样的调度决策应该发生在船舶完成锚地选择之后,因此对应弧的起点应该是 ,由于是决策某个时刻是否使用,所以终点应该属于 中那些同一位置但时间相差1的节点(锚地选择唯一性)。我们将这样弧归属到集合 。 - 子问题还需要决策某个船舶(进港)开始靠泊的时刻,这个调度决策可能发生在船舶选哪个航道什么时刻进港这个调度决策后,表达不适用锚地直接去泊位,也可能发生在某个时刻是否占用锚地这个调度决策之后,表达从锚地去泊位。所以这个调度决策对应的弧的起点可以是
中的元素,也可以是 中的元素。终点可以是船舶可以泊位以及可用时刻笛卡尔积中的某个元素。我们把这样的弧归结到集合 。 - 船舶还需要决策船舶(出港)选择哪个锚地
。这个调度决策对应弧的起点应该是船舶作业完成时刻之后,你会发现目前没有Node集合表达船舶在泊位作业完成时刻,所以新定义一个集合 。终点定义为 。对应的弧归属于 。 - 船舶还需要决策船舶(出港)某个时刻是否使用锚地
,定义 为可用锚地和时刻的笛卡尔积,起点和终点都应该是这个集合中的对应元素,把对应弧归属到集合 - 船舶还需要决策船舶(出港)选择哪个航道和时刻,这个调度决策可能发生在船舶在泊作业完成之后,也可能是在锚地停留之后。所以弧起点可以是
中的元素,也可以是 中的元素,终点定义为船舶可以航道和可到达航道外侧时刻的笛卡尔积 中的某个元素。把对应的弧归属到 。 - 最后整个调度结束,对应弧的起点应该是
中的某个元素,终点是 ,对应的弧归属到 。
至此,时空网络就构造出来了,从
最后一步,我们需要把目标函数中各个调度决策的系数映射成弧的成本。
上图是我们构造的时空网络中弧与调度决策的对应关系。从而,弧成本的定义实际上就是调度决策对应的代价。为了便于观察调度决策对应的代价,我们对子问题目标函数做一下改写,按照调度决策变量合并同类项,写成下面这种形式,我认为他文中那个转换是错误的,但思想没问题,是得做这步转换,按自己推的是这样::
总结一下时空网络构造技巧
- 时空网络中Node包含地点和时刻信息,Arc是Node到Node的有向弧,对应的是一条调度信息,进而Arc的选择就对应于某一个调度或分配决策。
- Node和Arc的定义可以从原问题包含哪些调度决策变量入手,思考对应的调度决策对应弧的起点应该是什么,终点应该是什么。且弧和顶点的构造必须保证对应调度决策是合法的,否则不应该生成相应的弧或节点。
- 将对偶值映射到弧成本上怎么实现?关键看对偶值对应的是哪个调度决策A,进而分析调度决策A又是对应的哪个弧B,从而对偶值就放到弧B上面。
代码
使用Colossus包求解#!/usr/bin/env python # -*- coding: utf-8 -*- # @Date : 2026-09-23 15:00:00 # @Author : syuansheng (Dalian Maritime University) """ 复现论文:Integrated vessel traffic scheduling and berth allocation with restricted channel widths 集成式船舶交通调度与泊位分配问题(VSP-BAP),采用列生成求解紧凑模型 P' 的线性松弛, 最后恢复整数性求解受限整数主问题得到可行整数解。 """ import sys sys.path.append('../../../') from Colossus import ToolBox, Environment, Engine from data import build_toy_instance from initial_cols import generateInitialCols from rmp_model import buildRestrictedMP from exact_sp import build_vessel_network, updateAndSolveSP_Exact # 算例(论文 3.3 节小规模算例,已知最优成本 22) dataloader = build_toy_instance() # 预构建每艘船舶的时间拓展网络(拓扑结构固定,仅弧成本随对偶变量变化) tenetworks = {i: build_vessel_network(i, dataloader) for i in dataloader.V} for i in dataloader.V: net = tenetworks[i] print('>>船舶 {0} 时间拓展网络:节点 {1} 个,弧 {2} 条'.format(i, len(net.node_lst), len(net.arc_lst))) fingerprint_lst = [(i,) for i in dataloader.V] # 把框架需要的工具搭建好 toolbox = ToolBox() toolbox.addTool('generateInitialCols', generateInitialCols) toolbox.addTool('buildRestrictedMP', buildRestrictedMP, dataloader=dataloader) toolbox.addTool('updateAndSolveSP_Exact', updateAndSolveSP_Exact, tenetwork=tenetworks, dataloader=dataloader) # 求解环境搭建好 # 决策变量 lam(服务方案列,需生成)与 u(拒绝服务,需恢复整数)在松弛主问题中被松弛为连续变量 environment = Environment('lam', ['lam', 'u'], ['BINARY', 'BINARY'], fingerprint_lst=fingerprint_lst) environment.heuristic_solve_on = False # 只使用精确定价 environment.col_file_output_on = False # 不输出 .col 文件 environment.result_file_output_on = True # 输出 .result 文件 environment.lp_file_ouput_on = False # 初始 RMP 约束为空,写 LP 文件异常,关闭 environment.ilp_file_output_on = False # 关闭不可行性文件输出 # 构造求解引擎并运行 engine = Engine('VSP_BAP_engine', environment, toolbox) engine.runEngine()
整体上来说,小规模算例上的求解还是非常快的,1s不到:
>>船舶 1 时间拓展网络:节点 198 个,弧 403 条
>>船舶 2 时间拓展网络:节点 166 个,弧 291 条
>>船舶 3 时间拓展网络:节点 166 个,弧 288 条
>>船舶 4 时间拓展网络:节点 180 个,弧 337 条
>>Engine配置信息
[1] Environment(variables_name='lam', relaxed_variables_name=['lam', 'u'], relaxed_variables_type=['BINARY', 'BINARY'], fingerprint_lst=[(1,), (2,), (3,), (4,)], barrier_method_on=True, exact_solve_on=True, heuristic_solve_on=False, mini_batch_percent=1.0, lp_file_ouput_on=False, ilp_file_output_on=False, col_file_output_on=False, result_file_output_on=True, dual_alpha=0, truncate_threshold=0.001, max_not_bertter=inf, max_time_limit=3600.0, rmp_log_to_console=False, current_reduced_cost_on=False, in_rmp_optimal_basis_when_on=False)
[2] ToolBox contains tools: buildRestrictedMP, generateInitialCols, updateAndSolveSP_Exact
>>Note:在调用runEngine前请先检查配置是否满足你的要求
>>预备:初始列构造完毕
Cardinal Optimizer v8.0.6. Build date Aug 7 2026
Copyright Cardinal Operations 2026. All Rights Reserved
>>限制主问题数学模型构造完毕,进入列生成迭代轮次
>>列生成迭代过程:
+-----+---------+---------+--------------+----------+------------+------------+------------+
|iter |time(s) |cols_num |rmp_obj |no_better |rmp_time(s) |sps_method |sps_time(s) |
+-----+---------+---------+--------------+----------+------------+------------+------------+
|1 |0.01 |0 |4,000.00 |0 |0.00 |exact |0.00 |
|2 |0.01 |4 |3,000.00 |0 |0.00 |exact |0.00 |
|3 |0.02 |8 |2,000.00 |0 |0.00 |exact |0.00 |
|4 |0.03 |12 |2,000.00 |1 |0.00 |exact |0.00 |
|5 |0.04 |16 |2,000.00 |2 |0.00 |exact |0.00 |
|6 |0.04 |19 |1,006.00 |0 |0.00 |exact |0.00 |
|7 |0.05 |22 |1,006.00 |1 |0.00 |exact |0.00 |
|8 |0.06 |25 |1,006.00 |2 |0.00 |exact |0.00 |
|9 |0.06 |29 |1,006.00 |3 |0.00 |exact |0.00 |
|10 |0.07 |33 |509.00 |0 |0.00 |exact |0.00 |
|11 |0.08 |37 |14.00 |0 |0.00 |exact |0.00 |
|12 |0.09 |40 |12.00 |0 |0.00 |exact |0.00 |
|13 |0.09 |42 |12.00 |1 |0.00 |exact |0.00 |
+-----+---------+---------+--------------+----------+------------+------------+------------+
>>列生成迭代终止原因:无法再找到改善列
>>列生成迭代轮次终止,下面将根据现有ColumnPool对主问题进行求解:
Setting parameter 'Logging' to 1
Setting parameter 'LpMethod' to -1
Model fingerprint: 5ced3ba1
Using Cardinal Optimizer v8.0.6 on Windows (25H2 Build 26200 - x86_64)
The CPU model is 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
Hardware has 4 physical cores and 8 logical cores. Using instruction set X86_AVX512_E1 (14)
Minimizing a MIP problem
The original problem has:
277 rows, 46 columns and 572 non-zero elements
46 binaries
Starting the MIP solver with 8 threads and 32 tasks
Presolving the problem
The presolved problem has:
19 rows, 29 columns and 137 non-zero elements
29 binaries
Problem info:
Range of matrix coefficients: [1e+00,1e+00]
Range of rhs coefficients: [1e+00,1e+00]
Range of bound coefficients: [1e+00,1e+00]
Range of cost coefficients: [1e+03,1e+03]
Density of cost: 100.0%
Nodes Active LPit/n IntInf BestBound BestSolution Gap Time
0 1 -- 0 -2.487200e+04 -- Inf 0.01s
H 0 1 -- 0 -2.487200e+04 4.000000e+03 116.08% 0.01s
H 0 1 -- 0 -2.487200e+04 2.006000e+03 108.07% 0.01s
H 0 1 -- 0 -2.487200e+04 1.006000e+03 104.04% 0.01s
* 0 1 -- 0 1.200000e+01 1.200000e+01 0.000% 0.01s
0 1 -- 0 1.200000e+01 1.200000e+01 0.000% 0.01s
1 0 0.0 0 1.200000e+01 1.200000e+01 0.000% 0.01s
1 0 0.0 0 1.200000e+01 1.200000e+01 0.000% 0.01s
Best solution : 12.000000000
Best bound : 12.000000000
Best gap : 0.0000%
Solve time : 0.01
Solve node : 1
MIP status : solved
Solution status : integer optimal (relative gap limit 0.0001)
Violations : absolute relative
bounds : 0 0
rows : 0 0
integrality : 0
>>已经将Engine=VSP_BAP_engine的求解结果写入VSP_BAP_engine.result文件中
[Finished in 422ms]









