船舶スケジューリングモデル

実務的な船舶スケジューリングモデルの定式化と解法

海運在庫配送問題

maritime inventory routing problemのベンチマーク問題例 MIRPLIB: https://mirplib.scl.gatech.edu/home

添字と集合

• \(t \in \mathcal{T}\): 期の集合（\(T = |\mathcal{T}|\)）
• \(v \in \mathcal{V}\) (vc \(\in \mathcal{VC}\)): 船舶の集合（船舶クラスの集合）
• \(j \in \mathcal{J}^P\) (\(r \in \mathcal{R}^P\)): 生産、すなわち積荷の港（地域）の集合
• \(j \in \mathcal{J}^C\) (\(r \in \mathcal{R}^C\)): 消費、すなわち荷降ろしの港（地域）の集合
• \(j \in \mathcal{J}\) (\(r \in \mathcal{R}\)): 全ての港（地域）の集合：\(\mathcal{J} = \mathcal{J}^P \cup \mathcal{J}^C\) および \(\mathcal{R} = \mathcal{R}^P \cup \mathcal{R}^C\)
• \(n \in \mathcal{N}\): 通常の点；港と期の組の集合：\(\mathcal{N} = \{n = (j, t) ; j \in \mathcal{J}, t \in \mathcal{T}\}\)
• \(n \in \mathcal{N}_{0, T+1}\): すべての点の集合（開始点 \(n_0\) と終了点 \(n_{T+1}\) を含む）
• \(a \in \mathcal{A}\): 全ての枝の集合
• \(a \in \mathcal{A}^v\): 船舶 \(v \in \mathcal{V}\) に関連する枝の集合
• \(a \in \mathcal{F}_{S_n}^v\): 点 \(n = (j, t) \in \mathcal{N}_{0, T+1}\) および船舶 \(v \in \mathcal{V}\) に対する後続点の集合 (forward star)
• \(a \in \mathcal{R}_{S_n}^v\): 点 \(n = (j, t) \in \mathcal{N}_{0, T+1}\) および船舶 \(v \in \mathcal{V}\) に対する先行点の集合 (reverse star)

データ

• \(B_j\): 港\(j \in \mathcal{J}\)のバース数（バースの上限）
• \(C_a^v\): 船舶\(v \in \mathcal{V}\)が枝\(a = ((j_1,t_1),(j_2,t_2)) \in A^v\)を通過する費用
• \(D_{j,t}^{\min}\)(\(D_{j,t}^{\max}\)): 期\(t \in \mathcal{T}\)に港\(j \in \mathcal{J}\)で生産/消費できる最小（最大）単位数
• \(\Delta_j\): 港\(j \in \mathcal{J}^P\)の場合は \(+1\)、港\(j \in \mathcal{J}^C\)の場合は \(-1\) を取るパラメータ
• \(F_{j,t}^{\min}\)(\(F_{j,t}^{\max}\)): 期\(t \in \mathcal{T}\)に単一の船舶から港\(j \in \mathcal{J}\)に積み込まれる/積み降ろされる製品の最小（最大）量
• \(Q^v\): 船舶\(v \in \mathcal{V}\)の容量
• \(R_{j,t}\): 期\(t \in \mathcal{T}\)に港\(j \in \mathcal{J}\)で荷降ろしされた製品の単位販売収益
• \(S_{j,t}^{\min}\)(\(S_{j,t}^{\max}\)): 期\(t \in \mathcal{T}\)に港\(j \in \mathcal{J}\)での在庫の下限（容量）
• \(s_{j,0}\): 港\(j \in \mathcal{J}\)の初期在庫
• \(s_0^v\): 船舶\(v \in \mathcal{V}\)の初期在庫

変数

• \(d_{j,t}\): 期間\(t \in \mathcal{T}\)に港\(j \in \mathcal{J}\)で生産/消費する量（連続変数）
• \(f_{j,t}^v\): 期間\(t \in \mathcal{T}\)に船舶\(v \in \mathcal{V}\)から港\(j \in \mathcal{J}\)に積み込まれる/積み降ろされる量（連続変数）
• \(s_{j,t}\): 期間\(t \in \mathcal{T}\)の終わりに港\(j \in \mathcal{J}\)で利用可能な在庫数（連続変数）
• \(s_{t}^v\): 期間\(t \in \mathcal{T}\)の終わりに船舶\(v \in \mathcal{V}\)で利用可能な在庫数（連続変数）
• \(x_{a}^v\): 船舶\(v \in \mathcal{V}\)が枝\(a\)を使用する場合に \(1\) となる\(0\)-\(1\)変数
• \(z_{j,t}^v\): 点\(n = (j,t) \in \mathcal{N}\)で船舶\(v \in \mathcal{V}\)が製品を積み込むまたは積み降ろすことができる場合に \(1\) となる\(0\)-\(1\)変数

定式化

max

\[ \sum_{n=(j,t) \in \mathcal{N}} \sum_{v \in \mathcal{V}} R_{j,t} f_{j,t}^v - \sum_{a \in \mathcal{A}^v} C_a^v x_a^v \]

s.t.

\[ \sum_{a \in \mathcal{F}_{S_n}^v} x_a^v - \sum_{a \in \mathcal{R}_{S_n}^v} x_a^v = \begin{cases} +1 & \text{if } n = n_0 \\ -1 & \text{if } n = n_{T+1} \\ 0 & \text{if } n \in \mathcal{N} \end{cases} \quad \forall n \in \mathcal{N}_0, T+1, \forall v \in \mathcal{V} \]

\[ s_{j,t} = s_{j,t-1} + \Delta_j d_{j,t} - \sum_{v \in \mathcal{V}} \Delta_j f_{j,t}^v \quad \forall n = (j,t) \in \mathcal{N} \]

\[ s_{j,t}^v = s_{j,t-1}^v + \sum_{\{n=(j,t) \in \mathcal{N}\}} \Delta_j f_{j,t}^v \quad \forall t \in \mathcal{T}, \forall v \in \mathcal{V} \]

\[ \sum_{v \in \mathcal{V}} z_{j,t}^v \leq B_j \quad \forall n = (j,t) \in \mathcal{N} \]

\[ z_{j,t}^v \leq \sum_{a \in \mathcal{R}_{S_n}^v} x_a^v \quad \forall n = (j,t) \in \mathcal{N}, \forall v \in \mathcal{V} \]

\[ F_{j,t}^{\text{min} v} z_{j,t}^v \leq f_{j,t}^v \leq F_{j,t}^{\text{max} v} z_{j,t}^v \quad \forall n = (j,t) \in \mathcal{N}, \forall v \in \mathcal{V} \]

\[ D_{j,t}^{\text{min} v} \leq d_{j,t} \leq D_{j,t}^{\text{max} v} \quad \forall n = (j,t) \in \mathcal{N}, \forall v \in \mathcal{V} \]

\[ S_{j,t}^{\text{min} v} \leq s_{j,t} \leq S_{j,t}^{\text{max} v} \quad \forall n = (j,t) \in \mathcal{N}, \forall v \in \mathcal{V} \]

\[ 0 \leq s_{j,t}^v \leq Q^v \quad \forall n = (j,t) \in \mathcal{N}, \forall v \in \mathcal{V} \]

\[ x_a^v \in \{0, 1\} \quad \forall v \in \mathcal{V}, \forall a \in \mathcal{A}^v \]

\[ z_{j,t}^v \in \{0, 1\} \quad \forall n = (j,t) \in \mathcal{N}, \forall v \in \mathcal{V} \]

在庫配送問題

在庫配車問題（IRP）は、ある製品を供給者から複数の顧客に一定の期間内に配送しなければならないという配送問題です。 各顧客は最大在庫レベルを設定します。供給者は各顧客の在庫を監視し、在庫切れが発生しないように補充ポリシーを決定します（ベンダー管理在庫ポリシー）。 供給者から顧客への出荷は、一定の容量を持つ車両によって行われます。

集合

• \(\mathcal{V}\): 全ての点（顧客と供給地点を表すデポ \(0\)）の集合
• \(\mathcal{V}'\): 顧客の集合
• \(\mathcal{T}\): 期の集合
• \(\mathcal{K}\): 運搬車の集合

パラメータ

• \(h_i\): 点 \(i\) に関連する在庫費用
• \(c_{ij}\): 点 \(i\) から点 \(j\) への移動費用
• \(r^t\): 期 \(t\) における供給量
• \(d_i^t\): 期 \(t\) における点 \(i\) の需要量
• \(C_i\): 点 \(i\) の最大容量
• \(Q\): 全体の最大供給量

変数

• \(I_i^t\): 期 \(t\) における点 \(i\) の在庫量
• \(q_i^t\): 期 \(t\) における点 \(i\) への配達量
• \(x_{ij}^t\): 期 \(t\) における点 \(i\) から点 \(j\) へ運搬車が移動するとき \(1\) の \(0\)-\(1\) 変数
• \(y_i^t\): 期 \(t\) における点 \(i\) を配送するか否かを表す \(0\)-\(1\) 変数

minimize

\[ \sum_{i \in \mathcal{V}} \sum_{t \in \mathcal{T}} h_i I_i^t + \sum_{i \in \mathcal{V}} \sum_{j \in \mathcal{V}, i < j} \sum_{t \in \mathcal{T}} c_{ij} x_{ij}^t \]

subject to

\[ I_0^t = I_0^{t-1} + r^t - \sum_{i \in \mathcal{V}'} q_i^t \quad t \in \mathcal{T} \]

\[ I_0^t \geq 0 \quad t \in \mathcal{T} \]

\[ I_i^t = I_i^{t-1} + q_i^t - d_i^t \quad i \in \mathcal{V}' \quad t \in \mathcal{T} \]

\[ I_i^t \geq 0 \quad i \in \mathcal{V}' \quad t \in \mathcal{T} \]

\[ q_i^t \geq C_i y_i^t - I_i^{t-1} \quad i \in \mathcal{V} \quad t \in \mathcal{T} \]

\[ q_i^t \leq C_i - I_i^{t-1} \quad i \in \mathcal{V} \quad t \in \mathcal{T} \]

\[ q_i^t \leq C_i y_i^t \quad i \in \mathcal{V} \quad t \in \mathcal{T} \]

\[ \sum_{i \in \mathcal{V}'} q_i^t \leq Q \quad t \in \mathcal{T} \]

\[ \sum_{i \in \mathcal{V}} q_i^t \leq Q y_0^t \quad t \in \mathcal{T} \]

\[ \sum_{j \in \mathcal{V}, j < i} x_{ji}^t + \sum_{j \in \mathcal{V}, j > i} x_{ij}^t = 2y_i^t \quad i \in \mathcal{V} \quad t \in \mathcal{T} \]

\[ \sum_{i \in \mathcal{V}_m} \sum_{j \in \mathcal{V}, j < i} x_{ji}^t + \sum_{i \in \mathcal{V}_m} \sum_{j \in \mathcal{V}, j > i} x_{ij}^t \leq \sum_{i \in \mathcal{V}} y_i^t - y_m^t \quad \mathcal{S} \subseteq \mathcal{V} \quad t \in \mathcal{T} \quad m \in \mathcal{S} \]

\[ q_i^t \geq 0 \quad i \in \mathcal{V} \quad t \in \mathcal{T} \]

\[ x_{ij}^t \in \{0, 1\} \quad i, j \in \mathcal{V}, i \neq j \quad t \in \mathcal{T} \]

\[ x_i^t \in \{0, 1, 2\} \quad i \in \mathcal{V} \quad t \in \mathcal{T} \]

\[ y_i^t \in \{0, 1\} \quad i \in \mathcal{V} \quad t \in \mathcal{T} \]

ベンチマーク問題例

Archettiらのインスタンス

データファイルの形式は以下の通りです：

最初の行には、顧客の数、計画期間の離散時間の数、および車両の輸送容量が記載されています。

2行目には供給者に関する情報が含まれています：

供給者のインデックス

x座標

y座標

在庫の初期レベル

計画期間の各離散時間で利用可能な製品の量

単位在庫コスト

次の行には、各顧客について以下の情報が含まれています：

顧客のインデックス

x座標

y座標

在庫の初期レベル

在庫の最大レベル

在庫の最小レベル

計画期間の各離散時間で吸収される量

単位在庫コスト

---

source

compute_dist

 compute_dist (xi, xj, yi, yj)

---

source

compute_distance_supplier

 compute_distance_supplier (x_coord_supplier, y_coord_supplier, x_coord,
                            y_coord)

---

source

compute_distance_matrix

 compute_distance_matrix (x_coord, y_coord)

---

source

read_input_irp

 read_input_irp (filename)

---

source

read_elem

 read_elem (filename)

folder = "../data/inventory_routing/"
instance_file ="abs1n5.dat"
n, horizon, capacity, start_level_supplier, production_rate_supplier, \
        holding_cost_supplier, start_level, max_level, demand_rate, holding_cost, \
        dist_matrix_data, dist_supplier_data = read_input_irp(folder+instance_file)

#距離行列の準備（n+1が供給地点）
for i, row in enumerate(dist_matrix_data[:]):
    row.append(dist_supplier_data[i])

dist_matrix_data.append(dist_supplier_data+[0])

print("顧客数=", n)
print("計画期間数=", horizon)
print("船の容量=", capacity)
print("供給地点の初期在庫レベル", start_level_supplier)
print("生産速度", production_rate_supplier)
print("供給地点の在庫費用=", holding_cost_supplier)
print("顧客の初期在庫レベル=", start_level)
print("最大在庫レベル=", max_level)
print("消費速度=", demand_rate)
print("在庫費用=", holding_cost)
print("顧客間の距離行列=", dist_matrix_data)
print("供給地点と顧客間の距離=", dist_supplier_data)

顧客数= 5
計画期間数= 3
船の容量= 289
供給地点の初期在庫レベル 510
生産速度 193
供給地点の在庫費用= 0.3
顧客の初期在庫レベル= [130, 70, 58, 48, 11]
最大在庫レベル= [195, 105, 116, 72, 22]
消費速度= [65, 35, 58, 24, 11]
在庫費用= [0.23, 0.32, 0.33, 0.23, 0.18]
顧客間の距離行列= [[0, 265, 102, 214, 226, 85], [265, 0, 366, 368, 238, 349], [102, 366, 0, 207, 302, 17], [214, 368, 207, 0, 431, 203], [226, 238, 302, 431, 0, 289], [85, 349, 17, 203, 289, 0]]
供給地点と顧客間の距離= [85, 349, 17, 203, 289]

数理最適化ソルバーのコード

供給者の時刻tにおける在庫レベルは、時刻 \(t-1\) のレベルに、時刻 \(t-1\) に利用可能な製品量を加え、 時刻 \(t-1\) に顧客へ出荷された総量を差し引いたものである。 供給者における在庫切れ制約は、各配送時刻tにおいて供給者の在庫レベルが時刻tに顧客に配送される総量を出荷するのに十分であることを保証するものである。

各顧客について、時刻tにおける在庫レベルは、時刻 \(t-1\) のレベルに、時刻 \(t-1\) に供給者から顧客に出荷された製品量を加え、 時刻 \(t-1\) に消費された製品量を差し引いたものである。顧客における在庫切れ制約は、各顧客の時刻tにおける在庫レベルが正であることを保証するものである。

容量制約は、時刻tに車両に積載される製品の総量が輸送容量を超えないことを保証するものである。

最大レベル制約は、供給者から各顧客に出荷される製品量が各顧客の最大在庫レベルを超えないことを保証するものである。

目的は、すべてのコスト、すなわち供給者の総在庫コスト、顧客の総在庫コスト、および総輸送コストの合計を最小化することである。

---

source

inventory_routing

 inventory_routing (n, horizon, capacity, start_level_supplier,
                    production_rate_supplier, holding_cost_supplier,
                    start_level, max_level, demand_rate, holding_cost,
                    dist_matrix_data, dist_supplier_data)

x, y, u, delivery, I = inventory_routing(n, horizon, capacity, start_level_supplier, production_rate_supplier, \
        holding_cost_supplier, start_level, max_level, demand_rate, holding_cost, \
        dist_matrix_data, dist_supplier_data)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 162 rows, 177 columns and 561 nonzeros
Model fingerprint: 0xd4baeafd
Variable types: 51 continuous, 126 integer (126 binary)
Coefficient statistics:
  Matrix range     [1e+00, 2e+02]
  Objective range  [2e-01, 4e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [4e+00, 7e+02]
Found heuristic solution: objective 2245.6700000
Presolve removed 22 rows and 28 columns
Presolve time: 0.00s
Presolved: 140 rows, 149 columns, 513 nonzeros
Variable types: 41 continuous, 108 integer (108 binary)

Root relaxation: objective 1.341727e+03, 71 iterations, 0.00 seconds (0.00 work units)

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

     0     0 1341.72667    0   19 2245.67000 1341.72667  40.3%     -    0s
H    0     0                    2218.3700000 1411.84000  36.4%     -    0s
     0     0 1660.00000    0   22 2218.37000 1660.00000  25.2%     -    0s
H    0     0                    1974.3000000 1660.00000  15.9%     -    0s
H    0     0                    1949.3000000 1660.00000  14.8%     -    0s
     0     0 1847.60666    0   18 1949.30000 1847.60666  5.22%     -    0s
     0     0 1849.19000    0   18 1949.30000 1849.19000  5.14%     -    0s
     0     0 1849.54000    0   14 1949.30000 1849.54000  5.12%     -    0s
     0     0 1849.80000    0   24 1949.30000 1849.80000  5.10%     -    0s
H    0     0                    1947.7000000 1849.80000  5.03%     -    0s
H    0     0                    1864.3000000 1849.80000  0.78%     -    0s

Cutting planes:
  Learned: 8
  Gomory: 7
  Implied bound: 6
  Clique: 1
  MIR: 13
  Flow cover: 2
  Flow path: 4
  Zero half: 5
  Mod-K: 4
  RLT: 1

Explored 1 nodes (150 simplex iterations) in 0.05 seconds (0.01 work units)
Thread count was 16 (of 16 available processors)

Solution count 5: 1864.3 1947.7 1949.3 ... 2218.37

Optimal solution found (tolerance 1.00e-04)
Best objective 1.864300000000e+03, best bound 1.864300000000e+03, gap 0.0000%
Opt.value = 1864.3