# Solving Tramp Ship Routing and Scheduling Problems

## Introduction

This work was conducted within the subject TIØ4500 - Managerial Economics and Operations Research, a central part of the Master of Science In Engineering in Industrial Economics and Technology Management at the Norwegian University of Science and Technology (NTNU). The subject TIØ4500 acts as a preparatory course for the final master thesis to be written in the spring of 2023. This work and the upcoming master thesis concern Tramp Ship Routing and Scheduling Problems with Flexible Cargo Quantities and Bunker Optimization and originated through discussion with industry partners Western Bulk and Maritime Optima. As far as I know, it is the first academic study studying a Tramp Ship Routing and Scheduling problem that combines Flexible Cargo Quantities and Bunker Optimization. This article summarizes the main findings. The full report is available here.

## Background

Maritime transportation is the backbone of international trade, dominating other modes of transportation in terms of total volume transported. Due to the economic effects of absolute advantages, comparative advantages, technical improvements, and economies of scale, more than 80% of the total volume traded internationally was carried on board vessels in 2019 (UNCTAD, 2022). For bulk operators in the shipping industry, a deep understanding of trade patterns and skills regarding securing transport contracts and cargo transport can make their operations financially lucrative.

In the studied problem, a dry bulk operator like Western Bulk is to deploy their fleet of hired vessels to maximize the profit for the entire fleet. To do so, they have to pick between a set of mandatory contracted cargoes and optional cargoes, potentially paying a cost for not servicing a contracted cargo. Profit is generated by the cargo transported. The variable costs include port, canal, and fuel purchase costs. As such, the dry bulk operator has to construct profit-maximizing routes and schedules for each vessel in their fleet. In the dry bulk industry, an operator can often choose the exact amount of cargo to transport from a predetermined More or Less in Owner's Option (MoLOO) interval. Thus, the operator needs to decide how much cargo to transport. As vessels may carry multiple cargoes simultaneously, the profit-maximizing amount of each cargo can be challenging to decide. Further, as the ships eventually run out of fuel, the operator must choose where to stop for fuel and the quantity of bunker (fuel) to purchase. Western Bulk relies on charting and operations managers to make such decisions. This study aims to provide these divisions with an additional tool to aid their decision-making process.

This article formulates a mathematical arc flow model that captures the complex operational environment of dry bulk operators such as Western Bulk. Results show that commercial solvers can provide optimal solutions to the arc flow formulation for smaller test instances. For larger test instances, finding an optimal solution becomes too time-consuming. A mathematical path flow model was formulated to overcome the limitations of the arc flow model. To find optimal solutions, the solver of the path flow model leverages an a priori column generation approach. This solution method generates the set of all feasible routes. Each route is then optimized to construct schedules with an optimal profit. The commercial path flow model solver selects the schedules maximizing the fleet-specific profit. The a priori column generation approach obtains optimal solutions for test instances of up to 30 cargoes, 10 vessels, and 10 bunker ports within 30 minutes.

## Assumptions

The design process of the model resulted in some relatively strict assumptions. Next semester's Master Thesis will extend the model to loosen these assumptions. The following lists some of the most notable assumptions:

• Freight rates, bunker prices, and speeds are fixed for the duration of the planning horizon and are assumed to be known.
• The fleet size is fixed for the duration of the planning horizon. The contracted and mandatory may be serviced by vessels chartered in the spot market. Only vessels in the fleet may service optional spot cargoes.
• Some vessels and ports may be incompatible due to berth water depths and loading/discharging gear restrictions.
• Cargoes are uncoupled such that there are no restrictions on the order of how cargo types are transported. In real life, a vessel's cargo holds must be cleaned at a cost after transporting cement, or the vessel is limited to transporting similarly "dirty" cargoes. Dirty cargoes are not modeled in this research.
• Vessels are not allowed to fill bunker and load/discharge cargo concurrently.

## Solution Methods

To model Western Bulk's operational environment, this research suggests a Mixed-Integer Linear Programming (MILP) arc flow formulation. The arc flow formulation models all of the fleet-specific constraints and the vessel-specific constraints together. Commercial solvers such as Gurobi are capable of solving such models to optimality. However, as the instance size grows, the solution time grows exponentially. The largest instance that the solver of the arc flow formulation was capable of solving within one hour consisted of 9 cargoes, two vessels, and four bunker nodes. The full arc flow formulation is presented at the end of this article

A Dantzig-Wolfe decomposition was applied to generate a path flow formulation of the original model so the limitations of the arc flow formulation could be overcome. The path flow formulation is considerably smaller than the arc flow formulation, as there are only two sets of constraints that are fleet-specific in the original arc flow model. The rest are vessel-specific and may be solved independently. Let a route be defined as a consecutive sequence of nodes a vessel will visit, e.g., {0, 23, 5, 15, 1, 11, 30}. Such a route consists of an artificial origin and destination node and potentially bunker option nodes, pickup nodes, and delivery nodes. Let a schedule be a route in which the optimal cargo to transport and the optimal bunker to purchase have been decided to optimize the vessel-specific schedule profit. Given all such optimal schedules and their corresponding profits, solving the path flow formulation yields the optimal assignments of schedules to the fleet of vessels that maximize the fleet-specific profit.

The solution method of the path flow formulation is referred to as a priori column generation and can be summarized as follows:

1. Generate all feasible routes for every vessel in the fleet
2. Optimize every route to generate their corresponding schedule with an optimal vessel-specific schedule profit.
3. Solve the path flow formulation to assign the best combination of schedules to the fleet of vessels.

The a priori column generation approach solved instances of 30 cargoes, 10 vessels, and 10 bunker ports within 30 minutes.

### Route Generation

All feasible routes for each vessel were generated by performing a recursive Depth-First-Search of every vessel graph. In line 19, the algorithm checks for feasibility when considering a potential move to a new node. If vessel-specific constraints such as the ones that model the cargo balance, the bunker balance, or the time balance are broken, the algorithm skips to the next node. If a new node is feasible, the algorithm updates the time, bunker, and cargo resources to reflect correct values at the next node, as shown in line 21. The algorithm recursively calls itself in line 22.

### Schedule Generation

A Linear Programming (LP) Problem must be solved to generate the optimal schedules for each route. It consists of all of the vessel-specific constraints. As the route is given, and there are no fleet-specific constraints, all binary variables in the original arc flow formulation are excluded. Thus, the vessel schedule optimization problem is an LP problem that can be solved in polynomial time. In contrast, the original MILP arc flow formulation is NP-Hard. As such, it is possible to solve a lot of schedule optimization problems in a shorter amount of time than one would spend solving a similar path flow formulation instance. The full formulation of the schedule generation is not presented in this article but is available in the report.

### Path Flow Formulation

Let $$\mathscr{R}_v$$ be a feasible sequence of node visits for vessel $$v$$ comprising a route.  $$\mathcal{R}_v$$ denotes the set of all feasible routes for vessel $$v$$. For each feasible route $$\mathscr{R}_v$$ for vessel $$v$$, the optimal schedule must be found. Each such schedule is associated with an optimal vessel-specific profit, denoted $$P_{sv}$$. Let $$\mathcal{S}_v$$ be the set of optimal schedules for vessel $$v$$ derived from transporting the optimal amount of cargo and purchasing the optimal amount of bunker.  Further, let the sets $$\mathcal{V}$$, $$\mathcal{N}$$, $$\mathcal{N}^C$$, and $$\mathcal{N}^O$$ denote the set of vessels, cargo-related nodes, contracted cargoes, and optional spot cargoes, respectively. As in the arc flow model, $$C^S_i$$ denotes the cost of servicing the cargo at pickup node $$i$$ with a spot ship. A new parameter, $$A_{isv}$$, is introduced and set equal to 1 if vessel $$v$$ carries cargo $$i$$ in schedule $$s$$, and 0 otherwise. The variable $$x_{sv}$$ defines a binary schedule variable equal to 1 if vessel $$v$$ is chosen to sail schedule $$s$$, and 0 otherwise. Finally, the binary variables $$y_i \in \mathcal{N}^O$$ and $$z_i \in \mathcal{N}^C$$ denote whether or not a spot cargo is serviced and whether or not a spot ship is used to service a Contract of Affreightment (CoA) cargo, respectively. Tables 1-3 summarize the notation used in the path flow formulation.

Table 1: Path Flow Sets
Set Notation    Set Description
$$\mathcal{V}$$ Set of vessels
$$\mathcal{N}$$ Set of all cargo-related nodes
$$\mathcal{N}^C$$ Set of pickup nodes for the mandatory contracted cargoes
$$\mathcal{N}^O$$ Set of pickup nodes for the optional spot cargoes
$$\mathcal{R}_v$$ Set of all feasible routes for vessel $$v$$
$$\mathcal{S}_v$$ Set of all optimal schedules for vessel $$v$$

Table 2: Path Flow Parameters
Parameter Notation    Parameter Domain    Parameter Description
$$P_{sv}$$ $$v \in \mathcal{V}, s \in \mathcal{S}_v$$ Vessel specific profit generated by performing
schedule $$s$$ for vessel $$v$$
$$C^{S}_{i}$$ $$i \in \mathcal{N}^P$$ Cost of servicing the cargo at node $$i$$ with
a spot ship
$$A_{isv}$$ $$i \in \mathcal{N}, v \in \mathcal{V}, s \in \mathcal{S}_v$$ 1 if vessel $$v$$ carries cargo $$i$$ in schedule $$s$$,
and 0 otherwise

Table 3: Path Flow Variables
Variable Notation    Variable Domain    Variable Description
$$x_{sv}$$ $$v \in \mathcal{V}, s \in \mathcal{S}_v$$ 1 if vessel $$v$$ is chosen to sail schedule $$s$$, 0 otherwise
$$y_{i}$$ $$i \in \mathcal{N}^O$$ 1 if the optional spot cargo is transported, 0 otherwise
$$z_{i}$$ $$i \in \mathcal{N}^P$$ 1 if cargo at node $$i$$ is serviced by vessel from the spot
market, 0 otherwise

A new set of constraints must be included to ensure that each ship is assigned exactly one schedule. The Objective Function (1) and Constraints (2)-(6) present the path flow reformulation to the arc flow model presented.

\begin{align}
\max & \quad \sum_{v \in \mathcal{V}} \sum_{s \in \mathcal{S}_v} P_{sv} x_{sv} - \sum_{i \in \mathcal{N}^C} C^S_i z_i \tag{1}\\
\text{s.t.} && \\
&\sum_{v \in \mathcal{V}} \sum_{s \in \mathcal{S}_v} A_{isv} x_{sv} + z_i = 1 && i \in \mathcal{N}^C \label{eq-master-spot-ship} \tag{2}\\
&\sum_{v \in \mathcal{V}} \sum_{s \in \mathcal{S}_v} A_{isv} x_{sv} - y_i = 0 && i \in \mathcal{N}^O \label{eq-master-spot-cargo} \tag{3}\\
&\sum_{s \in \mathcal{S}_v} x_{sv} = 1 && v \in \mathcal{V} \label{eq-master-assignment}\tag{4}\\
&x_{sv} \in \{0, 1\} && s \in \mathcal{S}_v, v \in \mathcal{V} \label{eq-master-binary}\tag{5}
\end{align}

The Objective Function (1) maximizes the fleet profit, which consists of the vessel-specific profit less the cost of hiring a spot ship to service a CoA cargo. Constraints (2) state that each contracted cargo is either transported by a vessel in the fleet or by a spot ship. Constraints (3) specify that optional cargoes may be serviced by a vessel in the fleet. Constraints (4) ensure that each vessel is assigned exactly one route. Finally, Constraints (5) define the domains of the binary schedule assignment variables.

### Test Instances and Results

The test instances were generated by randomly sampling a set of voyages performed by Western Bulk's Panamax and Supramax vessels in the Indo-Pacific region in 2022. All distances between ports were provided by Maritime Optima. Bunker prices were provided by Maritime Optima's BunkerEx interface. Spot ship costs were based on FearnPulse rates. Canal costs were based on the Suez Toll Calculator by Wilhelmsen. Loading/discharge rates were uniformly sampled from 7,500 to 30,000 MT/Day. Port costs were sampled from 40,000 to 100,000 USD. Vessel speeds were set to 14 knots. The minimum/maximum bunker capacity of a vessel was assumed to be 500/2,500 MT. Bunker consumption was assumed 25 MT/day during sailing and 2.5 MT/day in port.

Classes of test instances consist of five test instances. For example, C30V10B10 denotes the class of test instances that consist of 30 cargoes (10 CoA and 20 spot), 10 vessels and 10 bunker nodes. The arc flow and the path flow models were implemented in C++ leveraging the Gurobi Optimization, LLC (2022) C++ interface under an academic license. C++ was used to make route generation and schedule generation more efficient.

For all the generated test instances, the solver of the arc flow model was run for up to an hour. If the solver could not produce an optimal solution within one hour, the best primal and dual bounds were reported. The path flow solution method consists of route generation, schedule generation, and setting up and solving the path flow model. The combined solution method is referred to as the a Priori Column Generation (aPCG) approach. The aPCG approach was allowed to complete solving an ongoing test instance even though the total duration of the solution method exceeded an hour because the solution method would not otherwise provide a solution. However, if the solution time of a test instance exceeded an hour, larger-sized classes of test instances were not attempted to be solved.

Tables 4 and 6 summarize relevant characteristics of the models built from the generated test instances. Tables 5 and 7 display relevant results from solving the models constructed from the same test instances. In all tables, the Instance Class columns denote the classes of generated test instances, specifying the number of cargoes, vessels, and bunker nodes. Every class consists of five randomly generated test instances. The number of spot cargoes is twice that of CoA cargoes, and the More or Less in Owner's Option (MoLOO) flexibility limits were set to $$\pm$$ 10%. The Size columns signify the size of the test instances classes.

In Table 4, the CPUAF column states the solution times of the arc flow model in seconds averaged over the test instances in each class. The NVar, NBin, and NConstr columns display the average number of variables, binary variables, and constraints, respectively. Finally, the PH column gives the average number of days in the planning horizon for each test instance class. As the duration of the planning horizon is set to the largest randomly generated delivery node end-time of a test instance, each test instance has a different planning horizon. The column Solved in Table 5 shows the number of test instances within each class that were solved to optimality for the arc flow model within one hour. The BestBoundAF column displays the best dual bound found when solving the arc flow model, averaged over the test instances in each class. Similarly, the ObjValsAF column states the best primal bounds. The Gap column shows the average percentage difference between the best primal and dual bounds.

Table 4: Arc flow model characteristics averaged over the test instances in each instance class
Instance Class Size CPUAF NVar NBin NConstr PH
(seconds)       (Days)
C6V2B4 Small 86 346 232 1,253 69.7
C9V2B4 Small 153 469 330 1,780 75.8
C9V3B4 Small 2,890 647 443 2,435 72.1
C15V5B10 Medium 3,600 4,060 3,351 17,257 76.1
C30V5B10 Medium 3,600 9,921 8,666 43,489 76.9
C30V10B10 Large 3,600 20,002 17,510 87,926 81.5
C45V10B10 Large 3,600 32,900 29,360 146,380 80.1

Table 5: Arc flow model results averaged over the test instances in each instance class
Instance Class Size Solved BestBoundAF ObjValAF Gap
(USD) (USD) (%)
C6V2B4 Small 5/5 3.76e+06 3.76e+06 0.00
C9V2B4 Small 5/5 5.22e+06 5.22e+06 0.00
C9V3B4 Small 2/5 9.24e+06 6.75e+06 34.26
C15V5B10 Medium 0/5 2.23e+07 1.11e+07 104.39
C30V5B10 Medium 0/5 4.43e+07 1.28e+07 250.77
C30V10B10 Large 0/5 4.99e+07 2.76e+07 81.56
C45V10B10 Large 0/5 6.99e+07 2.62e+07 169.93

In Table 6, columns CPURG, CPUSG and CPUPF display the average run time in seconds for route generation, schedule generation, and solving the path flow model, respectively. The CPUaPCG column shows the total duration of the aPCG approach, including route generation, schedule generation, and setting up and solving the path flow model. The values are displayed in seconds and averaged over the generated test instances for each class. Experiments show that setting up the path flow model is relatively slow. Finally, column NR reports the average number of routes generated by the different classes of test instances.

In Table 7, the Solved column shows the number of test instances within each class that were solved to optimality by the aPCG approach within one hour. The NSS and NSC columns state the number of spot ships used and the number of spot cargoes transported in the optimal solution. Their values are averaged over the test instances in each class. Further, the ObjValaPCG column states the optimal solution found by the aPCG approach, averaged over the test instances in each class. Thus, if Table 7 included a Gap column, its values would all equate to zero. Finally, the $$\Delta$$ column reports the average percentage improvements of the optimal objective value found when solving the path flow model compared to the best objective value found when solving the arc flow model.

Table 6: Path flow model characteristics averaged over the test instances in each instance class
Instance Class Size CPURG CPUSG CPUPF CPUaPCG NR
(Seconds) (Seconds) (Seconds) (Seconds)
C6V2B4 Small 0.1 0.2 0.001 0.3 816
C9V2B4 Small 0.1 0.3 0.001 0.6 1370
C9V3B4 Small 0.1 0.3 0.002 0.6 1581
C15V5B10 Medium 12.8 23.2 0.055 56.3 90,194
C30V5B10 Medium 94.3 117.36 0.326 416.1 406,435
C30V10B10 Large 369.7 434.3 1.229 1,515.7 1,312,606
C45V10B10 Large 813.2 702.8 2.065 3,708.0 2,095,022

Table 7: Path flow model results averaged over the test instances in each instance class
Instance Class Size Solved NSS NSC ObjValaPCG $$\Delta$$
(USD) (%)
C6V2B4 Small 5/5 1 2 3.76e+06 0.00
C9V2B4 Small 5/5 1 2 5.22e+06 0.00
C9V3B4 Small 5/5 1 3 6.76e+06 0.15
C15V5B10 Medium 5/5 1 6 1.26e+07 13.51
C30V5B10 Medium 5/5 5 7 1.55e+07 21.09
C30V10B10 Large 5/5 0 14 3.14e+07 13.77
C45V10B10 Large 2/5 3 13 3.47e+07 32.44

In general, solving the path flow model is significantly faster than solving the arc flow model for the Small instance classes. For the C9V2B4 class, optimal solutions were obtained by solving the path flow model in about half a second on average, whereas solving the arc flow took about two and a half minutes. Tables 5 and 7 further show that the same objective values were found when solving the two models for the Small test instance classes. Moreover, the two test instances of the C9V3B4 class that were solved to optimality by the commercial arc flow model solver had run times of more than half an hour. In contrast, the aPCG approach solved all test instances of the C9V3B4 class to optimality with average run times of 0.6 seconds. The significant discrepancies in run time duration are attributable to the total number of feasible routes in the problems. With an average number of total routes being less than 1,600 for the Small classes, the complete enumeration of all feasible routes combined with solving the vessel schedule optimization problem for each route is an efficient approach.

Nevertheless, the run times of the aPCG approach also grow with the test instance size. As seen in column NR in Table 7, the number of feasible routes rapidly grows as the number of cargoes, vessels, and bunker nodes increases. For the Medium test instance classes, the number of feasible routes is, on average, in the range of a few hundred thousand. For these classes, the aPCG approach results in run times of a few minutes. However, for the C45V10B10 test instance class, the number of feasible routes is, on average, in the range of millions. Consequently, on average, the aPCG approach takes more than an hour to find the optimal solution. The C45V10B10 effectively shows the limitations of the aPCG approach. For this class of test instances, the exhaustive enumeration of all feasible routes combined with solving the vessel schedule optimization problem for each route demands considerable time. Another potential limitation of the aPCG approach concerns memory usage. Generating millions of feasible routes could quickly exhaust the available memory of any hardware used to run the aPCG approach. Although such problems were not encountered when solving the generated test instances in this report, it should be noted that the aPCG approach demands a significant amount of memory usage.

### The Arc Flow Formulation

The notation and formulation presented in this section follow the arc flow formulation of the Tramp Ship Routing and Scheduling Problem presented by Christiansen et al. (2007), the flexible cargo quantities extension formulated in Christiansen and Fagerholt (2014), and the bunker management extension devised by Vilhelmsen et al. (2014).

Let the fleet of vessels be represented by the set $$\mathcal{V}$$, indexed by $$v$$. Assume $$N$$ cargoes are available for transport, indexed by $$i$$. The sets of pickup and delivery nodes are given by $$\mathcal{N}^P = \{0, \ldots, N\}$$ and $$\mathcal{N}^D = \{N + 1, \ldots, 2N\}$$, respectively. These sets allow each cargo $$i$$ to be represented by a pickup node $$i \in \mathcal{N}^P$$ and a delivery node $$N + i \in \mathcal{N}^D$$. As such, the set of all cargo-related nodes is represented as $$\mathcal{N} = \mathcal{N}^P \cup \mathcal{N}^D$$. The set of pickup nodes $$\mathcal{N}^P$$ is partitioned into $$\mathcal{N}^P = \mathcal{N}^C \cup \mathcal{N}^O$$, where $$\mathcal{N}^C$$ and $$\mathcal{N}^O$$ denote the mandatory contracted and optional spot cargoes, respectively. A vessel has a node of origin $$o(v)$$, a geographical position at sea or in port. It is assigned an artificial destination node $$d(v)$$, which a vessel will travel to after its last planned discharge node service. Due to cargo capacity, depth, and cargo handling equipment restrictions, some vessels and nodes in $$\mathcal{N}$$ may be incompatible. Thus, the set of nodes a vessel $$v$$ may visit is represented as $$\mathcal{N}_v \subseteq \mathcal{N} \cup \{o(v), d(v)\}$$. $$(\mathcal{N}_v, \mathcal{A}_v)$$ defines a base network for each vessel $$v$$, where the set of arcs $$\mathcal{A}_v \subseteq \{(i, j) | i \in \mathcal{N}_v, j \in \mathcal{N}_v\}$$ represents all arcs traversable by vessel $$v$$ with respect to time and bunker consumption. The pickup and delivery nodes compatible with vessel $$v$$ are defined as $$\mathcal{N}^P_v =\mathcal{N}^P \cap \mathcal{N}_v$$ and $$\mathcal{N}^D_v = \mathcal{N}^D \cap \mathcal{N}_v$$ for convenience.

The next step is extending the base network $$(\mathcal{N}_v, \mathcal{A}_v)$$ to allow vessel visits at bunker option nodes. Let $$\mathcal{B}$$ be the set of bunker option nodes. Similarly, as for the cargo-related nodes, vessels and bunker option nodes may be incompatible. Thus, set $$\mathcal{B}_v \subseteq \mathcal{B}$$ denotes the bunker option nodes that vessel $$v$$ may visit. The set of base nodes $$\mathcal{N}_v$$ for vessel $$v$$ is extended to include the nodes in $$\mathcal{B}_v$$ such that $$\hat{\mathcal{N}}_v= \mathcal{N}_v \cup \mathcal{B}_v$$. Further, the set of base arcs $$\mathcal{A}_v$$ is extended by adding all arcs connecting nodes in $$\mathcal{N}_v \setminus d(v)$$ with nodes in $$\mathcal{B}_v$$ that are traversable by vessel $$v$$ with respect to time and bunker. The arcs from the destination node $$d(v)$$ are excluded as the next node to visit is unknown. As such, the extended set of arcs is defined as $$\hat{\mathcal{A}}_v= \mathcal{A}_v \cup \mathcal{A}^B_v$$, where $$\mathcal{A}^B_v$$ is the set of arcs for vessel $$v$$ connecting bunker option nodes $$\mathcal{B}_v$$ to the nodes in $$\mathcal{N}_v \setminus d(v)$$. The extended cargo-bunker network is thereby $$(\hat{\mathcal{N}}_v, \hat{\mathcal{A}}_v)$$.

Let $$T^{S}_{ijv}$$ be the sailing time from node $$i$$ to node $$j$$ for a given vessel $$v$$ and arc $$(i, j) \in \hat{\mathcal{A}}_v$$. If node $$j$$ is a bunker node, a fixed time is added to $$T^S_{ijv}$$, signifying the time spent bunkering. The time required to load or discharge one unit of cargo $$i$$ with vessel $$v$$ is defined as $$T^{Q}_{iv}$$. The variable voyage cost $$C_{ijv}$$ accounts for the costs of visiting node $$i$$ and the sailing costs from node $$i$$ to node $$j$$ for vessel $$v$$. It is important to note that the cost of the purchased bunker is not included in this parameter as the purchase of bunker will be modeled separately. Instead, let $$P_i$$ denote the cost of purchasing one unit of bunker available in bunker option node $$i \in \mathcal{B}$$. The bonus bunker remaining at the destination node $$d(v)$$ is represented as $$P$$ and calculated as the mean of all available bunker node prices. The charter cost of servicing a contracted cargo $$i$$ by a spot ship is denoted $$C^S_i$$. Further, the bunker consumption at sea is denoted $$B^S_{ijv}$$ for a vessel $$v$$ traversing the arc $$(i, j) \in \hat{\mathcal{A}}_v$$. The bunker consumed by vessel $$v$$ in node $$i$$ per unit of time is represented by $$B^P_{iv}$$. There is a unit revenue $$R_i$$ generated for transporting cargo $$i \in \mathcal{N}^P$$. For each node $$i \in \hat{\mathcal{A}}_v$$ visited by vessel $$v$$, there is a time window $$\left[ \underline{T}_{iv}, \overline{T}_{iv} \right]$$ for the pickup and delivery time, defined by the earliest time $$\underline{T}_{iv}$$ and the latest time $$\overline{T}_{iv}$$ of the visit. Additionally, the bunker level on board each vessel $$v$$ must be within a lower limit $$\underline{B}_v$$ and an upper limit $$\overline{B}_v$$ to ensure realistic and safe operation. The quantity of cargo that can be transported is flexible within interval $$\left[\underline{Q}_i, \overline{Q}_i \right]$$, where $$\underline{Q}_i$$ is the minimum quantity and $$\overline{Q}_i$$ is the maximum quantity that must be transported if cargo $$i$$ is serviced. Finally, the cargo carrying capacity of vessel $$v$$ is denoted $$K_v$$.

The binary variable $$x_{ijv}$$ is assigned the value 1 if vessel $$v$$ traverses the arc $$(i, j) \in \hat{\mathcal{A}}_v$$, and 0 otherwise. The variables $$t_{iv}$$ denote the time vessel $$v$$ begins service at node $$i$$. The cargo load on board a vessel $$v$$ when it leaves node $$i$$ is represented by the variable $$l^C_{iv}$$. Similarly, the variable $$l^B_{iv}$$ represents the bunker load on board vessel $$v$$ just after completing service at node $$i$$. The quantity of cargo $$i$$ transported by vessel $$v$$ is represented by $$q_{iv}$$, and the quantity of bunker purchased by vessel $$v$$ at bunker option node $$i$$ is given by the variable $$b_{iv}$$. To represent whether a contracted cargo is serviced by a spot ship, the binary variable $$z_i$$ is assigned the value 1, or 0 otherwise. Finally, the binary variable $$y_i$$ is assigned a value of 1 if the optional spot cargo $$i$$ is serviced, and 0 otherwise.

Tables 8 - 10 provide a summary of the sets, parameters, and variables introduced in this section, respectively.

Table 8: Model Sets
Set Notation Set Description
$$\mathcal{V}$$ Set of vessels
$$\mathcal{B}$$ Set of bunker option nodes
$$\mathcal{B}_v$$ Set of bunker option nodes that vessel $$v$$ may service
$$\mathcal{N}$$ Set of all cargo-related nodes
$$\mathcal{N}^D$$ Set of delivery nodes
$$\mathcal{N}^P$$ Set of pickup nodes
$$\mathcal{N}^C$$ Set of pickup nodes for the mandatory contracted cargoes
$$\mathcal{N}^O$$ Set of pickup nodes for the optional spot cargoes
$$\mathcal{N}_v$$ Set of nodes that can be serviced by vessel $$v$$, including the origin node and an
artificial destination node
$$\mathcal{N}^P_v$$ Set of pickup nodes that vessel $$v$$ may service
$$\mathcal{N}^D_v$$ Set of delivery nodes that vessel $$v$$ may service
$$\mathcal{A}_v$$ Set of all arcs traversable by vessel $$v$$
$$\hat{\mathcal{N}}_v$$ Extended set of nodes and bunker option nodes that vessel $$v$$ may service
$$\hat{\mathcal{A}}_v$$ Extended set of arcs between nodes and bunker option nodes that vessel $$v$$ may
traverse

Table 9: Model Parameters
Parameter Notation Parameter Domain Parameter Description
$$N$$   Number of available cargoes
$$T^{S}_{ijv}$$ $$v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v$$ Sailing time from node $$i$$ directly to node $$j$$ for
vessel $$v$$
$$T^{Q}_{iv}$$ $$v \in \mathcal{V}, i \in \hat{\mathcal{N}_v}$$ Time required to load or discharge one unit of
cargo at node $$i$$ with vessel $$v$$
$$C_{ijv}$$ $$v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v$$ Cost of servicing node $$i$$ and sailing directly from
node $$i$$ to node $$j$$ with vessel $$v$$, cost of
$$P_{i}$$ $$i \in \mathcal{B}$$ Price of purchasing one unit of bunker at the
bunker option $$i$$
$$P$$   Unit price of bonus bunker
$$C^{S}_{i}$$ $$i \in \mathcal{N}^P$$ Cost of servicing the cargo at node $$i$$ with
a spot ship
$$B^{S}_{ijv}$$ $$v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v$$ Total bunker consumption for vessel $$v$$ while
sailing directly from node $$i$$ to node $$j$$
$$B^{P}_{iv}$$ $$v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v$$ Port bunker consumption for vessel $$v$$ while in
node $$i$$
$$B^0_v$$ $$v \in \mathcal{V}$$ Initial bunker level on board vessel $$v$$
$$R_{i}$$ $$i \in \mathcal{N}^P$$ Revenue generated from transporting one unit
of cargo from node $$i$$
$$\overline{T}_{iv}$$ $$v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v$$ The latest time at which vessel $$v$$ may begin its
service at node $$i$$
$$\underline{T}_{iv}$$ $$v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v$$ The earliest time at which vessel $$v$$ may begin its
service at node $$i$$
$$\overline{B}_v$$ $$v \in \mathcal{V}$$ Maximum bunker level for vessel $$v$$
$$\underline{B}_v$$ $$v \in \mathcal{V}$$ Minimum bunker level for vessel $$v$$
$$\overline{Q}_i$$ $$i \in \mathcal{N}^P$$ Maximum quantity of the cargo at node $$i$$ to be
transported
$$\underline{Q}_i$$ $$i \in \mathcal{N}^P$$ Minimum quantity of the cargo at node $$i$$ to be
transported
$$K_{v}$$ $$v \in \mathcal{V}$$ Cargo carrying capacity of vessel $$v$$

Table 10: Model Variables
Variable Notation Variable Domain Variable Description
$$x_{ijv}$$ $$v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v$$ 1 if vessel $$v$$ services node $$i$$ just before node $$j$$,
0 otherwise
$$t_{iv}$$ $$v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v$$ Time of which vessel $$v$$ begins service at node $$i$$
$$b_{iv}$$ $$v \in \mathcal{V}, i \in \mathcal{B}_v$$ The bunker quantity purchased by vessel $$v$$ at node $$i$$
$$l^{C}_{iv}$$ $$v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v$$ Cargo load on board vessel $$v$$ just after completing
service at node $$i$$
$$l^{B}_{iv}$$ $$v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v$$ Bunker load on board vessel $$v$$ just after completing
service at node $$i$$
$$q_{iv}$$ $$v \in \mathcal{V}, i \in \mathcal{N}^P_v$$ Quantity of cargo $$i$$ that is transported by vessel $$v$$
$$y_{i}$$ $$i \in \mathcal{N}^O$$ 1 if the optional spot cargo is transported, 0 otherwise
$$z_{i}$$ $$i \in \mathcal{N}^P$$ 1 if cargo at node $$i$$ is serviced by vessel from
the spot market, 0 otherwise

\begin{split}
\max \quad & \overbrace{\sum_{v \in \mathcal{V}} \sum_{i \in \mathcal{N}^P} R_{i} q_{iv}}^{\text{Freight Rate} \times \text{Quantity}} - \overbrace{\sum_{v \in \mathcal{V}} \sum_{(i, j) \in \hat{\mathcal{A}}_v} C_{ijv} x_{ijv}}^{\text{Sailing Cost}} - \overbrace{\sum_{i \in \mathcal{N}^C} C^{S}_{i}z_{i}}^{\text{Spot Ship Cost}} \\
& - \underbrace{\sum_{v \in \mathcal{V}} \sum_{i \in \mathcal{B}_v} P_{i} b_{iv}}_{\text{Bunker Cost}} + \underbrace{\sum_{v \in \mathcal{V}} P \cdot (l^{B}_{d(v)v} - B^{0}_{v})}_{\text{Cost/Value of Surfeit/Excess Bunker}}\\
\text{S.T.} &
\end{split} \label{eq-objective-function} \tag{6}

\begin{align}
\sum_{v \in \mathcal{V}} \sum_{j \in \hat{\mathcal{N}}_v} x_{ijv} + z_{i} = 1 && i \in \mathcal{N}^C  \label{eq-contract-cargo}\tag{7}\\
\sum_{v \in \mathcal{V}}  \sum_{j \in \hat{\mathcal{N}}_v} x_{ijv} - y_{i} = 0 && i \in \mathcal{N}^O \label{eq-spot-cargo}\tag{8}\\
\sum_{j \in \hat{\mathcal{N}}_{v} } x_{o(v)jv} = 1 && v \in \mathcal{V} \label{eq-start-voyage}\tag{9}\\
\sum_{j \in \hat{\mathcal{N}}_v} x_{ijv} - \sum_{j \in \hat{\mathcal{N}}_v} x_{jiv} = 0 && v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v \setminus \{ o(v), d(v)  \} \label{eq-voyage-flow}\tag{10}\\
\sum_{i \in \hat{\mathcal{N}}_{v} }  x_{i d(v) v} = 1 && v \in \mathcal{V} \label{eq-end-voyage}\tag{11}\\
\sum_{j \in \hat{\mathcal{N}}_v} x_{ijv} - \sum_{j \in \hat{\mathcal{N}}_v} x_{N+i, jv} = 0 && v \in \mathcal{V}, i \in \mathcal{N}^P_v \label{eq-service-both-nodes} \tag{12}\\
x_{ijv}\left(t_{iv} + T^{Q}_{iv} q_{iv} + T^{S}_{ijv} - t_{jv}\right) \leq 0 && v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v  \label{eq-pickup-time-balance}\tag{13}\\
&& | i \in \mathcal{N}^P_v \nonumber\\[1em]
x_{ijv}\left(t_{iv} + T^{Q}_{iv} q_{i-N,v} + T^{S}_{ijv} - t_{jv}\right) \leq 0 && v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v  \label{eq-delivery-time-balance}\tag{14}\\
&& | i \in \mathcal{N}^D_v  \nonumber\\[1em]
x_{ijv}\left(t_{iv} + T^{S}_{ijv} - t_{jv} \right) \leq 0 && v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v  \label{eq-bunker-time-balance}\tag{15}\\
&& | i \notin \mathcal{N}^P_v \wedge i \notin \mathcal{N}^D_v \nonumber\\[1em]
t_{iv} + T^{Q}_{iv} q_{iv} + T^{S}_{i,N + i,v} - t_{N+i,v}  \leq 0 && v \in \mathcal{V}, i \in \mathcal{N}^{P}_{v} \label{eq-node-precedence}\tag{16}\\
\underline{T}_{iv} \leq t_{iv} \leq \overline{T}_{iv} && v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v \label{eq-time-limits} \tag{17}\\
x_{{ijv}}\left(l^{C}_{iv} + q_{jv} - l^{C}_{jv}\right)  = 0 && v \in \mathcal{V}, (i,j) \in \hat{\mathcal{A}}_v  \label{eq-pickup-cargo-balance}\tag{18}\\
&& | j \in \mathcal{N}^P_v \nonumber \\[1em]
x_{{i,N+j,v}}\left(l^{C}_{iv} - q_{jv} - l^{C}_{N+j,v}\right) = 0 && v \in \mathcal{V}, (i, N+j) \in \hat{\mathcal{A}}_v  \label{eq-delivery-cargo-balance}\tag{19}\\
&& | j \in \mathcal{N}^P_v \nonumber\\[1em]
x_{{ijv}}\left(l^{C}_{iv} - l^{C}_{jv}\right)= 0 && v \in \mathcal{V}, (i,j) \in \hat{\mathcal{A}}_v  \label{eq-bunker-cargo-balance}\tag{20}\\
&&  | j \in \mathcal{B}_v \nonumber\\[1em]
l^C_{o(v),v} = 0 && v \in \mathcal{V} \label{eq-initial-cargo}\tag{21}\\
q_{iv} \leq l^{C}_{iv} \leq \sum_{j \in \hat{\mathcal{N}}_v} K_v x_{ijv} && v \in \mathcal{V}, i \in \mathcal{N}^P_v \label{eq-load-pickup}\tag{22}\\
0 \leq l^{C}_{N+i, v} \leq \sum_{j \in \hat{\mathcal{N}}_v} \left(K_v - q_{iv}\right)  x_{N+i, jv}  && v \in \mathcal{V}, i \in \mathcal{N}^P_v \label{eq-load-discharge}\tag{23}\\
\sum_{j \in \hat{\mathcal{N}}_v} \underline{Q}_i x_{ijv} \leq q_{iv} \leq \sum_{j \in \hat{\mathcal{N}}_v}  \overline{Q}_i x_{ijv} &&  v \in \mathcal{V}, i \in \mathcal{N}^P_v \label{eq-cargo-carry-limits} \tag{24}\\
x_{ijv}\left(l^{B}_{iv} - B^S_{ijv} - B^P_{jv} T^{Q}_{jv}q_{jv} - l^B_{jv} \right) = 0  &&v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v \label{eq-bunker-pickup-balance}\tag{25}\\
&& | j \in \mathcal{N}^P_v \nonumber \\[1em]
x_{ijv}\left(l^{B}_{iv} - B^S_{ijv} - B^P_{iv} T^{Q}_{jv}q_{j-N,v} - l^B_{jv} \right) =  0  &&v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v \label{eq-bunker-delivery-balance}\tag{26}\\
&& | j \in \mathcal{N}^D_v \nonumber \\[1em]
x_{ijv}\left(l^{B}_{iv} - B^S_{ijv} + b_{jv} - l^B_{jv}\right) = 0  &&v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v \label{eq-bunker-bunker-balance}\tag{27}\\
&& | j \in \mathcal{B}_v \nonumber \\[1em]
x_{ijv}\left(l^{B}_{iv} - B^S_{ijv} - l^B_{jv} \right) = 0  &&v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v  \label{eq-bunker-od-balance}\tag{28}\\
&&| j \in \{d(v)\}\nonumber\\[1em]
l^{B}_{o(v)v} = B^0_v  &&v \in \mathcal{V} \label{eq-bunker-start} \tag{29}\\

\sum_{j \in \hat{\mathcal{N}}_v}  \left(\underline{B}_v +  B^S_{ijv} + B^P_{jv} T^{Q}_{jv} q_{jv}    \right) x_{ijv}  \leq  l^B_{iv} \leq \overline{B}_v \sum_{j \in \hat{\mathcal{N}}_v} x_{ijv} &&v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v \label{eq-bunker-cargo-limits}\tag{30}\\
b_{iv} \leq \left(  \overline{B}_v -  \underline{B}_v  \right) \sum_{j \in \hat{\mathcal{N}}_v} x_{ijv} && v \in \mathcal{V}, i \in \mathcal{B}_v \label{eq-bunker-load-limits}    \tag{31}\\
x_{ijv} \in \{0,1\} && v \in \mathcal{V}, (i, j) \in \hat{\mathcal{A}}_v \label{eq-take-cargo}\tag{32}\\
y_i \in \{0,1\} && i \in \mathcal{N}^O \label{eq-take-spot}\tag{33}\\
z_i \in \{0,1\}&& i \in \mathcal{N}^C\label{eq-take-spot-ship}\tag{34}\\
t_{iv} \in \mathbb{R}^{+} && v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v \label{eq-positive-time}\tag{35}\\
q_{iv} \in \mathbb{R}^{+} && v \in \mathcal{V}, i \in \mathcal{N}^{P}_v \label{eq-positive-moloo}\tag{36}\\
b_{iv} \in \mathbb{R}^{+} && v \in \mathcal{V}, i \in \mathcal{B} \label{eq-positive-bunker-purchase}\tag{37}\\
l^{C}_{iv} \in \mathbb{R}^{+} && v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v \label{eq-positive-cargo-load}\tag{38}\\
l^{B}_{iv} \in \mathbb{R}^{+} && v \in \mathcal{V}, i \in \hat{\mathcal{N}}_v \label{eq-positive-bunker-load}\tag{39}
\end{align}

The objective function (1) aims to maximize the total profit for the deployed fleet. The first term specifies the revenue generated for both contracted and optional cargoes as a function of the quantity transported. The second term specifies the variable voyage costs associated with servicing a cargo. The final term on the first line represents the charter costs associated with servicing a contracted cargo with a spot ship. Together, the first line describes the objective presented by Christiansen et al. (2014). The fourth term describes the cost associated with purchasing bunker, and the fifth term models the value of the bonus bunker remaining at the destination node. Note that if the bunker level at the destination is higher than at the origin, this term will positively contribute to the objective. Conversely, a lower bunker level at the destination than at the origin yields a negative contribution.

Constraints (7)-(12) denote the network flow constraints of the arc flow formulation. Constraints (7) ensure that all contract cargoes are transported by a vessel in the predefined fleet or by a spot ship. In contrast, Constraints (8) state that all spot cargoes must be transported by a vessel in the predefined fleet or not at all. Constraints (9) ensure that all vessels begin a schedule from their origin node and travel to a pickup node, a bunker option node, or directly to their destination node. Constraints (10) denote the flow conservation, while Constraints (11) ensure that all vessels travel to their destination node, either from a delivery node, a bunker option node, or their origin node. Together, Constraints (9) - (11) describe the flow on the sailing route used by vessel $$v$$. For each cargo $$i$$ and vessel $$v$$, Constraints (12) ensure that the same vessel services both the pickup and the cargo's corresponding delivery node.

Constraints (13)-(17) handle the temporal aspects of the problem. Constraints (13) state that if vessel $$v$$ travels directly from node $$i$$ to node $$j$$, then, the time at which a vessel $$v$$ begins loading of cargo in node $$j$$ must be greater than or equal to the time at which the vessel began loading of cargo in node $$i$$ plus the sum of a quantity dependent loading time at node $$i$$ and the sailing time from node $$i$$ to node $$j$$. Similarly, Constraints (14) specify the time progression for delivery nodes if vessel $$v$$ travel directly from node $$i$$ to node $$j$$. Note, however, that the cargo quantity variable $$q_{iv}$$ is only defined for pickup nodes, and as such the term that specifies the variable discharge time uses a different index. Constraints (15) describe a similar time flow for origin, bunker, and destination nodes as there is no cargo loading or discharging at these nodes. Constraints (16) are precedence constraints forcing the pickup node to be serviced before its corresponding delivery node. The inequality signs of Constraints (13) - (16) give vessels the option to wait before beginning their service at node $$j$$. Finally, Constraints (17) are the time window constraints of within a vessel must begin their service at node $$i$$.

Constraints (18)-(24) describe the restrictions related to the handling of the cargo transported. Constraints (18) state that if a vessel $$v$$ travels directly from node $$i$$ to a pickup node $$j$$, then, the cargo load on board the vessel just after completing service at node $$i$$ plus the cargo loaded at node $$j$$ must be equal to the cargo level on board the vessel just after completing service at node $$j$$. Constraints (19) specify a similar cargo flow if vessel $$v$$ travels directly from node $$i$$ to delivery node $$N+j$$. Further, Constraints (20) enforce that if a vessel $$v$$ travels to a bunker node $$j$$, then the cargo load on board the vessel stays the same as there is no cargo loading or discharging at bunker nodes. Note that the equality sign in Constraints (18) - (20) ensure that no cargo is misplaced during the shipment. Constraints (22) describe the cargo load capacity restrictions at pickup nodes, i.e., that the cargo loaded by vessel $$v$$ at pickup node $$i$$ must be less or equal to the cargo load on board just after the vessel completes its service of node $$i$$ which must be less or equal to the total cargo carrying capacity of the vessel. Constraints (23) specify the cargo discharge capacity restrictions at delivery nodes, i.e., that the cargo load on board vessel $$v$$ just after completing its service at delivery node $$N+i$$ must be greater or equal to zero and less than or equal to the cargo load capacity of the vessel less the cargo loaded at pickup node $$i$$. Finally, Constraints (24) ensure that the cargo loaded at pickup port $$i$$ is within the minimum and maximum cargo load limits. The summations in Constraints (22), (23), and (24) are introduced to strengthen the formulation of the model by enforcing only one of the $$x_{ijv}$$ variables involved in the summation to be non-zero.

Constraints (25) - (31) describe the bunker constraints of the model. Constraints (25) state that if vessel $$v$$ travels directly from node $$i$$ to pickup node $$j$$, then, the bunker load on board just after completing service at node $$i$$ less the bunker consumption of sailing from node $$i$$ to node $$j$$ less the bunker consumed while loading at node $$j$$ must be equal to the bunker load on board just after completing service at node $$j$$. Constraints (26) specify the same balance for the delivery nodes $$j$$, using the similar indexing scheme as before. Constraints (27) specify the bunker load balance for bunker nodes $$j$$ at which there is no cargo handling. There is, however, bunker loading at node $$j$$ which is added to the bunker load balance. Constraints (28) describe the bunker load balance during travel to the destination node as there is no bunkering nor bunker consumption in port. The quality signs in (25) - (28) ensure that no bunker is lost during travel. Constraints (29) specify the bunker level on board vessel $$v$$ at the beginning of the planning period. Constraints (30) state the bunker load on board vessel $$v$$ just after completing its service at node $$i$$ must be greater than or equal to the minimum bunker level of the vessel plus the bunker consumed during sailing and cargo handling in node $$j$$ and less than or equal to the total bunker capacity of the vessel. Constraints (31) ensure that the bunker purchased by vessel $$v$$ at bunker node $$i$$ is less than or equal to the difference between the minimum and maximum bunker load capacity. The summations in (30) and (31) are again strengthening the formulation of the model, enforcing only one of the $$x_{ijv}$$ variables involved in the summation to be non-zero.

Finally, Constraints (32) - (39) specify the variables' domains.  The $$x_{ijv}$$ variables in Constraints (32) are binary as they specify whether a vessel $$v$$ travels directly from node $$i$$ to node $$j$$, or not. The $$y_i$$ variables in Constraints (33) are binary as they indicate whether a spot cargo is serviced, or not. The $$z_i$$ variables in Constraints (34) are binary as they state whether a spot ship is used to transport a contracted cargo or not. The remaining variables specified by Constraints (35)-(39) can take fractional values and are required to be non-negative. Constraints (35) state that the time variables $$t_{iv}$$ can never take on values from before the planning period begins. Constraints (36) denote that a vessel $$v$$ may never transport negative cargo from node $$i$$. Constraints (37) specify that a vessel $$v$$ may never purchase a negative amount of bunker at bunker node $$i$$. Constraints (38) ensure that a vessel $$v$$ may never have a negative cargo load on board. Finally, Constraints (39) state that there may never be a negative amount of bunker on board vessel $$v$$ at node $$i$$.

### Written by

Simen Omholt-Jensen
Student / Part-time Data Science