1 Introduction

Since their introduction in the 1960s, shipping containers have become the main way to transport goods around the world. About 75% of the global trade volume is carried by sea, and half of that part is shipped in containers (Lee and Song 2017). Over the years, the vessels that are used to ship these containers between deep-sea ports all around the world have increased substantially and are now able to carry more than 10,000 containers. In general, the deep-sea port is not the final destination of a container and the container needs to be shipped further inland. The large number of containers that arrive at a container terminal on a single deep-sea vessel requires efficient transportation from the deep-sea port to the final destination, which is the so-called hinterland transportation.

For the hinterland transportation, three different modes of transportation could be used: barges, trains, and trucks. It is also possible to use a combination of these modes, for instance, shipping the container with a barge to an inland container terminal and using a truck to transport it from the inland terminal to its destination. The concept in which a container is allowed to be shipped on multiple modes of transportation is referred to as both multimodal and synchromodal transportation. The key difference between multimodal and synchromodal transportation is that in the latter the mode of transportation might still change when the container is already in transit (Dong et al. 2018). The European Commission and large ports in Europe, such as the port of Rotterdam, are aiming for a modal shift from trucks towards trains, and especially barges (European Commission 2011; Port of Rotterdam Authority 2011). The main advantage of barge and train transportation is that it is usually cheaper than using trucks. Furthermore, barge and train transportation is also more eco-friendly than truck transportation. Moreover, deep-sea ports are often located in urban areas with a considerable amount of road congestion, so not using trucks can help in reducing road congestion.

However, in the European Union, the modal shift has not yet taken place. In the last decade, the share of inland waterway, train, and truck transport have all remained constant in the European Union (Eurostat 2019). A major drawback of barge and train transportation is that the planning process is much more challenging than for truck transportation. First of all, barges and trains have longer transit times than trucks and as a result, the transportation plan needs to be made earlier. It often happens that not all the required information is available when the planning for barges or trucks needs to be made. Moreover, only if the number of shipped containers is large enough, the economies of scale of barges and trains ensure that these modes of transportation are cheaper than truck transportation. Finally, barges are also heavily influenced by the congestion at deep-sea ports. The growing volume of containers that are transshipped in combination with the increasing size of deep-sea vessels puts a large pressure on the operation of the deep-sea terminals. Hence, the number of containers that can be loaded and unloaded on a barge is limited. We will refer to the number of containers that can be (un)loaded at a terminal as the number of moves.

Unfortunately, this number of moves at the deep-sea terminals is often unknown, because of the delay of deep-sea vessels. These deep-sea vessels are delayed, because of, for instance, bad weather conditions or a late departure at a previous port. The service of deep-sea vessels has priority over that of barges and thus the delay of deep-sea vessels also influences the service of barges which has become rather unreliable. For example, in the port of Rotterdam, waiting times of more than 8 h are not uncommon. At the moment, barges spend about 30–40% of the time they are in the port of Rotterdam waiting to be served at container terminals (Port of Rotterdam Authority 2019). In this paper, we consider a planning problem faced by an inland terminal in the port of Amsterdam. In this problem, it has to be decided which containers to load on which barge. If many containers are loaded, the expected costs for late minute adjustments, that are incurred when the number of moves is insufficient, are high. On the other hand, if only a few export containers are shipped per barge, a large number will be shipped per truck, which results in higher transportation costs.

The contribution of this paper is threefold: first, we present a new operational planning problem faced by an inland terminal at the port of Amsterdam. This is the first problem in the literature in which the uncertain service of barges at deep-sea terminals is studied. In this paper, we formulate this problem as a two-stage stochastic problem with recourse (Birge and Louveaux 2011). The second contribution is that we give a Sample Average Approximation (SAA) method to solve this problem. This method uses Monte Carlo sampling to solve stochastic optimization problems. To the best of our knowledge, we are the first to apply SAA to the field of hinterland container transportation. This technique allows us to solve more realistic problems than the existing methods in the literature. Although the SAA method can produce (almost) optimal solutions for small instances, its running time for larger instances is larger. Therefore, our third contribution is that we propose a heuristic based on Stochastic Programming (SP).

The remainder of this paper is organized as follows. The relevant literature is discussed in Sect. 2. In Sect. 3, the problem is described in more detail and it is translated into a mathematical model in Sect. 4. Three different solution methods will be discussed in Sect. 5 and the numerical results for these three methods are given in Sect. 6. Finally, the conclusions will be presented in Sect. 7.

2 Literature review

Over the last years, there has been an increased interest in operational problems for multimodal and synchromodal transportation. In these operational problems, containers have to be assigned to existing barge, train, or truck services. We divide the literature in three types of assignment problems: offline assignment problem, online assignment problems and stochastic assignment problems. In offline assignment problems, the containers are assigned to the services if for every container the relevant information is known. This in contrast to online problems, in which each incoming container has to be assigned directly to a mode of transportation before having any information about container shipments that are booked later. In stochastic assignment problems, part of the information is uncertain at the moment the planning is made.

Let us first focus on the offline assignment problems. In Baykasoglu and Subulan (2016), the optimal offline assignment of containers to modes in a multimodal network is presented. They use a multiobjective problem formulation in which the transportation costs and the \({\text {CO}}_2\)-emissions are minimized, while the service level is maximized. Another offline multiobjective problem, that is solely focusing on import containers, is formulated in Zweers et al. (2019). In this problem, the utilization of a barge is maximized while minimizing the demurrage and storage costs and visiting as few terminals as possible. In Heggen et al. (2019), the planning for the long-haul transportation is combined with the planning for the transportation by truck between the inland terminal and the final destination. This integrated approach results in cost savings and better utilization of the capacity of the long-haul transportation mode. In the offline problem of Pérez Rivera and Mes (2016), the planner can decide to postpone a shipment, which does not give any direct costs. Nevertheless, the shipment has to be shipped after the planning horizon, and since the characteristics of the containers which arrive after the planning horizon are unknown, the future costs are stochastic.

Second, we discuss the online assignment problems. An online problem in which the transportation costs and penalties for being late are minimized is studied by van Riessen et al. (2016). They compute the optimal offline solutions and use them to construct a decision tree. This tree is used to decide online on how to assign incoming containers. Another online planning problem is studied by Mes and Iacob (2016). Their objective is to minimize the costs, delays and \({\text {CO}}_2\)-emissions of a shipment. For each incoming container, the k best paths in the network for these objectives are constructed. The planners can then select the best one. Finally, Wang et al. (2016) study the assignment of containers to modes of transportation from a revenue management perspective. In this work, the flexibility of accepting or refusing an incoming order is introduced. If a shipment is accepted, the service on which it will be shipped has to be decided.

To the best of our knowledge, there are only three papers in the literature that consider stochastic assignment problems (Pérez Rivera and Mes 2017; Zhang et al. 2016; Zuidwijk and Veenstra 2015). In these three papers, the unknown factor is the number of containers that needs to be transported. This uncertainty is caused by the fact that the arrival of a deep-sea vessel could be delayed or that the release of a container by customs might take longer. In Zuidwijk and Veenstra (2015), this problem is studied for the first time. They consider a problem in which there is only a single barge that visits one terminal. The moment the barge visits the terminal has to be decided, in which a trade-off is made between the number of containers that can be transported by the barge and the arrival time of the barge at the inland port. Optimal Pareto-frontiers are given for different information scenarios and with that, the value of information can be determined. Zhang et al. (2016) study a problem in which a hinterland operator needs to transport the containers on a deep-sea vessel that is approaching the port. The planner needs to decide on the amount of capacity for each of the different modes of transportation. Each mode of transportation incurs different costs but also has a different speed. Certain containers need to be shipped with a fast and expensive mode of transportation, because otherwise, they will be too late at their final destination, while other containers can go on the slow and cheap mode of transportation. The total demand for containers is known, but the number of containers that needs to go on each mode of transportation is unknown. Similar to Zuidwijk and Veenstra (2015) and Zhang et al. (2016) also study a theoretical problem that they can solve to optimality and look into different information scenarios. The problem described in Pérez Rivera and Mes (2017), is the most similar to our problem. In this problem, every day a barge makes a round-trip from a single inland terminal to a set of deep-sea terminals. They formulate a Markov Decision Problem to decide on both which terminals to visit with a barge and how many containers of each type should be delivered to or be picked up from a certain terminal. For each terminal that is visited by a barge, a penalty has to be paid. The uncertainty in this model lies in the fact that the demand for each of the container types is unknown.

The three problems studied by Pérez Rivera and Mes (2017), Zhang et al. (2016) and Zuidwijk and Veenstra (2015) assume that the costs for transporting a container on a barge are the same for all containers. This assumption mainly serves to make calculations easier, but it is often unrealistic in practice. To include different costs per container, we will use the technique called Sample Average Approximation (SAA), introduced by Kleywegt et al. (2001). In Shapiro et al. (2009), a theoretical overview of the SAA method is given and we refer to Kim et al. (2015) for an overview that is more orientated to the application of the SAA method. The SAA method has been applied to many different areas, but we focus on problems in transportation. In Verweij et al. (2003), the SAA method is used to solve the Stochastic Vehicle Routing Problem. Another transportation problem, the Stochastic Multiperiod Location Transportation Problem, is also solved using SAA in Klibi et al. (2010). In this problem, the design of a transportation network and the routing decisions over multiple planning periods inside this network are taken into account. As the SAA method is hard to solve for large instances of this problem, also heuristics are proposed to solve the problem hierarchically. Another application of the SAA that is related to our work is in Long et al. (2012) in which SAA is used to solve a problem in empty container repositioning.

Finally, the paper of Toktas et al. (2006), in which the Generalized Assignment Problem with stochastic capacities is studied, has to be mentioned. Our problem can be seen as a generalization of this problem. In Toktas et al. (2006), simple heuristics to deal with the uncertainty in the capacity constraints are compared. We will use the two methods that perform best in their study to our problem and discuss them in more detail in Sect. 5.3.

3 Problem formulation

We consider a problem in which a single inland terminal needs to export containers to multiple deep-sea terminals and import containers from the same deep-sea terminals. There is a fixed planning horizon in which these containers need to be transported. In that planning period, the inland terminal has a set of barges at its disposal. Besides barges, containers can also be transported on trucks. For both the import and export containers, we assume that every container has the following properties:

  • given size in Twenty feet Equivalent Unit (TEU);

  • fixed transportation costs for each barge and truck;

  • deep-sea terminal at which it needs to be picked up (import containers) or needs to be delivered (export containers).

TEU is a standardized size for containers and most containers have a size of 1 or 2 TEU. The costs for transporting a container are fixed for each barge and truck, because the inland terminal has a contract with barge operators and trucking companies to ship containers for a given price. It is important to note that the costs for each container and barge combination are potentially different, because the barges could be operated by different barge operators. Moreover, for some containers demurrage, storage and detention costs need to be paid. Demurrage costs have to be paid for import containers that are located longer than a certain period at a deep-sea terminal. On top of that, if an import container arrives early at the inland terminal it has to be stored inside the terminal, which results in storage costs for the inland terminal. Furthermore, an export container needs to be back at the deep-sea terminal after a given number of days, because otherwise, detention costs need to be paid per day it is late. Especially, the demurrage and detention costs are substantial and influence the transportation plan. Allowing for different transportation costs for different barges also gives us the flexibility to include the fact that some containers cannot be transported on specific barges because it has, for instance, not arrived at the deep-sea terminal. If the costs for transporting a container on a specific barge is set high enough, the container will never be assigned to that barge.

Considering the barges, we make the following assumptions. Every barge:

  • makes a round trip starting at the inland terminal, going to the deep-sea port and returning at the inland terminal;

  • has a capacity in TEU;

  • has a maximum number of terminals it can visit in a round trip;

  • has an uncertain number of moves at every terminal.

The barge operator decides the moment that a barge is leaving the inland terminal and also the moment it needs to leave the deep-sea port. In a large deep-sea port, there are too many terminals to be all visited by a single barge. That is why the number of terminals that can be visited by a barge is limited. The inland terminal can decide on the terminals that are visited, but the route of the barge is decided by the barge operator. The congestion at the deep-sea port causes deep-sea terminals to limit the number of moves that can be done for every barge. The inland terminal is allowed to use all these moves solely to unload containers or load containers, but also each possible combination of the two.

For truck transportation we make the following two assumptions:

  • each truck can only ship a single container;

  • there is an unlimited number of trucks.

Trucking a container is more expensive than shipping a container on a barge, because each truck can only transport a single container at a time. Since truck transportation is relatively expensive not many trucks will be used in general for the transportation of containers. Hence, the assumption that the number of available trucks is unlimited is not restricting us much from reality. Furthermore, the benefit of this assumption is that there always exists a feasible solution. We do not consider other modes of transportation, such as trains, in this problem, although they could be incorporated in our framework.

In general, inland terminals are not located in the direct neighborhood of a deep-sea terminal. Therefore, the transport plan is made several days in advance. However, at that time, the number of moves at a deep-sea terminal is either unknown, or it is unreliable. Hence, the number of moves is modeled as a stochastic variable. In this problem, we assume that the total number of containers that can be loaded and unloaded at a terminal is only revealed after all barges have left the inland terminal. At that time, all barges are already loaded with export containers. In case a barge is shipping more export containers for a certain terminal than the number of moves, the barge cannot unload all its containers in the allocated slot. For each container that cannot be unloaded, a recourse action is required, which could be unloading the container at another terminal or returning it to the inland terminal. In both cases, extra transportation costs are incurred and the planner at the inland terminal has to do extra work, which may be costly. The main problem we are facing is to determine which terminals to visit with each barge and which export containers to load on each barge. Below, we give an illustration of the problem.

3.1 Problem illustration

In Fig. 1, an illustration is given of the problem we consider and a possible solution. It is important to note that this is not necessarily the optimal solution. In Fig. 1a, the input data of the problem is given. For simplicity, we ignore the cost structure in this example. We consider an example with two barges and both barges can visit at most two of the three deep-sea terminals. Besides that, the capacity of both barges is eight containers. The containers at the inland terminal are grouped with respect to the deep-sea terminal to which they need to go. For instance, container E12 is the export container with number 2 that needs to go to the first deep-sea terminal. We assume that every container can be shipped on both of the barges. Every deep-sea terminal has for every barge a certain number of moves available. Nevertheless, at the beginning of the planning phase, this number is unknown and therefore it is represented by a question mark in Fig. 1a.

Fig. 1
figure 1

Illustration of the problem formulation

In Fig. 1b, it has been decided which deep-sea terminals are visited by which barge and which export containers are loaded on which barge. Barge 1 is visiting deep-sea terminals 1 and 2 and barge 2 is visiting deep-sea terminal 2 and 3. In Fig. 1b, it is also depicted which export containers are loaded on the two barges. Although there is capacity left on the second barge, some containers for deep-sea terminal 3 are shipped per truck (containers E33, E34, and E35). These containers are shipped per truck, because the number of moves at the deep-sea terminals is still unknown at this stage.

In Fig. 1c, the numbers of moves at the deep-sea terminals are revealed. The number of moves at a deep-sea terminal is only revealed if a barge is visiting the terminal. For instance, deep-sea terminal 1 is only visited by barge 1 and thus the number of moves for this barge is revealed (2), but the number of moves for the second barge is not relevant. Four export containers for deep-sea terminal 1 were loaded on barge 1, but the number of moves is only two, meaning that for two of the four containers a recourse action is required.

In Fig. 1d, it is shown, for each barge, how many moves are left for the import containers. The number of moves left is the same as the number of moves in Fig. 1c minus the number of export containers unloaded at that terminal. For instance, the number of moves for barge 2 at deep-sea terminal 2 was three, but since also three export containers were unloaded at the deep-sea terminal, no more moves are left for barge 2. As the first barge has four moves left at the deep-sea terminal 2, that barge can ship four of the five import containers located at the second deep-sea terminal. The fifth container has to be shipped per truck, similar to the five containers located at deep-sea terminal 1 because no moves are left at that terminal for barge 1. The second barge has seven moves left at the third deep-sea terminal, but only three import containers are located at that terminal. In retrospect, it would have been better to load more export containers for terminal 3 on the second barge, because the number of moves was sufficiently large. The assignment of the import containers that is done in Fig. 1d is relatively easy compared to the assignment of the export containers in Fig. 1b because the number of moves has been revealed. Hence, the problem has become completely deterministic.

4 Mathematical model

In this section, the problem described in Sect. 3 is formulated as a two-stage stochastic problem. We consider a problem in which there are m barges and we use j as an index for the barges. Each barge j has a capacity of \(u_j\) TEU and can visit at most \(N_j\) terminals. Furthermore, there are \(n^e\) types of export containers and \(n^i\) types of import containers. We use t as an index both for the type of import and export containers. To a type of container belong all containers with the same characteristics. For an import (export) container these characteristics are: the size in TEU of a container \(w^i_t\) (\(w^e_t\)), the cost of transporting a container of type t on barge \(j=1, \ldots , m\) is \(c^i_{tj}\) (\(c^e_{tj}\)) and per truck it is \(c^i_{t0}\) (\(c^e_{t0}\)). Moreover, all containers in a type are located at the same terminal. The total number of import container and export containers of type t is given by, respectively, \(d^i_t\) and \(d^e_t\). In total there are l terminals and we use the index k to refer to a terminal. The set of import and export containers at terminal k is, respectively, \({\mathcal {I}}(k)\) and \({\mathcal {E}}(k)\). The number of moves by barge j at terminal k is given by the stochastic variable \(\varPhi _{jk}\). We use \(\phi _{jk}\) to refer to a realization of that stochastic variable. Finally, if more than \(\phi _{jk}\) containers from the sets \({\mathcal {I}}(k)\) and \({\mathcal {E}}(k)\) are assigned to barge j, then a cost of \(q_{jk}\) per container above the threshold \(\phi _{jk}\) is incurred. The notation for the input parameters is summarized in Table 10 in “Appendix 1”.

The remainder of this section is organized as follows: first, we formulate in Sect. 4.1 a model in which we assume that the number of moves at a terminal is deterministic. In Sect. 4.2, the two-stage stochastic problem with recourse is formulated.

4.1 Deterministic model

We first assume that the number of moves at terminal k for barge j is the deterministic value \(\phi _{jk}\). We use three different types of decision variables: \(X^e_{tj}\), \(X^i_{tj}\) and \(Y_{jk}\). The variable \(Y_{jk}\) is a binary variable indicating that terminal k is visited by barge j. The variable \(X^e_{tj}\) indicates the number of export containers of type t that are transported on barge j, if \(j= 1, \ldots ,m;\) and if \(j=0\), it gives the number of export containers of type t on a truck. The variables \(X^i _{tj}\) represent the same, but then for import containers. This deterministic problem can be modeled as an Integer Linear Program (ILP) in the following way:

$$\begin{aligned} \min \quad \sum _{t =1}^{n^i}\sum _{j=0}^{m} c^i_{tj}X^i_{tj} + \sum _{t=1}^{n^e} \sum _{j=0}^m c^e_{tj}X^e_{tj} \end{aligned}$$
(1)

subject to:

$$\begin{aligned}&\sum _{j=0}^mX^i_{tj} = d^i_t \qquad t =1, \ldots , n^i \end{aligned}$$
(2)
$$\begin{aligned}&\sum _{j=0}^mX^e_{tj} = d^e_t \qquad t =1, \ldots , n^e\end{aligned}$$
(3)
$$\begin{aligned}&\sum _{t=1}^{n^i} w^i_tX^i_{tj}\le u_j \qquad j =1, \ldots , m \end{aligned}$$
(4)
$$\begin{aligned}&\sum _{t =1}^{n^e} w^e_tX^e_{tj}\le u_j \qquad j =1, \ldots , m\end{aligned}$$
(5)
$$\begin{aligned}&\sum _{t \in {\mathcal {I}}(k)}X^i_{tj} + \sum _{t \in {\mathcal {E}}(k) } X^e_{tj} \le \phi _{jk} \qquad j=1, \ldots , m,\quad k = 1, \ldots , l \end{aligned}$$
(6)
$$\begin{aligned}&\sum _{k =1}^l Y_{jk} \le N_j \qquad j=1, \ldots , m \end{aligned}$$
(7)
$$\begin{aligned}&X^i_{tj} \le d^i_t Y_{jk} \qquad j=1, \ldots , m, \quad k =1, \ldots , l,\quad t =1, \ldots , n^i \end{aligned}$$
(8)
$$\begin{aligned}&X^e_{tj} \le d^e_t Y_{jk} \qquad j=1, \ldots , m, \quad k =1, \ldots , l,\quad t =1, \ldots , n^e\end{aligned}$$
(9)
$$\begin{aligned}&X^i_{tj}\in {\mathbb {N}}_0 \qquad t =1, \ldots n^i,\quad j =1, \ldots , m\end{aligned}$$
(10)
$$\begin{aligned}&X^e_{tj} \in {\mathbb {N}}_0\qquad t =1, \ldots n^e,\quad j =1, \ldots , m\end{aligned}$$
(11)
$$\begin{aligned}&Y_{jk} \in \{0,1\} \qquad j =1, \ldots , m, \quad k = 1, \ldots , l. \end{aligned}$$
(12)

The objective function (1) consists of two sums. The first sum represents the total costs for barge and truck transportation of the import containers. The second sum is the same as the first sum, except for the fact that it represents the total costs for the export containers. Constraints (2) and (3) ensure that for each type of import, respectively export container, the entire demand is transported. Constraints (4) and (5) ensure that the total number of import and export containers assigned to a barge does not violate the capacity of a barge. The number of containers on a barge that can be transshipped at a terminal is limited by constraint (6). Furthermore, in constraint (7) the number of terminals that can be visited by a barge is limited. Constraints (8) and (9) couple the X-variables and Y-variables. Constraints (10) and (11) ensure that the number of export and import containers on a truck and barge is non-negative. The Y-variables are ensured to be binary in constraint (12).

The capacity of the barge is considered for all import and export containers [constraints (4) and (5)]. However, there is no constraint to enforce that the number of export containers unloaded at, or after, the kth terminal and the number of loaded import containers before visiting the kth terminals is not exceeding the capacity of the barge. Nevertheless, (4) and (5) indirectly imply this type of constraint, because, if these two constraints are satisfied, there exists a route visiting the terminals such that the capacity constraint is never violated in the entire trip. To see that, let us denote the difference between the export containers unloaded at terminal k by barge j and the number of import containers loaded at terminal k by barge j as \(\delta _{jk}:=\sum _{t \in {\mathcal {E}}(k)}X_{tj}^e - \sum _{t\in {\mathcal {I}}(k)}X_{tj}^i\). If barge j visits the terminals in non-increasing order of \(\delta _{jk}\), then total the number of containers on barge j will decrease after each terminal visit as long as \(\delta _{jk}>0\) and thus the capacity will never be exceeded. If \(\delta _{jk}\) is negative for a terminal, more import containers are loaded than export containers are unloaded for that terminal, which could potentially lead to a violation of the capacity of the barge. Nevertheless, for all following terminals that are visited by the barge j, the total number of unloaded export containers is smaller than the import containers that are loaded at those terminals. Since all import containers satisfy the capacity constraint, it is impossible that, after a visit to a terminal, the total remaining export containers in combination with the already loaded import containers are violating the capacity constraint of a barge.

4.2 Two-stage stochastic problem with recourse

To incorporate the uncertainty in the number of moves for a barge at a terminal, we model the problem as a two-stage stochastic problem with recourse. For an introduction to stochastic problems with recourse, we refer to Birge and Louveaux (2011). The general idea of a two-stage stochastic problem is that the decision variables can be divided into two groups: the first-stage and the second-stage decision variables. The first-stage decision variables have to be decided before any of the realizations of the stochastic variables are known, whereas the value for the second-stage decision variables can be determined when the realizations of the stochastic variables are known. In our problem, the first-stage decisions are the decisions about which terminals to visit by which barge and which export containers to ship on which barge or truck. The second-stage decisions are the decisions concerning the transportation of import containers.

In a recourse problem, it is assumed that besides the second-stage decision variables, also recourse actions can be taken after the realization of the random variables. Usually, these actions have high costs and are only chosen to ‘repair’ the first-stage decisions in case the realization of the random variables is ‘bad’. In our problem, the recourse action is the opportunity to unload export containers at another terminal and use trucks to ship them to the original destination. This action is chosen if, for a certain terminal, more export containers are on a barge than the number of moves at that terminal.

In the ILP-formulation of (1)–(12) above, the realization \(\phi _{jk}\) was assumed to be deterministic and known before any decision is made. Nevertheless, in the two-stage stochastic problem, the realization \(\phi _{jk}\) is only known when the decisions for the \(X^i\)-variables need to be made. The decisions for \(X^e\) and Y need to be made in such a way that the expected costs in the second stage are minimized. In other words, the two-stage stochastic problem can be formulated as follows:

$$\begin{aligned} \min \qquad \sum _{t =1}^{n^e}\sum _{j=0}^{m} c^e_{tj}X^e_{tj} + {\mathbb {E}}\left[ Q(X^e,Y,\varPhi )\right] \end{aligned}$$
(13)

subject to:

$$\begin{aligned}&{\text {Constraints: }} (3), (5), (7), (9),(11) \, {\text {and}}\,(12), \end{aligned}$$
(14)

in which \(Q(X^e,Y, \phi )\) is defined as:

$$\begin{aligned} Q(X^e,Y,\phi ):= \min \sum _{t =1}^{n^i}\sum _{j=0}^{m} c^i_{tj}X^i_{tj} + \sum _{j=1}^m\sum _{k = 1}^l q_{jk}Z_{jk} \end{aligned}$$
(15)

subject to:

$$\begin{aligned}&\sum _{j=1}^m\sum _{t=1}^{n^i} X^i_{tj} +\sum _{j=1}^m\sum _{t=1}^{n^e} X^e_{tj}-Z_{jk}\le \phi _{jk} \qquad j =1, \ldots , m \quad k =1, \ldots , l \end{aligned}$$
(16)
$$\begin{aligned}&Z_{jk} \in {\mathbb {N}}_0 \qquad j =1, \ldots , m \quad k =1, \ldots , l \end{aligned}$$
(17)
$$\begin{aligned}&{\text {Constraints: }} (2), (4),(8)\, {\text {and}}\, (10). \end{aligned}$$
(18)

The problem in (13)–(14) is the first stage problem and the second stage problem is given in (15)–(18). The function \(Q(X^e,Y, \phi )\) represents the optimal costs for the import containers given the first-stage decision variables \(X^e\) and Y, and realization \(\phi\). In the formulation (15)–(18), the variable \(Z_{jk}\) indicates the number of export containers on barge j that could not be unloaded at terminal k because the maximum number of moves has been reached at that terminal. For these containers, a recourse action is required. The second stage is a deterministic assignment problem that can be solved relatively easily to optimality. Nevertheless, computing \({\mathbb {E}}[Q(X^e,Y,\varPhi )]\) is computationally very expensive because the number of realizations can be enormous and it is not possible to derive a closed-form expression for \(Q(X^e,Y,\phi )\). Therefore, in the following section, we describe three different methods to solve this two-stage stochastic problem.

5 Solution methods

In this section, we describe five different methods to solve the problem (13)–(18). The first method that we use to solve our problem is the SAA method. The solution of this method is based on the realizations of the random variable \(\varPhi _{jk}\). The solution obtained by the SAA method converges to the optimal solution if the number of realizations increases. Nevertheless, the running time of this method also increases with a growing number of realizations. If all information is known long before the planning process has to finish, then it is not problematic if the running time of the SAA method. However, in many situations, the required information is only available a short amount of time before the planning has to be made. Hence, also a faster method than the SAA method is required. We use Stochastic Programming (SP) to derive such a method and thus, we refer to this method as the SP-based method. The idea behind this method is the following: if we simplify our problem enough, it can be solved to optimality using stochastic programming. In the SP-based method, the characteristics of this optimal solution are used to derive a fast solution for the original problem. Although the solution SP-based method is based on the solution for a simplified problem, the idea is that this simplified problem resembles enough of the original problem, such that the SP-based solution is close to the optimal solution. To compare the performance of the SP-based method, we also describe in this section three fast methods that have been proposed before in the literature. These methods are applied with success to other stochastic assignment problems, but are less tailored towards our problem.

In all methods, the stochastic variable \(\varPhi _{jk}\) will be replaced by a deterministic value, such that we again can formulate the problem as an ILP as illustrated in Sect. 4.1. If the variables \(Y_{jk}\) are fixed for \(j=1, \ldots , m\) and \(k= 1, \ldots , l\), then this ILP is equivalent to a variant of the Generalized Assignment Problem (GAP). In Benders and van Nunen (1983), it has been shown that for a linear relaxation of the GAP, the number of fractional variables is smaller than the number of capacity constraints that hold with equality. For each barge, there are at most \(N_j\) capacity constraints of the type of constraint (6), and one capacity constraint of type (4) and one of type (5). Hence, the number of fractional variables in the LP-relaxation will be relatively small. In Zweers et al. (2019), it is shown that for a similar problem the objective function with only fractional \(X^e\)-variables, \(X^i\)-variables, and Z-variables is close to the integral objective function. Moreover, the computation time for the former problem is much shorter than the latter. Therefore, in all solutions methods, we drop the integrality constraints for the \(X^e\)-variables, \(X^i\)-variables, and Z-variables.

Once the first-stage decision variables are fixed and the realization of \(\phi\) is known, determining the optimal values for the second-stage variables is a simple assignment problem. Therefore, we will only describe methods to find solutions for the first-stage variables. Both the terminals to visit by each barge and the transport decision of each export container are first-stage decision variables. Nevertheless, from the decision which export containers to transport on which barge follows directly which terminals are visited by a barge. Hence, only reporting the transportation plan of the export containers is sufficient as a first-stage decision. The remainder of this section is organized as follows: in Sect. 5.1, we give the SAA method and the SP-based method is given in Sect. 5.2. Afterward, in Sect. 5.3, we give three methods that perform well in the study of Toktas et al. (2006).

5.1 Sample average approximation method

The main idea of SAA is that in many stochastic optimization problems the expectation in the objective function is hard to compute exactly, but that for a given realization of the random variable the problem is easier to solve. This is also the case in our problem, in which the ILP given in (1)–(12) is (relatively) easy to solve and the stochastic problem (13)–(18) is hard to solve. The goal of the SAA method is to approximate \({\mathbb {E}}\left[ Q(X^e,Y,\varPhi )\right]\) by using a set of N vectors of realizations \(\bar{\phi }_N := (\phi ^1,\phi ^2, \ldots , \phi ^N)\). We refer to a vector of realizations \(\phi ^n\) as a scenario and use the index n to refer to a scenario. A scenario \(\phi ^n\) consists of one realization \(\phi _{jk}\) of \(\varPhi _{jk}\) for every \(j=1, \ldots , m\) and \(k= 1, \ldots , l\). For each scenario, the value for \(Q(X^e,Y, \phi ^n)\) can be evaluated, so \({\mathbb {E}}\left[ Q(X^e,Y,\varPhi )\right]\) can be approximated by \(\frac{1}{N} \sum _{n=1}^N Q(X^e,Y,\phi ^n)\).

In the SAA-method, for each scenario n, the assignment of the import containers depends on the value \(\phi ^n\). Hence, a decision variable \(X^{in}_{tj}\) is needed that indicates the number of import containers of type t transported on barge j for scenario n. On top of that, a variable \(Z_{jk}^n\) is needed to indicate the number of containers on barge j for terminal k for which a recourse action is required in scenario n. Nevertheless, the export containers and the set of visited terminals should be the same for every realization, thus the variables \(X^e_{tj}\) and \(Y_{jk}\) do not depend on the realization of \(\phi ^n\). Using these variables, the SAA method for our problem can be formulated as the following ILP:

$$\begin{aligned} \min \sum _{t=1}^{n^e} \sum _{j=0}^m c^e_{tj}X^{e}_{tj}+ \frac{1}{N} \sum _{n=1}^N \left( \sum _{t =1}^{n^i} \sum _{j=0}^m c^i_{tj}X^{in}_{tj} + \sum _{j =1}^m\sum _{k=1}^lq_{jk}Z^{n}_{jk}\right) \end{aligned}$$
(19)

subject to:

$$\begin{aligned}&\sum _{j=0}^mX^{in}_{tj} = d^i_t \qquad t = 1, \ldots n^i \quad n =1 , \ldots , N \end{aligned}$$
(20)
$$\begin{aligned}&\sum _{t =1}^{n^i} w^i_tX^{in}_{tj}\le u_j \qquad j = 1, \ldots , m\quad n =1, \ldots , N\end{aligned}$$
(21)
$$\begin{aligned}&X^{in}_{tj} \le d^i_t Y_{jk} \qquad \qquad j= 1, \ldots , m \qquad k = 1, \ldots , l \quad t = 1, \ldots , n^i \quad n =1, \ldots , N\end{aligned}$$
(22)
$$\begin{aligned}&\sum _{t \in {\mathcal {I}}(k)}X^{in}_{tj} + \sum _{t \in {\mathcal {E}}(k) } X^{e}_{tj} - Z^{n}_{jk} \le \phi _{jk}^n \qquad j =1, \ldots , m \quad k =1, \ldots , l\quad n=1, \ldots , N \end{aligned}$$
(23)
$$\begin{aligned}&{\text {Constraints:}}\, (3),(5), (7),(9), (11)\end{aligned}$$
(24)
$$\begin{aligned}&X^{in}_{tj}\in {\mathbb {N}}_0 \qquad t =1, \ldots , n^i \quad j=1, \ldots , m\quad n = 1, \ldots , N \end{aligned}$$
(25)
$$\begin{aligned}&Z^{n}_{jk} \in {\mathbb {N}}_0\qquad j= 1, \ldots , m \quad k = 1, \ldots , m. \end{aligned}$$
(26)

Let us denote the optimal solution of the problem (19)–(26) by the quadruple \((X^{e^*},X^{in^*}, Y^*, Z^{n^*})\). It is important to realize that this solution is optimal for the objective function (19), but it is not necessarily the optimal solution for the objective function (13). We will refer to the optimal first stage decisions of the problem (19)–(26) as \(\hat{x}_N\) and to the value of the objective function (19) under \(\hat{x}_N\) and \(\bar{\phi }_N\) as \(\hat{z}_N(\hat{x}_N,\bar{\phi }_N)\). Let us denote the solution produced by the SAA method by \(x^{SAA}\). One could take for \(x^{SAA}\) the value of \(\hat{x}_N\), but in order to get a better solution than \(\hat{x}_N\) the problem (19)–(26) can be solved multiple times. We could solve the problem (19)–(26) for M different sets of N samples, which we denote by \(\bar{\phi }_N^1,\bar{\phi }_N^2, \ldots , \bar{\phi }_N^M\). This will result in M solutions, \(\hat{x}_N^1, \hat{x}_N^2, \ldots , \hat{x}_N^M\). From these M solutions, only one can be executed in practice. To select the best of these M solutions, we use a different sample \(\bar{\phi }_{N'}\) consisting of \(N'\) scenarios to calculate \(\hat{z}_{N'}(\hat{x}_N^m,\bar{\phi }_{N'})\), for \(m = 1, 2, \ldots ,M\). One could see the original M samples of size N as a training set and this new set with a sample of size \(N'\) as a validation set. It is common practice in SAA to select the solution \(x^{SAA}\) that has the lowest objective function under the validation set. In other words,

$$\begin{aligned} x^{SAA} := {{\,\mathrm{arg min}\,}}\left\{ \hat{z}_{N'}\left( \hat{x}_N,\bar{\phi }_{N'}\right) : \hat{x}_{N} \in \left\{ \hat{x}_N^1, \hat{x}_N^2, \ldots , \hat{x}_N^M\right\} \right\} . \end{aligned}$$

An advantage of the SAA method is that can also be used to calculate a lower bound and an upper bound for the optimal solution and these bound can then be used to calculate an optimality gap. We first show how a lower bound for the optimal solution can be obtained. As the solution \(\hat{x}_N^m\) was determined by using the sample \(\bar{\phi }_N^m\), the value \(\hat{z}_N(\hat{x}_N^m,\bar{\phi }_N^m)\) is negatively biased estimator for the optimal value of the real objective function (13), which we denote by \(z^*\). In other words, \(\hat{z}_N(\hat{x}_N^m,\bar{\phi }_N^m)\) is a lower bound for the optimal solution and the following inequality holds (Mak et al. 1999):

$$\begin{aligned} {\mathbb {E}}\left[ \hat{z}_N(\hat{x}_N^m, \bar{\phi }_N^m)\right] \le z^*. \end{aligned}$$

If we take the average over all M values for \(\hat{z}_N(\hat{x}_N^m, \bar{\phi }_N^m)\) for \(m= 1, \ldots , M\), we get \(\bar{t}^M_N:= \frac{1}{M}\sum _{m=1}^M \hat{z}_N(\hat{x}_N^m, \bar{\phi }_N^m)\), which is in expectation equal to \(\hat{z}_N(\hat{x}_N^m, \bar{\phi }_N^m)\). Therefore, we use \(\bar{t}^M_N\) as an estimation for the lower bound of the objective function of the problem in (13)–(18).

Let us now focus on the upper bound for the optimal solution. As the solution \(\hat{x}_N^m\) is a feasible solution for (13)–(18), the value \(\hat{z}_{N'}(\hat{x}_N^m,\bar{\phi }_{N'})\) is an upper bound for the optimal objective function. Similarly to \(\bar{t}^M_N\), we define \(\bar{v}^M_{N'}\) as the average over all M objective functions given the scenarios in the validation set, i.e, \(\bar{v}^M_{N'} := \frac{1}{M}\sum _{m=1}^M \hat{z}_{N'}(\hat{x}_N^m,\bar{\phi }_{N'})\). All in all, we have that \({\mathbb {E}}[\bar{v}^M_{N'}]\ge z^*\). Therefore, an estimation for the optimality gap is given by \(\bar{v}_{N'}^M - \bar{t}_N^M\). Since both the \(\bar{t}_{N'}^M\) and \(\bar{v}_{N}^M\) are estimates, the Central Limit Theorem can be used to take accuracy into account. An estimator of the optimality gap that takes accuracy into account is (Kleywegt et al. 2001):

$$\begin{aligned} \bar{v}_{N'}^M - \bar{t}_{N'}^M + z_{\alpha }\frac{\bar{S}^M}{\sqrt{M}}, \end{aligned}$$
(27)

in which \(z_\alpha :=\varPhi ^{-1}(1-\alpha )\), where \(\varPhi (z)\) is the cumulative distribution function of the standard normal distribution and \(\frac{\bar{S}^M}{\sqrt{M}}\) is given by:

$$\begin{aligned} \frac{\bar{S}^M}{\sqrt{M}}:= \sqrt{\frac{1}{M(M-1)} \sum _{m=1}^M\left( \left( \hat{z}_{N'}(\hat{x}_N^m,\bar{\phi }_{N'})-\hat{z}_{N}(\hat{x}_N^m,\bar{\phi }^m_{N})\right) -\left( \bar{v}_{N'}^M - \bar{t}_N^M\right) \right) ^2}. \end{aligned}$$
(28)

All in all, the SAA method consists of two phases: a training and a validation phase. In the training phase, M different problems with N scenarios are solved and in the validation phase, the quality of these solutions is investigated for \(N'\) different scenarios. It should be noted that the validation phase is less computationally intensive than the training phase. In the training phase, a decision for both the first-stage and the second-stage variables needs to be made. To do so, all scenarios in the training set have to be considered simultaneously. On the other hand, in the validation phase, the first-stage decision variables and the realization of the stochastic variables are fixed and only a decision for the second-stage decision variables is needed. Consequently, each of the \(N'\) scenarios can be investigated independently.

5.2 Stochastic programming method

In this section, we explain how the SP-based method is derived. We explain in Sect. 5.2.1 how to get an optimal solution for the simplified problem and in Sect. 5.2.2 this solution will be used to get a solution for the problem (13)–(18).

5.2.1 Optimal solution simplified problem

Let us consider a situation with only one barge and one terminal. The number of export containers that has to be delivered to that terminal is denoted by \(d^e\), and the number of import containers that need to be collected from that terminal is denoted by \(d^i\). Moreover, we assume that the barge transportation costs for all import and export containers are \(c^i\) and \(c^e\), respectively. Furthermore, the costs of transporting an import or an export container with a truck are denoted by \(c^i_t\) and \(c^e_t\), respectively. The number of moves at the terminal is the stochastic variable \(\varPhi\), distributed according to the cumulative distribution function \(F(\cdot )\). Finally, the costs of a recourse action for an export container are q.

Let us use x to denote the decision variable on the number of export containers on the barge. We assume that this decision variable is continuous. From the decision x, the number of export containers on a truck, \(d^e-x\), follows directly. Furthermore, we assume that the variable \(\varPhi\) is that restrictive that we can ignore the capacity of a barge. If x export containers are transported and \(\phi\) is the realization of the number of moves, then no import containers can be shipped by barge if \(x\ge \phi\). In case \(x<\phi\), at most \(\phi -x\) import containers can be shipped per barge. Hence, given x, the expected number of import containers per barge is equal to \({\mathbb {E}}[\min \{d^i,\max \{0,\varPhi -x\}\}]\). Furthermore, the expected number of export containers for which we need to perform a recourse action is given by \({\mathbb {E}}[\max \{0, x-\varPhi \}]\). In Table 11, that is given in “Appendix 1”, we give an overview of all notation that is described above.

For \(0\le x\le d^e\), the expected total costs T(x) can be given by the following expression:

$$\begin{aligned} T(x)&= c^ex + c^e_t(d^e-x) + c^i {\mathbb {E}}[\min \{d^i,\max \{0,\varPhi -x\}\}] \nonumber \\&\quad + c^i_t (d^i - {\mathbb {E}}[\min \{d^i,\max \{0,\varPhi -x\}\}])+q {\mathbb {E}}[\max \{0, x - \varPhi \}]. \end{aligned}$$
(29)

Equation (29) consists of five terms: the first two terms are deterministic because they represent, respectively, the barge and truck costs for the export containers. The third part of the sum gives the expected barge costs for the import containers and the expected truck costs for the import containers are given by the fourth part. Finally, the last term in Eq. (29) corresponds to the recourse costs for the export containers. In “Appendix 2”, we show that Eq. (29) can be simplified into:

$$\begin{aligned} T(x)&= c^ex + c^e_t(d^e-x) + c^i \left( d^i - \int _{x}^{d^i+x} F(t) dt \right) \\&\quad + c^i_t \left( \int _{x}^{d^i+x} F(t) dt \right) + q \int _{0}^{x}F(t)dt. \end{aligned}$$

The optimization problem can now be given as:

$$\begin{aligned}&\min _x \qquad T(x) \end{aligned}$$
(30)
$$\begin{aligned}&{\hbox {subject to:}}\quad g_1(x) =x -d^e\le 0\end{aligned}$$
(31)
$$\begin{aligned}&g_2(x) = - \,x\le 0. \end{aligned}$$
(32)

This is a problem with simple recourse and thus the expected recourse costs are given by a function that is convex in x (Birge and Louveaux 2011). In other words, we can use the Karush–Kuhn–Tucker (KKT) conditions, to calculate the optimal number of export containers on a barge. The KKT conditions are given by:

$$\begin{aligned}&\frac{dT(x)}{dx} = -\,\lambda _1 \frac{d g_1(x)}{dx}-\lambda _2 \frac{d g_2(x)}{dx}\\&\lambda _1 g_1(x) = 0\\&\lambda _2 g_2(x) = 0\\&\lambda _1,\lambda _2 \ge 0\\&g_1(x),g_2(x)\le 0, \end{aligned}$$

which can be reformulated as:

$$\begin{aligned}&c^e -c_t^e +\left( c_t^i-c^i\right) F\left( d^i+x\right) + \left( c^i -c^i_t +q \right) F(x) =-\lambda _1+\lambda _2\nonumber \\&\lambda _1 \left( x-d^e\right) = 0\nonumber \\&\lambda _2 (-x) = 0\nonumber \\&\lambda \ge 0\nonumber \\&(x-d^e) \le 0\nonumber \\&-x\le 0. \end{aligned}$$
(33)

The system of equalities (33) above has four different types of solutions. Namely, (1) \(\lambda _1>0\) and \(\lambda _2>0\), (2) \(\lambda _1 = 0\) and \(\lambda _2>0\), (3) \(\lambda _1>0\) and \(\lambda _2=0\), and (iv) \(\lambda _1=\lambda _2=0\). A feasible solution for (33) with both \(\lambda _1>0\) and \(\lambda _2>0\) is only possible in the trivial case in which \(d^e=0\). If \(\lambda _1=0\) and \(\lambda _2>0\), then the third equality gives us that \(x^* =0\). Equivalently, if \(\lambda _1>0\) and \(\lambda _2 =0\), then because of the second equality we know that \(x^*\) should be \(d^e\). Finally, if \(\lambda _1=\lambda _2=0\), then the first equality implies that the optimal \(x^*\) should satisfy:

$$\begin{aligned} c^e -c_t^e +\left( c_t^i-c^i\right) F\left( d^i+x^*\right) + \left( c^i -c^i_t +q \right) F(x^*) =0. \end{aligned}$$
(34)

Although it is not possible to derive an analytical expression for \(x^*\), the equality of Eq. (34) can be solved using, for instance, the Newton–Raphson method. Since the cumulative distribution function \(F(\cdot )\) is a non-decreasing function, we know that the expression in Eq. (34), is also non-increasing in x. Consequently, there is at most one value \(x^*\ge 0\) for which the equality in Eq. (34) holds. Hence, there exists an easy method to find the solution for the system of equalities in (33). If

$$\begin{aligned} \frac{dT(0)}{dx}&=c^e -c_t^e +\left( c_t^i-c^i\right) F\left( d^i\right) + \left( c^i -c^i_t +q \right) F(0) \\&=c^e -c_t^e +\left( c_t^i-c^i\right) F\left( d^i\right) \ge 0, \end{aligned}$$

then \(x^*=0\). Otherwise, the solution to the equality in Eq. (34) has to be calculated. If that solution is smaller than \(d^e\), then it is the optimal solution, otherwise \(x^* = d^e\).

5.2.2 Heuristic stochastic programming method

In the previous section, we have given the optimal solution for a simplified problem. There are two aspects in which the problem of the previous section was easier than the problem described in Sect. 4. First of all, the assumption that there is only a single barge visiting a single terminal and second, that the transportation costs are the same for all containers. Nevertheless, the solution to (30)–(32) can be used as a basis to find a solution for the problem described in Sect. 4. The goal of the SP-based method is to find a constraint for the number of export containers on a barge for a specific terminal.

Let us first focus on terminal k and barge j. Let us assume that we have an estimate for the number of import and export containers on terminal k that could be transported on barge j (i.e., \(\delta ^i_{jk}\) and \(\delta ^e_{jk}\)), an estimate for the costs of shipping an import and an export container from terminal k on barge j (i.e., \(\gamma ^i_{jk}\) and \(\gamma ^e_{jk}\)) and an estimate for the transportation costs per truck for the import (i.e., \(\gamma ^i_{0_{k}}\)) and export containers (i.e., \(\gamma ^e_{0_{k}}\)) on terminal k. Furthermore, we assume that barge j is the only barge visiting terminal k. Hence, if containers from terminal k are not on barge j, they are transported by truck. Finally, we denote the sum of all export containers for terminal k that are on barge j by \(X^e_{jk}\). Using these assumptions, the total transportation costs for all containers on terminal k can be approximated in the same way as in the optimization problem (30)–(32). Therefore, we can use an equality similar to the equality in Eq. (34), to find a solution for \(X^e_{jk}\), namely

$$\begin{aligned} \gamma _{jk}^e - \gamma _{0_k}^e + \left( \gamma _{0_k}^i - \gamma _{jk}^i\right) F_{jk}(\delta ^i_{jk}+X^e_{jk}) + (\gamma _{jk}^i - \gamma _{0_k}^i + q_{jk}) F_{jk}(X^e_{jk})=0. \end{aligned}$$
(35)

Besides, deciding how many export containers for terminal k will be on barge j, also the specific containers that will be loaded on this barge have to be selected. These two decisions cannot be taken independently. If only one barge is visiting a terminal and only a single export container for a terminal is transported on that barge j, then it is easy to determine which container is the best to transport on barge j. If a container is not transported on barge j, then it is transported by truck. The container that benefits most from being transported on barge j instead of the truck is the container for which the difference between the truck transportation costs and the barge transportation costs is the largest, i.e., \(\tau _j = {{\,\mathrm{arg max}\,}}_{t \in {\mathcal {E}}(k)}\{c^e_{t0}-c^e_{tj}\}\). Therefore, if \(X^e_{jk}=1\), the estimates \(\gamma _{jk}^e\) and \(\gamma _{0_k}^e\) are not needed and can be replaced by \(c^e_{\tau _j j}\) and \(c^e_{\tau _j 0}\). Similarly, if \(X^e_{jk}=2\), the container for which the difference between the truck and barge transportation costs is second largest is chosen. Consequently, it becomes clear that the estimators \(\gamma _{jk}^e\) or \(\gamma _{0_k}^e\) are not needed but that the real value \(c^e_{t0}\) and \(c^e_{tj}\) can be used in Eq. (35).

figure a

Similar to the solution method described in Sect. 5.2.1, we can increase the number of export containers for a terminal that is assigned to a barge until the equality in (35) is satisfied. In Algorithm 1, this procedure is formalized. This algorithm is used to derive a constraint for the maximum number of export containers \(M_{jk}\) that can be transported with barge j to terminal k. When determining \(M_{jk}\), we assumed that barge j was the only barge visiting terminal k. Consequently, all containers can be selected for transportation on barge j. It might be that a container on terminal k has extremely high truck transportation costs, and is therefore used in calculating the \(M_{jk}\)-value for every barge j. Nevertheless, that container can only be transported on one barge. Hence, even if there are multiple barges visiting terminal k, the number of export containers on barge j to terminal k is always less than \(M_{jk}\). With the constraint from Algorithm 1, we get the following deterministic ILP formulation:

$$\begin{aligned} \min \quad \sum _{t =1}^{n^i} \sum _{j=0}^m c^i_{tj}X^i_{tj} + \sum _{t=1}^{n^e} \sum _{j=0}^m c^e_{tj}X^e_{tj} \end{aligned}$$
(36)

subject to:

$$\begin{aligned}&\sum _{t\in {\mathcal {E}}(k)} X^e_{tj}\le M_{jk} \qquad j =1, \ldots , m \quad k =1, \ldots , l \end{aligned}$$
(37)
$$\begin{aligned}&\sum _{t\in {\mathcal {I}}(k)}^{n^i} X^i_{tj} +\sum _{t\in {\mathcal {E}}(k)} X^e_{tj}\le {\mathbb {E}}[\varPhi _{jk}] \qquad j= 1, \ldots , m \quad k = 1, \ldots , l \end{aligned}$$
(38)
$$\begin{aligned}&{\text {Constraints }} (2){-}(5)\,{\text {and}}\,(7){-}(12). \end{aligned}$$
(39)

If we only add constraint (37) to the ILP, then the import containers would not be restricted by the number of moves at a terminal. As a consequence, the solution of the ILP would probably overestimate the number of import containers that can be shipped on a barge. Hence, constraint (38) is added to limit the number of import containers shipped on a barge.

5.3 Other methods

In this section, we give three methods that are given in Toktas et al. (2006) to solve assignment problems with stochastic capacity constraints. We explain how these methods can be used to solve the problem (13)–(18). In Sect. 5.3.1, we give the expectational method, the risk-averse trimmed mean method is explained in Sect. 5.3.2 and the comparative performance evaluation method is described in Sect. 5.3.3.

5.3.1 Expectational method

In this section, the expectational method is described. In this method, we replace the stochastic variable \(\varPhi _{jk}\) by its expectation. So the first-stage objective function is modified into:

$$\begin{aligned} \min \, \sum _{t =1}^{n^e}\sum _{j=0}^{m} c^e_{tj}X^e_{tj} + Q(X^e,Y,{\mathbb {E}}\left[ \varPhi \right] ). \end{aligned}$$

Since the expectation is a deterministic value, the problem can be formulated as an ILP, similar to the ILP given in (1)–(12). The only difference is that constraint (6) is replaced by:

$$\begin{aligned} \sum _{t\in {\mathcal {I}}(k)} X^i_{tj} +\sum _{t\in {\mathcal {E}}(k)} X^e_{tj}\le {\mathbb {E}}[\varPhi _{jk}] \qquad j =1, \ldots , m \quad k =1, \ldots , l. \end{aligned}$$
(40)

This expectational method is rather naive since it only uses the expectation and no other information of the distribution. Nevertheless, it is a good benchmark to see what one can gain by incorporating more knowledge of the distribution.

5.3.2 Risk-averse trimmed mean method

The risk-averse trimmed mean method is almost identical to the expecataional method. The only difference is that instead of the expectation a more conservative estimate for the number of moves is used, which is called the risk-averse trimmed mean. This risk-averse trimmed mean is defined as follows:

$$\begin{aligned} {\mathbb {A}}_\rho \left[ \varPhi _{jk}\right] := {\mathbb {E}}\left[ \varPhi _{jk}|F(\phi _{jk})\le \rho \right] . \end{aligned}$$

The interpretation of the risk-averse trimmed mean is that the expectation of the values less than or equal to the \(\rho\)-quantile is taken. In the risk-averse trimmed mean method, the constraint (6) in the ILP (1)-12) is replaced by:

$$\begin{aligned} \sum _{t\in {\mathcal {I}}(k)} X^i_{tj} +\sum _{t\in {\mathcal {E}}(k)} X^e_{tj}\le {\mathbb {A}}_\rho [\varPhi _{jk}] \qquad j =1, \ldots , m \quad k =1, \ldots , l. \end{aligned}$$

The idea behind the risk-averse trimmed mean method is that taking a more conservative estimate for the number of moves than the expectation, results in fewer recourse costs. The parameter \(\rho\) can be used to indicate the risk one is willing to take. If \(\rho = 0\), then \({\mathbb {A}}_{\rho }=0\) and if \(\rho =1\), then \({\mathbb {A}}_{\rho }[\varPhi _{jk}]={\mathbb {E}}[\varPhi _{jk}]\).

5.3.3 Comparative performance evaluation method

The Comparative Permorance Evaluation (CPE) method is, similar to the SAA method, a sample based method. Let us have a vector of N scenarios \(\bar{\phi }_N:=\left( \phi ^1, \phi ^2, \ldots , \phi ^N\right)\). For a single scenario \(\phi ^n\), the deterministic problem (1)–(12) can (relatively) easily be solved, which gives a solution \((X^{en},Y^n)\). Using that solution, the value \({\mathbb {E}}\left[ Q(X^{en},Y^n,\varPhi )\right]\) can be estimated by: \(\frac{1}{N}\sum _{n'=1}^NQ\left( X^{en},Y^n,\phi ^{n'}\right)\). The idea of the CPE method is to select out of the N solutions \((X^{en}, Y^n)\) that are based on a single scenario, the solution that has the lowest estimated costs over all scenario in \(\bar{\phi }_N\). In other words,

$$\begin{aligned} x^{CPE} :={\mathop {{\hbox {arg min}}}\limits _{n = 1, \ldots , N}} \left\{ \sum _{t=1}^{n^e}\sum _{j=0}^m c_{tj}^eX^{en}_{{tj}} +\frac{1}{N}\sum _{n'=1}^NQ\left( X^{en},Y^n,\phi ^{n'}\right) \right\} . \end{aligned}$$

6 Numerical results

The quality of the solution methods described in the previous section is investigated in this section using numerical experiments. In Sect. 6.1, we describe three different terminal types that are used in the experiments. A terminal type reflects the distribution for the number of moves at a deep-sea terminal. Next, in Sect. 6.2, small instances are solved to investigate the convergence rate of the SAA method to the optimal solution and the quality of the solution produced by the other methods. Large instances are solved in Sect. 6.3. In this section, we focus on the scalability of the SAA method and compare the outcome of the best parameters for the SAA with the SP-based method and the methods discussed in Sect. 5.3. Finally, in Sect. 6.4, we look at the quality of the methods if the number of moves is correlated between barges.

6.1 Terminal types

In practice, the congestion is not the same at each deep-sea terminal. Moreover, the way the terminals deal with congestion also differs per terminal. Therefore, we consider three different types of terminals: predictable, unpredictable and open-closed terminals. At a predictable terminal, the available number of moves for each barge does not vary much. This terminal might not be really congested or it is conservative in the number of moves it allocates to barges and with that, it can always have a rather constant number of moves. We model the predictable terminal with a Poisson distribution with expectation \(\lambda\). On the other hand, at an unpredictable terminal, the number of moves is sometimes very limited and sometimes very large. Based on a limited amount of data for the moves at container terminals, we conjecture that at least some terminals are of this kind. We model these terminals with a geometric distribution with parameter p and the support \(k = 0,1, \ldots\).

Finally, we consider the open-closed terminals, which are either closed or they are open. If the terminal is closed not a single container can be loaded or unloaded, but in case the terminal is open, then its moves are rather predictable. We model the open-closed terminals with a bimodal distribution that takes with probability \(\frac{1}{3}\) the value 0 and with probability \(\frac{2}{3}\), the value is distributed according to a Poisson distribution with expectation \(\mu\). The difference between an unpredictable terminal and an open-closed terminal is that the first one is accepting barges also if the terminal has very limited time available for loading and unloading. On the other hand, the open-closed terminal accepts a barge visit only if it has enough time to handle a large number of containers.

To investigate the consequences of the level of variability for the number of moves, all distributions will have the same expectation, namely \({\mathbb {E}}[\varPhi ]\). For the three different distributions, the variance of the number of moves can be expressed in terms of this \({\mathbb {E}}[\varPhi ]\). The variance of the Poisson distribution is \({\mathbb {E}}[\varPhi ]\), for the bimodal distribution it is \(\frac{1}{2}\left( {\mathbb {E}}[\varPhi ]\right) ^2+{\mathbb {E}}[\varPhi ]\) and the variance of the geometric distribution is \(\left( {\mathbb {E}}[\varPhi ]\right) ^2 + {\mathbb {E}}[\varPhi ]\). Thus, the variance of the geometric distribution is about twice the variance of the bimodal distribution.

All these three distributions are discrete, because not a fractional number of containers can be loaded. Nevertheless, the methods described in Sect. 5 also apply to continuous probability distributions. Furthermore, in this paper, we assume that all the moves for all terminals and all barges in a port have the same distribution. However, the three solution methods described in Sect. 5 can all be applied to an instance with unique distributions for every barge and terminal combination.

6.2 Small instances

First, we solve ten relatively small instances to investigate the convergence of the SAA solution to the optimal solution and to evaluate the quality of the other methods. Both these small instances and the large instances that are solved in Sect. 6.3 are based on real data from the inland terminal at the port of Amsterdam. The small instances consist of 50 export and 50 import containers that can be transported on two barges with a capacity of 50 TEU. Every container has a size of either 1 or 2 TEU and the costs for transporting an import or export container on a barge are uniformly distributed between 25 and 100, whereas the costs of transporting a container per truck are uniformly distributed between 150 and 200. The costs for a recourse action for an export container are 300. The containers are distributed over five different deep-sea terminals and each barge is allowed to visit three of them. At each terminal, the expected number of moves is 10 for all terminal types.

The number of scenarios in the training set, denoted by N, is varied to evaluate the consequence of the size of the training set. We have used twenty different runs of the SAA method, i.e., \(M=20\), and we have used \(N'=5000\), in other words, the validation set consists of 5000 samples. Out of the twenty different solutions we select the solution with the lowest objective function for the validation set. This solution is a negatively biased estimate for the value of the SAA solution and thus, another validation set of 5000 samples is used to obtain an unbiased estimate for the objective function. In Table 1, the average of this objective function and the optimality gap for the SAA solution are given for the three different terminal types and different numbers of scenarios in the training set (N). The optimality gap in this table is calculated by dividing the outcome of Eq. (27) by the lower bound of the optimal solution \(\bar{t}^M_{N'}\). We use for this optimality gap and all remaining optimality gaps \(\alpha =0.025\). For each combination of N and terminal type, the average over the optimality gaps of the ten instances is given.

The most important conclusion from Table 1 is that for \(N=1000\), the average optimality gap for all three terminal types is small, namely 0.3% for the predictable terminals, 0.2% for the open-closed terminals and 0.0% for the unpredictable terminals. Hence, if 1000 scenarios are used, the SAA method can find a solution that is (close to) optimal. Furthermore, it can be seen that for \(N=10\), the predictable terminal has the smallest optimality gap. This observation can be explained by the fact that the Poisson distribution has the smallest variance. Hence, the scenarios in the training and validation are more similar than for the unpredictable and open-closed terminals. For those two terminal types, a sharp decrease in the optimality gap is seen for an increasing size of the training set.

Table 1 Average objective function and optimality gap of the SAA solution for the small instances with increasing numbers of scenarios

It becomes clear from Table 1 that the predictable terminal has significantly lower costs than the two other terminal types. In other words, the inland container terminal would greatly benefit from the deep-sea terminals having a more reliable number of moves. On the other hand, the difference in the objective between the unpredictable and the open-closed terminal types is not that large. Finally, for all three terminal types, the objective function for \(N=1000\) is about 6% better than the solution found for \(N=10\). The running time for the training phase of the SAA increases if N gets larger. Nevertheless, for the small instances even for \(N=1000\), the training phase was able to finish within 10 min. The running time of the validation phase only depends on M and \(N'\) and is independent of N. If \(M=20\) and \(N'=5000\), the running time of the validation phase is about 12 min.

Now that we have shown that the solution of the SAA method converges to the optimal solution, it is also interesting to look into the quality of the solutions for the SP-based method and the other three methods. For the SP-based method, three different types of parameters need to be defined, namely \(\gamma ^i_{jk}\), \(\gamma ^i_{0k}\) and \(\delta ^i_{jk}\). For setting the value of \(\delta ^i_{jk}\), we count how many import containers are at terminal k available for transportation by barge j. Thereafter, we divide this value by the average number of times a terminal is visited by a barge. In these small instances, the two barges visit both three terminals. Since there are in total five terminals, the average number of times a terminal is visited by a barge is \(\frac{6}{5}\). It is important to realize that in Algorithm 1 only the difference between \(\gamma _{jk}^i\) and \(\gamma _{0k}^i\) is used. We have decided to perform some numerical experiments to decide upon the best parameter setting for \(\gamma _{0k}^i-\gamma _{jk}^i\). In Table 2, the average objective function and, between brackets, its standard deviation are given for all three terminal types and \(\gamma _{0k}^i -\gamma _{jk}^i\) is 50, 100, and 150. We see that for all three terminal types the best results are obtained if \(\gamma _{0k}^i -\gamma _{jk}^i\) is set to 50. So in the remaining of this section, we use these parameter settings for the SP-based methods. For the risk-averse trimmed mean method, we have to decide on the value for the parameter \(\rho\). In Toktas et al. (2006), it is shown that the best solutions are obtained if \(\rho = 0.8\), so we have decided to also use this value.

Table 2 Average objective function and its standard deviation for the 10 small instances for the SP-based method using different settings for \(\gamma _{0k}^i-\gamma _{jk}^i\)
Table 3 Average objective function and its standard deviation for the 10 small instances for the five solution methods and the three different terminal types

In Table 3, the objective function for the expectational method (EXP), the risk-averse trimmed mean method (RAT), the CPE method and the SP-based method are compared with the SAA solution. For all methods, we give the average objective function over the ten instances and the standard deviation of the objective function is given between brackets. The SAA solution in this table is the solution obtained by using \(N=1000\) and we have seen above that this solution is close to the optimal solution. The first conclusion to draw from Table 3 is that the SP-based method is the method that produces the best solutions and the worst solutions are from the CPE method. Moreover, the difference between the objective functions for the expectational method and the SAA method is significant and thus taking the uncertainty for the number of moves is beneficial. For the predictable terminals, the risk-averse trimmed mean is too conservative and it is better to choose the expectation, whereas for the unpredictable terminals the variance of the number of moves is higher and the risk-averse trimmed mean is performing better than the expectation. For the open-closed terminals, the difference between the objective function of the expectational method and the SAA solution is with 7.6% only slightly larger than the difference between the objective function of the SP-based method and the SAA solution (7.1%). The fact that these two methods produce solutions that have a rather similar objective function is mainly a coincidence because the solutions themselves are quite different. The expectational method ships more export containers per barge and with that also fewer export containers per truck than the SP-based method. Given the specific parameters that we have used, the extra recourse costs incurred by the expectational method are comparable with the costs it saves by using fewer trucks.

It should be noted that the standard deviations of the objective functions are substantial, so the conclusion drawn above should be made with some reservations. Moreover, it is remarkable to see that the SAA method has the largest standard deviation from all methods. For three instances, the SAA method finds solutions that result in a much lower objective function than for the other instances. For these three instances, the other methods do not find these good solutions and thus the objective function of these methods has a lower standard deviation.

The quality of the SP-based method is comparable to the SAA method for \(N=10\). Nevertheless, the running time of the SP-based method is negligible compared to the SAA method. Only Algorithm 1 and a single ILP have to be solved to obtain a solution for the SP-based method, which is all done in about a second, whereas for the SAA method with \(N=10\) the training phase takes 3–5 s and the validation phase 12 min.

6.3 Large instances

For the small instances, we have shown that the SAA method converges to the optimal solution and that the SP-based method is the best heuristic method. In this section, we investigate the performance of the five methods for larger instances. We consider ten instances that consist of 750 import containers, 750 export containers, and four barges with a capacity of 250 TEU. Similar to the small instances, the containers have a size 1 or 2 TEU and the characteristics of the costs are the same as for the small instances. On top of that, we have added the condition that with probability \(\frac{1}{4}\) a container cannot be transported on a barge. With this condition, we model the situation in which containers are not available for transportation on a certain barge because they have not arrived at the deep-sea port yet or have to be at the customer earlier than the arrival of the barge at the inland terminal. The containers are transshipped via ten deep-sea terminals and each barge is only allowed to visit five of them. The expectation of the number of moves at a terminal is set to 75.

Table 4 Average running time in seconds of the training phase of the SAA method for the large instance, for the three terminal types and different numbers of scenarios

For the small instances, the running time of the SAA method is small enough to solve the problem almost to optimality within a reasonable time. Nevertheless, for the large instances, the running times become too big for large values of M, N and \(N'\). In Table 4, the running times of the training phase are given when \(M=1\) and for different values of N. For predictable terminals, the SAA method takes much longer than for the two other terminal types. Moreover, the running time for the predictable and unpredictable terminals when \(N=50\) is about twenty-five times larger than for \(N=10\). The running time for the open-closed terminals when \(N=50\) is even about forty times larger than for \(N=10\). For the predictable and open-closed terminals, the running time is about five times larger for \(N=100\) than for \(N=50\). Since we want to compute a solution in 3–4 h, we have to decide not to use \(N=100\) for the predictable terminals. We believe that the running time for the predictable terminals is much larger because the variety in the number of moves at a terminal is much lower and thus the scenarios are more similar. Consequently, the value of the solutions for visiting different terminals are close to each other and as a result, the branch-and-bound method in the ILP solver needs longer to find the optimal solution.

The running time of the validation phase does not vary much for the different terminal types and is linear in \(N'\). Solving the underlying ILP for a single scenario takes about 0.06 s. Hence, if we denote the running time from Table 4 by r, the total running time for the SAA method with parameters is M, N and \(N'\) is \(Mr + 0.06MN'\) s. We use this formula to create different parameter settings, for which we expect that the SAA method can be solved within 3–4 h. For the unpredictable and open-closed terminal types, we have created six different parameter settings. The average value and standard deviation of the objective functions, the optimality gaps and the running times for these six different parameter settings are given in Table 5. For predictable terminals, only four parameter settings are considered because the running time for this terminal type is much larger. The average value and standard deviation of the objective functions, the optimality gaps and the running times for the predictable terminal type are given in Table 6. Similar as for the small instances, we use a set consisting of 5000 scenarios to calculate the objective function in Tables 5 and 6.

Table 5 The average objective function for the ten large instances for the upredictable and open-closed terminals and different parameter settings of the SAA method
Table 6 The average objective function for the ten large instances for the predictable terminal type and different parameter settings of the SAA method

Based on the results from Table 5, the parameters \(M=20\), \(N=50\), and \(N'=1000\) gives the best results for the unpredictable and the open-closed terminal types. For both types of terminals, the objective function for \(N=50\) is substantially lower than for \(N=10\). Increasing the number of SAA runs only decreases the objective function slightly. The difference between \(N=50\) and \(N=100\) is not as big as between \(N=10\) and \(N=50\). In Table 1, we have already seen that for the small instances the biggest improvement was also made by N going from 10 to 50. Moreover, for \(N=100\) only four runs of the SAA algorithm could be performed. At first it might be surprising that the optimality gap for both terminals for \(N=100\), \(M = 4\), \(N'=5000\) is smaller than for \(N=50\), \(M=20\), \(N'=1000\), but that the value of the objective function is larger. Nevertheless, this can be explained by the fact that the lower bound \((\bar{t}^M_{N'})\) for larger N is stronger. If we would use the lower bound for \(N=100\) for \(N=50\), the optimality gap for the latter will be smaller than for the former.

The best parameters for the predictable terminals are, according to Table 6, \(N=10\), \(M=50\) and \(N'=1500\). The objective function for \(N=50\), \(M=2\) and \(N=5000\) is only slightly worse, if more runs of the SAA method had been possible, the objective function would, most probably, have been better for \(N=50\). Compared with the unpredictable and open-closed terminal types a single run of the training phase of the SAA method has a longer running time. However, we see in Table 6 that the solution produced by the SAA method for the predictable terminals has an optimality gap of only 0.2%, which is much smaller than the optimality gaps for the unpredictable and open-closed terminals.

Table 7 Average objective function and its standard deviation for the 10 large instances for the SP-based method using different settings for \(\gamma _{0k}^i-\gamma _{jk}^i\)
Table 8 Average objective function over ten large instances for the five solution methods and for the three different terminal types

In Table 7, the consequences of different values for \(\gamma _{0k}^i-\gamma _{jk}^i\) for the SP-based method are investigated. For the large instances, the best value for the difference between the two gamma values is 150, in contrast to the small instances for which 50 gave the best results. However, it should be noted that for all three terminal types the difference between the three different settings for the gamma value is extremely small.

In Table 8, the SAA solutions are compared, in a similar fashion as Table 3, to the solutions of the expectational method, the risk-averse trimmed mean method, the CPE method, and the SP-based method. Also for the large instances, the SP-based method has on average the smallest gap with the SAA solution. However, in contrast to the small instances, the SP-based method is not for all three terminal types the best: for the open-closed terminal, the risk-averse trimmed mean method is performing better. Nevertheless, the risk-averse trimmed mean method gives for the predictable terminal type a solution that is 10% worse than the SAA method. Hence, the SP-based method is more robust for different terminal types than the risk-averse trimmed mean method. For the predictable terminal types, we see that the expectational method, the CPE method and the SP-based method all produce solutions that are close to the best SAA solution that is found. Therefore, we may conclude that it might not be that beneficial to include the stochasticity of the number of moves into account: only using the expectation already produces good results for the predictable terminals. That is mainly due to the fact that if the moves are Poisson distributed, then it will hardly happen that the number of moves is exceeding the number of export containers loaded on the barge and thus few recourse costs have to be paid. On the other hand, for the unpredictable and open-closed terminal types, the expectational method results in much recourse costs and the risk-averse trimmed mean method gives lower cost, because fewer export containers are loaded using this method.

For unpredictable and open-closed terminals, one should keep in mind that the SAA solutions had still an optimality gap of a few percentages. Hence, for the methods given in Table 8, the difference with the optimal solution is likely to be bigger than the reported difference with the SAA solution. Another thing to keep in mind is that the running time of the SP-based method is only a couple of seconds. Hence, the SP-based method is a good scalable alternative for the SAA method. A final observation to be made is that, similar to the small examples, a predictable terminal gives by far the lowest costs. However, although for the small instances the unpredictable and open-closed terminals had almost the same value for the objective function, for the large instances the objective function for the unpredictable terminals is much lower than the open-closed terminals. A possible explanation for this could be that the expectation for the number of moves for the small and large instances differ. For the open-closed terminals, the number of moves is either zero or it is Poisson distributed. Since the expectation for the number of moves for the small instances is lower than for the large instances, the Poisson distribution for the small instances also has a lower expectation. Consequently, the difference between being open or closed is less for the small instances than for the large instances. Although the variance for the unpredictable terminals is higher, the realizations are more evenly distributed over the support of the probability distributions. Hence, it is possible to have a better trade-off between the recourse costs and the transportation costs.

6.4 Correlated instances

In this section, we investigate the performances of the five methods when the number of moves is correlated. The instances are the same as used for the small instance in Sect. 6.2, but the difference is that there is a positive correlation for the number of moves at terminal k between the two different barges. Let \(\varPsi _{1k}\) and \(\varPsi _{2k}\) be the correlated number of moves at terminal k for these two barges. Given the two uncorrelated random variables \(\varPhi _{1k}\) and \(\varPhi _{2k}\) and a value \(\alpha \in [0,1]\), the variables \(\varPsi _{1k}\) and \(\varPsi _{2k}\) are defined as follows:

$$\begin{aligned} \varPsi _{1k}&= \varPhi _{1k} \quad k =1, \ldots , l\\ \varPsi _{2k}&= \lceil \alpha \varPhi _{1k}\rceil + \lfloor (1-\alpha )\varPhi _{2k}\rfloor \quad k = 1, \ldots , l. \end{aligned}$$

The value for \(\alpha\) is a parameter to determine the amount of correlation: if \(\alpha\) is zero, then \(\varPsi _{1k}\) and \(\varPsi _{2k}\) are uncorrelated and if \(\alpha =1\), then the variables \(\varPsi _{1k}\) and \(\varPsi _{2k}\) always have the same value. In this paper, we assume that \(\varPhi _{1k}\) and \(\varPhi _{2k}\) are drawn from the same distribution. Consequently, the variables \(\varPsi _{1k}\) and \(\varPsi _{2k}\) are uniquely determined by a value for \(\alpha\) and a distribution for \(\varPhi _{jk}\) for \(j=1,2\). Moreover, the expectation of all four random variables \(\varPhi _{1k}\), \(\varPhi _{2k}\), \(\varPsi _{1k}\) and \(\varPsi _{2k}\) is the same. It is trivially to see that this statement is true for the first three random variables and for \(\varPsi _{2k}\) we have:

$$\begin{aligned} {\mathbb {E}}\left[ \varPsi _{2k}\right]&= {\mathbb {E}}\left[ \lceil \alpha \varPhi _{1k}\rceil + \lfloor (1-\alpha )\varPhi _{2k}\rfloor \right] \\&= {\mathbb {E}}\left[ \lceil \alpha \varPhi _{1k}\rceil \right] +{\mathbb {E}}\left[ \varPhi _{2k}\right] - {\mathbb {E}}\left[ \lceil \alpha \varPhi _{2k}\rceil \right] \\&= {\mathbb {E}}\left[ \varPhi _{2k}\right] . \end{aligned}$$

The distribution of \(\varPsi _{1k}\) is the same as used for \(\varPhi _{1k}\), but \(\varPsi _{2k}\) follows a different distribution. In “Appendix 3”, it is shown that the probability density function of \(\varPsi _{2k}\) can be expressed in terms of the distribution function \(F(\cdot )\) for \(\varPhi _{jk}\), namely:

$$\begin{aligned} {\mathbb {P}}\left( \varPsi _{2k}= \psi \right) =\sum _{\phi =0}^{\psi } \left( F\left( \frac{\phi }{\alpha }\right) - F\left( \frac{\phi -1}{\alpha }\right) \right) \left( F\left( \left\lceil \frac{\psi -\phi +\alpha }{1-\alpha }\right\rceil \right) - F\left( \left\lceil \frac{\psi -\phi +\alpha -1}{1-\alpha }\right\rceil \right) \right) , \end{aligned}$$

for \(\psi = 0, 1, \ldots\) and \(k = 1, 2, \ldots , l.\) Knowing the probability distribution for \(\varPsi _{1k}\) and \(\varPsi _{2k}\), the expectational method, the risk-averse trimmed mean method, and the SP-based method can be applied to the correlated instances. Furthermore, it is also possible to generate correlated samples and thus the SAA and CPE methods can also be used for the correlated instance.

In Table 9, the results for the five different solution methods for the predictable and unpredictable terminals and \(\alpha =0.25, 0.5\) and 0.75 are given. We have chosen not to use the open-closed distribution because in the correlated distribution the realizations do longer follow the idea of being either 0 or having a high value. The main conclusion for the correlated samples is that the SP-based method is also the method that is the closest to the SAA method. Nevertheless, for the unpredictable terminals, the risk-averse trimmed mean method is performing almost as good. Furthermore, the value of \(\alpha\) does not have a big influence on the performance of the methods. The only difference between the small instances and the correlated instances is the distribution of the number of moves. If we compare the results of Table 9 with the results for the small instances in Table 3, we see that the value of the objective for the solution of the SAA method is slightly larger for the correlated instances than for the small instances. However, for the other methods, the value of the objective function for the correlated instances is about the same as for the small instances or even slightly better. A possible explanation could be that in the SAA method, it is not explicitly defined that the visits for a terminal by the two barges are correlated. The method can only learn that from the data, which is also the case for the CPE method. However, the other three methods can use the correlated distribution function.

Table 9 Average objective function over the 10 correlated instances for the five solution methods and the predictable and unpredictable terminal types

7 Conclusion

In this paper, we have studied a problem motivated by real-life practice: an inland terminal in the port of Amsterdam needs to decide how to ship containers from and towards congested deep-sea terminals. In this problem, the number of containers that can be loaded and unloaded at a deep-sea terminal is unknown when the export containers are loaded on the barges. We have modeled this problem as a two-stage stochastic program with recourse. We have presented an SAA method that can solve small instances almost to optimality. Nevertheless, the SAA method is not scalable and thus, for larger instances, it takes too long to produce almost optimal solutions. As in practice fast solutions are required, we have also developed a fast heuristic method. The idea behind this method is that stochastic programming can be used to find the optimal solution for a simplified problem. The characteristics of this optimal solution to the simplified problem are used in the SP-based method that we have presented to solve the original problem. We have compared the results of this SP-based method with three methods for general stochastic assignment problems: the expectational method, the risk-averse trimmed mean method, and the comparative performance evaluation method. We have tested these methods for three different terminal types: predictable, unpredictable and open-closed terminals. Moreover, we have also used three different types of instances: instances with a small number of containers and an uncorrelated number of moves, instances with a large number of containers and an uncorrelated number of moves and instances with a small number of containers and a correlated number of moves.

The SAA method produces for the small instances almost optimal solutions. For the larger instances and certain terminal types, the solution that is produced by the SAA method is about 2–4% from the optimal solution. The SP-based method performs almost always better than the other three methods and for the large instances, it can compute in a couple of seconds a solution that is only less than 1% worse than the solution of the SAA method, which requires a couple of hours of computing time. We have also seen that the performance of the five methods is not much different if the number of moves at a terminal is correlated between all the barges. All in all, if the planning is allowed to take a few hours, the SAA method is the best method to use, but the SP-based method is shown to be a good alternative for a faster solution. Concerning the different terminal types, we conclude that predictable terminals result in the lowest cost. For the small instances, the difference between the unpredictable and open-closed terminals are not that large, but for the large instances, we conclude that the inland terminal has lower costs if the terminal is unpredictable than if the terminal is of the type open-closed.

The SAA method was not able to solve the large instances to optimality. In this paper, a simple implementation of the SAA method was used. In further research, it could be investigated if a more advanced decomposition to solve the SAA method will result in faster solutions for the SAA method. At the moment, the SAA method takes up to a couple of hours to solve the problem, whereas the SP-based method only uses a couple of seconds. However, in practice, one might be willing to wait a couple of minutes for a good solution. Hence, a direction of further research could be to find another method in which the running time and the solution quality are in between the SAA method and the SP-based method. In the current formulation of the problem, we have made some simplifying assumptions which could also be relaxed in further research. Including train transport as a mode of transportation that is cheaper than truck transportation, but more expensive barge transportation could be an interesting option. Moreover, including the route a barge has to sail to visit the terminals is also an option. In that case, one also has to make sure that the capacity of the barge is not violating between two deep-sea terminals. Moreover, solving the resulting ILP formulation is expected to take longer. Finally, one could also include that the number of moves for the first barge is revealed before the export containers for the other barges has to be loaded. The resulting problem would be a multi-stage stochastic problem.