Algorithm 3

A simplified version of the algorithm can be referred to Data-Driven Path-Based Ship Scheduling Optimization Algorithm.pdf

输入

1. 输入-港口信息表.csv

   港口代码	港口卸货效率	港口装货效率	港口中文	港口英文	经度	纬度

2. 输入-港口限制表.csv

   港口代码	船舶类型	进港限制	出港限制

3. 输入-港间航行表.csv

   开始港口	结束港口	航行时间（小时）

4. 输入-货盘.csv

   货盘ID	客户ID	客户	排期计划编号	装货港口代码	装港最大载重	装货效率	卸货港口代码	卸货港最大载重	卸货效率	卸货港经度	卸货港维度	纸浆装货量（吨）	纸浆体积	装货航区	卸货航区	受载开始	受载截止	等级客户	允许拼货	允许拆货

5. 输入-运力.csv

   mmsi	船舶名称	运力释放港口	运力释放时间	载重吨	包装舱容	船舶长度	船舶宽度	船舶高度	船舶类型

路径生成

该部分的主要功能是在给定船舶及其他相关条件，即港间平均航行时间（或船舶航速和港间距）、港口装卸货效率等数据，为船舶生成所有可行的路径。

该部分操作如下，

首先基于货盘信息，为每一个货盘的装货港口和卸货港口建立节点，即若有$m$个货盘，那么建立$2m$个节点，每个装货节点包括货盘序号、货盘ID、客户ID、港口英文、港口代码、港口进出港吃水限制、港口长宽高限制、纸浆装货量（吨）、纸浆体积、客户优先级、受载开始、受载截止、卸货航区、装（卸）货效率、延误时间、相邻节点等信息，卸货节点根据实际需求录入装货节点中的部分信息，然后建立节点与节点之间的连接（即弧）形成运输网络，

两个节点之间的连接有四种可能：

1. 装货港口-装货港口：根据受载期尝试互连，即两个装货港货物的受载期若有重合，则考虑互联，否则让受载期较早的装货节点连接受载期较晚的装货节点
2. 装货港口-卸货港口：全连接
3. 卸货港口-装货港口：鉴于实际业务对于避免绕航并且同一艘船在同一港口不应多次停靠等的要求，船舶不被允许在卸过货后重新回到装货港口装货
4. 卸货港口-卸货港口：全连接

注：只考虑同一航区之间节点的连接情况，不同航区之间不连通

在完成装卸货港口节点图的建立后，读入船舶数据，为每条船的运力释放港口构建节点，并让该节点与所有装货港口节点连接。若有$n$条船，则目前总共$2m+n$个节点。

路径生成的过程可以视作让每条船在图中进行广度优先搜索。我们使用堆栈来实现该过程，每次从堆栈中pop出储存该船当前所在港口current_port、挂靠过的港口visited_ports、当前时间current_time、当前运力使用率（载重使用率和仓容使用率）current_cap_util、当前纸浆载重current_cargo_weight、当前纸浆体积current_cargo_volume、运输过的纸浆编号以及延误时间assigned_cargos、当前所载纸浆编号cargo_indices_on_ship、空放时长empty_sailing_duration、挂靠过的港口名visited_ports_names等等信息的元组，

遍历所有和当前港口相连的港口，检查

• 是否挂靠过该港口（业务要求同一港口不应多次停靠，但鉴于不同节点可能对应着同一港口，因此每次检查是否挂靠过该港口时需要豁免相邻挂靠港口是同一港口的情况）
• 船舶长宽高应满足港口限制
• 如果该节点是卸货节点，船上需装有该节点需要卸载的纸浆
• 根据港口货盘信息，计算抵港装卸货后船舶的预计载重以及纸浆体积等数据，和港口吃水限制进行比对
• 计算预计抵港时间，该时间点不应晚于受载期结束时间加上最大允许延误时间

若以上任意条件不满足则continue（在诸多约束条件的限制下，预计运输网络的规模可控），否则允许抵港，并将更新后的相关信息存入堆栈内，然后不断重复该过程，直到堆栈为空

路径选择

1. 决策变量

   设$a_0=0$，$a_1$为船$1$的候选路径数，$a_2$为船$2$的候选路径数，……，$a_n$为船$n$的候选路径数，

   决策变量$x$可以分为三部分，即$x=\begin{bmatrix} x_1 \ x_2 \ x_3 \end{bmatrix}$

   $x_1$为一个维度为$\sum_{i=1}^n a_i \times 1$的列向量，向量第$1+\sum_{i=0}^{j-1} a_i$到$\sum_{i=0}^{j} a_i$的部分对应船$j$各条候选路径，

   对于$k \in [1+\sum_{i=0}^{j-1} a_i, \sum_{i=0}^{j} a_i]\cap \mathbb{Z}$，如果船$j$选择它的候选路径中的第$k-\sum_{i=0}^{j-1}a_i$条，$x_{k}=1$，$x_{k}=1$$x_k = 0$

   $x_2$为一个维度为$m \times 1$的列向量，该向量为松弛变量，用于实现对货盘违约的惩罚， $x_{2i}\in{0,1}\quad \forall x_{2i} \in x_2$

   $x_3$是一个标量，该向量同为松弛变量，用于实现对未能满足集装箱需求的惩罚， $x_{3}\in{0,1}$

2. 目标函数

   对于每条路径，我们都有包括航次天数、载重使用率、仓容使用率、空放时长、挂靠港口数量、货盘分配情况以及延误时间在内的一系列数据。我们将每个数据点从所有路径中提取出来，并将它们组合成一个向量。

   例如，我们将所有路径的航次天数提取出来，形成一个向量；我们将所有路径的载重使用率提取出来，形成另一个向量，以此类推。

   这样，我们就可以得到航次天ship_route_sailing_days、载重使用率ship_load_utilization、仓容使用率ship_capacity_utilization、空放时长empty_sailing_days、挂靠港数num_ports_of_call、是否分配货盘ship_sailing_penalty、延误时间delay_array等向量，每个向量都代表了一种特定的数据类型。各向量维度为$1\times\sum_{i=1}^n a_i$，每个向量第$1+\sum_{i=0}^{j-1} a_i$到$\sum_{i=0}^{j} a_i$的部分对应船$j$各条候选路径。将这7个向量线性按一定系数线性组合后得到$c_{0}^\top$，再在$c_{0}^\top$尾部添加违约的惩罚系数（维度为$m\times1$的向量，具体大小和每个货盘的客户优先级有关）和未能满足集装箱需求的惩罚系数，得到$c^\top$

3. 限制条件

   1. 为实现每条船只能选择一条路径

      构造一$n\times\sum_{i=1}^{n}a_{i}$矩阵route_choice，$a_{i}$为船$i$的候选路径数（候选路径中包含空路径），该矩阵第$j$行第$1+\sum_{i=0}^{j-1} a_i$到$\sum_{i=0}^{j} a_i$的位置为1，其他为0，令该矩阵乘决策变量列向量等于$\mathbb{1}$向量，即route_choice @ x_1 == np.ones(route_choice.shape[0])

   2. 为实现集装箱业务需求

      根据集装箱业务对于不同时间段ship_preparation_period对于抵达对应港口船只数量的要求，我们用1和0来表示所有候选路径的航行结束时间是否在要求的时间窗口内，形成一$t\times\sum_{i=1}^{n}a_{i}$矩阵container_liner_requirements，其中$t$为时间窗口数，$a_i$为船$i$的候选路径数，矩阵的每行分别对应ship_preparation_period中一个时间窗口的数据。

      为了表示各条路径所属的船型，我们构造了一个$\sum_{i=1}^{n}a_{i}\times1$向量ship_size，其中$a_i$为船$i$的候选路径数，该向量被分为$n$个部分，每一部分对应一艘船的所有候选路径。具体来说，对于船$j$，其对应的部分是向量的第$1+\sum_{i=0}^{j-1} a_i$到$\sum_{i=0}^{j} a_i$个元素。我们用$10^{2k}$来表示各种类型的船，其中$k\in\mathbb{N}$。例如，如果第一艘船是小船，那么向量的第一部分全是1；如果第二艘船是大船，那么向量的第二部分全是100。

      我们选择使用$10^{2k}$来表示各种类型的船，是为了确保在计算在规定时间段内抵达对应港口船只数量时，不会出现两种不同的组合方式得到相同的结果。具体来说，如果我们总共需要$x$条小船，$y$条大船，不妨将小船、大船的对应数据设置为$1$和$10^2$，那么我们选择完对应路径后，对应路径中符合时间窗口部分的船型数据和应该是$x+10^2 y$。由于每种类型的船的数量一般不会超过100条，因此这个值是唯一的，不会有第二种组合方式得到相同的结果。这种方法使我们能够有效地表示各个路径所属的船型，并避免在求解过程中出现歧义。

      最后，我们将container_liner_requirements矩阵的每一行与ship_size向量进行按元素乘法，并将结果重新赋给container_liner_requirements矩阵。此时的container_liner_requirements矩阵维度仍为$t\times\sum_{i=1}^{n}a_{i}$，其中0代表该候选路径的航行结束时间不在要求的时间窗口内，非0值则代表某个特定型号的船的该条路径在要求的时间窗口内。这样，我们就可以综合考虑各候选路径的航行结束时间和船型。

      例如，若集装箱业务需求为在月中月末一共准备至少五条大船五条，我们可以令container_liner_requirements.sum(axis=0) @ x_1 >= 501

      注：若之后新增多种船型需求而非简单区分大小船，可以用类似方法替代，比如用$10^0$、$10^2$、$10^4$等依次代表62K船型、64K船型、68K船型等，另外若对于不同时间窗口满足情况有不同偏好，也可在container_liner_requirements.sum(axis=0)之前先对container_liner_requirements不同行乘不同系数之后再计算列和

   3. 为实现每个货盘最多只能被一条船运送一次

      根据候选路径，提取出每条路径服务过的货盘序号，然后根据所有候选路径的服务过的货盘信息构造一$m\times\sum_{i=1}^{n}a_{i}$矩阵delivery_matrix，$m$是货盘数，$a_{i}$为船$i$的候选路径数。该矩阵第$j$行第$1+\sum_{i=0}^{n-1} a_i$到$\sum_{i=0}^{n} a_i$列用1和0来分别在选择船$k$的各条候选路径时，货盘$j$是否会被船$k$服务，鉴于货盘并非一定会被运送，delivery_matrix @ x_1 <= np.ones(delivery_matrix.shape[0])

   注：实际操作时我引入了松弛变量来将不等式化为等式，并把多个限制条件矩阵合并为一个大的矩阵，通过np.vstack、np.hstack、np.eye、np.zeros等函数实现

基于实际运行情况进行的调整

1. 是否卸货

   若纸浆可以拆货为两部分，可以考虑在导入货盘信息时将所有（或部分）可拆的货盘拆分为两个节点，但存在的问题是暂时没有较好的方案得知将纸浆按何等比例拆分能获得较好效果，且在拆货较多时，图中节点数会增加不少，程序运行时间或显著增加但最后履约率可能只有少量提升

2. 延误时间的处理策略
   1. 第一种策略是直接设定一个较大的最大延误天数。这种策略的局限性在于，由于允许的延误天数较大，可能会导致可行的船舶航线排列组合数量增加，从而使得程序运行时间延长。这种策略的优势在于，由于所有货盘都允许较大的延误天数，会导致有更多的订单能够在延期后的时间内履约，从而提高履约率。
   2. 第二种策略是逐步增加允许的最大延误天数。具体步骤如下：
      1. 在初始阶段，我们设定一个较小的最大延误天数，这样可以减少各港口排列组合形成的可行航线数量。
      2. 在每次运行后，我们将未履约的订单的允许最大延误天数增加1，并锁定那些结果较为理想的船货匹配。这样在后续最大延误天数变大时，由于锁定了一部分船和货，虽然允许的最大延误天数增加了，但各港口排列组合形成的可行航线数量会显著降低，从而可以在一定程度上控制程序运行时间。
      3. 我们将重复上述操作，直到所有货物都被锁定或者允许的最大延误天数达到5为止。这种策略旨在通过动态调整延误时间和优化船货匹配来平衡效率和运行时间。这种策略的不足之处在于，由于逐步增加延误时间，会导致考虑并未考虑到所有的可行解，因此最优解的履约率或许不如第一种方案。

3. 卸货港口之间的抵港顺序

   在实际操作中，或许可以考虑限制船舶必须从南到北卸货，或者先卸重货后卸轻货，给定卸货顺序的限制条件会一定程度上减少代码运行时间

排期格式标准化

Internship Certificate