Container shipping service selection and cargo routing with transshipment limits

Accepted Manuscript

Container Shipping Service Selection and Cargo Routing with Transshipment Limits Anantaram Balakrishnan , Christian Vad Karsten PII: DOI: Reference:

S0377-2217(17)30459-9 10.1016/j.ejor.2017.05.031 EOR 14459

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

17 February 2016 11 May 2017 13 May 2017

Please cite this article as: Anantaram Balakrishnan , Christian Vad Karsten , Container Shipping Service Selection and Cargo Routing with Transshipment Limits, European Journal of Operational Research (2017), doi: 10.1016/j.ejor.2017.05.031

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Highlights  Paper proposes new model for liner shipping service planning with limited container transshipments.  Limiting the number of times a container is transshipped is important for customers.  Optimization model jointly decides shipping services, demand selection, and cargo routing.  Solution procedure combines problem reduction, model strengthening, and heuristic.  Results for realistic liner shipping problems show effectiveness of solution approach.

1

ACCEPTED MANUSCRIPT

Container Shipping Service Selection and Cargo Routing with Transshipment Limits Anantaram Balakrishnan McCombs School of Business, University of Texas at Austin, Austin, Texas 78712, USA

Christian Vad Karsten

CR IP T

[email protected]

Department of Management Engineering, Technical University of Denmark, Kgs. Lyngby, Denmark [email protected]

AN US

Abstract

We address the tactical planning problem faced by container liner shipping companies to select a set of sailing services from a given pool of candidate services and route available cargo over the chosen services so as to maximize profit. One of the distinctive features of our model is that it incorporates limits on the number of transshipments for each container, a common service requirement in practice. These limits can

M

vary by shipment attributes such as origin and destination, and cargo priority. We propose a new stageindexed multi-commodity flow model that is based on an augmented network containing links (representing sub-paths) between every pair of ports visited by a candidate service. This sub-path

ED

structure, together with our approach of indexing the flow variables by transportation stage, enables the model to accurately capture transshipment costs and enforce transshipment limits.

To reduce the

PT

computational time to solve this problem, we develop preprocessing steps that exploit network structure to eliminate variables, describe valid inequalities to strengthen the model’s linear programming

CE

relaxation, and propose an optimization-based heuristic algorithm to generate good initial solutions. We report successful computational results for realistic problem instances from a benchmark suite of liner

AC

shipping problems, solved using a standard solver applied to our reduced and strengthened model.

Key words: OR in maritime industry, Decision support, Transportation, Networks

2

ACCEPTED MANUSCRIPT

1. Introduction Maritime transportation is a vital component of the modern global trading system. The share of goods transported globally on container ships has increased steadily over the past decades and is expected to grow further due to the economic and environmental advantages of ocean transport compared to other modes of transportation. Since container ships are very expensive to acquire and operate, liner shipping

CR IP T

companies need to utilize their assets effectively by judiciously choosing their sailing routes and deciding which demands to meet in order to maximize profit while ensuring adequate service to customers. This paper develops and solves an optimization model to support tactical decisions regarding which set of services to operate and how to route container flows on the chosen services, taking into account the demand patterns, associated revenues, available ships and their capacities, and operating costs. An important feature of our model is its ability to limit the number of transshipments for each container, a

AN US

common practice among liner shipping companies to assure good service to customers (Brouer et al. 2014a). The problem we address, which we call the Liner Service Planning (LSP) problem, is defined as follows. We are given the anticipated demand between various ports that the liner shipping company serves, and a candidate pool of sailing services. Each service is defined over a cyclic route spanning multiple ports, with an associated assignment of ships and service frequency. The LSP problem entails

M

selecting a subset of services from the candidate pool, deciding how much of each demand to serve, and routing the selected demand, subject to transshipment limits and vessel capacities, so as to maximize

ED

profit.

In this paper, we develop a new model for liner service planning, discuss modeling and methodological

PT

enhancements to effectively solve this large-scale optimization problem, and computationally validate the model and solution approach. We propose a novel multi-commodity model based on flows along sub-

CE

paths, indexed by transportation stage, to capture the transshipment constraints and costs for each commodity. The model can readily incorporate practical container routing issues such as cabotage rules, regional policies, and embargoes. For this model, we outline a problem reduction procedure that exploits

AC

the network structure and routing requirements to eliminate or combine some variables. We then describe valid inequalities to increase the linear programming lower bounds, and develop an optimization-based heuristic procedure to generate good initial solutions. By reducing the size of the formulation, raising the lower bounds, and generating good initial upper bounds, these techniques reduce the computational time needed to solve the problem optimally (or near-optimally, with performance guarantees) using solvers that apply branch-and-bound. To demonstrate the effectiveness of our solution procedure, we present computational results for realistic problem instances based on the Liner-Lib benchmark problems (Brouer

3

ACCEPTED MANUSCRIPT

et al. 2014a).

Next, we provide some background on maritime container transport and related

optimization models to motivate and position our work. The International Maritime Organization estimates that 90% of global trade is carried by sea, with container ships transporting around 60 percent of the seaborne goods by value (World Shipping Council 2009).

Modern cost and energy efficient container vessels can carry almost 20,000 twenty-foot

equivalent units (TEU) of cargo. Since these ships cost more than 100 million dollars per vessel,

CR IP T

operating a global ocean transport network requires enormous capital investments. With such large investments, it is necessary to ensure high utilization of ships to be able to offer low shipping cost and remain competitive. Global liner shipping networks divide their geographical coverage regions into major trade lanes that follow the North-South, East-West, and intra-regional trade patterns of the world. Within each trade lane, a carrier may operate multiple services, each consisting of a cyclic sailing route

AN US

that visits a given subset of ports. A service consists of several vessels of similar size that are scheduled to ply the route at regular intervals (e.g., one week apart) so that every port on the route receives periodic (e.g., weekly) service. Brouer et al. (2014a) provide a broad introduction to the liner shipping context and operations. An effective shipping network will exploit transshipment opportunities at intermediate ports to serve demand between many different origin-destination pairs while maintaining high fleet utilization. At the same time, to meet customer requirements, liner shipping companies must also limit the number of

M

transshipments of each container. Generally, customers prefer to have no more than two or three transshipments per container in order to reduce the risk of damage or loss, missed connections, and long Some types of cargo (e.g., hazardous goods or high value items) may need to be

ED

layover times.

transported without any transshipments whereas others may require less stringent limits on the number of transshipments. Therefore, the service planning model must have the ability to impose demand-specific

CE

and regulations.

PT

transshipment limits. Container routing decisions may also be constrained by other operational policies

For the liner shipping context, researchers have proposed various optimization models to support decisions at the strategic, tactical, and operational levels. The long-term strategic issues relate to the

AC

acquisition and deployment of highly capital intensive assets, including decisions regarding the markets to be served, fleet size and mix, and the general configuration of the network. At the lowest level, shortterm operational decisions focus on day-to-day issues such as loading containers on a ship, responding to disruptions, and repositioning empty containers. This paper focuses on medium-term tactical decisions, namely, which services to operate, what available vessels to assign to these services, which demands to meet, and how to route these demands on the network. These decisions are reviewed periodically to

4

ACCEPTED MANUSCRIPT

determine whether to add new services or modify existing services in response to changing demand patterns. Traditionally, liner shipping companies have relied on experienced planners to manually select candidate routes and associated services based on the company’s service strategy. However, choosing the services for one trade lane or region at a time, as is often the case with manual planning, causes ripple effects

CR IP T

throughout the network, leading to suboptimal overall performance in terms of profitability and resource usage. Our model seeks to optimize these service selection decisions from a network-wide perspective, taking into account the potential demand between various ports and the available ships. Instead of designing routes and services from scratch, as many strategic network design models in the literature do, our model takes as input the candidate set of routes and services that planners have identified, and focuses on selecting the best subset of these services for the company to operate so as to maximize profits from

AN US

container transport. This approach of service selection (versus service design) at the tactical level permits incorporating cargo routing considerations such as transshipment limits and is also appealing to practitioners because it permits planners to specify which services to consider. These candidate services may include those that the company has previously operated as well as potential new services driven by

M

strategic considerations or obtained by applying a network design model to identify promising routes. This paper’s contributions to the literature span modeling, methodology, and computational validation.

ED

We present a new mathematical formulation for service selection and demand routing with transshipment limits, a problem that has not been adequately addressed in the liner shipping optimization literature. By representing the decision variables as flows along sub-paths (rather than flows on the arcs of the physical

PT

network, i.e., sailing edges, or on end-to-end paths) and by indexing these flows by transshipment stage, our model accounts for container transshipments and costs, and requires only polynomial number (in

CE

terms of network size, number of commodities and candidate services, and maximum number of transshipments) variables and constraints. So, the model does not require specialized solution methods such as column generation. To effectively solve this NP-hard problem, we propose problem reduction

AC

techniques, model strengthening, and an optimization-based heuristic approach to generate good initial upper bounds. Computational results using our reduced and strengthened model (solved by a standard solver) for realistic problem instances confirm that our enhancements significantly reduce computational time and yield good solutions.

The rest of this paper is organized as follows: Section 2 reviews the literature on related models in liner shipping.

Section 3 formally defines the LSP problem, introduces the model, discusses problem

5

ACCEPTED MANUSCRIPT

complexity, and outlines additional features that the model can incorporate and the merits of our modeling approach. Section 4 presents our modeling and methodological enhancements. Section 5 reports on our computational results, and Section 6 offers concluding remarks.

2. Literature Review Christiansen, Fagerholt, and Ronen (2004) and Christiansen et al. (2013) provide comprehensive reviews

CR IP T

of papers in the area of maritime transportation over past few decades. For liner shipping, researchers have addressed problems at the three planning levels – strategic, tactical, and operational. Agarwal and Ergun (2008) and Meng et al. (2014) discuss hierarchical planning frameworks for liner shipping, and Kjeldsen (2011) develops a classification scheme for routing and scheduling problems in liner shipping. Since the literature on optimizing liner shipping operations is extensive and several recent papers review this literature, we focus our review on selected papers related to our work, namely those that consider

AN US

network design, fleet deployment, and cargo routing, with emphasis on those that consider transshipment costs or constraints. As our review indicates, most models do not consider limits on the number of transshipments and, unlike our work, most papers do not focus on obtaining solutions that are optimal or provably close to optimal.

M

Network design refers to decisions regarding the configuration of cyclic routes and deployment of available ships for a given service frequency (e.g., weekly) taking into account container transport demand patterns. Because the problem entails developing sailing routes from scratch rather than selecting

ED

services from a candidate list, it is quite complex even without additional cargo routing restrictions such as limited transshipments. Meng et al. (2014) review the literature on optimization models and methods

PT

for network design, and conclude that no approach has been successful in designing networks of realistic size while also accounting for industry constraints, service requirements, and routing decisions. Within the broad category of network design models, researchers have studied many variants and special cases

CE

such as designing a single route (e.g., Plum et al. 2014), modeling networks with special structure (e.g., hub and feeder network, Meng and Wang 2011), and transporting cargo without transshipment (Wang and

AC

Meng, 2014). For general networks, Agarwal and Ergun (2008) formulate a network design problem with weekly service frequency which they call the simultaneous ship scheduling and cargo routing model. The model does not consider transshipment costs (or limits) since, as the paper notes, it is difficult to account for such costs if the network is not known. The authors propose and test three solution methods – a greedy heuristic, a column generation algorithm, and methods based on Benders decomposition. Computational results are reported for randomly generated networks with up to 20 ports, but sparse traffic patterns.

Mulder and Dekker (2014) develop a heuristic method for network design (without

transshipment costs) based on a genetic algorithm, and report results for a Europe-Asia network. 6

ACCEPTED MANUSCRIPT

Recently, researchers have developed models that do account for transshipment costs; however, most papers either employ heuristic methods or are only able to solve small problem instances to optimality. Álvarez (2009) considers the problem of joint routing and deployment of a fleet of container vessels on cyclic routes, and incorporates cargo transshipment costs (but, the model does not capture these costs accurately when a route visits the same port more than once). The author solves the model using tabu

CR IP T

search and a standard MIP solver. Brouer et al. (2014a) extend the model of Álvarez (2009) to correctly account for transshipment costs for butterfly networks, and introduce a service frequency restriction. A heuristic route generation method is able to produce service networks of reasonable size but with varying quality. Reinhardt and Pisinger (2012) propose an exact branch-and-cut algorithm for network design, without frequency restrictions, but their method can only solve fairly small problems. Plum, Pisinger, and

AN US

Sigurd (2014) develop a compact formulation of the network design problem based on service flows and incorporating transshipment costs. The authors report computational results using CPLEX for two small Liner-Lib instances (Brouer et al. 2014a), Baltic and West Africa (WAF), neither of which could be solved to optimality. For the model presented in Brouer et al. (2014a), Brouer, Desaulniers, and Pisinger (2014b) develop a matheuristic to perform incremental network optimization, i.e., improve the shipping routes taking into account the cargo flows. Although these papers consider transshipment costs, they do

M

not incorporate any restrictions, such as transshipment limits, on the container routes.

ED

The literature on models that combine network design or service selection with transshipment limits or other cargo routing restrictions is very limited. Meng and Wang (2011) consider the problem of selecting services over a set of hub and feeder ports, with each feeder port assigned to exactly one hub port and no

PT

transshipment permitted when transporting cargo between hubs. Hence, the possible container paths are limited and can be readily enumerated. The authors describe a mixed-integer programming model to

CE

minimize cost to meet all demand, and present computational results for an Asia–Europe–Oceania network containing six hubs at which transshipments occur. Note that due to the special structure of this hub and feeder network and the policy of not permitting transshipments for inter-hub movements, cargo

AC

between any origin-destination pair is transshipped at most twice, and so explicit transshipment limits are not needed. At lower levels in the planning hierarchy, Meng and Wang (2012) consider fleet deployment and time-constrained container routing for a given set of sailing routes. Wang et al. (2012) propose and solve an integer program for the operational problem of selecting container paths for a single OD pair, subject to transit time and cabotage restrictions. Similarly, Song and Dong (2012) solve the operational problem of joint cargo routing and empty container repositioning, accounting for demurrage and inventory cost of empty containers. Both these papers address cost minimization problems assuming that 7

ACCEPTED MANUSCRIPT

all demand must be satisfied. Wang, Meng, and Lee (2016) address the problem of deciding what portion of transit-time-sensitive demand to accept and load on a given set of services, show that this problem is polynomially solvable, and report computational results using actual data. Recently, Karsten et al. (2016) developed a column generation procedure to solve a time-constrained multicommodity flow problem that can be used to route available cargo over a given set of services. This corresponding pricing subproblem requires solving a resource constrained shortest path problem (that can also be adapted to handle

network design, and report computational results.

CR IP T

transshipment limits). They embed this procedure in Brouer, Desaulniers, and Pisinger’s matheuristic for

We conclude the literature review by noting that, although most prior models for liner shipping do not consider transshipment limits, the broader network optimization literature has considered analogous

AN US

constraints for some general classes of problems. Balakrishnan and Altinkemer (1992) were among the first researchers to propose and model hop constraints (i.e., limits on the number of hops or arcs on which a commodity is routed) in the context of fixed-charge network design. Since then, other authors have studied models with hop constraints for special cases such as shortest paths, spanning trees, and Steiner trees (e.g., Dahl and Gouveia 2004, Dahl et al. 2006, and Gouveia and Magnanti 2003). These models are not directly applicable to liner shipping since they do not incorporate features such as service selection,

M

vessel deployment, and sailing edge capacities.

ED

3. LSP Problem: Definition and Model Formulation In this section, we formally define the LSP problem, present a multi-commodity flow formulation defined

PT

over an augmented network whose links represent sub-paths, discuss the problem’s computational complexity, outline additional features that the model can accommodate, and contrast our modeling

CE

approach with conventional arc flow and path flow models. 3.1 Problem statement and network representation

AC

Given a pool of candidate services, the number of vessels available in each class, and the estimated demand (containers) between various ports, the LSP problem entails selecting the best subset of services from the candidate pool, deciding how much of each origin-to-destination demand to serve, and routing these demands on the chosen services. In this context, a service (also known as rotation) refers to a specified sailing route, consisting of a sequence of sailing legs, with an associated number of ships of a particular class deployed on that route, service frequency, and ship speed. The routing decisions must satisfy limits on the number of transshipments for each demand, and the capacity on each sailing leg of every chosen service. The goal of the LSP model is to maximize profit, obtained by subtracting the costs 8

ACCEPTED MANUSCRIPT

for operating ships and container loading/unloading/transshipment from the revenue obtained by meeting customer demands. (Equivalently, by assigning an opportunity cost or penalty per unit of unmet demand, we can frame the problem as a cost minimization model; see discussion in Section 3.2.) The liner shipping transportation network consists of a set of ports N, and a set of candidate services R, with sailing edges e  Er in each service r  R; a sailing edge represents the portion of a ship’s itinerary

CR IP T

between two successive ports of call on the service route. Let V denote the set of vessel classes or types needed for the candidate services. The firm has M v available vessels of each class v  V.

For each

service r, we are given the required number mrv of vessels of class v  V, the capacity ge (in terms of number of containers) of each sailing edge e  Er of that service, and the fixed cost fr for selecting the service. This cost includes the amortized cost of the vessels assigned to the service, fuel and other

AN US

operating costs for sailing and idling, canal costs, and port call costs. In practice, the available capacity of a sailing leg may differ from the capacity of the vessels assigned to the corresponding service if, for instance, there is a vessel sharing agreement (VSA) for the service. Hence, instead of assuming the same capacity on each sailing edge of a service, we permit the capacity g e to vary by edge. As inputs to the model, we are given the demand, i.e., container traffic available to be transported, between various origin-

M

destination ports. We associate a commodity k  K with each such demand. This commodity originates at port sk and needs to be transported to destination port tk . If needed, we can further distinguish

ED

between containers having the same origin and destination but different characteristics, e.g., different limits on number of transshipments, by defining multiple commodities with different attributes between the same origin-destination pair. Let d k denote the number of containers of commodity k available for

PT

transport. We permit splitting the flows of a commodity, i.e., we can route containers of a commodity along multiple paths from the corresponding origin to destination as long as each path satisfies the limits

CE

on the number of transshipments. The firm is not required to transport all dk units of demand for commodity k, but incurs a cost or penalty of qk per container of unmet demand. This cost may represent,

AC

for instance, the opportunity loss (e.g., lost revenue, loss of goodwill) for not transporting the container. One of the distinctive and key features of our model is the way we define the underlying network for multi-commodity flows. Specifically, to incorporate the limits on the number of transshipments and capture their costs, we introduce an augmented multi-commodity flow network based on sub-paths. We define a sub-path of a commodity’s route as the portion of this route in which the container travels on a single service. So, if the container is loaded onto a ship for service r at port i and unloaded from the ship at port j, we denote the segment from i to j on this service as the sub-path . We include the service 9

ACCEPTED MANUSCRIPT

index r in the notation for a sub-path since there may be several different services that can transport the container from i to j. Note that, if a service r makes nr port calls, we can have as many as nr (nr  1) sub-paths associated with that service (see Figure 1). But, as we discuss later, operational

policies may preclude using some of these sub-paths. Our approach of modeling container movements as flows over sub-paths lies between two conventional modeling strategies – as flows over the arcs of the

CR IP T

physical network, i.e., over sailing edges, or flows over origin-to-destination paths. Later (in Section 3.6) we discuss the significant advantages of using our intermediate approach compared to the arc-flow and path-flow models.

Our augmented network contains one node for each port and one link for every sub-path of each service. Figure 1 shows an original service, say service r, covering five ports, and the augmented network with

AN US

links corresponding to the sub-paths of this service (in this and subsequent figures, the sailing direction is counterclockwise). Each (directed) sub-path consists of a sequence of underlying sailing edges of the service. For instance, the sub-path from i4 to i1 on this service r contains the sailing edges (i4,

CE

PT

ED

M

i5) and (i5, i1).

AC

Figure 1: Sub-paths corresponding to a service (counterclockwise) that visits five ports Figure 2 shows the augmented network for a liner shipping system with three services, each of which visits five ports and shares one port in common with another service. Containers transported on one service can be transshipped to another service at a common port. (In general, services may have more than one common port, creating many alternate paths for container flows.) As in Figure 1, each service induces a complete subgraph (shown with light, dashed, or dotted lines in Figure 2) of the augmented network, containing one link connecting every pair of ports visited by the service. Together, the services 10

ACCEPTED MANUSCRIPT

in the network span 13 ports. The links shown in bold in Figure 2 represent one possible sequence of services to transport a container from its origin a to its destination d. Under this routing plan, the container uses portions (sub-paths) of all three services, and is transshipped twice, once at each

M

AN US

CR IP T

intermediate ports b and c.

ED

Figure 2: Augmented network for example with 3 services spanning 13 ports In our model formulation, the flow of a commodity k on a sub-path will represent the containers

PT

of commodity k that are loaded on service r at port i (which may be the origin or an intermediate port at which the commodity is transshipped from another inbound service) and unloaded (and possibly

CE

transshipped) at port j. The unit cost (per container) of this sub-path includes the variable costs of the sailing edges embedded in this sub-path as well as the costs for container handling (loading, unloading) at ports i and j. We next elaborate on these latter costs. The firm incurs costs for loading a container at its

AC

origin, and unloading it at the destination.

At intermediate ports where transshipments occur, the

container must first be unloaded from one ship and then loaded onto another. These container handling costs vary by port. In practice, the cost of a transshipment (unloading and re-loading) is lower (by around 15%, based on discussions with planners at a liner shipping company) than the sum of costs for one unload and one load operation (some ports have a larger cost difference).

Our model, therefore,

distinguishes between the cost of transshipping from one service to another at an intermediate port from the cost of loading/unloading at the origin/destination of the commodity by making the sub-path costs

11

ACCEPTED MANUSCRIPT

commodity-dependent as follows. Let  skk , tkk , and  kj respectively denote the cost of loading a container of commodity k at its origin port sk , unloading the container at its destination tk , and transshipping the container at an intermediate port j. Then, if cijrk denotes the cost for routing one container of commodity k on sub-path , we set this cost equal to the sum of the following cost components: (i) the total cost to carry (transport) the container on all the sailing edges e  Er of service r that belong to the sub-path;

CR IP T

(ii) the transshipment cost  kj at port j if this port is an intermediate port for commodity k, or the destination unloading cost  tkk if j  tk ; and, (iii) the origin loading cost  skk if i  sk .

Each sub-path corresponds to a transportation stage, which we define as the portion of a container’s itinerary (from port i to port j) in which a container is carried on one service r. Thus, as

AN US

illustrated in Figure 2, a container’s route consists of a sequence of transportation stages (henceforth abbreviated as “stages”) from origin to destination, with transshipment between stages. The number of transshipments for the container is one less than the number of sub-paths on which it travels. Hence, we can enforce the transshipment limits by limiting the number of stages on a container’s route. Let hk denote the maximum permitted number of stages for commodity k. Define Ar as the set of sub-paths of

M

service r; this set contains sub-paths for every pair of ports i and j that service r visits. Later, we will discuss how to prune this set by excluding certain sub-paths for each commodity based on problem

ED

and network structure, operating policies, or regulations, thereby reducing the number of decision variables in the problem formulation. Table 1 summarizes the notation for the sets and parameters in our model.

CE

PT

Sets and indices N Set of ports (nodes) in the shipping network; i, j  N K Set of commodities; k  K R Set of candidate services; r  R sk , tk Origin and destination for commodity k  K

SEijr

Set of sailing edges in service r; e  E r Set of sub-paths on service r Set of sailing edges of service r in sub-path

V

Set of vessel classes; v  V

AC

Er Ar

Parameters Fixed cost for selecting service r fr Cost per container if commodity k is routed on sub-path cijrk dk hk

Number of containers available to be transported (demand) for commodity k Maximum permitted number of transportation stages for commodity k 12

ACCEPTED MANUSCRIPT

nr v r

m

Mv qk

ge

Number of port calls in service r Number of vessels of class v needed for service r Number of available vessels of class v Penalty for not meeting one unit (container) of demand for commodity k Capacity (number of containers) of sailing edge e of service r

CR IP T

The use of sub-paths as the basis for our augmented network is new, although it has some similarity to the segments introduced by Meng and Wang (2011) and the links defined by Bell et al. (2011). In these papers, the segments and links are simpler than general sub-paths or not service-specific, whereas in our model each sub-path has an associated service and can have an arbitrary number of embedded sailing edges, providing more modeling flexibility. Note that our sub-path structure is not limited to simple

AN US

routes that visit each port exactly once; it also extends to complex routes such as a butterfly or conveyor belt services (Brouer et al. 2014a) that may visit a port more than once. We refer to such ports as multivisit ports. Suppose a sub-path goes through a multi-visit port l more than once. The solution can choose between of assigning containers to this sub-path (without intermediate transshipments), or use two shorter sub-paths, one from i to l and the other from l to j; in the latter case, containers must be unloaded and then re-loaded (on to the same service) at port l. The model’s decision on whether to use

M

the single longer sub-path or the two shorter sub-paths depends on the cost tradeoffs (i.e., cost of additional sailing edges in the longer sub-path versus additional transshipment cost for the two sub-paths)

3.2 Model formulation

ED

and the capacity usage on the sailing edges of the longer sub-path.

PT

In our multi-commodity model defined over the augmented network, instead of adding explicit constraints on the number of transshipments, we define flow variables for each possible transportation stage of every

CE

commodity. We index these stages consecutively, starting with stage 1 at the container’s origin. In general, stage h represents the hth sub-path or service on a container’s route after it leaves the origin. In this scheme, h = 1 corresponds to the container’s first sub-path when it leaves the origin. If the port, say, j

AC

at which this sub-path ends is an intermediate port for the container, then h = 2 corresponds to the second sub-path from port j on which the container travels, and so on. For each commodity k, we index its subpath flow decision variables by transportation stage h, for h = 1, 2, …, hk (recall that hk is the maximum number of stages on which the commodity can routed). These stage-indexed flow variables are then linked across stages using appropriate flow conservation constraints. By defining flow variables for only hk stages, any origin-to-destination path can contain no more than hk stages, as required by the

transshipment limits. This approach also has the considerable advantage of permitting the solution to split 13

ACCEPTED MANUSCRIPT

flows of a commodity across multiple origin-to-destination paths while ensuring that each of these paths satisfies the transshipment limits. We next formally define the decision variables, and formulate the LSP problem. Our LSP model uses three sets of decision variables representing respectively the decisions regarding which services to select, how to route each commodity, and how much demand of each commodity to

CR IP T

forgo. We define these decision variables as follows: xr = 1 if service r is used, and 0 otherwise, for all r  R;

uijrhk = flow of commodity k using sub-path of service r as the hth stage, for r  R,  i, j, r  Ar ,

and h = 1, 2, …, hk ; and,

zk = amount of demand (number of containers) not met for commodity k  K.

AN US

We can then formulate the LSP problem as the following mixed-integer program, denoted as [LSP]: hk

[LSP]

 fr xr   cijrk uijrhk   qk zk

Minimize

rR

kK rR h1

 

u1skk jr  zk  d k

 

uijrhk  

rR  sk , j ,r Ar

rR i:i , j ,r Ar

hk

rR l: j ,l ,r Ar



kK h1 i , j ,r Ar :eSEijr

m x

uijrhk  ge xr

 Mv

PT

rR

v r r

uijrhk  0

u hjlr1,k  0

ED





M

subject to

(1)

kK

k  K ,

(2)

k  K , j  N \{sk , tk }, h  1,..., hk  1,

(3)

r  R, e  E r ,

(4)

v V ,

(5)

k  K , r  R,  i, j, r  Ar , h  1,..., hk , (6) k  K , and

(7)

xr {0,1}

r  R .

(8)

AC

CE

zk  0

The objective function (1) minimizes the total cost, consisting of the fixed costs for the selected services, the cost of transporting goods on each sub-path, and the penalties for unmet demand. Observe that, instead of formulating the problem as a profit maximization model, we have framed it as a cost minimization problem by considering the penalty or opportunity cost for unmet demand. Alternatively, if qk represents the revenue per unit for transporting commodity k on the firm’s network, we can express

the objective function as:

14

ACCEPTED MANUSCRIPT

Maximize

hk  k hk  q ( d  z )  f x     r r  cijr uijr  . k k k kK kK rR h1  rR 

This objective is equivalent to the minimization objective (1). We use the minimization form in our implementation and when reporting the computational results.

CR IP T

Constraints (2) assign the flow of each commodity k to the sub-paths  sk , j, r  incident from the origin port sk for this commodity (in stage h = 1), and specify that the total flow on these sub-paths together with the unmet demand (variable zk ) must equal the commodity’s demand. Constraints (3) are the flow balance constraints at each intermediate port j  sk , tk requiring that the total flow of a commodity entering the port via sub-paths corresponding to stage h must be equal the total flow leaving that port on

AN US

sub-paths for stage (h + 1), for h = 1, 2, …, hk  1 . These constraints, together with constraints (2), ensure that all the flows of commodity k (on one or more paths from the origin) enter the destination node within hk or fewer stages, as required, and the total flow entering the destination is the total flow leaving the origin node, excluding the unmet demand. Constraints (4) serve to both impose the capacity of each sailing edge and to ensure that we assign flows to a sub-path only if the corresponding service r is selected. The left-hand side of this constraint is the total flow on all the sub-paths that contain sailing

M

edge e. The right-hand side specifies that flow can be positive only if the corresponding service r is chosen (i.e., xr  1 ) and must not exceed the capacity of edge e. Constraints (5) ensure that the total

ED

number of vessels of each class needed to operate the chosen services does not exceed the liner company’s available fleet of that vessel class. Constraints (6) to (8) are the non-negativity and binary

PT

constraints.

CE

If nˆ (| N |) denotes the maximum number of ports covered by a service and hˆ is the maximum number of transshipments permitted across all commodities, then formulation [LSP] has O(nˆ 2 | R || K | hˆ)

AC

variables, of which |R| variables are binary, and O(| K || N | hˆ  nˆ | R |  | V |) constraints. Note that, if sub-path costs vary by transportation stage, we can capture these costs by replacing the coefficient cijrk in the objective function with a stage-dependent unit cost cijrhk . The model can also accommodate other considerations using additional side constraints.

For instance, suppose the firm

wants to ensure that the transshipment cost per container of commodity k (averaged over all containers of that commodity that are transported by the company) does not exceed a pre-specified value, say, k. To satisfy this requirement, we add the constraint:

15

ACCEPTED MANUSCRIPT

hk





rR h1 i , j ,r Ar , j tk



 u1skk jr  ,  rR sk , j ,r Ar , j tk 

 kj uijrhk   k  



where  kj is the cost of transshipping a container of commodity k at intermediate port j.

These

observations illustrate the versatility of our LSP model formulation.

CR IP T

3.3 Problem complexity Agarwal and Ergun (2008) and Brouer et al. (2014a) have shown that the liner shipping network design problem is NP-hard, as are some special cases such as the problem in which the candidate rotations are given but the chosen rotations must visit all the ports. In the following proposition we show that the LSP problem is also NP-hard.

AN US

Proposition: The LSP problem is NP-hard.

Proof: To establish this result, we transform the knapsack problem to a special case of the LSP problem in polynomial time. Given n items, i = 1, 2, …, n, each with a value  i (or profit) and weight i , and a knapsack with a given capacity b, the knapsack problem requires selecting a subset of items with maximum total value that will fit within the knapsack. Given any instance of the knapsack problem, we

M

create an equivalent LSP problem instance having n commodities, one corresponding to each item, defined over a network with (n + 1) nodes (ports), one for each item i together with an additional node, say, node 0 that serves as the common destination for all commodities. There are b available vessels, all

ED

of the same type, and n potential services to choose from. For i = 1, 2, …, n, service i plies between node i and node 0, and requires i vessels. Commodity i has port i as its origin and port 0 as its destination.

PT

Its demand d i can be fully met using the i vessels of service i, and each unit of this demand yields a profit of  i / di . Observe that, by construction, if service i is not selected then none of the demand for

CE

commodity i can be satisfied; conversely, if this service i is selected, then we can transported all available units of commodity i. The optimal solution to the equivalent LSP problem selects a subset of services (items), assigns the needed vessels from the available fleet to each of these services, and serves all of the

AC

demand for the corresponding commodities. Hence, if commodity i is served in this solution, the corresponding knapsack solution includes item i in the subset of chosen items. Since the knapsack problem is NP-hard (Nemhauser and Wolsey 1988), so is the LSP problem.



3.4 Pruning the set of sub-path flow variables We can reduce the number of flow variables uijrhk in model [LSP] by eliminating, a priori, flows of commodity k on certain sub-paths for specific stages h. To capture such reductions, we define the 16

ACCEPTED MANUSCRIPT

set Brhk  Ar as the set of sub-paths of service r that are eligible as the hth transportation stage for commodity k, for h = 1, 2, …, hk .

Hence, instead of defining variables uijrhk for all sub-paths

 i, j, r  Ar , as we did when introducing formulation [LSP], we only need to define these variables for

the sub-paths in Brhk for each k  K, r  R, and h = 1, 2, …, hk . (Correspondingly, we replace the summation over  i, j, r  Ar in the left-hand sides of constraints (2), (3), and (4) with the summation

CR IP T

over the smaller set  i, j, r  Brhk for the appropriate values of h.) We next discuss some natural ways to reduce the number of flow variables based on the properties of commodity routes in the optimal solution. Later, we discuss how trade and operational policies permit further reduction in the sets Brhk (in Section 3.5), and propose additional strategies for variable elimination based on the topology of the

AN US

augmented network (in Section 4.1).

Since sub-path costs are positive, in any optimal LSP solution, the route for every container must be a simple path, from its origin to destination, in the augmented network, i.e., this path will not visit the same transshipment port more than once. Consequently, we can eliminate some sub-paths from the sets Brhk in the following ways:

For each commodity k, service r, and stage h, we do not include in Brhk any sub-path that

M



ends at the origin (i.e., with j  sk ) or starts at the destination (i.e., with i  tk ). Sub-paths that start at the origin sk for commodity k can only be used at stage h = 1. So, for h =

ED



2, …, hk , we can eliminate from Brhk any sub-path with i  sk . Sub-paths for the highest permitted stage hk must end at the destination tk of commodity k. So,

PT



for h  hk , we can eliminate from Brhk k any sub-path with j  tk .

CE

Applying these rules to define the sets Brhk , and including in the LSP model only the flow variables uijrhk corresponding to sub-paths in Brhk reduces the size of the model but does not change its linear

AC

programming lower bound. Next, we discuss some practical issues – regulations and operational policies – that our model can readily accommodate.

3.5 Incorporating trade and operational policies Due to the trade policies of various countries, liner shipping companies must follow certain cabotage rules such as restrictions on the transshipment and internal transport of goods within certain countries. There are also rules for hazardous goods and various embargoes between countries. In the standard

17

ACCEPTED MANUSCRIPT

network flow formulation for cargo routing, one possible way to account for cabotage rules is to add them as resource constraints in the shortest path calculations when solving the routing sub-problems. For our model, on the other hand, we can simply remove sub-paths that are not permitted for a commodity during a preprocessing step. Consider, for instance, a cabotage rule (applicable in some countries) that prohibits transshipment within the country for any container that originates in or is destined to a port in the country, and permits no more than one transshipment in the country for other shipments. Figure 3 shows an

CR IP T

example with five ports, two of which are in one country (ports i1 and i5 , shaded in the figure) and the other three in another country. Then, if the above cabotage rule applies in both countries, we can eliminate the corresponding sub-paths for such commodities, as shown in the right-hand side of Figure 3. For our model, we delete any sub-path of a service r that violates cabotage rules from the set Brhk of feasible sub-paths for commodity k for every stage h, and omit the corresponding sub-path flow variables

AN US

from the model. So, accounting for such rules reduces the model size rather than complicating it, while

PT

ED

M

also making the model more realistic.

Figure 3: Eliminating sub-paths that violate maritime cabotage rules

CE

Next, we discuss how to adapt the model to meet certain operational policies of the liner shipping firm. Networks that are optimized to minimize cost and increase network utilization often assign containers to detours that have unused capacity (Karsten et al. 2015). In practice, there are some sub-paths that a

AC

planner would rule out a priori. For instance, the planner may specify that a container that must be moved on a service from port i to port j in the same region (e.g., the Americas), must not be assigned to a service that visits another distant region (e.g., Asia) enroute from i to j. This rule ensures that containers do not follow circuitous routes. To incorporate this regional policy, we simply remove from the sets Brhk (and omit corresponding flow variables) any sub-paths that start and end in the same region but include an intermediate port in another region.

18

ACCEPTED MANUSCRIPT

3.6 Alternative modeling approaches We conclude the discussion of LSP problem formulation by outlining why our modeling approach, using flows on sub-paths, to formulate the LSP problem is superior to models based on either “arc” flows or “path” flows. (See Brouer, Pisinger, and Spoorendonk (2011) for a discussion of arc flow and path flow models for a problem that involves allocating and routing cargo on a given service network.)

CR IP T

A natural way to model movement of goods (e.g., containers) is to represent them as flows on the arcs of the physical (service) network, i.e., flows on the individual sailing edges of each service. We will refer to this approach as an “arc flow” model. (Note our formulation is different from an arc flow model because its stage-indexed flow variables are defined over links of the augmented network which represent subpaths, and not the arcs of the physical network.) Since the total number of sailing edges is smaller than

AN US

the number of sub-paths in our augmented network, the arc flow model would require fewer flow variables than sub-path-based formulation [LSP]. However, imposing the transshipment limits in the arc flow model is quite difficult, and requires many additional binary variables and constraints, particularly if flow splitting is permitted. Specifically, to capture these limits, the arc flow model needs additional transshipment indicator (binary) variables for each commodity, at every transshipment port, and for every pair of services that visit this port, as well as binary variables to trace which sailing edges are used to

M

transport a commodity. The model also requires numerous additional constraints including forcing constraints that relate the binary variables to the flow values. The linear programming relaxation of this

ED

model can be quite weak particularly because of the “big-M” coefficients in the forcing constraints. Since our formulation [LSP] does not require these additional variables and constraints, it is better suited for

PT

solving the problem using exact methods such as branch-and-bound.

At the other extreme, we can formulate the LSP problem using path flow variables. In its complete form,

CE

this formulation requires one flow variable for each commodity and every feasible origin-to-destination path (satisfying the limit on the number of transshipments) for that commodity. The number of origin-to-

AC

destination paths increases exponentially with the number of permitted transshipments. So, a common strategy is to solve such models using a column generation approach that iteratively identifies promising paths by solving a pricing subproblem, given the solution (dual prices) to the linear programming relaxation of the current master problem. For the LSP problem, this subproblem is itself NP-hard if the maximum number of transshipments is not fixed. Moreover, since the LSP problem also requires selecting the set of services from the candidate set, the master problem needs to include forcing constraints to ensure that, whenever one or more units of a commodity are routed on a newly generated path, all the embedded services (or sub-paths) are selected in the master problem’s solution. Hence, the 19

ACCEPTED MANUSCRIPT

number of constraints in the master problem also grows as new columns are added. Finally, due to integrality of the service selection variables, we need to apply a branch-and-price approach to find the optimal solution. These requirements make the path flow model unappealing from the perspective of effectively solving the problem.

In contrast, since our sub-path based LSP model has a polynomial number of variables, we can apply

CR IP T

standard solvers, and do not require customized decomposition algorithms (such as branch-and-price with associated algorithms to solve the subproblems). Moreover, we expect its linear programming relaxation to be tight. Thus, our sub-path-based model has several advantages over the two extreme alternatives of arc flow and path flow models. It achieves a good tradeoff between model size, computational time, and implementation difficulty, as validated by our successful computational experience (Section 5).

AN US

4. Modeling and Methodological Enhancements

This section discusses three types of enhancements – problem reduction, model strengthening, and an optimization-based heuristic approach – to improve computational performance for solving the LSP

M

model to optimality (or close to optimality).

4.1 Problem reduction

Due to the structure of the service network and after incorporating the previous operational policies (from

ED

Section 3.5), for each commodity, some of the flow variables in the model [LSP] may correspond to subpaths that cannot serve as a particular stage h in any feasible origin-to-destination path for that

PT

commodity. We can eliminate such variables a priori (by removing the sub-paths from Brhk ), thereby reducing the problem size. As an illustration, suppose a commodity k must be transported with less than

CE

three transshipments from its origin to destination. Let N k and N k respectively denote the set of ports nodes that are adjacent to the origin sk and to the destination tk of commodity k in the augmented

AC

network. Then, we can eliminate from Brhk all the sub-paths with either i  N k or j  N k (or both) since these sub-paths cannot be part of any feasible route from sk to tk . More generally, based on the topology of the augmented network, we can determine if a particular sub-path can serve as the hth transportation stage on any origin-to-destination route for commodity k. If this sub-path is not eligible for stage h, we can eliminate it from the set Brhk . As another variable reduction strategy, when an intermediate port has only one incoming (or outgoing) sub-path for a commodity k at a port j, we can merge this sub-path with each of its predecessor (or 20

ACCEPTED MANUSCRIPT

successor) sub-paths, and replace the individual sub-path flow variables with composite variables representing flows on the merged sub-paths. To illustrate this method, consider a port j that has only one incoming sub-path , entering this port, that can serve as the hth stage for a commodity k. Let AO hj 1,k denote the set of all outgoing sub-paths  j, l , r '  , leaving port j, that can serve as the (h+1)st

stage for commodity k. In this case, we can replace the original flow variables on the incoming sub-path

CR IP T

and outgoing sub-paths  j, l , r '  with combined flow variables that bypass the intermediate node j; further, we can delete the flow conservation constraint for the hth stage at node j. Specifically, we perform the following local transformations at node j to get an equivalent model: 

for each outgoing sub-path  j, l , r '  AOhj 1,k , corresponding to stage (h + 1) for commodity k, define a “composite” flow variable wijlhk,r ,r ' which denotes the flow of containers from port i to

AN US

port l via port j, using service r followed by service r'. This flow variable corresponds to movement over two consecutive stages, h and (h+1), for commodity k. Assign a cost (coefficient in the objective function (1)) of (cijrk  c kjlr ' ) to this variable, and replace the original flow variable u hjl,1,rjlk with wijlhk,r ,r ' in the constraints of the model formulation;

replace the original flow variable uijrh ,k with the sum model formulation; and,

 j ,l ,r 'AOhj 1,k

wijlhk,r ,r ' in the constraints of the

delete the flow conservation constraint (3) at node j corresponding to commodity k and stage h.

ED





M



Figure 4 pictorially illustrates this transformation. In this figure, intermediate port j has only one

PT

incoming sub-path for commodity k at stage h, permitting us to merge this sub-path with the succeeding

AC

CE

(outgoing) sub-paths.

Figure 4: Combining variables for a commodity k at a node j with in-degree of one

With the above transformation, the reduced model is equivalent to the original model (i.e., for each feasible solution to one model, the other model has an equal cost feasible solution). An analogous transformation applies when there is only one sub-path leaving an intermediate port j at some 21

ACCEPTED MANUSCRIPT

stage h; in this case, we can combine the appropriate incoming sub-path flow variables and the outgoing sub-path flow variables into composite variables. Moreover, this process extends to more complex situations (e.g., when a node j has only one incoming and one outgoing sub-path, in which case we can combine flow variables across three stages).

These techniques reduce the size of the model by

eliminating some variables and constraints.

CR IP T

4.2 Strengthening the LSP model

We refer to the formulation [LSP] after applying the variable reduction steps of Sections 3.4, 3.5, and 4.1 as the Base model for the LSP problem. We next discuss some valid inequalities to strengthen this model, i.e., increase the value of its linear programming (LP) relaxation; in turn, this lower bound improvement permits solvers that apply branch-and-bound to prune more branches and reduce computational time. The

AN US

inequalities we propose are essentially disaggregate forcing constraints that combine the linkage between the flow variables and service selection variables together with the capacity restrictions on individual sailing edges.

The validity of the inequalities stem from the fact that the flow of a commodity k on a sub-path cannot exceed the commodity’s demand d k and the capacity of that sub-path. Further, this flow is

M

permitted only if the corresponding service r is selected.

These requirements, together with the

observation that the capacity of a sub-path is the smallest edge capacity over all the sailing edges in the

ED

sub-path, imply that the following set of forcing constraints are valid for the LSP model: uijrhk  min{dk ,min eSEijr ge }xr

k  K , r  R,  i, j, r  Brhk , h  1,2,..., hk ,

(9)

PT

where SEijr denotes the subset of sailing edges in service r that correspond to sub-path . These constraints are analogous to disaggregate forcing constraints that have been proposed previously to

CE

strengthen the formulation of facility location and fixed charge network design problems (e.g., Balakrishnan, Magnanti, and Wong 1989, Aardal, Pochet, and Wolsey 1995). However, for the LSP problem, these constraints are too numerous to include in the Base model (since we require one constraint

AC

for every commodity, sub-path, and stage). Instead, we next propose two other families of disaggregate forcing constraints that are both tighter than constraints (9) and fewer in number. Disaggregate (commodity-edge) constraints We can strengthen inequality (9) by aggregating mutually exclusive flows on the original sailing edges. For a given commodity and service, we consider the total flow on all sub-paths at any stage h that have a common sailing edge e. This total flow must not exceed the capacity of edge e since, with positive

22

ACCEPTED MANUSCRIPT

transportation costs, no optimal solution will route a commodity multiple times over the same edge. Further, the total flow must be less than the commodity’s demand, implying that the following inequalities are valid: hk

 

h1 i , j ,r :eSEijr

uijrhk  min{d k , ge }xr k  K , r  R, e  E r .

(10)

CR IP T

These inequalities are stronger than constraints (9). To see this, suppose the sailing edge e '  Eijr has the lowest capacity among all sailing edges in sub-path . Then, the right-hand side coefficient (of variable xr ) in constraint (10) corresponding to edge e' is the same as the right-hand side coefficient in constraint (9); however, the left-hand side of constraint (10) sums the flows over multiple sub-paths and stages whereas the left-hand side of constraint (9) contains only one of these flow variables. Hence, constraints (10) are at least as tight as constraints (9). Moreover, constraints (10) require only one

for every sub-path of the service. Disaggregate (commodity-service) constraints

AN US

inequality per edge of each service (for each commodity), whereas constraints (9) require one inequality

We can also aggregate the flows on the sub-paths associated with each service r to obtain another class of

M

valid inequalities that tighten the Base model and are not implied by inequalities (10). This class of inequalities is motivated by the observation that the flow of a commodity k on a service r at stage h must not exceed demand dk if the service is chosen. Hence, the following constraint is valid: i , j ,r Brhk

uijrhk  min{dk ,max i , j ,r Bhk min eSEijr ge }xr

ED



r

k  K , r  R, h  2,..., hk  1 .

(11)

PT

For h = 2, …, hk  1 , constraint (11) strengthens formulation [LSP] (for h = 1 and h  hk , this inequality is implied by our previous constraints).

CE

Our implementation strengthens the Base model [LSP] by adding to it an effective subset of constraints from the inequality classes (10) and (11) before initiating the branch-and-bound procedure (Section 5.4

AC

provides more details of our approach). 4.3 LP-based heuristic procedure To obtain good solutions for the LSP problem, we apply a heuristic method that iteratively rounds (up or down) fractional values for the service selection ( xr ) variables in the solution to the LP relaxation of the problem. At each iteration, we select the highest or lowest fractional value among all the fractional xr values in the current LP solution, round this value to one or zero respectively, and re-solve the LP. If

23

ACCEPTED MANUSCRIPT

rounding the variable to one violates the fleet availability constraint, we set it to zero. The procedure stops when all the service selection variables have integer values.

The following pseudo-code

summarizes this procedure. Solve LP relaxation of LSP model

fix integer elements of x to current values of 0 or 1 pick the smallest fractional element xr ' , if xr '  0 then fix xr ' to 0 else pick the largest fractional element xr " ,

CR IP T

while x is fractional

if feasible in terms of fleet availability then fix xr " to 1

AN US

else fix xr " to 0 re-solve LP relaxation endwhile

The threshold value  for rounding down or up can be adjusted. During our computational tests, we

M

observed that only services with a relatively high fractional value (e.g., around 0.7) are included in the optimal solution. Further, some services have a high initial fractional value but are not necessarily

ED

included in final solution. Therefore, our implementation first rounds down low x-values (thereby discarding some services) before rounding up x-variables with high fractional values. By eliminating unattractive services first, this approach retains the flexibility of using available ships for later choices

PT

instead of committing them early for services that are included. For our test problem instances (see Section 5), the LP-based heuristic yielded solutions that are close to

CE

optimality. In contrast, the initial upper bounds for the LSP problem generated by solvers such as Gurobi and CPLEX are often quite poor (e.g., the initial upper bound may select no service, thus incurring high

AC

penalty costs for not meeting demand). Since our heuristic procedure identifies good solutions, we can use it to generate the initial upper bound to warm-start any branch-and-bound procedure. The heuristic solution can also serve as an interim recommendation to the planner if solving the problem optimally takes a long time. In the next section, we report computational results when the LSP model is applied to practical problem instances, and assess the benefits of our enhancements such as adding valid inequalities and applying the LP-based heuristic procedure. 24

ACCEPTED MANUSCRIPT

5. Computational Results We implemented the model and solution method in C++, using the Boost Graph Library to handle the graph construction and preprocessing, and solved the LP relaxation and mixed-integer programs using Gurobi 6.0. The tests were performed on a computer with an Intel Xeon CPU X5550 2.67GHz and 24 GB RAM. When solving the mixed-integer program, we set a CPU time limit of 12 hours (43,200

CR IP T

seconds). In our later summary of computational results (Tables 3 and 4), if the branch-and-bound procedure terminates within the time limit with the optimal solution (i.e., with a gap of 10–8 or less), we record the final gap as zero, and show the actual computational time. Otherwise, if the procedure stopped because of the time limit, we report the final gap (between the best upper and lower bounds) at termination.

AN US

5.1 Test problems

We tested problem instances based on the data for four common container shipping sectors – Baltic Sea (Baltic), West Africa (WAF), Mediterranean (Med), and Pacific (Pac) – provided in the Liner-Lib benchmark suite (www.linerlib.org, Brouer et al. 2014a). For each of these sectors, we have data on the origin and destination ports of each commodity, forecasted demand (containers), and freight rate.

M

Information for each port in the sector includes name, longitude/latitude, country, geographical region, cabotage rules, unit load/unload cost, unit transshipment cost, fixed port call costs, and variable port call

ED

costs that can vary with vessel capacity. The fleet consists of different vessel classes with varying TEU capacities; for the test problem instances, these capacities are 900 TEU, 1,600 TEU, 2,400 TEU, and 4,800 TEU. The number of available vessels varies by problem instance, with smaller instances having

PT

fewer vessels than larger instances. For each vessel, the data includes its capacity, bunker consumption for sailing and idling, fees for traversing the Suez and Panama canals, and the time-charter rate. We

CE

calculate the cost parameters needed for our model (e.g., cost for each service and unit penalty for not meeting demand for each commodity) using the approach described in Brouer et al. (2014a).

All

AC

demands can be transshipped at most twice. To generate candidate services, we applied the matheuristic described in Brouer et al. (2014b).

This

method permits constructing complex routes that visit the same port multiple times; our test problems contained several such routes. For each service, the vessel class assigned to this service takes into account the draft requirements (e.g., larger vessels cannot visit small ports) and canal limitations (larger vessels will have to use an alternative and usually longer path) of the ports on the route. By applying the heuristic method multiple times, we obtained different sets of services that we added to the candidate

25

ACCEPTED MANUSCRIPT

pool. Table 1 summarizes the dimensions of our problem instances. The Baltic sector is a relatively small problem that is easy to solve. No. of Ports

Baltic Sea (Baltic) West Africa (WAF) Mediterranean (Med) Pacific (Pac)

12 19 39 45

No. of Commodities 22 38 369 722

No. of Vessels (vessel classes) 6 (2) 42 (2) 28 (3) 100 (4)

No. of Services 7 24 22 31

CR IP T

Sector

Table 1: Characteristics of test problems

5.2 Model dimensions

Table 2 summarizes the sizes of the problem formulation for the different problem instances. Column 2

AN US

shows the number of sub-paths before and after applying the cabotage rules (as discussed in Section 3.5). As we might expect, cabotage rules eliminate more sub-paths for larger problems since these problems cover more countries. Column 3 shows the number of commodity flow variables in formulation [LSP] after applying the variable elimination methods of Section 3.4 as well as the cabotage rules and regional policies discussed in Section 3.5. Regional policies only apply (and yield reductions) to the Pacific sector

M

problem instance because the routes in this problem span ports in both the Asian east coast and the west coast of the Americas.

Baltic WAF Med Pac

No. of sub-paths No. of commodity flow variables After Original Initial * After reduction cabotage rules

ED

Problem

CE

PT

122 422 1,218 2,910

*

121 404 1,013 2,347

1,566 9,347 139,987 791,221

1,329 8,727 133,220 780,845

No. of constraints 233 920 11,990 32,875

Table 2: Problem reduction and formulation sizes

after applying the variable elimination strategies discussed in Sections 3.4 and 3.5

AC

Column 4 shows the number of flow variables remaining after we apply the reduction method described in Section 4.1 that determines which sub-paths are eligible for each stage of a commodity and eliminates the flow variables for sub-paths that are not eligible. For the problem scenarios we tested, only modest reduction was possible by combining sub-path flow variables into composite variables (Section 4.1) since the augmented network has few nodes with in-degree or out-degree of one (as required for this reduction). Column 5 shows the number of constraints in the model. All of the following computational results are based on the reduced model after applying all of our preprocessing strategies. 26

ACCEPTED MANUSCRIPT

5.3 Results for the Base model The first four columns in Table 3 show the results for the Base model, before adding our valid inequalities (10) and (11), and for the Strong model after adding these inequalities. For these computational runs, we did not apply our heuristic algorithm to warm-start the branch-and-bound procedure (the improvement

Base Model

Problem

Strong Model

Baltic

143

Final Gap* (%) 0

B&B nodes 0

CPU Time (sec.) 0.21

# Valid ineq. retained 245

WAF

18

0

1,722

3.93

1,288

Med

61

0

169

543

6,939

Pac

82

11

512

> 43,200

18,165

Initial LB (%)

B&B nodes

53

Final Gap * (%) 0

0

CPU Time (sec.) 0.13

9

0

1,012

3.31

7

0

46

304

35

0

319

23,678

AN US

Initial LB (%)

CR IP T

due to warm-start is discussed in Section 5.5).

Table 3: Computational Results for base model and after model strengthening (without warm-start)

M

Initial LB % = (Best Final Upper Bound – Initial LP value)/Best Final Upper Bound Final Gap % = (Best Final Upper Bound – Final Lower Bound)/Best Final Upper Bound # Valid inequalities = Number of valid inequalities (10), (11) included in Strong model * gaps of 0 indicate that the branch-and-bound procedure terminated before the time limit, when the final gap was 10–8 or less

ED

For the Base model, the Initial lower bound (LB) gap, defined as the (Best Final Upper Bound – Initial LP value)/Best Final Upper Bound) exceeds 50% for three out of the four problems, where the Best Final Upper Bound is the value of the best solution at termination (in all problem instances, the Strong model

PT

gave the best final solution). Further, using the Base model, the Pac scenario is not solved to optimality

CE

within the time limit of 12 hours (but the gap reduces from 82% to a final value of 11% at termination). 5.4 Effect of strengthening the model Table 3 shows that the Base model’s LP is weak, resulting in a relatively high Initial LB gap. The last

AC

five columns of Table 3 show the improved results using the Strong model. The valid inequalities (10) and (11) can be added a priori before applying branch-and-bound or added as cutting planes (when violated) at intermediate stages of the solution procedure. For our test problems, we found that an effective strategy is to generate all the valid inequalities a priori, but include only a selected subset in the model when solving the integer program. Specifically, we first add all the inequalities to formulation [LSP] and solve its LP relaxation. Then, we retain in the model only those inequalities that are tight (or nearly tight, with a slack not exceeding 0.1) and remove the other inequalities (that are not tight) before 27

ACCEPTED MANUSCRIPT

applying the branch-and-bound procedure. As the results in Table 3 demonstrate, including the subset of tight valid inequalities (10) and (11) significantly improves the initial lower bound of the formulation for all instances. For the larger scenarios, the strong model’s initial gap is considerably lower, e.g., in Med, the initial gap decreases from 61% to 7%.

This improvement in initial lower bound significantly

accelerates the branch-and-bound procedure. The total computational time for the Strong model is up to 70% lower than that for the Base model, and the difficult Pac problem scenario can now be solved to

CR IP T

optimality well before the branch-and-bound termination time limit. Also, with the valid inequalities, the procedure needs to explore fewer branch-and-bound nodes before terminating. 5.5 Heuristic performance

We next discuss the performance our LP-based iterative rounding heuristic described in Section 4.3 (since the Baltic problem is easy to solve, we did not apply the heuristic to this problem). The first column of

AN US

Table 4 show the quality of the heuristic solutions, measured as the gap between the heuristic upper bound and the final upper bound when the branch-and-bound procedure terminates (as a fraction of the final upper bound). Columns 2 and 3 show the number of rounding iterations for the heuristic, and its computational time.

M

ED

Problem

Heuristic performance No. of Heuristic Heuristic rounding solution gap runtime iterations (%) (sec.) 8 2.9 0.15 12 3.4 38 11 3.5 301

PT

WAF Med Pac

MIP with Heuristic warm-start MIP Final B&B MIP Gap (%) nodes runtime (sec.) 0 869 2.12 0 32 237 0 352 22,633

Table 4: Computational results for LP-based rounding heuristic.

CE

Heuristic gap % = (Heuristic Upper Bound – Final Upper Bound)/Final Upper Bound MIP Final Gap % = (Final Upper Bound – Final Lower Bound)/Final Upper Bound

As the results in this table demonstrate, the heuristic is very quick and produces good solutions that are

AC

within 3.5 % of the final upper bound obtained after branch-and-bound. Only around ten rounding iterations (iterative solution of the LP relaxation) are needed, and the heuristic’s computational time is a small fraction of the time needed to for the branch-and-bound procedure. Although solvers such as Gurobi and CPLEX can generate initial upper bounds using built-in methods, for the LSP problem, these initial upper bounds are very poor and have costs that are often orders of magnitude higher than the optimal value. (However, the solvers improve this solution quickly during the

28

ACCEPTED MANUSCRIPT

branch-and-bound process.) The last four columns show the effect of providing the heuristic solution as the initial upper bound to warm-start the branch-and-bound procedure with the Strong model. Compared to the results without this warm-start method (in Table 3, for the Strong model), the initial upper bound improves very significantly with warm-start. The overall running time (including the time to find the heuristic solution) also decreases, although only modestly. In general, the best solution is found quickly when starting with a good solution, but most of the computational time is spent for proving optimality of

CR IP T

this solution.

Since the heuristic procedure yields solutions that are within 3.5% from optimal and does not require much computational time, just applying this heuristic approach may be useful when planners require quick solution times to analyze different scenarios and perform sensitivity analysis.

AN US

Table 5 summarizes the characteristics of the best solutions found. The availability of ships permits selecting only around half (or less) of the candidate services, but the optimal solution is able to cover a significant portion of the available demand (up to 98% for the largest problem instance). Thus, the model’s solution obtains high revenue despite the limits on number of transshipments.

ED

Baltic WAF Med Pac

No. of services selected out of candidate services 3/7 10/24 7/22 18/31

M

Problem

% of available demand transported 87.9 % 96.9 % 90.5 % 98.1 %

6. Conclusions

PT

Table 5: Solution characteristics of best found solutions using the strengthened model.

CE

Our Liner Service Planning model provides a new approach to support tactical decisions by focusing on selecting services from a specified set of candidate services and explicitly considering service assurance

AC

in terms of maximum number of transshipments. Incorporating transshipment limits for each demand gives realistic cargo routing solutions. Our modeling approach facilitates including many operational requirements such as cabotage rules and regional routing policies to ensure that the solutions are implementable in practice. Moreover, such rules and policies also reduce the model size, thereby improving computational performance. Unlike some prior models, our model can accurately account for the cost of transshipments and loading/unloading operations. The valid inequalities that we developed significantly improve the lower bound and reduce computational time, permitting us to solve even the largest problem to optimality within the time allotted. Our LP-based rounding heuristic was able to find 29

ACCEPTED MANUSCRIPT

good quality solutions in relatively short time, making it suitable both as a stand-alone procedure and as a means to warm-start the exact solution method. Future research directions include developing other modeling approaches or further methodological enhancements to solve much larger instances, and other optimization-based approaches including path-based decomposition or heuristic strategies.

Acknowledgements

CR IP T

The authors wish to thank Professors Stefan Ropke and David Pisinger for valuable comments. Christian Vad Karsten was supported in part by The Danish Maritime Fund under the Competitive Liner Shipping Network Design project.

References

AN US

Aardal, K., Pochet, Y., Wolsey, L.A. (1995) Capacitated facility location: Valid inequalities and facets. Mathematics of Operations Research, 20, 562-582. Agarwal, R., Ergun, Ö. (2008) Ship scheduling and network design for cargo routing in liner shipping. Transportation Science, 42(2), 175-196.

M

Álvarez, J.F. (2009) Joint routing and deployment of a fleet of container vessels. Maritime Economics & Logistics, 11(2), 186-208. Balakrishnan, A., Altinkemer, K. (1992) Using a hop-constrained model to generate alternative communication network design. ORSA Journal on Computing, 4(2), 192-205.

ED

Balakrishnan, A., Magnanti, T.L., Wong, R.T. (1989) A dual-ascent procedure for large-scale uncapacitated network design. Operations Research, 37, 716-740.

PT

Bell, M.G., Liu, X., Angeloudis, P., Fonzone, A., Hosseinloo, S.H. (2011) A frequency-based maritime container assignment model. Transportation Research Part B: Methodological, 45(8), 1152-1161.

CE

Brouer, B.D., Alvarez, J.F., Plum, C.E., Pisinger, D., Sigurd, M.M. (2014a) A base integer programming model and benchmark suite for liner-shipping network design. Transportation Science, 48(2), 281312.

AC

Brouer, B.D., Desaulniers, G., Pisinger, D. (2014b) A matheuristic for the liner shipping network design problem. Transportation Research Part E: Logistics and Transportation Review, 72, 42-59. Brouer, B.D., Pisinger, D., Spoorendonk, S. (2011) Liner Shipping Cargo Allocation with Repositioning of Empty Containers. INFOR, 49 (2), 109–124. Christiansen, M., Fagerholt, K., Nygreen, B., Ronen, D. (2013) Ship routing and scheduling in the new millennium. European Journal of Operational Research, 228(3), 467-483. Christiansen, M., Fagerholt, K., Ronen, D. (2004) Ship routing and scheduling: Status and perspectives. Transportation Science, 38(1), 1-18. Dahl, G., Gouveia, L. (2004) On the directed hop-constrained shortest path problem. Operations 30

ACCEPTED MANUSCRIPT

Research Letters, 32, 15-22. Dahl, G., Gouveia, L., Requejo, C. (2006) On formulations and methods for the hop-constrained minimum spanning tree problem. In Handbook of Optimization in Telecommunications. Springer. 493515. Gouveia, L., Magnanti, T.L. (2003) Network flow models for designing diameter-constrained minimumspanning and Steiner trees. Networks, 41, 159-173.

CR IP T

Karsten, C. V., Brouer, B. D., Desaulniers, G., & Pisinger, D. (2016) Time constrained liner shipping network design. Transportation Research Part E: Logistics and Transportation Review, in press. Karsten, C.V., Pisinger, D., Ropke, S., Brouer, B.D. (2015) The time constrained multi-commodity network flow problem and its application to liner shipping network design. Transportation Research Part E: Logistics and Transportation Review, 76, 122-138. Kjeldsen, K.H. (2011) Classification of ship routing and scheduling problems in liner shipping. INFOR: Information Systems and Operational Research, 49(2), 139-152.

AN US

Meng, Q., Wang, S. (2011) Liner shipping service network design with empty container repositioning. Transportation Research Part E: Logistics and Transportation Review, 47(5), 695-708. Meng, Q., Wang, S. (2012) Liner ship fleet deployment with week-dependent container shipment demand. European Journal of Operational Research, 222(2), 241-252. Meng, Q., Wang, S., Andersson, H., Thun, K. (2014) Containership routing and scheduling in liner shipping: overview and future research directions. Transportation Science, 48(2), 265-280.

M

Mulder, J., Dekker, R. (2014) Methods for strategic liner shipping network design. European Journal of Operational Research, 235(2), 367-377.

ED

Nemhauser, G.L., Wolsey, L.A. (1988) Integer and Combinatorial Optimization. Wiley-Interscience. Plum, C. E., Pisinger, D., Salazar-González, J. J., Sigurd, M. M. (2014) Single liner shipping service design. Computers & Operations Research, 45, 1-6.

PT

Plum, C. E., Pisinger, D., Sigurd, M. M. (2014) A service flow model for the liner shipping network design problem. European Journal of Operational Research, 235(2), 378-386.

CE

Reinhardt, L.B., Pisinger, D. (2012) A branch and cut algorithm for the container shipping network design problem. Flexible Services and Manufacturing Journal, 24(3), 349-374.

AC

Song, D.P., Dong, J.X. (2012) Cargo routing and empty container repositioning in multiple shipping service routes. Transportation Research Part B, 46 (10), 1556-1575. Wang, S., Meng, Q. (2014) Liner shipping network design with deadlines. Computers & Operations Research, 41, 140-149. Wang, S., Meng, Q., Lee, C.-Y. (2016) Liner container assignment model with transit-time-sensitive container shipment demand and its applications. Transportation Research Part B, 90, 135-155. Wang, S., Meng, Q., Sun, Z. (2013) Container routing in liner shipping. Transportation Research Part E, 49, 1-7.

31

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

World Shipping Council. (2009) The liner shipping industry and carbon emissions policy. http://www.worldshipping.org/pdf/liner_shipping_co2emissions_policy_september.pdf

32