Solving Chartering and Tramp Shipping 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.

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.

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 N_R 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.