Planning Annual LNG Deliveries with Transshipment

Abstract

The introduction of transshipment ports in the liquefied natural gas (LNG) supply chain in recent years offers additional flexibility, but also challenges to the planning of the annual delivery program. We present a new variant of the LNG-annual delivery program (ADP) planning problem by considering transshipment as well as time-dependent sailing times. We present a continuous time formulation for the LNG-ADP problem and propose a rolling horizon heuristic to solve the problem. Both the model and heuristic were used to solve a case inspired by the Yamal LNG project. The computational results show that the heuristic provides good solutions within a relatively short amount of time, especially compared to the exact solution methods. However, there is a trade-off between computational time and solution quality when designing the rolling horizon heuristic. The results also show the impact storage capacity at the transshipment port has on the total cost.

1. Introduction

The recent years have witnessed an increasing demand for natural gas. This is partly due to natural gas being considered a cleaner source of energy when compared to other fossil fuels. Especially Asian countries such as China, India, and South Korea are expected to further increase their consumption of natural gas [1]. However, the major natural gas producers are Russia, Iran, Qatar and the US. For such long distances between the major producers and customers that cross continents, it is more economic and efficient to liquefy the natural gas and transport it with specialized liquefied natural gas (LNG) carriers. According to [2], the global LNG imports have been increasing since 2013, reaching 313.8 MT in 2018 and this trend is expected to continue. Supply of LNG is also increasing, for example in Russia, where the annual production capacity of the Yamal LNG project will increase from 16.5 million tonnes per annum (MTPA) to approximately 17.5 MTPA once the fourth train is ready for production [3]. In addition, the Arctic-2 LNG project with an annual production capacity of 19.8 MTPA is currently under construction [4]. This expansion of the markets brings forth greater needs for transporting LNG.

A typical LNG supply chain is shown in Figure 1, where the focus of this paper is highlighted by the dashed box. After being extracted, the natural gas is liquefied and stored in tanks at the production terminals. It is then transported by LNG carriers to the destination terminals, where it is stored and eventually regasified. Finally, it is delivered to the end users.

In practice, the LNG producer enters into long-term contracts, which specify monthly and/or annual delivery quantities required at the destination ports [5]. It is the producer’s responsibility to deliver LNG to the customers, while also being in charge of production and inventory management. The producer creates an annual delivery program (ADP) for each customer, which specifies the loading and unloading port for each shipment of LNG as well as departure and arrival dates for a period of 12 to 18 months [6]. Usually, the objective of the ADP is to serve the long-term contracts at minimum cost. When planning the ADP, inventory management has to be considered as it is of vital importance in the LNG industry. The inventory levels must be maintained within given limits. If the inventory level is too high or too low, the production may be forced to shut down, which may cause severe losses and hence is undesirable. In addition to inventory limits, the planned ADP must feasible with regard to berth availability at the different ports. It is common to commit most, but not all, of the production capacity to medium- and long-term contracts [7]. The Russian liquefaction plants of the Yamal LNG and Sakhalin Energy projects for example, have entered into medium- and long-term contracts accounting for approximately 90% and 95% of the annual production, respectively [2]. The excess production of LNG is usually sold to different spot markets, providing an opportunity for the producer to gain extra revenue or release pressure on the production inventory. Choosing which spot market to serve and when is usually also part of the ADP planning process.

We present a LNG-ADP problem inspired by the Yamal LNG project located on the Yamal peninsula in Northern Russia. To transport LNG from Sabetta to the customers, specially designed ice-going LNG carriers are needed due to the presence of sea ice. These ice-going LNG carriers are designed to operate in the arctic waters of the Northern Sea Route, but are outperformed in open water by conventional LNG carriers. As the eastbound route from Sabetta to Asia might be unavailable during winter due to environmental conditions, the Asian customers are served through a transshipment port in Europe. By using the transshipment port, the producer avoids longer than necessary voyages in open water with ice-going LNG carriers, as conventional LNG carriers can take over the LNG and deliver it to customers. The direct and transshipment routes differ in distances, and consequently sailing time as well as costs. The producer then needs to decide not only on the departure dates, but also through which route to serve the customers. The introduction of a transshipment port into the network offers additional flexibility, but it also makes the planning of the annual delivery program even more challenging.

We refer to the problem as LNG-ADP with transshipment. In our problem, each load of LNG is only allowed to be transshipped at most once. The problem has emerged in recent years with the Yamal LNG project entering production and is a new variant of the tactical LNG-ADP (e.g., [5,6]). Our work contributes to research in LNG-ADP planning in the following aspects: (1) we enrich the LNG-ADP planning problem by introducing the option of using transshipment, (2) we present a continuous-time formulation for LNG-ADP problem and (3) we develop a rolling horizon heuristic (RHH) for the problem, which provides good solutions. The computational results clearly indicate that there is a trade-off between solution quality and required computational time. The results also show that the inventory capacity at the transshipment port has a huge influence on the total cost.

The remainder of this paper is organized as follows. Section 2 reviews the relevant literature for LNG-ADP problems and related LNG inventory routing problems (LNG-IRP). A detailed description of the problem is presented in Section 3, while Section 4 introduces the mathematical formulation. Section 5 describes the RHH followed by a computational study and discussion of results in Section 6. Finally, Section 7 concludes the paper.

2. Literature Review

As the LNG-ADP problem combines ship routing and scheduling with inventory management, it belongs to the category of maritime inventory routing problems (MIRP). MIRPs have been extensively studied in the literature and we refer the interested reader to [8,9,10] for more general surveys on maritime inventory routing problems.

The features of LNG-ADP problems are compared to more general MIRPs in [5]. The authors comment that the network structure of LNG-ADP is often simple with only one production port. In addition, the problem usually only considers full shiploads, implying that a voyage only contains one loading port and one unloading port. On the other hand, LNG-ADP can include multiple types of LNG products and the size of the LNG-ADP problem is often large due to the large size of fleet as well as the long planning horizon. Moreover, including contract management, i.e., the producer’s decision regarding to which spot market to sell the excess LNG in, also complicates the problem. They present a branch-and-cut algorithm for the problem as well as several types of valid inequalities. The computational study confirms that the problem is very complex, but some of the presented valid inequalities are very effective.

One of the first papers specifically addressing the LNG-ADP problem is [6] and the authors propose a construction and improvement (C&I) scheme for solving the problem which combines a rolling horizon scheme with a mixed integer programming (MIP) based heuristic. The rolling horizon heuristic solves a subproblem with a shortened horizon, then fixes the decision variables in chronological order and rolls forward, until it finally reaches the end of the planning horizon. In the subsequent improvement stage, only a subset of binary variables are created to reduce the solution space and achieve a better computational performance. The same problem is also solved in [7] using a different C&I heuristic. In the construction stage, a set of solutions is generated by means of multi-start approach based on a greedy principle where unfulfilled contracts and less flexible carriers are assigned first. During the improvement stage, two different heuristics are applied to improve the solution. The first one is based on local search, while the second is same MIP heuristic as the one presented in [6]. The heuristics can be used independent of each other, but using the MIP heuristics after the local search produces better results. A similar problem is also solved in [11]. The authors apply the Dantzig–Wolfe decomposition framework treating routing decisions for individual carrier as subproblem and employ both branch-and-bound and local search to solve the subproblems.

The LNG-ADP problem is also addressed in [12]. The authors propose a compact formulation of the model which takes advantage of the homogeneous fleet. The compact model is solved both using exact methods and a fix-and-relax heuristic. For the fix-and-relax heuristic, the whole planning horizon is first divided into several periods. The integer requirements are then relaxed in all periods except the first one and the problem is solved. The integer variables of the first period are then locked while the integer requirement are reintroduced for the second period. This continues until a feasible solution for the entire planning horizon has been found. The authors report that the fix-and-relax heuristic reduces run time by up to 95%.

A comprehensive tactical LNG-ADP problem is presented in [13] where the decisions are not limited to delivery planning and inventory management, but extended to the production rate at the production terminal. The authors also allow split-cargoes and include penalties for under- and over-deliveries. They design a vessel routing heuristic to solve the proposed model. The contracts are ranked. Then, one by one, vessels are assigned to the contracts. In a following step, the departure and arrival times are fixed before cargo volumes to the different contracts are decided. By using multiple ranking criteria, the solutions can be diversified.

The LNG-ADP problem is closely related to the LNG-IRP. Both LNG-ADP and LNG-IRP include ship routing and scheduling as well as inventory management, but the network structure of LNG-IRP is often more complex with more than one production port and allowing for less-than-shipload deliveries [5]. Less-than-shipload deliveries are commonly modeled by dividing the ship into several compartments, where one or more compartments can be delivered at different ports. Nevertheless, the two problems often share some common modeling and solution techniques. The first LNG-IRP problem in the literature was introduced by [14]. In addition to deciding delivery quantities, inventory management as well as routing and scheduling for a heterogeneous fleet, the authors also include the daily production rate as decision variable. The model allows for partial unloading in the model and aims at maximizing the total profit while satisfying vessel capacity and inventory limits. The fix-and-relax heuristic developed in [15] for solving a LNG-IRP can also be applied to the model presented in [12] for the LNG-ADP. The LNG-IRP models presented in [16,17] both arise from the same real-world case. According to [16], the model can also be used to develop an ADP. In [17], the problem is solved using a large-neighborhood-search (LNS) framework, where all binary variables are fixed except those concerning two chosen carriers. In [16], the LNS framework is extended into a C&I heuristic. A greedy construction stage assigns the vessel with the highest capacity to the port with the highest demand, iterating the production port and time period. In the first improvement phase they fix the routing decisions and adjust the departure time within a time window, followed by a two-carrier local search. In [18], the construction stage uses both a rolling time scheme with a round-trip-voyage subproblem and a greedy randomized adaptive search procedure. In addition to the time window and two-carrier local search, the authors add a third local search called singleton search where one solution element, e.g., one carrier or one terminal, is freed while other solution variables remain fixed.

A feature that distinguishes LNG transportation from other maritime cargo is the presence of boil-off. Due to the low boil point, LNG evaporates (boils-off) everyday and the boil-off gas is commonly used as fuel for the carrier [16]. Most of the literature above mention boil-off, but only [14,16,18] explicitly include it in the model. In addition, sailing times are considered to be constant throughout the year in most of the literature except for [16]. For an ADP problem with a planning horizon of around one year, seasonal variations in speed can be expected, which may have an impact of the quality of solutions. Using constant sailing time may cause scheduling issues in practice, and hence seasonal sailing time should be addressed in the model.

LNG-ADP problems are often formulated with discrete time, see e.g., [5,6,7,12,13], whereas continues time formulations are common in MIRPs, see [19,20,21]. The majority of LNG-IRP and almost all LNG-ADP are solved with heuristics, usually consisting of one or more elements such as rolling horizons, greedy construction heuristics and local search improvement heuristics.

To the best of our knowledge, transshipments have not been previously considered in LNG-ADP problems. Our paper extends the existing literature in this respect. Including transshipment of LNG in the LNG-ADP model allows a more detailed modeling of the physical flow of LNG, e.g., to capture re-exports arising from Russia’s liquefaction plants on the Yamal peninsula. In addition, we provide new a continuous time formulation for the LNG-ADP problem and show how to solve it using a rolling horizon heuristic.

3. Problem Description

When developing the ADP, the LNG producer aims at minimizing the total cost of serving long-term delivery contracts. LNG not delivered to long-term customers may be sold in a spot market for additional revenues. The objective is then to minimize the sum of transportation costs plus penalties for deviating from specified delivery volumes in the long-term contracts minus revenues from spot market sales. The decisions include the departure times for the LNG carriers from both the production port and transshipment port as well as their destination, the amount to deliver each month to the long-term customers, through which route to serve them, and whether to serve a spot market. The proposed plan has to be be feasible regarding availability of the routes, inventory requirements and port availability.

We consider one production port and one transshipment port and multiple customers that are either long-term customers or spot markets. Depending on how a customer can be reached from the production port, we distinguish between three types of customers. Figure 2 illustrates the network structure considered in this paper. Long-term customers and spot market of type 1 can only be visited through the transshipment port. Long-term customers and the spot market of type 2 may be visited either directly from the production port or via the transshipment port. However, the direct route is only seasonally available (represented by the dashed line in Figure 2). Long-term customers of type 3 are always visited directly from the production port. Note that each type of customer may contain one or more customers.

The producer operates a fleet consisting of two different types of LNG carriers. The two types differ in size and cost as well as operation area. Carriers of type A are dedicated to load LNG at the production port and transport it to the transshipment port or customer ports of type 2 or 3. Carriers of type B on the other hand operate between the transshipment port and customer ports of type 1 or 2. Only full shiploads are considered for both types of carriers. However, as a consequence of boil-off, the amount unloaded is smaller than the amount loaded. Due to seasonal variations in weather conditions, the speed of a vessel and thus the sailing time between two ports and varies depending on the starting time of a voyage. We assume that at the beginning of the planning horizon, all carriers of type A are located at the production port while all carriers of type B are at the transshipment port.

The required delivery volumes of LNG are specified in long-term contracts, both on a monthly and an annual basis. LNG contracts include some delivery flexibility, i.e., they allow for small deviations from the specified volumes, usually $\pm 10\%$ of the contracted volume [5]. On a monthly level, these deviations (both under- and over-deliveries) are penalized at fixed rates per m${}^3$. This fixed monthly deviation penalty for under-deliveries is chosen higher than the revenue from the spot market to ensure that the long-term contracts are prioritized over sales to the spot market. To provide an incentive for the producer to match annual demand, we penalize annual deviations using a piecewise linear penalty function. Smaller deviations, here defined as deviations less than one shipload, have a lower penalty per m${}^3$ than large deviations, i.e., deviations exceeding one shipload. Note that the annual penalty rate for large under-deliveries is chosen higher than the monthly penalty rate to even out monthly deviations (see also [7]). Figure 3 illustrates the annual penalty function. There are no delivery requirements for the spot markets.

The producer is also responsible for the inventory management of the storage tanks linked to the liquefaction terminal at the production port and the storage tanks at the transshipment port. The inventory level at all storage tanks must be maintained within given limits. There are no inventory requirements for the delivery ports.

In addition, the ADP is subject to port operation constraints. At all ports, only one LNG carrier can be served at a time. In addition, there is a berth time requirement between two consecutive visits at the same port. This berth time accounts for example for the time it takes to load or unload the LNG carrier.

4. Model Formulation

The problem of determining the optimal ADP is formulated as a mixed integer linear programming (MILP) problem. Before presenting the mathematical formulation, we first present the notation used in the model.

4.1. Notation

Let us introduce the following notation:

Sets
$G^C$	Set of delivery periods to customers
K	Set of departure periods determining sailing speed of the LNG carriers
$K^W$	Set of departure periods when the direct route to some customers is closed, $K^W\subset K$.
P	Set of ports.
$P_v^A$	Set of unloading ports that can be visited by carrier v.
$P^{AP}$	Set of ports that can be visited from the production port, $P^{AP}=\left\{P_v^A\right|v\in V^A\}$.
$P^{AT}$	Set of ports that can be visited from the transshipment port, $P^{AT}=\left\{P_v^A\right|v\in V^B\}$.
$P^C$	Set of customer ports, $P^C\subset P$.
$P^D$	Set of ports that can be visited both directly and via transshipment, $P^D\subset P^C$.
$P^L$	Set of long-term customers, $P^L\subset P^C$.
$P^P$	Set of production ports, $P^P\subset P$.
$P^S$	Set of spot markets, $P^S\subset P^C$.
$P^T$	Set of transshipment ports, $P^T\subset P$.
$S^A$	Set of possible nodes $(i,m)$ which indicates the m-th port call at port i.
V	Set of LNG carriers.
$V^A$	Set of LNG carriers of type A, $V^A\subset V$. These carriers only load at the production port.
$V^B$	Set of LNG carriers of type B, $V^B\subset V$. These carriers only load at the transshipment port.
Indices
g	Periods index for delivery, $g\in G^C$.
$i,j$	Port index, $i,j\in P$.
k	Period index for departure time, $k\in K$.
$m,n$	Port call index.
v	Carrier index, $v\in V$.
Parameters
$B_{ijv}$	Boil-off for a round trip from port i to port j with carrier v.
$C_{ijv}$	Cost of a round trip from port i to port j with carrier v.
$D_{jg}$	Demand (in m${}^3$) at customer j in delivery period g, $j\in P^L$.
$L^v$	Loading capacity for carriers of type v, $v\in V$.
$o_v$	Initial position of carrier v, $v\in V$.
$P_i$	Production rate at port i, $i\in P^P$.
$Q_j^{+}$	Penalty per m${}^3$ for annual over-delivery at customer j exceeding one shipload, $j\in P^L$.
$Q_j^{-}$	Penalty per m${}^3$ for annual over-delivery at customer j below one shipload, $j\in P^L$.
$Q_{jg}$	Penalty per m${}^3$ for over-delivery in delivery period g at customer j, $j\in P^L$.
$R_{jv}$	Revenue from selling one shipload LNG on carrier v to spot market j.
$s_{i0}$	Initial storage at port i, $i\in P^P\cup P^T$.
${\overline{S}}_i$	Upper bound on inventory level at port i.
${\underline{S}}_i$	Lower bound on inventory level at port i.
$T_{end}$	End of planning horizon.
$T_{eoy}$	End of year, $T_{eoy}<T_{end}$.
$T_g$	Start of delivery period g, $g\in G^C$.
$T_k$	Start of departure period k, $k\in K$.
$T_i^O$	Minimum operational time required by port i between two consecutive visits.
$T_{ijk}^v$	Sailing time of carrier v for traveling from port i to j when starting in departure period k.
$U_j^{+}$	Penalty per m${}^3$ for annual under-delivery at customer j exceeding one-shipload, $j\in P^L$.
$U_j^{-}$	Penalty per m${}^3$ for annual under-delivery at customer j below one-shipload, $j\in P^L$.
$U_{jg}$	Penalty per m${}^3$ for under-delivery in delivery period g at customer j, $j\in P^L$.
$ϵ$	Sufficiently small number.
${\overline{\mu }}_i$	Maximum number of port calls at port i, $i\in P$.
Decision variables
$q_j^{+}$	Over-delivery (in m${}^3$) at customer j exceeding one shipload over the entire planning horizon,
	$j\in P^L$.
$q_j^{-}$	Over-delivery (in m${}^3$) at customer j below one shipload over the entire planning horizon, $j\in P^L$.
$q_{jg}$	Over-delivery (in m${}^3$) at customer j in delivery period g, $j\in P^L$.
$s_{im}$	Inventory level at port i at the beginning of the $m-th$ port call, $(i,m)\in S^A$.
$t_{im}$	Start of voyage from node $(i,m)$, $(i,m)\in S^A$.
$u_j^{+}$	Under-delivery (in m${}^3$) at customer j exceeding one shipload over the entire planning horizon,
	$j\in P^L$.
$u_j^{-}$	Under-delivery (in m${}^3$) at customer j below one shipload over the entire planning horizon, $j\in P^L$.
$u_{jg}$	Under-delivery (in m${}^3$) at customer j in delivery period g, $j\in P^L$.
$w_{imv}$	1 if node $(i,m)$ is visited by carrier v, 0 otherwise, $(i,m)\in S^A,v\in V$.
$y_{im}$	1 if node $(i,m)$ is visited by any carrier, 0 otherwise, $(i,m)\in S^A$.
$z_{imjnvg}$	1 if carrier v starts the voyage from unloading node $(j,n)$ back to loading node $(i,m)$ in delivery
	period g, 0 otherwise, $i\in P^P\cup P^T,j\in P_v^A$.
${\alpha }_{imjnvk}$	1 if a round trip of carrier v from loading node $(i,m)$ to unloading node $(j,n)$ starts in departure
	period k, 0 otherwise, $i\in P^P\cup P^T,j\in P_v^A$.
${\beta }_{imk}$	1 if start of voyage from node $(i,m)$ is greater or equal to the start time of departure period $T_k$,
	0 otherwise, $i\in P^P\cup P^T$.
${\gamma }_{jng}$	1 if start of voyage from node $(j,n)$ is greater or equal to the start time of delivery period $T_g$,
	0 otherwise, $j\in P^L$.

4.2. Model Formulation

In this section, we present the mathematical formulation for the LNG-ADP problem with transshipment. The LNG-ADP with transshipment involves the routing and scheduling of LNG carriers, inventory management at production and transshipment ports as well as satisfying customer demand. In addition, the sailing time of the LNG carriers depends on the departure time due to seasonal variations in weather conditions. We have chosen to group the constraints according to the different aspects of the model.

4.2.1. Objective Function

The objective is to minimize the sum of transportation costs plus penalties for deviating from the long-term delivery commitments minus the revenues from sales to the spot markets. The first term is the transportation cost. The second and third term are the penalties for over- and under-deliveries in each delivery period and over the entire planning horizon, respectively. The last term is the revenue from spot sales.

$$\begin{array}{c}min{\sum }_{(i,m)\in S^A}{\sum }_{(j,n)\in S^A}{\sum }_{v\in V}{\sum }_{k\in K}{C}_{ijv}{\alpha }_{imjnvk}+\hfill \\ {\sum }_{j\in P^L}{\sum }_{g\in G^C}(Q_{jg}^{+}{q}_{jg}+U_{jg}^{-}{u}_{jg})+{\sum }_{j\in P^L}(Q_j^{+}{q}_j^{+}+Q_j^{-}{q}_j^{-}+U_j^{+}{u}_j^{+}+U_j^{-}{u}_j^{-})-\hfill \\ {\sum }_{(i,m)\in S^A}{\sum }_{\left\{(j,n)|j\in P^S\right\}}{\sum }_{v\in V}{\sum }_{k\in K}{R}_{jv}{\alpha }_{imjnvk}\hfill \end{array}$$

(1)

4.2.2. Routing and Symmetry Breaking Constraints

Constraints in Equations (2)–(4) are the flow conservation constraints for production port, transshipment port and customer ports, respectively. On the left hand side, $w_{imv}$ indicates if a port call $(i,m)$ is carried by v. The right hand side (RHS) of Equation (2) sum up the number of departure voyages from the production port by v over all departure periods. For transshipment ports in Equation (3), a port call can be an arrival of carrier type A or a departure of carrier type B. Equation (4) takes care of the arriving vessels in customer ports. Equation (5) links port calls to the carrier sailing. Equation (6) prevents direct voyages sailed by carrier type A when the direct route is unavailable. Equation (7) and Equation (8) are symmetry breaking constraints. Equation (7) makes sure the m-th port call is only possible if the previous port call is realized, while Equation (8) ensures that carrier v is used more often than $v+1$. Equation (9) makes sure that if there are multiple carriers sailing between the same loading port and unloading port, then the carrier that departs earlier also arrives earlier.

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle w_{imv}=\sum _{(j,n)\in S^A}\sum _{k\in K}{\alpha }_{imjnvk}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P,(i,m)\in S^A,v\in V^A,}\end{array}$$

(2)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle w_{imv}=\sum _{\left\{(j,n)|j\in P^P\right\}}\sum _{k\in K}{\alpha }_{jnimvk}+\sum _{(j,n^{\prime })\in S^A}\sum _{k\in K}{\alpha }_{imj{n}^{{}^{\prime }}vk}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^T,(i,m)\in S^A,v\in V,}\end{array}$$

(3)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle w_{imv}=\sum _{\left\{(j,n)|j=o_v\right\}}\sum _{k\in K}{\alpha }_{jnimvk}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^C,(i,m)\in S^A,v\in V,}\end{array}$$

(4)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{\left\{v\in V|i\in \left\{o_v\right\}\cup P_v^A\right\}}{w}_{imv}=y_{im}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}(i,m)\in S^A,}\end{array}$$

(5)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{\left\{(i,m)\in S^A|i\in P^P\right\}}\sum _{n\le {\overline{\mu }}_j}\sum _{v\in V^A}{\alpha }_{imjnvk}=0}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}j\in P^D,k\in K^W,}\end{array}$$

(6)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle y_{i,m-1}-y_{im}\ge 0}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}m>1,(i,m)\in S^A,}\end{array}$$

(7)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{(i,m)\in S^A}{w}_{imv}\ge \sum _{(i,m)\in S^A}{w}_{imv+1}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}v,v+1\in V^{A}orv,v+1\in V^B,}\end{array}$$

(8)

$$\begin{array}{cc}\multicolumn{1}{c}{(n-1)(1-\sum _{v\in V}\sum _{k\in K}{\alpha }_{imjnvk})\ge \sum _{m<l<{\overline{\mu }}_i}\sum _{b<n}\sum _{v\in V}\sum _{k\in K}{\alpha }_{iljbvk}}& \multicolumn{1}{c}{}\\ & \multicolumn{1}{c}{(i,m)\in S^A,(j,n)\in S^A,}\end{array}$$

(9)

4.2.3. Constraints for Grouping Departures and Contracted Deliveries

Equations (10)–(12) link a voyage to a departure period and ensure that the correct sailing time is used. Note that the continuous time formulation requires a maximum number of port calls for each port and even port calls without actual voyages are assigned a departure time. Inequalities in Equation (10) ensure that only real voyages are linked to the correct departure period. Similarly, the deliveries at long-term customers are handled by Equations (13)–(16). With Equations (13) and (14), all contracted deliveries are grouped based on the start time of the voyage leaving the delivery port. Note that port operations (unloading in case of a delivery port) mark the start of the voyage. Equation (15) ensures that a departure from production or transshipment port is linked to an arrival at the customer port. Equation (16) links a voyage to the corresponding delivery period.

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{(j,n)\in S^A}\sum _{v\in V}{\alpha }_{imjnvk}\le {\beta }_{imk}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P\cup P^T,(i,m)\in S^A,k\in K,}\end{array}$$

(10)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{k\in K}{\beta }_{imk}=1}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P\cup P^T,(i,m)\in S^A,}\end{array}$$

(11)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{k\in K}{T}_k{\beta }_{imk}\le t_{im}\le \sum _{k\in K}{T}_{k+1}{\beta }_{imk}-ϵ}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P\cup P^T,(i,m)\in S^A,}\end{array}$$

(12)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{g\in G}{\gamma }_{jng}=1}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}j\in P^L,(j,n)\in S^A,}\end{array}$$

(13)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{g\in G}{T}_g{\gamma }_{jng}\le t_{jn}\le \sum _{g\in G}{T}_{g+1}{\gamma }_{jng}-ϵ}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}j\in P^L,(j,n)\in S^A,}\end{array}$$

(14)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{g\in G}{z}_{imjnvg}=\sum _{k\in K}{\alpha }_{imjnvk}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}j\in P^L,(i,m),(j,n)\in S^A,v\in V,}\end{array}$$

(15)

$$\begin{array}{cc}\hfill \sum _{(i,m)\in S^A}\sum _{v\in V}{z}_{imjnvg}\le {\gamma }_{jng}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}j\in P^L,(j,n)\in S^A,g\in G^C,}\end{array}$$

(16)

4.2.4. Sailing Time Constraints

Equation (17) requires that, if $t_{im}$ and $t_{jn}$ belong to the same round trip voyage, the difference between the start time of the voyage at port i and the start time of the voyage at port j has to be greater or equal to port operation time at port i and the sailing time from port i to port j. Equation (18) states that after arriving at an unloading port, a carrier must have enough time to finish the port operations and return to the loading port before it starts a new voyage. Equation (17), together with Equation (18) ensure that a carrier can complete a round trip before starting a new voyage. Equation (19) ensures that time between two consecutive visits at the same port is longer than the minimum port operation time. As carriers operate all year round, at the end of year, carriers may start voyages for delivery in the coming year. Equation (20) ensures that the latest voyage from the production port starts before the end of year, while Equation (21) allows type B carriers to delivery cargo until the end of the planning horizon.

$$\begin{array}{c}\hfill {\displaystyle t_{im}+\sum _{v\in V}\sum _{k\in K}{T}_{ijk}^v{\alpha }_{imjnvk}+T_i^O{y}_{im}\le t_{jn}+T_{end}(1-\sum _{v\in V}\sum _{k\in K}{\alpha }_{imjnvk})}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ \multicolumn{1}{c}{i\in P^P,j\in P^{AP}ori\in P^T,j\in P^{AT};(i,m),(j,n)\in S^A,}\end{array}$$

(17)

$$\begin{array}{c}\multicolumn{1}{c}{{\displaystyle t_{jn}+\sum _{m<m^{{}^{\prime }}}\sum _{k\in K}{T}_{ijk}^v{\alpha }_{imjnvk}+T_j^O{y}_{jn}\le t_{i{m}^{{}^{\prime }}}+T_{end}(2-w_{i{m}^{{}^{\prime }}v}-\sum _{m<m^{{}^{\prime }}}\sum _{k\in K}{T}_{ijk}^v{\alpha }_{imjnvk})}}\\ \multicolumn{1}{c}{i\in P^P,v\in V^{A}ori\in P^T,v\in V^B;j\in P_v^A,(i,m^{{}^{\prime }}),(j,n)\in S^A,}\end{array}$$

(18)

$$\begin{array}{cc}\multicolumn{1}{c}{t_{im}-t_{i,m-1}\ge T_i^O{y}_{i,m-1}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}m>1,(i,m)\in S^A,}\end{array}$$

(19)

$$\begin{array}{cc}\multicolumn{1}{c}{t_{im}\le T_{eoy}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P,(i,m)\in S^A,}\end{array}$$

(20)

$$\begin{array}{cc}\multicolumn{1}{c}{t_{im}\le T_{end}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^T\cup P^C,(i,m)\in S^A,}\end{array}$$

(21)

4.2.5. Inventory Constraints

The inventory level at the production port is controlled by Equations (22)–(26). Equation (22) calculates the inventory for the first port call. Equation (23) records the inventory between two consecutive visits considering the volume loaded on the LNG carrier and produced since the last port calls. Equation (24) ensures that the inventory level does not fall below the lower inventory limit after loading. Equation (25) makes sure that the inventory level after the last port call plus the production between the last port call and the end of year does not exceed the upper inventory limit. As no departures are allowed from the production port after the end of year, no limits on the inventory level at the production port are imposed afterwards. Equation (26) ensures that the inventory never exceeds the upper inventory limit. The inventory constraints for the transshipment port are similar to the production port, but simpler. Equation (27) sets the inventory level before the first port call equal to the initial inventory level. Equation (28) handles the inventory change between two consecutive port calls. Depending on whether the previous port call is an LNG delivery or an LNG pickup, the new inventory level either increases by the amount of LNG unloaded (boil-off deducted) or decreases by the amount loaded. As there is no production at the transshipment port, the inventory level before port call $(i,m)$ is equal to the inventory level after the port call $(i,m-1)$. The inventory level at the transshipment port has to be within the inventory limits at all time, see Equation (29).

$$\begin{array}{cc}\multicolumn{1}{c}{s_{i1}=s_{i0}+P_i{t}_{i,1}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P,}\end{array}$$

(22)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle s_{im}=s_{i,m-1}-\sum _{v\in V^A}{L}^v{w}_{i,m-1,v}+P_i(t_{i,m}-t_{i,m-1})}}& \hfill \\ & \multicolumn{1}{c}{i\in P^P,m>1,(i,m)\in S^A,}\end{array}$$

(23)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle s_{im}-\sum _{v\in V^A}{L}^v{w}_{imv}+P_i{T}_i^O{y}_{im}\ge \underline{S_i}}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P,(i,m)\in S^A,}\end{array}$$

(24)

$$\begin{array}{cc}\hfill s_{i,{\overline{\mu }}_i}-\sum _{v\in V^A}{L}^v{w}_{i,{\overline{\mu }}_i,v}+P_i(T_{eoy}-t_{i\overline{{\mu }_i}})\le \overline{S_i}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P,}\end{array}$$

(25)

$$\begin{array}{cc}\multicolumn{1}{c}{s_{im}\le \overline{S_i}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^P,(i,m)\in S^A,}\end{array}$$

(26)

$$\begin{array}{cc}\multicolumn{1}{c}{s_{i1}=s_{i0}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^T,}\end{array}$$

(27)

$$\begin{array}{c}\multicolumn{1}{c}{{\displaystyle s_{im}=s_{i,m-1}+\sum _{(j,n)\in S^A|j=o_v}\sum _{v\in V^A}\sum _{k\in K}(L^v-B_{jiv}){\alpha }_{jni,m-1,vk}-\sum _{v\in V^B}{L}^v{w}_{i,m-1,v}}}\\ \multicolumn{1}{c}{i\in P^T,m>1,(i,m)\in S^A,}\end{array}$$

(28)

$$\begin{array}{cc}\multicolumn{1}{c}{\underline{S_i}\le s_{im}\le \overline{S_i}}& \multicolumn{1}{c}{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in P^T,(i,m)\in S^A,}\end{array}$$

(29)

4.2.6. Contract Management and Non-Negativity Constraints

Equation (30) determines the deviation in delivered volumes from demand in delivery period g for each long-term customer. Similarly, Equation (31), determines the total deviation from demand over the entire planning horizon. As annual deviations below one shipload and exceeding one shipload are penalized differently, Equations (32) and (33) help assign the correct volumes to the different deviation variables. Equations (34) and (35) are the non-negativity constraints for the continuous variables and the binary requirements for the binary variables, respectively. Please note that the indices have been omitted.

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{(i,m)\in S^A}\sum _{n\le {\overline{\mu }}_j}\sum _{v\in V}(L^v-B_{ijv})z_{imjnvg}+u_{jg}-q_{jg}=D_{jg}}}& j\in P^L,g\in G_C,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\multicolumn{1}{c}{{\displaystyle \sum _{(i,m)\in S^A}\sum _{n\le {\overline{\mu }}_j}\sum _{v\in V}\sum _{g\in G_C}(L^v-B_{ijv})z_{imjnvg}+u_j^{+}+u_j^{-}-q_j^{+}-q_j^{-}=\sum _{g\in G_C}{D}_{jg}}}& \multicolumn{1}{c}{j\in P^L,}\end{array}$$

(31)

$$\begin{array}{cc}\multicolumn{1}{c}{q_j^{-}\le L^v}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}v\in V^A,j\in P^L,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\multicolumn{1}{c}{u_j^{-}\le L^v}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}v\in V^A,j\in P^L,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\multicolumn{1}{c}{q,q^{+},q^{-},s,t,u,u^{+},u^{-}\ge 0,}& \hfill \end{array}$$

(34)

$$\begin{array}{cc}\multicolumn{1}{c}{\alpha ,\beta ,\gamma ,w,y,z\in \{0,1\}.}& \hfill \end{array}$$

(35)

5. Solution Method

As pointed out in [5], one of the challenges for LNG-ADP problem is the problem size, which is mainly due to the long planning horizon. A common attempt to reduce problem size has been to shorten the planning horizon and then solve consecutive problems in a rolling horizon framework until a feasible solution for the entire planning horizon is available, see e.g., [6,7,22].

Most of the papers on LNG-ADP, and especially the ones using RHH, use a discrete-time formulation, which provides a very strong link between decisions and time. For continuous time formulations, the use of RHH is less common as the link between decisions and time is weaker. Binary ship routing decisions indicate the sequence of port calls, but the departure times are decision variables as well. It is therefore difficult to know in advance which departure (or arrival) will be within the shortened planning horizon. Notable exceptions with respect to the use of RHH in inventory routing problems are [23,24], but none of the papers consider the LNG-ADP problem.

In Section 5.1, we present the main idea of RHH and address the differences between discrete time formulations and continuous time formulations in more detail. Our RHH implementation for solving the problem presented in this paper is then discussed in Section 5.2.

5.1. The Rolling Horizon Framework

In a rolling horizon framework, the whole planning horizon is divided into overlapping subhorizons. Each subhorizon is usually split into three different periods: The first period is the locked period where most if not all decisions are already taken and can no longer be changed. The central period is the second period where no decisions are restricted or relaxed and thirdly, the forecast period where certain relaxations apply. The problem is then solved consecutively for all subhorizons, such that the decisions made in the central period of iteration k become part of the locked period in iteration $k+1$, whereas the forecast period of iteration k becomes the central period in iteration $k+1$. The subhorizon of iteration $k+1$ is then extended by a new forecast period, see also Figure 4.

When using a discrete time formulation, the time index clearly indicates which of the three periods a decision variable belongs to. This makes it very easy to restrict or relax a decision. For LNG-ADP with discrete time formulations, it is common to relax the binary requirements on the shipping variables during the forecast period (see [6,7]). Once the problem has been solved in iteration k, the variables belonging to the central period are fixed in iteration $k+1$ as part of the locked period in the latter iteration. For the variables that belonged to the forecast period in iteration k, but that are now part of the central period of iteration $k+1$, binary requirements are reintroduced and the problem is solved again. This procedure continues until the entire planning horizon is covered.

For a continuous time formulation, it is not possible to a-priori relax the binary requirements of the shipping variables as the continuous departure time variable determines whether a particular departure belongs to the central period or the forecast period. This makes it more difficult to reduce the number of binary variables beyond the shorted planning horizon in order to achieve a computationally tractable subproblem. The use of a fix-and-relax scheme (see [12,15]) also becomes more difficult with a continuous time formulation.

5.2. RHH for Continuous Time Formulations

Reducing the planning horizon reduces the problem size of our LNG-ADP. However, this alone is not sufficient to guarantee a computationally tractable subproblem. We therefore propose a different approach for the forecast period: the forecast period is chosen as short as possible, but long enough to allow for LNG carriers that leave the production port in the central period to deliver LNG to the customer port before the end of the planning horizon. The length of the forecast period is set to be equal to one delivery period g. As there are no departures from the production port during the forecast period, we also relax the constraints on the upper inventory limit at the production port during this period. Note that departures from the transshipment port are possible during the forecast period, but also these carriers have to reach the customer port before the end of the planning horizon.

The set of decision variables to lock in the next iteration is determined based on the value of time variable indicating the start of a voyage, $t_{im}$. If $t_{im}$ is smaller than the end of the central period, all corresponding decisions, i.e., carrier departures and/or arrivals as well as inventory levels, become part of the locked period in the next iteration. Note that departures from the transshipment port during the forecast period are not locked.

To further reduce the number of binary variables in the subproblems, we dynamically adjust the maximum number of port calls for each iteration. For the production port, total production over the planning horizon of the subproblem can be calculated, providing an estimate for the number of port calls at the production port. For long-term customers, the demand profile over the planning horizon provides this estimate. The maximum number of port calls at the transshipment port can be estimated by doubling the port calls at production port, implying that all voyages will use the transshipment port. Note that the transshipment port has both loading and unloading port calls. Port calls to the spot markets can be estimated by considering the difference between the port calls at the production port and the long-term customer ports. Each of these estimates are slightly increased (usually by one or two port calls) to provide some flexibility for the optimal solution to deviate from these numbers. This way, the total number of port calls increases from iteration to iteration. However, port calls from previous iterations have been locked, so the number of binary decision variables is fairly stable across iterations and the subproblems remain tractable in each iteration.

6. Computational Study

This section presents the computational study of a case inspired by the Yamal LNG project. We first introduce the input data, before we present and discuss the results. The model was implemented and solved using FICO Xpress 8.5. All computations were carried out on Lenovo NextScale nx360 M5 computers with two 3.4 GHz Intel 6-core E5-2643 CPUs and 512 GB RAM running a Linux operating system.

6.1. Input Data

6.1.1. Port Information

We consider long-term customers and spot markets in both Europe and Asia. The production port is located in Sabetta on the Yamal peninsula and the transshipment port is Zeebrugge in Belgium.

Table 1. Port information.

Port Number	Type	Direct Route	via Transshipment	Direct Route Distance (nm)	Distance from Transshipment (nm)
1	Production	-	-	-	-
2	Long-term	Seasonal	Yes	5959.7	11,067.2
3	Long-term	Seasonal	Yes	6054	11,160.6
4	Spot	Seasonal	Yes	5348.9	11,058.2
5	Long-term	Seasonal	Yes	10,268.2	10,955.9
6	Long-term	Seasonal	Yes	6010.9	10,862.2
7	Long-term	Yes	No	2784.9	-
8	Long-term	No	Yes	-	2080.2
9	Spot	No	Yes	-	741.7
10	Transshipment	Yes	-	2661.8	-

gives an overview of the ports considered in the problem and their distances from the production and transshipment port. The second column indicates the port type. The third and fourth column show whether or not a port can be reached from the production port (direct route) and the transshipment port respectively. Only Asian customers are subject to seasonal availability of the direct route. The last two columns provide the distance (in nautical miles) of a port from the production port and the transshipment port, respectively. Note that the distance in the last column does not include the distance from the production port to the transshipment port. For all distances, we refer to [26].

The production rate is estimated for a single train with an annual production capacity of 5.8 MTPA, and it is assumed to be constant throughout the year at 34,279 m${}^3$/day. Port operations at all ports take 1 day to complete.

The spot market in Asia is located in Japan while for Europe it is in Spain. Both Chinese and Indian ports represent the Asian long-term customers whereas ports in the UK and France hold the long-term contracts in Europe. The demand profiles for the Asian and European long-term customers are shown in Figure 5 and Figure 6, respectively. There is some seasonal variation in demand for one of the Asian customers and both of the European customers. About 95% of the total annual production is allocated to long-term customers. Of the long-term demand, 54% is allocated to customers in Asia while the remaining 46% are allocated to customers in Europe.

The price at the Asian spot market is taken from [27]. The average Northeast Asia spot price has been $9.78 per million British Thermal Units (MMBtu) or approximately $226 per m${}^3$. The European spot market price varies significantly from season to season and year to year. According to the data provided by [28], the price climbed from $5.001/MMBtu in March 2017 to $9.522/MMBtu in September 2018 before it dropped to $3.633/MMBtu by January 2020. We use a price of $6.93/MMBtu for the European spot market.

6.1.2. Carrier Information

The fleet in this problem consists of 15 carriers of type A and 25 carriers of type B. Carriers of type A are owned by the producer while carriers of type B are time-chartered. The loading capacity for a carrier of type A is 172,600 m${}^3$ while it is set to about 161,000 m${}^3$ for a carrier of type B [2].

LNG carriers of type A are ice-going vessels with a design speed of 19.5 knots in open water [29]. As conditions along the Northern Sea Route can be challenging during winter, we assume that carriers of type A operate at 18 knots during summer and 16 knots during winter, resulting in different sailing times depending on season. Note that the eastbound route from Sabetta to Asia is closed during winter. Carriers of type B are open water vessels. These vessels operate at higher speeds than carriers of type A and are capable of a maximum speed of 21 knots [30]. We assume an average speed of 19 knots during summer and 18 knots during winter. The sailing time for both seasons is calculated based on the information above and average sailing time for the whole year is used to estimate boil-off and operational cost.

For the operational cost of carriers of type A, we consider only the crew costs of $9222/day, whereas we consider a daily charter rate of $60,000/day for carriers of type B, see [31]. Carriers sailing from Zeebrugge to Asia use the Suez canal and canal fees apply. Canal fees are assumed to be $350,000 for a round trip voyage [32].

We assume that type A carriers use boil-off as fuel while type B carriers use fuel oil. The fuel oil cost is included in the charter rate. We distinguish between natural boil-off which is due to the heat exchange with the surroundings and forced boil-off to produce fuel. The natural boil-off rate is set to 0.11% [33], and we estimate forced boil-off for fuel to be 520 m${}^3$/day according to [31,34]. Additional fuel costs are not considered.

6.1.3. Deviation Penalties

Deviations from customer demand are penalized according to the following principles: First, we assume that customers are more willing to accept over-deliveries than under-deliveries. Over-deliveries are therefore penalized less than under-deliveries. Moreover, the annual penalty should be higher than the monthly penalty. By doing so, we allow some flexibility for deviating from monthly demand, but provide an incentive to fulfill the annual delivery commitment. The penalty for under-deliveries is set higher than the price at spot market, to prevent selling in the spot market at the expense of the long-term customers. Minor deviation are here considered to be within one shipload and are penalized less than major deviations.

Based on the principles discussed above, we set a monthly penalty for over-delivery at $45/m${}^3$ and at $330/m${}^3$ for under-delivery. The annual penalty for over-delivery below one shipload is $10/m${}^3$ and $500/m${}^3$ for over-delivery exceeding one shipload. For under-delivery below one shipload, the penalty is set at $14/m${}^3$ and at $700/m${}^3$ for all under-delivery that exceeds one shipload.

6.1.4. Planning Horizon and Initial Values

The whole planning horizon is illustrated in Figure 7. The planning horizon consists of 12 months plus one additional month before and after the 12 months period, 14 months in total. The first month of the planning horizon is used to initialize the problem. At the beginning of month 1, all carriers of type A are located at the production port and all carriers of type B are available at the transshipment port. Initial inventory levels are set to 50% of the available storage capacity at the production port and 0 for the transshipment port. There is no customer demand during month 1, but departures can be scheduled to deliver LNG to customers in month 2 and/or build inventory at the transshipment port to serve customers from there.

During the last month of the planning horizon, we no longer allow any departures from the production port, but must still satisfy customer demand. This way, we prevent the mathematical model from ceasing operations early and ensure a link to operations beyond the end of the central 12 months period. Without the additional 14th month, the model might choose to no longer plan any departures once customer demand for the 12 central months are satisfied.

We consider 6 months of winter followed by 6 months of summer. The direct route to Asian customers is only available during the summer season.

6.2. Case Study

6.2.1. Overview of the Cases

We set up two cases for the computational study. The two cases differ in the storage capacity at the transshipment. We considered a high capacity in the first case (H) and a low capacity in the second case (L). We first tried to solve the two full cases using FICO Xpress and applied then our RHH with different lengths of the central period. The full cases have of a planning horizon of 14 months in total, as described in Section 6.1.4. We set the maximum run time for these full cases to be 24 h. We refer to them as HF for the full case with high capacity and LF for the full case with low capacity, respectively. For the RHH, we implemented three instances per case with the central period being set to be 1 month, 2 months and 3 months respectively. In all of these instances, we set the forecast period to be 1 month. We refer to the instances using the length of the planning horizon, i.e., H2 to H4 and L2 to L4.

Table 2. Overview of test instances.

Instance	Storage Capacity at Transshipment	Planning Horizon for RHH (months)
HF	High	13+1
H2	High	1+1
H3	High	2+1
H4	High	3+1
LF	Low	13+1
L2	Low	1+1
L3	Low	2+1
L4	Low	3+1

provides an overview of the main characteristics of the different instances, i.e., the storage capacity at the transshipment port and the planning horizon for the RHH.

All RHH instances are initialized with the same first iteration. In the first iteration we solve a problem with a 4 months planning horizon where the last month is the forecast period. We then fix the decision variables for the first two months. This ensures that all instances have the same starting point. Note that the last iteration may have a shorter planning horizon for some instances as the entire planning horizon is limited to 14 months.

For each iteration of the RHH, we set a maximum run time of 4 h. Solutions within 1% of optimality are accepted to speed up the solution process.

6.2.2. Computational Results

No feasible solution is produced for the full cases, HF and LF, within the 24 h run time limit. Not even the linear relaxation of the problems is solved to optimality. As such, no lower bound for the full problem is available. The poor performance from the full cases emphasizes the need for heuristic methods.

The computational results for all instances are summarized in

Table 3. Computational results.

Instance	Best Solution	Best Bound Last Iteration	Gap Last Iteration	Not Optimal/Total Iterations	Total Run Time (h)
HF	Not solved	-	-	-	>24
H2	542,340,683	460,334,852	15.12%	1/12	4.5
H3	224,109,041	224,091,574	0%	0/7	1.3
H4	217,427,397	202,055,902	7.07%	2/5	11.4
LF	Not solved	-	-	-	>24
L2	Infeasible	-	-	-	-
L3	Infeasible	-	-	-	-
L4	1,075,488,728	1,075,402,397	0%	1/5	5.1

. For each instance we provide the objective function value of the best solution found by the RHH, the best bound and gap reported from the last iteration of the RHH, the number of iterations that have been solved and how many of them did not solve to optimality as well as the run time for the RHH. Please note that the gap reported in the fourth column is not the optimality gap with respect to the original problem, but the optimality gap in the last iteration given the decisions taken in previous iterations.

We see from Table 3 that the objective function value of the best solution decreases as the planning horizon of the RHH instance increases. Considering instances H2 and H3, the total cost of the best solution of H3 is approximately 60% cheaper than the best solution produced by the H2 instance. It is worth pointing out that the last iteration of H2 is not solved to optimality, but is stopped after 4 h with an optimality gap of 15.12%. Nevertheless, the best solution of H3 is about 50% lower than the best bound for the H2 instance.

The performance of H3 is also better than H2 with regard to the computational time. H3 solves 7 iterations within 1.3 h whereas H2 solves 12 iterations in 4.5 h. However, the run time of H2 is caused by the last iteration that cannot be solved to optimality within 4 h. The previous 11 iterations are actually solved within 0.5 h, resulting in a less time per iteration compared to H3 where 7 iterations are solved within 1.3 h. Still, extending the central period from 1 month to 2 months results in both a considerable reduction in cost and run time.

Comparing H3 to H4, the trade-off between solution quality and run time becomes more apparent. On one hand, the objective function value of the best solution of H4 is about 3% better than the best solution of H3. On the other hand, it takes almost 10 times longer to find this solution. Even when subtracting 8 h for the two iterations that have not been solved to optimality, the average run time of one optimally solved iteration of H4 is roughly equal to the run time of the entire instance H3. Concluding whether it is worth to extend the central period by an additional month is here more difficult and will in a practical setting most likely depend on whether the decision maker prefers faster solutions over better solutions.

The pattern of longer planning horizons leading to better solutions can also be observed for the low storage capacity case. In fact, both L2 and L3 become infeasible in the last iteration. Only with a planning horizon of 4 months does the RHH produce a feasible solution for the entire planning horizon. Please note that the last iteration of L4 is solved to optimality, but iteration 3 was stopped after 4 h as it could not be solved to optimality within that time.

Iteration 3 is also one of the iterations that could not be solved to optimality for instance H4. The reason this iteration is computationally demanding is that this is the only iteration that includes months from both the winter season and the summer season. As a consequence the model has to consider different travel times for winter and summer, greatly increasing the number of binary variables ${\alpha }_{imjnvk}$. In addition, the direct routes become available, adding flexibility to the network and further increasing the solution space.

Overall, the performance of the instances with 3-month central period is the best among the three different instances we tested for the two cases. Using even longer planning horizons might further improve the solution, but it might no longer be possible to solve the different iterations to optimality as these become more complex.

6.2.3. Differences in Solution Structure

When comparing the objective function value of the L4 solution to the objective function value of the solution of H4, we see that the L4 solution is almost five times as expensive. The only difference between the two cases is the storage capacity at the transshipment port and must therefore cause this increase in objective function value. The number of port calls at the different ports are summarized in

Table 4. Port calls for the best solutions. For the transshipment port (port 10) only the number of loading port calls is presented.

Instance	Port 1	Port 2	Port 3	Port 4	Port 5	Port 6	Port 7	Port 8	Port 9	Port 10
H4	79	6	5	7	18	5	17	19	3	43
L4	76	6	5	5	19	6	19	15	1	29

. Ports 1 and 10 are the productions port and transshipment port, respectively. Ports 4 and 9 are the two spot markets with the remaining ports representing the long-term customers.

When comparing the number of port calls, the two solutions appear to be similar. With respect to the long-term customers in Asia (ports 2, 3, 5 and 6), both solutions deliver the same number of shiploads to ports 2 and 3. The L4 solution delivers one additional shipload to each port 5 and 6. For the European customers however, H4 delivers fewer shiploads to port 7, whereas L4 delivers considerably fewer shiploads to port 8. The differences are even more pronounced when considering the volumes delivered to the long-term customers. See Figure 8 for a comparison of annual demand of long-term customers and deliveries to these customers for the two solutions.

Figure 8 shows that the L4 solution incurs both a higher over-delivery penalty in port 7 and a higher under-delivery penalty in port 8 compared to the H4 solution. Based on the input data used in these instances, some of the L4 solution’s over-delivery to port 7 should have been sent to port 8 or one of the sport markets (ports 4 and 9) instead. The main difference between port 7 and ports 4, 8 and 9 is that port 7 is directly accessible from the production port, whereas the other three ports are only accessible through the transshipment port 10. It is also worth noting that the delivery amount to port 5 is almost the same for both solutions, despite L4 sending an addition shipload. This is due to the fact that port 5 in the L4 solution is mainly served directly from the productions port with ice-going type A carriers (using LNG as fuel), whereas the H4 solution predominately supplies port 5 from the transshipment port with type B LNG carriers (operating on fuel oil). Looking closer at the transshipment port (port 10) in Table 4, we see that the L4 solution only has about 2/3 of the port calls of the H4 solution.

The main reason for the fewer port calls at the transshipment port is the reduced storage capacity. Due to the different loading capacity of the carriers type A and type B, the volume of LNG unloaded by type A carrier is greater than the loading capacity of the type B carrier. The difference remains in storage at the transshipment port. With low storage capacity at the transshipment port, there is a moment when the inventory level at the transshipment port is too low for a full shipload of a type B carrier, but also too high to accommodate another shipload from a type A carrier. The optimization model therefore stops using the transshipment port in solution L4 around day 320 and only uses direct shipments from the production port to the customers. Figure 9 shows the inventory level at the transshipment port for both solutions.

The lower storage capacity is also the reason that instances L2 and L3 become infeasible: In those instances, the usage of the transshipment port stops even earlier than in L4. As some of the direct routes become unavailable in winter, only one port can be served from the production port while the liquefaction plant continues to produce LNG. Eventually, there is no more storage capacity left at the production site and the problem becomes infeasible.

With the higher storage capacity, instance H4 avoids the problem of infeasibility, while being able to produce much cheaper solutions. These results highlight the importance of a systems perspective when designing and dimensioning an LNG supply chain with transshipment, as underdimensioned storage facilities can lead to large increases in operational cost, which reflects current challenges with the existing infrastructure.

In reality, the transshipment capacity is closer to the lower capacity case. In fact, the transshipment capacity in Zeebrugge is close to one shipload [35] and thus even lower than in our low storage capacity case. One real-world solution to this problem has been to bypass the transshipment bottleneck at Zeebrugge by using ship-to-ship transshipments [36] and alternative transshipment ports [37].

The presented model and solution method can generate a feasible annual delivery program for a LNG producer, explicitly considering the transshipment of LNG. The analysis carried out in this paper also emphasizes the trade-off between the cost of a solution and the run time needed to find this solution, providing some insights that can be useful for setting up the RHH. The model can also be used to analyze different network designs, for example to identify bottlenecks in the transportation network. In particular, the consequences of different storage capacities at transshipment ports can be studied. Still, the model should be further improved, e.g., by including ship-to-ship transshipments and partial shiploads to better represent the real-world operations of LNG transportation. This will be subject of future research.

7. Conclusions

We present a LNG-ADP planning problem with the option of using a transshipment port, inspired by the Yamal LNG project. In our problem, some customers cannot be served directly from the production port for parts of the year and are therefore served through a transshipment port. We provide a continuous time formulation for the problem and propose a rolling horizon heuristic to solve the problem. The heuristic is able to solve realistic problem instances with acceptable run time.

The computational study also highlights the influence of the length of the planning horizon on the quality of the solution. The longer the planning horizon of the subproblems, the better the solution. Still, there is a trade-off to consider as the improved solution quality is paid for by considerable increases in run time. We also show that underdimensioning the storage capacity at the transshipment port can cause massive increases in operational cost.

This paper provides some initial insight in the problem of shipping LNG through transshipment ports. However, more research is required to further improve the understanding of the system, e.g., through modeling partial loading and unloading at the different ports.