1 Introduction

Ever since man first set sails across the oceans, transporting goods between continents has been an important part of the economy. Today, maritime transportation is the backbone of the global economy and United Nations Conference on Trade And Development (UNCTAD) has in its review of maritime transportation shown that the volumes, counted in tonnes loaded, of international maritime trade has almost doubled since 2000, from slightly above 6 ⋅ 109 tonnes loaded in 2000 to almost 12 ⋅ 109 in 2018 [18]. Many types of goods are transported and many different vessels are used to do it, but maritime transportation is traditionally divided into three segments: industrial shipping, tramp shipping, and liner shipping, [9]. An industrial operator controls its own fleet of vessels and strives to minimize the costs of transporting its own cargoes. An operator within tramp shipping has similarities with taxi services, as the vessels follow the cargoes that become available in the market. Liner shipping resembles bus line operations since the vessels follow published schedules and itineraries. UNCTAD among others divided the goods transported into four categories; tanker trade, main bulk, other dry cargo, and containers [18]. The growth of global containerized trade has been around 5.8% yearly since 2000 and has as of 2018 reached 152 million 20-ft equivalent units compared with 63 millions in 2000. Most of the containerized trade is within liner shipping and the economic importance of liner shipping is therefore obvious.

Many different planning problems within liner shipping have been analyzed using operations research. Meng et al. [11] outline some of these, including fleet size and mix and network design on the strategic level, fleet deployment and speed optimization on the tactical level, and cargo routing and rescheduling on the operational level. For more comprehensive overviews of planning problems within liner shipping, the reader is referred to [5], [17], [10], and [6]. A liner shipping network consists of a number of cyclic routes, called services, visiting two or more ports each. Each service is assigned a given vessel class and sailed with a given frequency. Demand is given between pairs of ports, and can be transshipped at one or more intermediate ports. This means that the vessel picking up the demand may not be the vessel delivering it, but instead that it can be handled by any service, even if neither origin nor destination ports are part of that service. Designing a liner shipping network is a notoriously hard problem called the liner shipping network design problem (LSNDP); see [4] for an introduction to LSNDP. The practical aspects that can be included when defining the LSNDP are numerous and there are therefore many versions of the problem. Common aspects are, according to [6], transit time constraints, transshipment costs, rejected demands, speed optimization, and the structure of the services. The methods used to solve the LSNDP include both heuristics and exact methods, but most of them are based on a mathematical formulation of the problem and mathematical programming techniques. Here, we review the exact methods that have been proposed; for a detailed description of solution methods for the LSNDP, both heuristics and exact methods, the reader is referred to [6].

The seminal paper by [1] is the first that presents an exact solution method for the LSNDP. The problem is formulated over a time-space network, where a weekly frequency is assumed. Transit time is not included, and transshipment costs are calculated once the network has been designed and are not included in the design phase. Demand can be rejected, and there is a revenue for fulfilled demand. Complex services are possible as long as the visits to one port are not on the same day. In [13], the authors notice that disregarding transshipment costs is a weakness and propose a four-index formulation to remedy this. Butterfly services, i.e., services where one port can be visited twice, are handled by enumerating the order of the arcs of a service. Transit time is not included, and each service is served by a single vessel. All demand must be served. In [12], the authors notice that butterfly services may be too restrictive and propose a formulation allowing for multiple visits to multiple ports on each service. Weekly frequencies on the services are assumed and demand can be rejected. An analysis of the effect of different structures of the services is conducted by [16]. Simple services, i.e., no ports are allowed to be visited more than once on each service, are compared with butterfly services and even more complex services. The complex structure is modeled using a two-layer network, where each port is represented in both layers and the connection between the layers is only between nodes corresponding to the same port. This allows for many different structures such as butterfly and pendulum where all ports except two are visited twice; see Fig. 1. Weekly frequencies are assumed and all demand must be served. In a feeder network, all demand either originates from or is destined for a single port. The design problem on these networks is studied by [14]. Transshipment is not allowed and each demand is either accepted or rejected. In [2], the authors present a new compact formulation for the same version of the problem as studied by [16]. New valid inequalities are introduced to strengthen the formulation. One of their main findings is that the gain from using a two-layer formulation compared with a one-layer formulation is on average 2–4% at the cost of a significant increase in solution time.

Fig. 1
figure 1

Two examples of complex service structures. Left: A butterfly service where one port is visited twice. Right: A pendulum service where all ports except the ends are visited twice. The numbers by the arcs show the order of traversal

Full size image

This paper presents an improved formulation and a new branch-and-price algorithm for the LSNDP. The work in [16] is extended with new valid inequalities and a new solution method for the subproblem is developed. The rest of the paper is organized as follows. In Section 2, a problem description is given and the master problem is formulated together with valid inequalities. Section 3 outlines the proposed branch-and-price algorithm and Section 4 provides a computational study. Finally, Section 5 gives some concluding remarks.

2 Problem Description and Mathematical Formulation

Following the classification of liner shipping network design problems in [6], the problem studied here does not contain transit time constraints, speed optimization, or the opportunity to reject demand. It includes transshipment costs and the possibility to use complex route structures. We now give a description of the problem, the notation, and the mathematical formulation.

The liner shipping network design problem studied in this paper consists of a set of ports \(\mathcal {P}\). There are weekly demands, Dij, to be transported from port i to port j and a given fleet of vessels to perform the transportation. All demands must be fulfilled; hence, it is not possible to reject demand. The vessels are divided into vessel classes, \(\mathcal {V}\), and there is an upper bound on the number of vessels, Yv, available from each class v. The vessels are assigned to services, \(\mathcal {S}\), where each service represents a given sequence of ports. The duration of a service is bounded, but the services can be complex. For any combination of vessel class and service, we let Avs be the number of vessels needed to keep a given (typically weekly) frequency on that service.

Since it is possible to split a demand between different services and also to transship, there are multiple ways to load and unload goods along a service. There is no bound on the transit time for the demands, and they can be transshipped as many times as needed. We associate a set of delivery patterns \(\mathcal {W}_{vs}\) with each possible combination of service s and vessel class v. For a given delivery pattern \(w \in \mathcal {W}_{vs}\), let \(Q_{ijvsw}^{L}\) be the number of units with final destination j that are loaded onto service s of vessel class v at port i per week. This includes units that are loaded at their port of origin, as well as units that have been transshipped from another service. Similarly, we define \(Q_{ijvsw}^{U}\) as the number of units with final destination j that are unloaded from service s of vessel class v and delivery pattern w at port i per week. To minimize the number of delivery patterns in \(\mathcal {W}_{vs}\), we may use the property that every convex combination of feasible delivery patterns is a feasible delivery pattern; and thus, we only need to explicitly represent the extreme delivery patterns in order to express all of them. In [7] and [15], the authors introduce deliver patterns in a similar way in problems with split deliveries. For each combination of service s, vessel class v, and delivery pattern w, we denote Cvsw as the total cost of using service s of vessel class v and delivery pattern w. The cost includes a fixed cost for operating the vessel, sailing costs, port fees, and costs for loading and unloading goods.

The problem then consists of selecting a set of services, and assign a vessel class and delivery pattern to each selected service, such that all demands are transported from their origin port to their destination port. To model this, we introduce variables pvs which is 1 if service s of vessel class v is used and 0 otherwise, and yvsw representing the number of times service s of vessel class v with delivery pattern w is used. Using these variables, and the notation defined above, a mathematical model of the problem may be defined as follows:

$$ \min z = \underset{v \in \mathcal{V}}{\sum} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{vs}}{\sum} C_{vsw} y_{vsw} $$
(1)

Subject to:

$$ \begin{array}{@{}rcl@{}} \underset{v \in \mathcal{V}}{\sum} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{vs}}{\sum} (Q_{ijvsw}^{L} - Q_{ijvsw}^{U} ) y_{vsw} &=& D_{ij} \qquad\quad i, j \in \mathcal{P}, i \neq j \end{array} $$
(2)
$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{vs}}{\sum} A_{vs} y_{vsw} &\leq& Y_{v} \quad \qquad v \in \mathcal{V} \end{array} $$
(3)
$$ \begin{array}{@{}rcl@{}} \underset{w \in \mathcal{W}_{vs}}{\sum} y_{vsw} - p_{vs} &=& 0 \qquad\qquad v \in \mathcal{V}, s \in \mathcal{S} \end{array} $$
(4)
$$ \begin{array}{@{}rcl@{}} p_{vs} & \in& \{ 0, 1 \} \qquad v \in \mathcal{V}, s \in \mathcal{S} \end{array} $$
(5)
$$ \begin{array}{@{}rcl@{}} y_{vsw} &\geq& 0 \qquad\qquad v \in \mathcal{V}, s \in \mathcal{S}, w \in \mathcal{W}_{vs} \end{array} $$
(6)

The objective function (1) minimizes the total cost associated with the services used. The flow balance constraints (2) make sure that the flow of goods is preserved. Since the destination of the goods unloaded when using a service is known, transshipment costs can be handled in the generation of the service and delivery pattern and can therefore be omitted from the master problem. Constraints (3) ensure that there are enough vessels of each class for the services used. The convexity constraints (4) state that different delivery patterns for the same service can be weighted together, and constraints (5) are binary restrictions. The non-negativity of the variables is expressed in constraints (6).

To illustrate the formulation and the variables, assume that we have a small network of six ports; see Fig. 2. Three services are shown; the dashed service, s1, is operated by a vessel with capacity 12 from class v1 while the solid service s2 and the dotted service s3 are both operated by a vessel with capacity 7 from class v2. There is one demand defined, D15 = 10, i.e., 10 units are requested from port 1 to port 5. The solution we illustrate is to transport 10 unit from port 1 using service s1; 7 units are transshipped at port 3 and continue to the destination at port 5 using service s2. The last 3 units are transshipped at port 4 and continue to the destination using service s3. Only allowing extreme delivery patterns, i.e., patterns that cannot be expressed as convex combinations of other patterns, we define two patterns w1 and w2 for service s1. The non-zero coefficients of these patterns are \(Q^{L}_{15v_{1}s_{1}w_{1}} = 10\), \(Q^{U}_{35v_{1}s_{1}w_{1}} = 10\) and \(Q^{L}_{15v_{1}s_{1}w_{2}} = 10\), \(Q^{U}_{45v_{1}s_{1}w_{2}} = 10\). Correspondingly, we get \(Q^{L}_{35v_{2}s_{2}w_{3}} = 7\) for the pattern w3 for service s2. Note that the coefficient for unloading at the destination is excluded. Since the extreme delivery pattern w4 for service s3 with coefficient \(Q^{L}_{45v_{2}s_{3}w_{4}} = 7\) handles too many units, we also need to define a pattern w5 for service s3 with no non-zero coefficients. Using the defined variables, the solution is now expressed in Table 1.

Fig. 2
figure 2

A small network with six ports operated by three service

Full size image
Table 1 A possible solution for D15 = 10 in Fig. 2
Full size table

2.1 Strengthening the Formulation

In [2], the authors successfully applied the method of expressing the convex hull of a low-dimensional polyhedron as a convex combination of its extreme points to strengthen their formulation of the LSNDP. They specifically targeted the feasible region defined by capacity constraints over subsets of ports and called this an inner representation of the polyhedron. We have adopted this idea and adjusted it for the master problem. The set of subsets of \(\mathcal {P}\) is denoted \(\mathcal {U}\), and let \(\mathcal {P}^{S}_{u}\) be the subset of ports under consideration and \(\overline {\mathcal {P}}^{S}_{u}\) the complement of \(\mathcal {P}^{S}_{u}\). The least quantity of goods that must be transported either in or out of \(\mathcal {P}^{S}_{u}\) is \({\max \limits } \{ \underset {i \in \mathcal {P}^{S}_{u}}{\sum } \underset {j \in \overline {\mathcal {P}}^{S}_{u}}{\sum } D_{ij}, {\sum }_{i \in \overline {\mathcal {P}}^{S}_{u}} {\sum }_{j \in \mathcal {P}^{S}_{u}} D_{ij} \}\). If the number of vessel classes, \(\vert \mathcal {V} \vert \), is small, we can generate all points \({\chi _{u}^{S}} = (\chi _{u1}^{S}, \ldots , \chi _{u \vert \mathcal {V} \vert }^{S})\) where \(\chi _{uv}^{S}\) is the number of times a vessel of class v enters \(\mathcal {P}^{S}_{u}\) that fulfill:

$$ \underset{v \in \mathcal{V}}{\sum} Q_{v} \chi_{uv}^{S} \geq \max \left\{ \underset{i \in \mathcal{P}^{S}_{u}}{\sum} \underset{j \in \overline{\mathcal{P}}^{S}_{u}}{\sum} D_{ij}, \sum\limits_{i \in \overline{\mathcal{P}}^{S}_{u}} \sum\limits_{j \in \mathcal{P}^{S}_{u}} D_{ij} \right\} $$
(7)

and is part of the convex hull of the feasible region defined by constraint (7) intersected with \(\mathcal {X} = \{ \chi _{uv}^{S} \in \mathbb {N} \cup \{0\}, v \in \mathcal {V} \}\). Define \({{\Lambda }^{S}_{u}}\) to be the number of such points for subset \(\mathcal {P}^{S}_{u}\) and let \(\chi _{luv}^{S}\) be the value \(\chi _{uv}^{S}\) for point l. Since the sequence of ports visited is known for each service, we can define \(G_{us}^{S}\) to be the number of times service s enters subset \(\mathcal {P}^{S}_{u}\). With \(\mathcal {P}^{S}_{u} = \{1, 3, 5\}\) and service \(s = 1 \rightarrow 2 \rightarrow 4 \rightarrow 5 \rightarrow 3 \rightarrow 6 \rightarrow 1\) we get \(G_{us}^{S} = 2\) since s enters \(\mathcal {P}^{S}_{u}\) twice, \(6 \rightarrow 1\) and \(4 \rightarrow 5\). With this, the constraints presented in [2] can be adjusted to:

$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{vs}}{\sum} G_{us}^{S} y_{vsw} &\geq& \sum\limits_{l = 1}^{{{\Lambda}^{S}_{u}}} \chi_{luv}^{S} \lambda_{lu}^{S} \qquad v \in \mathcal{V}, u \in \mathcal{U} \end{array} $$
(8)
$$ \begin{array}{@{}rcl@{}} \sum\limits_{l = 1}^{{{\Lambda}^{S}_{u}}} \lambda_{lu}^{S} &\geq& 1 \qquad\qquad\qquad u \in \mathcal{U} \end{array} $$
(9)
$$ \begin{array}{@{}rcl@{}} \lambda_{lu}^{S} &\geq& 0 \qquad\qquad\qquad u \in \mathcal{U}, l = 1, \ldots, {{\Lambda}^{S}_{u}} \end{array} $$
(10)

where \(\lambda _{lu}^S\) is the fraction of point l of subset u. Constraints (8)–(10) can be added to the master problem (1)–(6) to strengthen it.

3 Solution Method

A major challenge with the mathematical model presented in Section 2 is that it contains one y-variable, henceforth referred to as a column, for each combination of vessel class, service, and delivery pattern. The total number of columns thus grows very fast with the number of ports and demands and it becomes cumbersome to explicitly generate all columns, even for problems with just a few ports and demands. To circumvent this problem, we have developed a solution method based on branch-and-price [3] to effectively solve this problem. Branch-and-price is a branch-and-bound (B&B) tree search algorithm, where the linear program (LP) solved to obtain dual bounds for each node of the B&B-tree is solved using column generation [8].

Column generation is a method to solve linear programs that starts by solving a restricted version of the problem, containing only a subset of the columns (variables). This problem is usually referred to as the restricted master problem (RMP), and solving the RMP gives a pessimistic bound on the optimal value of the LP. To improve this pessimistic bound, one or more subproblems are solved, which finds the column with the lowest reduced cost (assuming that the LP is a minimization problem) based on the dual information from the solution of the RMP. If one or more negative reduced cost columns are found, they are added to the RMP, which is then re-solved to obtain an improved pessimistic bound and a new dual solution, and the previous step is repeated. If the subproblems cannot obtain any negative reduced cost columns, then we have proved that the current solution of the RMP is in fact optimal for the LP relaxation of the full problem in the current node, and we have obtained a valid dual bound for that node. The solution is checked for integrality, and if it is not integral, new nodes are created using an appropriate branching strategy.

In the following, we describe first how to price (calculate the reduced cost of) a column in Section 3.1, before giving a mathematical formulation of the subproblems that are solved to find the minimum reduced cost columns in Section 3.2. Then, in Section 3.3, we describe two heuristic strategies to find negative reduced cost columns quickly, before discussing the branching strategies, and their impact on the subproblem in Section 3.4.

3.1 Pricing of a Column

To define the reduced cost of a column, we introduce \(\beta _{ij}^{*}\) as the dual value from the flow conservation constraints (2), \(\gamma _{v}^{*}\) as the dual value from constraints (3), and \(\alpha _{uv}^{*}\) as the dual value from constraints (8). The dual values from the other constraints do not affect the reduced cost of a column and are therefore omitted here. There are port fees related to vessels visiting the ports and handling costs related to loading and unloading at the ports. The cost structure for the vessels is a fixed cost for each vessel used and distance-dependent sailing costs. We denote the costs in the following way; \({C_{v}^{F}}\) is the weekly fixed cost for operating a vessel of class v, \({C^{K}_{i}}\) is the cost of loading/unloading one unit of goods at port i, and Cijv is the sailing cost between ports i and j including the port fee at port j for vessel class v. Service s is represented as a cyclic list of port visits \(\pi _{1} \rightarrow \pi _{2} \rightarrow {\cdots } \rightarrow \pi _{\rho _{s}}\), where ρs is the number of port visits of service s. To simplify the notation required to describe the last leg of a cycle back to the starting point, we define \(\pi _{\rho _{s}+1} = \pi _{1}\). The reduced cost of a column representing service s with delivery pattern w for vessel class v can now be stated:

$$ \begin{array}{@{}rcl@{}} \overline{C}_{vsw} &=& ({C_{v}^{F}} + \gamma_{v}^{*})A_{vs} + \sum\limits_{i = 1}^{\rho_{s}} C_{\pi_{i} \pi_{i+1} v} + \underset{u \in \mathcal{U}}{\sum} \alpha_{uv}^{*} G_{us} + \\ &+& \sum\limits_{i = 1}^{\rho_{s}} \underset{j \in \mathcal{P}}{\sum} \left( (C_{\pi_{i}}^{K} + \beta_{\pi_{i}j}^{*}) Q_{\pi_{i} jvsw}^{L} + (C_{\pi_{i}}^{K} - \beta_{\pi_{i}j}^{*}) Q_{\pi_{i} jvsw}^{U} \right) \end{array} $$
(11)

3.2 Subproblem

The subproblems used to generate columns for the RMP can be stated for each vessel class. Each subproblem is formulated over a two-layer graph \(\mathcal {G} = (\mathcal {N}, \mathcal {A})\), where each physical port is represented by exactly one node in each layer. Each layer is a complete subgraph and the layers are only connected by arcs between the nodes representing the same port. Figure 3 shows the ports (the left part) and the corresponding graph (the middle part). Note that all arcs between nodes in the same layer are omitted to increase readability. All simple cycles in the graph are defined as feasible services. With two layers, each port can be visited at most twice on a service, but not all services fulfilling this criterion are feasible. The right part of Fig. 3 shows an example of a service that visits no port more than twice that cannot be represented as a simple cycle in the graph.

Fig. 3
figure 3

The left part shows the location of the ports. The middle figure shows the two-layer graph. Each port is duplicated and the graph consists of two layers (slightly shifted and marked gray and white respectively). The double lines between the nodes representing the same port are the only connections between the layers. Each layer is a complete graph, but these arcs are omitted to increase readability. The right part is an example of a route that visits no port more than twice, but that cannot be represented as a simple cycle in the graph in the middle. The numbers by the arcs is the order of arc traversal

Full size image

The graph has two layers and we define the nodes of the graph as \(\mathcal {N} = \{ 1,\ldots , 2 \vert \mathcal {P} \vert \}\). The nodes are divided into an upper layer consisting of nodes \(1,\ldots ,\vert \mathcal {P} \vert \) and a lower layer consisting of nodes \(\vert \mathcal {P} \vert + 1, \ldots , 2\vert \mathcal {P} \vert \). Each port i is associated with both node i and \(\vert \mathcal {P} \vert + i\). The set of arcs \(\mathcal {A}\) is the union of all arcs connecting nodes in the upper layer, all arcs connecting nodes in the lower layer, and arcs connecting the layers, i.e., \((i, \vert \mathcal {P} \vert + i)\) and \((\vert \mathcal {P} \vert + i, i)\) for all ports i. Hence, it is only possible to change layers by using arcs connecting nodes representing the same port. We also define \(\mathcal {A}^{S}_{u}\) to be all arcs (i,j) for which \(i \in \mathcal {P}^{S}_{u}\) and \(j \in \overline {\mathcal {P}}^{S}_{u}\), i.e., all arcs leaving the subset \(\mathcal {P}^{S}_{u}\) of ports. We introduce Qv for the capacity of a vessel of class v and Tijv for the sailing time in weeks between port i and j for a vessel of class v. The maximum length of a service, in weeks, is restricted by \(\overline {T}\).

Both the construction of the service and the amount to load and unload in the ports are decisions in the subproblem. We therefore introduce xij as a binary variable stating if arc (i,j) is part of the service or not, \(q_{ik}^{L}\) as the amount of goods with final destination k that is loaded at node i and \(q_{ik}^{U}\) as the amount of goods with final destination k that is unloaded at node i. We also need the variable a for the number of vessels deployed on the service. To handle the capacity of the vessels, we define fijk for the amount of goods traveling along arc (i,j) with final destination k. Subtours, i.e., two or more disjoint cycles, are elusive in this formulation. We eliminate subtours by letting one of the visited nodes in the service be an imaginary depot. The flow of an imaginary product along the service is forced to increase for each port visit, and can only decrease at the depot. A subtour without the depot can therefore not exist. To handle this in the formulation, we introduce the binary variable wi stating if node i is the imaginary depot of the service or not, and zij as the flow of an imaginary product from node i to j. Since each port is represented by two different nodes, we use p(i) to denote the port associated with node i. The subproblem of the LSNDP for vessel class v can now be formulated:

$$ \begin{array}{@{}rcl@{}} \min z_{SP}(v) &=& ({C_{v}^{F}} + \gamma_{v}^{*})a + \underset{(i,j) \in \mathcal{A}}{\sum} C_{ijv} x_{ij} + \underset{u \in \mathcal{U}}{\sum} {\sum}_{(i,j) \in \mathcal{A}^{S}_{u}} \alpha_{uv} x_{ij} \\ &+& \underset{i \in \mathcal{N}}{\sum} \underset{k \in \mathcal{P}}{\sum} (C^{K}_{p(i)} + \beta_{p(i)k}^{*}) q_{ik}^{L} + \underset{i \in \mathcal{N}}{\sum} \underset{k \in \mathcal{P}}{\sum} (C^{K}_{p(i)} - \beta_{p(i)k}^{*}) q_{ik}^{U} \end{array} $$
(12)

Subject to:

$$ \begin{array}{@{}rcl@{}} \underset{(j,i) \in \mathcal{A}}{\sum} x_{ji} - \underset{(i,j) \in \mathcal{A}}{\sum} x_{ij} &=& 0 \qquad\qquad\qquad i \in \mathcal{N} \end{array} $$
(13)
$$ \begin{array}{@{}rcl@{}} \underset{(j,i) \in \mathcal{A}}{\sum} x_{ji} &\leq& 1 \qquad\qquad\qquad i \in \mathcal{N} \end{array} $$
(14)
$$ \begin{array}{@{}rcl@{}} \underset{i \in \mathcal{N}}{\sum} w_{i} &=& 1 \end{array} $$
(15)
$$ \begin{array}{@{}rcl@{}} z_{ij} - (\vert \mathcal{N} \vert - 1) x_{ij} &\leq& 0 \qquad\qquad\qquad (i, j) \in \mathcal{A} \end{array} $$
(16)
$$ \begin{array}{@{}rcl@{}} \underset{(i,j) \in \mathcal{A}}{\sum} z_{ij} - \underset{(j,i) \in \mathcal{A}}{\sum} z_{ji} + \vert \mathcal{N} \vert w_{i} &\geq& \underset{(i,j) \in \mathcal{A}}{\sum} x_{ij} \qquad i \in \mathcal{N} \end{array} $$
(17)
$$ \begin{array}{@{}rcl@{}} \underset{(i,j) \in \mathcal{A}}{\sum} f_{ijk} - \underset{(j,i) \in \mathcal{A}}{\sum} f_{jik} + q_{jk}^{L} - q_{jk}^{U} &=& 0 \qquad\qquad\qquad j \in \mathcal{N}, k \in \mathcal{P} \end{array} $$
(18)
$$ \begin{array}{@{}rcl@{}} \underset{k \in \mathcal{P}}{\sum} f_{ijk} - Q_{v} x_{ij} &\leq& 0 \qquad\qquad\qquad (i, j) \in \mathcal{A} \end{array} $$
(19)
$$ \begin{array}{@{}rcl@{}} f_{ijk} - H_{vk} x_{ij} &\leq& 0 \qquad\qquad\qquad (i, j) \in \mathcal{A}, k \in \mathcal{P} \end{array} $$
(20)
$$ \begin{array}{@{}rcl@{}} \underset{(i,j) \in \mathcal{A}}{\sum}T_{p(i)p(j)v}x_{ij} &\leq& a \end{array} $$
(21)
$$ \begin{array}{@{}rcl@{}} a &\in& \{ 0,1,{\ldots} , \overline{T} \} \end{array} $$
(22)
$$ \begin{array}{@{}rcl@{}} x_{ij} &\in& \lbrace 0, 1 \rbrace \qquad\qquad (i,j) \in \mathcal{A} \end{array} $$
(23)
$$ \begin{array}{@{}rcl@{}} 0 \leq f_{ijk} &\leq& \sum\limits_{l \in \mathcal{N}} D_{lk} \!\!\!\!\qquad\qquad (i, j) \in \mathcal{A}, k \in \mathcal{P} \end{array} $$
(24)
$$ \begin{array}{@{}rcl@{}} q_{ik}^{L}, q_{ik}^{U} &\geq& 0 \qquad\qquad\qquad i \in \mathcal{N}, k \in \mathcal{P} \end{array} $$
(25)

The objective function (12) is the reduced cost of the column that is created. Constraints (13) are balance constraints and force all services to be closed loops while constraints (14) ensure that each node is visited at most once. Constraints (15)–(17) eliminate subtours. Constraint (15) states that exactly one node is the imaginary depot. Constraints (16) connect the flow on the imaginary product with the route and constraints (17) handle the flow balance, if node i is not the depot, the outgoing flow is one larger than the incoming, while the constraints are redundant if node i is the depot. Constraints (18) balance the flow of cargo entering and leaving a node, and capacity restrictions are enforce by constraints (19) and (20). We have here introduced \(H_{vk} = {\min \limits } \left \{ Q_{v} , {\sum }_{l \in \mathcal {N}} D_{lk} \right \}\) to increase readability. That the number of vessels employed on the service must be enough to guarantee a weekly frequency is stated in constraints (21). Constraints (22) declare the number of deployed vessel to be integer and set the upper bound on the duration of the services. The other variables are declared in constraints (23)–(25).

If the optimal objective value to a subproblem is negative, the solution is used to create a new column that is added to the RMP. The cost of the new column is calculated as:

$$ C_{vsw} = {C_{v}^{F}} a^{*} + \underset{(i,j) \in \mathcal{A}}{\sum} C_{ijv} x_{ij}^{*} + \underset{i \in \mathcal{N}}{\sum} \underset{k \in \mathcal{P}}{\sum} C^{K}_{p(i)} (q_{ik}^{L*} + q_{ik}^{U*} ) $$
(26)

The route defined by \(x_{ij}^{*}\) is transformed to a cyclic list of port visits, and a check is made to see if the service already exists. Since the cost of a column consists of both sailing costs and costs for loading and unloading, an already existing service with a new delivery pattern can be optimal in the subproblem. In this case, the index of the existing service is used and a new delivery pattern is associated with this service. The loading and unloading parameters \(Q^{L}_{ikvsw}\) and \(Q_{ikvsw}^{U}\) are directly taken from \(q_{ik}^{L*}\) and \(q_{ik}^{U*}\), respectively, and the number of vessels deployed on the service, Avs is assigned a∗. Finally, the parameter Gus is given the value \({\sum }_{(i,j) \in \mathcal {A}^{S}_{u}} x_{ij}^{*}\).

Henceforth, the subproblem defined by objective function (12) and constraints (13)–(25) is called the full subproblem. To increase the number of columns created with each time the full subproblem is solved, we supply the MIP solver with a callback routine to store all encountered solutions during the search which have a negative reduced cost. A new column is created for each of these solutions, so that a single iteration of solving the full subproblem may yield multiple new columns in the master problem.

3.3 Acceleration Strategies

A major success criterion for applying branch-and-price is that the structure of the subproblem(s) can be exploited. In most applications, the subproblems are solved using some form of dynamic programming, e.g., by using labeling algorithms. However, constructing cyclic sequences of nodes and assigning vessels and cargo flows to these sequences (which is the natural subproblem for the LSNDP) is not well suited for dynamic programming, as there are multiple dependencies between ports that are not necessarily adjacent in the sequence and no defining state at any of the nodes. This makes it difficult to design an algorithm adhering to Bellmann’s optimality principle.

As it is not easy to design an effective exact solution method to solve the subproblems, it is important to design good heuristics to find most of the negative reduced cost columns quickly, thus minimizing the number of times we have to solve the subproblem as a MIP, which is known to be time consuming for large instances of the problem. We have therefore implemented the following two heuristics to obtain negative reduced cost columns from a given service:

  • Improve delivery patterns There is (possibly) a large number of feasible delivery patterns associated with every service. The integral variables in the subproblem decide the port sequence (xij) and the number of vessels needed (a), i.e., the design of the service. Therefore, if we focus on a specific service, we only need to solve a linear program to find the best delivery pattern for that service. This can be done for each service individually.

  • Local search We do a local search on an already existing service by either removing or adding a port visit. If port p is removed in the sequence \(i \rightarrow j \rightarrow p \rightarrow k\), it becomes \(i \rightarrow j \rightarrow k\). If j and k represent the same port, one of them is removed as well. If port p is added to the sequence \(i \rightarrow j \rightarrow k\) after j, it becomes \(i \rightarrow j \rightarrow p \rightarrow k\). This is not done if j or k represents the same port as p. After a port is added to or removed from a service, it is verified that the service represents a simple cycle in the two-layer graph. If the check is affirmative, a linear program is solved to find the best deliver pattern and the corresponding reduced cost.

Even though these heuristics are computationally fast, preliminary testing showed that it was too time consuming to run these heuristics for all services that are a part of at least one column in the RMP. Instead, we partition \(\mathcal {S}\) into three disjoint subsets; \(\mathcal {S}^{B+}\) is the set of all services that are part of a column with a positive value in the current basis, \(\mathcal {S}^{B0}\) is the set of all services associated with a column that has a reduced cost of zero and also a value of zero in the current basis, and \(\mathcal {S}^{N}\) is the set of all services that are not part of any column in the basis.

The heuristics are run in a hierarchical fashion, and if at least one negative reduced cost column is found in one step, the algorithm returns it to the RMP, which is then re-solved. The hierarchy of these heuristics is as follows:

  1. 1.

    Improve the delivery patterns for all services in \(\mathcal {S}^{B+}\)

  2. 2.

    Improve the delivery patterns for all services in \(\mathcal {S}^{B0}\)

  3. 3.

    Perform local search on all services in \(\mathcal {S}^{B+}\)

  4. 4.

    Perform local search on all services in \(\mathcal {S}^{B0}\)

If none of these four steps returns a negative reduced cost column, we solve the mathematical model presented in Section 3.2 as a MIP using a commercial solver. However, since we only need a negative reduced cost column, we set a time limit on the MIP, and abort the solver if at least one negative reduced cost column has been found by that time, or as soon as one is found after this limit. Furthermore, we solve each MIP subproblem sequentially, and return to the RMP as soon as one of the MIPs has found a negative reduced cost column.

3.4 Branching

When a node has been solved to optimality, i.e., when no column with a negative reduced cost is found, it is checked with respect to the integrality restrictions (4). If the check is affirmative, we have a feasible integer solution; otherwise, we have to branch. Designing a branching scheme is intricate since it affects the structure of the subproblem. On one hand, the scheme should reach integer solutions fast, while on the other the changes to the subproblem should be small. Since the subproblems here are solved using a MIP solver, we allow for the structure to change by introducing branching constraints. We have adopted the branching scheme from [16], where a four-stage procedure is proposed. The procedure is hierarchical and one stage must be integral before branching on the next.

First, we branch on the number of vessels used of each vessel class. We define \(\hat {a}_{v} = \underset {s \in \mathcal {S}}{\sum } \underset {w \in \mathcal {W}_{vs}}{\sum } A_{vs} y_{vsw}\) and \(v^{\prime } = \arg \min \limits _{v \in \mathcal {V}} \vert \hat {a}_{v} - \lfloor \hat {a}_{v} \rfloor - 0.5 \vert \). If \(\hat {a}_{v^{\prime }}\) is fractional, we create the branching constraints:

$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} A_{v's} y_{v'sw} &\leq& \lfloor \hat{a}_{v^{\prime}} \rfloor \end{array} $$
(27)
$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} A_{v's} y_{v'sw} &\geq& \lfloor \hat{a}_{v^{\prime}} \rfloor + 1 \end{array} $$
(28)

and use them to define the new subproblems. The dual values from these branching constraints are multiplied with the variable a in the objective function of the subproblem of vessel class \(v^{\prime }\).

Second, we branch on the number of visits to each port by vessels from a given class. We define \(\hat {a}_{iv} = \underset {s \in \mathcal {S}}{\sum } \underset {w \in \mathcal {W}_{vs}}{\sum } B_{is} y_{vsw}\) where Bis is the number of visits to port i by service s. We let \((i^{\prime },v^{\prime }) = \arg \min \limits _{(i,v) \in \mathcal {P} \times \mathcal {V}} \vert \hat {a}_{iv} - \lfloor \hat {a}_{iv} \rfloor - 0.5 \vert \). If \(\hat {a}_{i^{\prime }v^{\prime }}\) is fractional, we create the branching constraints:

$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} B_{i's} y_{v'sw} &\leq& \lfloor \hat{a}_{i^{\prime}v^{\prime}} \rfloor \end{array} $$
(29)
$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} B_{i's} y_{v'sw} &\geq& \lfloor \hat{a}_{i^{\prime}v^{\prime}} \rfloor + 1 \end{array} $$
(30)

and use them to define the new subproblems. The dual value from the branching constraints (29) and (30) is added to all arc variables \(x_{i^{\prime }j}\) and \(x_{\vert \mathcal {P} \vert + i^{\prime },j}\) except \(x_{i^{\prime },\vert \mathcal {P} \vert + i^{\prime }}\) and \(x_{\vert \mathcal {P} \vert + i^{\prime },i^{\prime }}\) in the objective function of the subproblem of vessel class \(v^{\prime }\) to encourage/discourage more visits to port \(i^{\prime }\).

Third, we branch on the total number of times vessels of class v sails from port i to port j. We define \(\hat {a}_{ijv} = \underset {s \in \mathcal {S}}{\sum } \underset {w \in \mathcal {W}_{vs}}{\sum } B_{ijs} y_{vsw}\), where Bijs is the number of times service s visits ports i, and j in sequence and \((i^{\prime },j^{\prime },v^{\prime }) = \arg \min \limits _{(i,j,v) \in \mathcal {P} \times \mathcal {P} \times \mathcal {V}} \vert \hat {a}_{ijv} - \lfloor \hat {a}_{ijv} \rfloor - 0.5 \vert \). If \(\hat {a}_{i^{\prime }j^{\prime }v^{\prime }}\) is fractional, we create the branching constraints:

$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} B_{i^{\prime}j's} y_{v'sw} &\leq& \lfloor \hat{a}_{i^{\prime}j^{\prime}v^{\prime}} \rfloor \end{array} $$
(31)
$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} B_{i^{\prime}j's} y_{v'sw} &\geq& \lfloor \hat{a}_{i^{\prime}j^{\prime}v^{\prime}} \rfloor + 1 \end{array} $$
(32)

and use them to define the new subproblems. The dual value from these branching constraints is added to arc variables \(x_{i^{\prime }j^{\prime }}\) and \(x_{\vert \mathcal {P} \vert + i^{\prime },\vert \mathcal {P} \vert + j^{\prime }}\) in the objective function of the subproblem of vessel class \(v^{\prime }\).

Branching using the first three strategies is often enough to ensure an integer solution, but there are cases when the solution is fractional even though none of the above strategies indicates a candidate for branching. Figure 4 shows a solution with four services used by vessels from class v. To ease the presentation, we introduce \(y_{s} = \underset {w \in \mathcal {W}_{vs}}{\sum } y_{vsw}\). With y1 = y2 = y3 = y4 = 0.5 and assuming that this gives an integer number of vessels used, we see that none of the strategies above gives a candidate for branching.

Fig. 4
figure 4

A solution showing that branching on individual arcs is not sufficient

Full size image

For the solution not to be integral, there must be fractional y-variables, and among these there must be at least one pair of services with at least one arc in common and at least one arc only present in one service. To have a branching strategy targeting pairs of arcs will therefore detect a candidate for branching and is thus enough to ensure an integer solution. Because of the cyclic nature of the services, the pair can always be constructed such that the tail of one arc is the head of the other arc. In Fig. 4, we see that the arcs (l,i) and (i,k) are not a candidate for branching since both s1 and s2 include these arcs. If we instead choose (i,k) and (k,j), only s1 includes both arcs.

The fourth strategy is therefore to branch on pairs of arcs having one port in common. We define \(\hat {a}_{ikjv} = \underset {s \in \mathcal {S}}{\sum } \underset {w \in \mathcal {W}_{vs}}{\sum } B_{ikjs} y_{vsw}\) where Bikjs is the number of times service s visits ports i, k, and j in sequence. Also, we define \((i^{\prime },k^{\prime }, j^{\prime },v^{\prime }) = \arg \min \limits _{(i,k,j,v) \in \mathcal {P} \times \mathcal {P} \times \mathcal {P} \times \mathcal {V}} \vert \hat {a}_{ikjv} - \lfloor \hat {a}_{ikjv} \rfloor - 0.5 \vert \). If \(\hat {a}_{i^{\prime }k^{\prime }j^{\prime }v^{\prime }}\) is fractional, we create the branching constraints:

$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} B_{i^{\prime}k^{\prime}j's} y_{v'sw} &\leq& \lfloor \hat{a}_{i^{\prime}k^{\prime}j^{\prime}v^{\prime}} \rfloor \end{array} $$
(33)
$$ \begin{array}{@{}rcl@{}} \underset{s \in \mathcal{S}}{\sum} \underset{w \in \mathcal{W}_{v's}}{\sum} B_{i^{\prime}k^{\prime}j's} y_{v'sw} &\geq& \lfloor \hat{a}_{i^{\prime}k^{\prime}j^{\prime}v^{\prime}} \rfloor + 1 \end{array} $$
(34)

and use them to define the new subproblems for vessel class \(v^{\prime }\). The dual value from these branching constraints cannot be directly added to the variables in the master problem. Instead, we introduce an integer variable \(v_{i^{\prime }k^{\prime }j^{\prime }}\) denoting the number of times the ports \(i^{\prime }\), \(k^{\prime }\), and \(j^{\prime }\) are visited in sequence and add it to the objective function with the dual values from constraint (33) or (34). The following constraints can be added to the subproblem to link the v variables with the x variables:

$$ \begin{array}{@{}rcl@{}} v_{i^{\prime}k^{\prime}j^{\prime}} &\leq& x_{i^{\prime}k^{\prime}} + x_{\vert \mathcal{P} \vert + i^{\prime}, \vert \mathcal{P} \vert + k^{\prime}} \end{array} $$
(35)
$$ \begin{array}{@{}rcl@{}} v_{i^{\prime}k^{\prime}j^{\prime}} &\leq& x_{k^{\prime}j^{\prime}} + x_{\vert \mathcal{P} \vert + k^{\prime}, \vert \mathcal{P} \vert + j^{\prime}} \end{array} $$
(36)
$$ \begin{array}{@{}rcl@{}} v_{i^{\prime}k^{\prime}j^{\prime}} &\leq& 2 - x_{k^{\prime}, \vert \mathcal{P} \vert + k^{\prime}} - x_{\vert \mathcal{P} \vert + k^{\prime},k^{\prime}} \end{array} $$
(37)
$$ \begin{array}{@{}rcl@{}} v_{i^{\prime}k^{\prime}j^{\prime}} &\geq& x_{i^{\prime}k^{\prime}} + x_{k^{\prime}, \vert \mathcal{P} \vert + k^{\prime}} + x_{\vert \mathcal{P} \vert + k^{\prime},\vert \mathcal{P} \vert + j^{\prime}} + x_{\vert \mathcal{P} \vert + i^{\prime},\vert \mathcal{P} \vert + k^{\prime}} + x_{\vert \mathcal{P} \vert + k^{\prime},k^{\prime}} + x_{k^{\prime}j^{\prime}} - 2 \end{array} $$
(38)
$$ \begin{array}{@{}rcl@{}} v_{i^{\prime}k^{\prime}j^{\prime}} &\geq& x_{i^{\prime}k^{\prime}} + x_{k^{\prime}j^{\prime}} - 1 \end{array} $$
(39)
$$ \begin{array}{@{}rcl@{}} v_{i^{\prime}k^{\prime}j^{\prime}} &\geq& x_{\vert \mathcal{P} \vert + i^{\prime},\vert \mathcal{P} \vert + k^{\prime}} + x_{\vert \mathcal{P} \vert + k^{\prime},\vert \mathcal{P} \vert + j^{\prime}} - 1 \end{array} $$
(40)

Since the dual variables corresponding to constraints (33) and (34) are non-positive and non-negative, respectively, only constraints (35), (36), and (37) are needed for for the subproblem created by constraints (33), while constraints (38), (39), and (40) are needed for the subproblem created by constraints (34).

4 Computational Study

The proposed branch-and-price algorithm was implemented in Java SE 8, and all mathematical models are solved using Gurobi Optimizer 9.0.1. The computational study was performed on a computer with 2 x 3.5GHz Intel Xeon Gold 6144 CPU 8 core CPUs and 384 GB RAM. The instances on which the algorithm is tested are presented in Section 4.1. Section 4.2 presents results from preliminary testing where different algorithmic choices are tested and evaluated. The results using this algorithm are compared with best known results and this comparison is made in Section 4.3.

4.1 Test Instances

The instances used in this computational study are taken from [2]. These instances are generated to resemble two common network structures: hub-and-spoke networks and feeder networks (see Fig. 5). The hub-and-spoke networks are usually intercontinental networks where most demand is between different continents. In the instances in [2], the hub-and-spoke networks have two regions with ports. Most of the demands are interregional but there are also demands between ports in the same region. The number of ports in each region is the same. General data about the instances are given in Table 2; for more specific information, see [2].

Fig. 5
figure 5

Left: Hub-and-spoke network with two regions, one hub in each region, one interregional service (solid line), three feeder services in the left region (dashed lines), and two feeder services in the right (dashed lines). Right: Feeder network with three services, the black node is the hub

Full size image
Table 2 General data about the instances
Full size table

Each instance is characterized by the structure of the network, the number of ports, and the number of demands. The instances are named S_P_D, where S = F means a feeder network and S = H means a hub-and-spoke network. P is the number of ports and D is the number of demands. Note that in a feeder network, D = 2(P − 1) since all demands either originate or are destined for the hub. This means that the instance H_6_10 is a hub-and-spoke network with six ports and ten demands.

4.2 Preliminary Testing

A subset of six representative instances are used for the preliminary testing. We first test the proposed strengthening of the master problem, constraints (8)–(10). The inner representation of the convex hull of the feasible region defined by capacity constraints over subsets of ports was added to the formulation presented in [2] and the computational study showed that this had a substantial effect on the objective value of the linear relaxation as well as the solution time. All subsets of cardinality up to four, i.e.\(\vert \mathcal {P}^{S}_{u} \vert \leq 4, u \in \mathcal {U}\), are generated and the corresponding constraints are added a priori to the master problem. Table 3 summarizes the results from this first test. We have compared the linear relaxation of four different models; the basic model of [2], Basic, the same model with the inner representation, Best, which was found to give the strongest relaxation in [2], the master problem proposed here, (1)–(6), and the same model strengthened with the inner representation, (1)–(6), (8)–(10). The table shows the relative gap between the linear relaxation and the best solution and is defined as \(G_{LP} = 100 \cdot (\overline {z}_{IP} - z_{LP} ) / (\overline {z}_{IP} - z_{LP}^{BASIC} )\), where zLP is the objective value of the linear relaxation, \(z_{LP}^{BASIC}\) is the objective value of linear relaxation of the basic model in [2], and \(\overline {z}_{IP}\) is the best objective value found in [2].

Table 3 A comparison of the relative gaps for the different linear relaxations
Full size table

First, we note that the linear relaxation of the master problem is stronger than the linear relaxation of the formulation with the inner representation as presented in [2] on the three largest instances. We also see a clear improvement in the objective value of the master problem with the inner representation, the mean remaining gap decreases from 39.1 to 18.2%. Based on this preliminary testing, we decide to continue the computational study with the formulation including the proposed strengthening of the master problem. This means that we use the models (1)–(6), and (8)–(10) unless otherwise stated.

Finding the column with the most negative reduced cost is not necessary in each iteration, instead we only need one column with negative reduced cost for the RMP to produce a new dual solution. Hopefully, the full subproblem only needs to be solved once for every vessel class in every node, to verify that there are no more columns with negative reduced cost and thus that the node is solved to optimality. Table 4 summarizes the results from including no acceleration strategies, No acc., including the improvement of the delivery patterns, Improve, and including both strategies, Both acc., with respect to the number of columns generated, while Table 5 focuses on the time spent.

Table 4 Results of the proposed acceleration strategies with respect to the number of columns generated
Full size table
Table 5 Results of the proposed acceleration strategies with respect to the time spent
Full size table

Table 4 shows the number of columns generated by the different strategies in the three different settings tested. CIP, CLP, and CLS are the numbers of columns found by the full subproblem, by the improve delivery pattern strategy, and by the local search strategy respectively. IIP is the number of times the full subproblem is called. When no acceleration strategies are used, the full subproblem generates on average 1.8 columns per call, while this number decreases to 0.4 when both strategies are used. This is a clear indication that the full subproblem is mainly verifying optimality and that the acceleration strategies find many columns. Ideally, CIP = 0 when acceleration strategies are introduced. We do not reach this, the only exception is F_4_6, but there is a substantial decrease in the number of columns generated by the full subproblem and the number of times the full subproblem is solved.

The effect of this decrease on the solution time is shown in Table 5. Here, TIP, TLP, and TLS are time spent in the full subproblem, in the improve delivery pattern strategy, and in the local search strategy, respectively. T is the total time spent in the algorithm. We see a dramatic decrease in solution time by including the acceleration strategies for the instances that are solved.

Finally, Table 6 summarizes the preliminary testing and shows the gap, G, solution time, T, and number of searched nodes, N, from the testing of the acceleration strategies. We define the relative gap between the lower bound after 1 h and the best solution found in [2] as \(G = 100 \cdot (\overline {z}_{IP} - z_{LB} ) / (\overline {z}_{IP} - z_{LP}^{BASIC} )\), where zLB is the lower bound after 1 h. The lower bound after 1 h and the solution time of the model in [2] with the inner representation is added as a comparison. Note that for instance F_12_22, the root node is not solved within 1 h. No method is clearly better than the other, [2] solves instance F_8_14 while the relative gap on instance H_8_10 is much higher. We see a clear positive effect when comparing the acceleration strategies; using both strategies gives much shorter solution times on the instances that are solved and a much higher number of nodes, and thus a better relative gap, on the instances that are not solved.

Table 6 Summarized results of the proposed acceleration strategies
Full size table

4.3 Comparison and Results

One of the main purposes with the Dantzig-Wolfe reformulation is to improve the dual bound. In Section 4.2, we saw that the lower bounds at the root node produced with the proposed formulation are clearly better than the corresponding bounds presented in [2]. Here, we present a more comprehensive study where all instances are included, and the computational time is set to 1 h. In Table 7, we present a comparison between results from the algorithm proposed here and results by the method proposed in [2]. The table shows the name of the instance, the root node gap calculated as \(G_{LP} = 100 \cdot (\overline {z}_{IP} - z_{LP} ) / \overline {z}_{IP}\), the final gap calculated as \(G = 100 \cdot (\overline {z}_{IP} - z_{LB} ) / \overline {z}_{IP}\), and the computational time for both algorithms. Note that the same upper bound is used for both algorithms. Instances that are solved to proven optimality are marked with a 0 in the G column.

Table 7 Comparison of the gaps between the method presented in [2] and the proposed algorithm
Full size table

We see that the method proposed in [2] solves 12 instances while our method solves eight. On the other hand, the average relative gap at the root node is almost halved, 17.3 compared with 9.1, in favor of our method and the average relative gap after 1 h is also smaller, 6.4 compared with 5.3. This improvement is especially clear on the larger instances, \(\vert \mathcal {P} \vert \geq 10\), where the relative gaps are 22.9 and 17.0 compared with 13.0 and 11.7. Note that instance F_12_22 is run until the root node is solved, the reported result is the root node solution and the time to solve the root node. This has been marked with an asterisk in the table. In [2], the authors report the gap of F_12_22 after 10 h to be 21.6, clearly worse than the root node gap reported here.

5 Conclusions

It is clear that the liner shipping network design problem is a notoriously hard problem. The transshipment possibilities decouple the loading and unloading decisions, and clearly weaken the formulation. Many other routing problems use dynamic programming for solving the pricing problem, but this is not a suitable method here. A reason for this is that you do not have a node with a well-defined state from where to start the algorithm. Often a depot is defined and the fact that the vehicle is either empty when leaving the depot or returning to it is used to start the algorithm. Another reason is transshipment, which creates a large number of loading and unloading possibilities in each node.

We find that capacity cuts defined as an inner representation of the convex hull of the feasible region from [2] significantly strengthen the master problem. The resulting LP relaxation of our model is considerably stronger than that of [2], giving an improved initial dual bound at the cost of computational time needed to solve the root node of the branch-and-bound tree. We see that closing the gap is challenging, and we are not able to reach optimality within 1 h for instances larger than six ports and ten demands.

There is still room for improvements within the approach taken here, and results could likely be improved by examining more strategies in terms of branching, branch-and-bound tree traversal, improved handling of columns in the master problem, and better heuristic methods for finding columns. However, there are limits to what can be done with algorithm engineering, and significant time is spent in the relatively complex full subproblem. This needs to be solved at least once in each fully explored branch-and-bound node to prove optimality, and makes the method impractical for finding optimal solutions to instances of any significant size. It is clear that heuristic methods are needed to find good primal solutions, but in terms of dual bounds we have found an improvement over previously published results on the same instances, especially for larger instances.