Stochastic Empty Container Repositioning Problem with CO2 Emission Considerations for an Intermodal Transportation System

As one of main challenge for carriers, empty container repositioning is subject to various uncertain factors in practice, which causes more operation costs. At the same time, the movements of empty containers can result in air pollution because of the CO2 emission, which has a negative impact on sustainable development. To incorporate environmental and stochastic characteristics of container shipping, in this paper, an empty container repositioning problem, taking into account CO2 emission, stochastic demand, and supply, is introduced in a sea–rail intermodal transportation system. This problem is formulated as a chance-constrained nonlinear integer programming model minimising the expected value of total weighted cost. A sample average approximation method is applied to convert this model into its deterministic equivalents, which is then solved by the proposed two-phase tabu search algorithm. A numerical example is studied to conclude that the stochastic demand and supply lead to more repositioning and CO2 emission-related cost. Total cost, inventory cost, and leasing cost increase with the variabilities of uncertain parameters. We also found that the total cost and other component costs are strongly dependent on the weights of repositioning cost and CO2 emission-related cost. Additionally, the sensitivity analysis is conducted on unit leasing cost.


Introduction
Due to the overwhelming growth of international trade in recent years, international freight transportation develops rapidly. Compared to unimodal freight transportation (e.g., railway, road, and maritime transportation), intermodal freight transportation is more competitive and advantageous with respect to security, flexibility, and productivity [1]. A significant amount of containers being transshipped are empty. As a result, intermodal freight transportation has been playing a more and more important role in global freight transportation [2]. Owing to different economic needs in different regions, the volumes of import and export trade in many regions are also quite different, which leads to the imbalance of global trade. For instance, the container flow from Eastern Asia to North America reaches a volume of 17.9 million TEUs, while the volume of the opposite direction (i.e., from North to Eastern Asia) was only 8.2 million TEUs, in 2017 [3]. Consequently, intermodal freight industries have to reposition empty containers, from the regions with empty container accumulation to the regions with empty container shortage. However, in practice, the surplus empty containers sometimes are not sufficient to fulfil the customer's demand. Thus, additional empty containers need to be leased from leasing companies, which is a more responsive way compared with repositioning. As a result, intermodal freight industries should make a highly efficient operation plan with minimal expenditure, incorporating the policies of leasing and purchasing containers to minimise total cost. A hybrid genetic algorithm was employed to solve the proposed model. Different from the literature above, some papers allow for ECRP with other decision-making problems, e.g., cargo routing problem, service network design problem, and ship deployment problem. For instance, Brouer et al. [17] addressed the cargo routing problem and ECRP simultaneously, and formulated it as a multi-commodity flow problem with interbalancing constraints, which was solved by the delayed column generation algorithm. The problem was further studied by Song and Dong [19] on a maritime transportation network with multiple service routes, deployed vessels, and regular voyages. Regarding service network design problem, Shintani et al. [14] studied the ECRP combined with a container liner service network design, which was formulated as a two-stage problem solved by a heuristic based on genetic algorithm. This problem was later investigated by Meng and Wang [18], who proposed a mixed integer linear programming model to formulate this problem. In their study, the concept of segment was firstly introduced to facilitate the problem formulation. A realistic case of Asia-Europe-Oceania shipping was carried out to demonstrate the proposed model, and concluded that the total operation cost could be reduced significantly by considering empty container repositioning. The ECRP raised by these two papers is limited to global ocean shipping. By contrast, An et al. [21] presented the service network design for inland water liner transportation with ECRP. By incorporating the specific features of inland water transportation, a nonlinear mixed integer programming model was constructed, and subsequently solved by a hybrid algorithm integrating genetic algorithm and integer programming method. Furthermore, Song and Dong [20] proposed a three-stage solution methodology, where the first, second, and third stage were used to deal with the service network design problem, ECRP, and ship deployment problem, respectively.
Compared with the literature about unimodal ECRP, some studies are concerned with the ECRP in an intermodal transportation system. Specifically, Choong et al. [35] analysed the effect of planning horizon length on mode selection of ECRP for intermodal transportation network, and indicated that the slower and cheaper transportation modes are more encouraged under the condition that the planning horizon is long enough. Meng et al. [36] developed a linear mixed integer programming model for intermodal liner service network design, which accounted for laden container routing, ECRP, and ship route selection. A combination of intermodal service network design problem and ECRP was studied by Braekers et al. [37] in a barge hinterland system, and discussed from different perspectives, i.e., the perspective of barge operators and the perspective of shipping lines. In addition to quantitative analysis on ECRP, Kolar et al. [38] described a research agenda of ECRP based on a qualitative analysis by conducting interviews with representatives of ocean carriers, and indicated some insightful findings.
All the literature mentioned above is studied under deterministic circumstances. However, generally, empty container repositioning is subjected to various stochastic factors in practice. For instance, the number of empty containers available cannot be forecasted precisely, because of stochastic demand and supply. For empty container demand, it is difficult to obtain exact future volume, due to a high dynamic environment and low precision forecast [2]. Similarly, the accurate future volume of empty container supply is also hard to estimate, due to the uncertain returning time of containers from customers [22]. Hence, both demand and supply should be described as uncertain parameters when addressing ECRP. To incorporate the stochastic nature of empty container repositioning, many efforts have been made to explore stochastic ECRP models (e.g., [17][18][19]28,[31][32][33][34]). Crainic et al. [23] focused on the ECRP that involved stochastic demand and supply modelled as random variables in the context of the land transportation of international maritime shipping. Moreover, Li et al. [24] addressed the ECRP in one port, where the amounts of import and export containers were treated as independent and identically distributed variables. They extended the results to the long-run average case and the infinite-horizon problem with a discounted factor. This problem was further generalized by Li et al. [25] to take multi-ports into consideration. The allocation plan of empty containers was derived by using a heuristic algorithm with ε-optimal policies. Lam et al. [27] presented a dynamic programming model for stochastic ECRP, and introduced an approximate approach based on temporal difference learning. Song and Zhang [30] constructed a fluid flow model for ECRP concerning stochastic demand, which is simulated as a two-state Markov chain. In order to examine the impact of different stochastic demand patterns on ECRP, Song and Dong [31] analysed the performances under three types of different demand distributions (i.e., normal distribution, uniform distribution, and exponential distribution). Lee et al. [32] focused on the joint container fleet sizing problem, and ECRP with the consideration of stochastic demand based on the feedback mechanism of inventory information, while Dong et al. [33] compared the results of different empty container repositioning policies under uncertain demand.
In addition to stochastic demand and supply, there are some papers concentrating on stochastic capacity. For example, Long et al. [2] set up a two-stage stochastic programming model involving uncertain demand, supply, ship weight capacity, and ship space capacity. Stage one determined the operation plan for repositioning empty containers with deterministic parameters, while stage two adjusted the decisions derived from stage one based on the probability distribution of random variables. Sample average approximation and progressive hedging strategy were applied to solve this model. A two-stage stochastic network model for dynamic ECRP was established by Cheung and Chen [22] in an ocean transportation network, to minimise the expected total cost by taking into account stochastic demand, supply, and capacity. Francesco et al. [29] addressed the ECRP in a maritime transportation system concerning stochastic demand, supply, residual transportation capacity of vessels, and loading and unloading capacity.
Additionally, for unimodal ECRP, Song [26] examined the ECRP from the perspective of the optimal stationary policy in two-terminal shuttle service system with uncertain demand. By exploring the structural characteristics of objective function, a Markov decision process method was used to obtain the results. Song and Dong [28] extended this model to the ECRP in a cyclic route under stochastic and dynamic conditions, which analysed the impact of demand patterns, fleet sizes and, in particular, degrees of uncertainty. For stochastic intermodal ECRP, Dong and Song [39] dealt with the container fleet sizing problem, combined with ECRP involving stochastic demand and inland transportation time, and revealed the relationship between the optimal container fleet size and the inland transportation time. Xie et al. [40] emphasised the coordination and sharing issue of empty container management with uncertain demand and supply between a liner company and a rail company. The optimal delivery policy, and its change with initial inventory state, were investigated under the centralised model, whereas a bilateral buy-back contract was designed in the decentralised model.
Although the literature on ECRP is substantial and extensive, most literature mainly focuses on the operation cost of repositioning empty containers, rather than environmental performance, e.g., CO 2 emissions. So far, the research on ECRP with CO 2 emissions is still scarce, both practically and theoretically. To the best of our knowledge, only Song and Xu [41], and Bernat et al. [34] studied CO 2 emissions whilst repositioning empty containers. Song and Xu [41] provided an operational activity-based method to calculate CO 2 emissions in relation to ECRP, and revealed the relationship among CO 2 emissions, empty container repositioning policies, and handling capacity. More recently, Bernat et al. [34] presented new periodic review policies integrating CO 2 emissions, repair option, street-turns, as well as stochastic parameters. A simulation model with metaheuristic approach was then developed to assess the optimality of these policies.
Empty container repositioning is an extremely important issue in the container shipping industry [11]. A highly efficient repositioning plan can effectively reduce the amount of repositioned and leased empty containers, and shorten the total repositioning distance from the nodes with surplus empty containers to the nodes with deficit empty containers, which leads to a reduction of total repositioning cost for container carriers. However, uncertainty in container shipping may have a negative effect on empty container movements [15], which causes the inefficiency, and even infeasibility, of the current repositioning plan. Thus, a reliable and robust repositioning plan is required to tackle these unpredictable factors in the container transportation system. At the same time, a well-organised repositioning plan can also reduce fuel consumption, traffic congestion, and emissions, due to the decrease of empty container movements. Accordingly, empty container repositioning not only has an economic impact on container carriers, but also an environmental influence on the sustainable development of containerised transportation.
To incorporate environmental and stochastic characteristics of container shipping, this paper proposes a chance-constrained programming model for ECRP in an intermodal transportation system, which allows for, explicitly, CO 2 emissions and uncertainty simultaneously. Stochastic demand and supply, in the model, are treated as random variables. The purpose of this paper is to determine how many empty containers should be repositioned, which route should be used to transport empty containers, as well as where and how many empty containers should be leased in order to satisfy the demands, which are referred to as distribution plan, routing plan, and leasing plan, respectively. In addition, we conduct a numerical example to evaluate the results obtained under deterministic and stochastic situations, as well as investigate the impact of model parameters and random variables on the total cost and component costs. The main contributions of this paper are three-fold. First and most significantly, it enriches the literature on ECRP by incorporating CO 2 emissions, and stochastic demand and supply. The uncertainty of ECRP facilitates us to construct a chance-constrained programming model, achieving a minimal expected total cost. Second, we use sample average approximation (SAA) to convert this model into its deterministic equivalents, and develop a two-phase tabu search algorithm to solve the associated SAA problem. Third, we examine the effects of unit leasing cost, stochastic demand and supply, as well as weights on costs and CO 2 emissions.
The rest of this paper is organised as follows. Section 2 formulates stochastic ECRP with CO 2 emissions as a chance-constrained programming model. Section 3 presents a hybrid solution methodology integrating a two-phase tabu search algorithm and SAA method. Section 4 illustrates the computational results. The conclusions and future research directions are given in Section 5.

Mathematical Model for Stochastic ECRP with CO 2 Emissions
The aim of this section is to present a chance-constrained mixed integer programming model for stochastic ECRP with CO 2 emissions, minimising the expected total cost. Section 2.1 introduces notations and assumptions, whilst Section 2.2 describes the mathematical model. Table 2 lists notations to be used in the following sections.

Index and set
V the set of nodes, V = V rail ∪ V port V rail the set of railway stations V port the set of ports V s rail the set of railway stations with surplus empty containers, V s rail ⊆ V rail V d rail the set of railway stations with deficit empty containers, V d rail ⊆ V rail V s port the set of ports with surplus empty containers, V s port ⊆ V port V d port the set of ports with deficit empty containers, V d port ⊆ V port A the set of arcs, A = A rail ∪ A ship A rail the set of railway arcs, including the arcs between two railway stations and the arcs between one railway station and one port A ship the set of ship arcs, i.e., the arcs between two ports The mathematic model for stochastic empty container repositioning is constructed reliant on the following assumptions.
Assumption 1: Compared with other transportation modes, road transportation is unsuitable and costly for transporting containers over long distance, due to its working hour restrictions for truck drivers, infrastructure deficits, and relatively higher operation costs. In addition, because of the pollution caused by the burning process of fuel, it is not an eco-friendly transportation mode. Therefore, only two transportation modes are taken into account, i.e., railway transportation and maritime transportation. Thus, the intermodal network we discuss in this problem only includes two types of nodes, i.e., railway stations and ports.
Assumption 2: The stochastic demand and supply of empty containers at any railway station and port are independent of time periods. That is, the variabilities of stochastic demand and supply are kept the same during the planning horizon.
Assumption 3: CO 2 emission is only incurred by transporting empty containers from the nodes with surplus empty containers to the nodes with deficit empty containers. In other words, it is assumed that leasing empty containers will not contribute to CO 2 emission. Assumption 4: There is no restriction on the amount of leased empty containers at any railway station and port at any time period. When surplus empty containers cannot satisfy the demand, sufficient empty containers are leased to fulfil the shortage in a complementary way. In addition, the leased empty containers are not required to be returned in the planning horizon.
Assumption 5: In order to simplify this model and reduce computational burden, the ECRP in this paper is addressed, regardless of multiple types of containers. Therefore, we only consider the standard containers, i.e., twenty-foot equivalent unit (TEU).
Assumption 6: All empty containers can be used without considering damage, repairs, and discards during the planning horizon. This assumption can be relaxed by adding the percentages of defective empty containers into the model if necessary. Assumption 7: The scope of this study is limited to the regional intermodal network, including railway hinterland network of ports and the ship routes among these ports. Additionally, the train services are all assumed to be direct, to ensure that inland transportation time is short enough compared with the one time period. Similarly, the travel times of all ship routes are much smaller than one time period. Consequently, we assume that all repositioned empty containers can be transported from their origins to their destinations within one time period.
Assumption 8: All demands of empty containers can be satisfied at each time period by repositioning and leasing empty containers, which is consistent with Assumption 4.
Assumption 9: For the stochastic demands at each railway station and port, they follow different uniform distributions, respectively, with the same variability but different expected values. The stochastic supplies are in the same manner as illustrated in Section 4.2.

A Chance-Constrained Nonlinear Integer Programming Model
The objective function of the model is to minimise the expected value of total weighted cost, which is the weighted sum of repositioning cost and CO 2 emission-related cost. Thus, the total weighted cost at time period t can be defined by Equation (1), where Q 1 X t and Q 2 X t are repositioning cost and CO 2 emission-related cost respectively, and X t is the vector of decision variables at time period t.
According to the notations as presented in Section 2.1, Q 1 X t and Q 2 X t are defined as below.
The repositioning cost is represented by Equation (2), which includes transportation cost, handling cost (i.e., the sum of loading cost and unloading cost), inventory cost, and leasing cost, while CO 2 emission-related cost is expressed by Equation (3).
Let ξ be a random vector including all random variables at time periods, i.e., The proposed stochastic empty container repositioning problem can be formulated as a chance-constrained nonlinear integer programming model.
subject to Pr Pr Pr Constraints (6)-(7) impose that the maximum number of repositioned empty containers from surplus railway stations and ports cannot exceed their available empty containers with a possibility of at least ε 1 . Constraints (8)- (9) enforce that all demands of empty containers must be satisfied by repositioning and/or leasing, with a possibility of at least ε 2 at deficit railway stations and ports.
(ii) Flow conservation constraints Constraints (10)-(13) guarantee the flow conservation of empty containers stored at all railway stations and ports. For example, the total number of empty containers stored at surplus railway station i at time t, equals the sum of empty containers stored at time t -1, and those supplied minus the demand of empty containers, and those repositioned to deficit railway stations and ports at time period t.
∑ p∈V s port ∑ q∈V d port x t pq y r,t pq δ r a ≤ cap r , ∀a ∈ A ship , ∀r ∈ R ship , ∀t = 1, 2, . . . , T, Constraint (14) limits the total number of empty containers transported by a train on arc a, while Constraint (15) restricts the arc flow of empty containers on ship route r. Constraints (16)-(19) state that the amount of empty containers loaded and unloaded is within the handling capacity at railway stations and ports. Constraints (20)-(21) ensure that the number of empty containers stored cannot exceed the inventory capacity at railway stations and ports.
(iv) Route selection constraints ∑ r∈R ship y r,t pq = 1, ∀p ∈ V s port , ∀q ∈ V d port , ∀t = 1, 2, . . . , T, Constraints (22)- (24) ensure that only one route can be used to reposition empty containers from one node with surplus empty containers to another with deficit empty containers, that is, the repositioned empty containers with the same origin and destination cannot be split over different routes.

Solution Methodology
In this section, we propose a heuristic solution methodology to solve the mathematical model for stochastic ECRP with CO 2 emissions in Section 2. Specifically, the SAA method is adopted to convert the chance-constrained mixed integer programming model into its deterministic equivalents in Section 3.1, while a two-phase tabu search algorithm is developed in Section 3.2.

SAA Method
The chance-constrained programming model was first proposed by Charnes and Cooper [42], and is then further studied in the domain of stochastic optimisation. But, now, it is still a challenging problem to solve. So far, there are some stochastic optimisation methods for solving the stochastic programming model, e.g., stochastic approximation, stochastic gradient descent, finite-difference stochastic approximation, simultaneous perturbation stochastic approximation, and scenario optimisation. However, these methods all have difficulties in calculating the values of uncertain functions accurately, especially for the approximation of chance constraints. SAA is a widely used approach proposed by Kleywegt et al. [43] for solving stochastic programming model based on Monte Carlo simulation. Compared with other stochastic programming methods, SAA can effectively approximate the expected total cost and two types of chance constraints regarding stochastic demand and supply in our model [44]. Therefore, we employ the SAA method to solve the stochastic ECRP with CO 2 emissions. In detail, we first generate an independent identically distributed sample ξ 1 , ξ 2 , . . . , ξ N , according to the probability distributions of random variables, which consists of N realisations of the random vector ξ, i.e., Then, the expected value of objective function E T ∑ t=1 Q X t , ξ can be approximated by sample [44]. Similarly, the possibilities that chance constraints are satisfied can be also approximated [45]. Take chance Constraint (6), for example, we first define the indicator function of (0, ∞), 1 (0,∞) : R → R , i.e., and then Equation (6) can be rewritten as According to law of large numbers theory, the value of the right-hand side will converge to the value of left-hand side, with the possibility of 1, as sample size N increases to infinity. Thus, the chance-constrained programming model for ECRP P0 can be converted into the corresponding deterministic SAA problem P1.
subject to Constraints (10)- (27), and the following Constraints (32)- (35): To evaluate the quality of optimal solution derived from SAA method, the optimality gap (i.e., the difference between lower bound and the objective value of optimal solution) need to be computed. In order to get a lower bound for true problem P0, we generate M independent identical samples with N realisations of random vector ξ. As a result, we can acquire M corresponding SAA problems. By solving them, we obtain M candidate solutions X 1 , X 2 , . . . , X M , and M associated objective values . . , f M N . Thus, as suggested by Verweij et al. [44], the lower bound f N can be estimated by Besides, as stated by Luedtke and Ahmed [46], a posteriori analysis should be conducted to see whether the chance constraints are satisfied for each candidate solution. Therefore, a large independent sample ξ 1 , ξ 2 , . . . , ξ N , comprising N realisations of random vector ξ, needs to be generated to check the feasibility of candidate solutions and select the feasible ones. Then, we recalculate the objective values of all feasible solutions denoted asX by and choose the one with the minimal total weighted cost as the optimal solution X * . Accordingly, the optimality gap can be calculated in Equation (38) as follows.

Two-Phase Tabu Search Algorithm
As a mature approach, tabu search algorithm has been extensively applied to transportation operation research, e.g., stop-skipping scheduling problem [47] and vehicle routing problem [48]. Therefore, a two-phase tabu search algorithm is designed to solve the deterministic SAA problem for stochastic ECRP with CO 2 emissions, which is composed of a distribution phase and routing phase. The procedure of two-phase tabu search algorithm, incorporating SAA method, is illustrated as a flowchart in Figure 1.

Distribution Phase
In the distribution phase, a tabu search algorithm is developed to determine how many empty containers are required to reposition from surplus nodes to deficit nodes, as well as where and how many empty containers to be leased at deficit nodes, i.e., distribution plan and leasing plan. The main elements are presented, next, in detail.

(i) Neighbourhood structure
In the procedure of tabu search, neighbourhood structure is a key factor for the performance and efficiency of this algorithm, which transforms the current solution into one of its neighbour solutions in the search space [49]. The neighbourhood structure of distribution phase is obtained by considering add or drop moves based on the current solution. Add moves are defined as the moves that increase the amount of repositioned empty containers between one surplus node and one deficit node, while drop moves are described as the moves that decrease the number of repositioned empty containers between two nodes.
For each decision variable of distribution plan t ij x , a random number is generated from uniform distribution on the interval (0, 1) to compare with the predetermined selection possibility, which indicates whether an add or drop move is performed for the variable. Then, at each iteration of distribution phase, an add/drop move sample is introduced to represent the move status of all distribution variables. The move status is sorted into three categories: drop move, no move, and add move, corresponding to values −1, 0, and 1, respectively, which the add/drop move sample can take, and the variability range of add or drop moves follows a uniform distribution between 1 TEU and 5 TEUs. Accordingly, all neighbour solutions are obtained, based on the add/drop move sample.
For the leasing plan, the corresponding decision variables, t i z and t p z (i.e., the amount of leasing containers at each deficit railway station or port), can be calculated by the difference between the number of repositioned and deficit empty containers. Considering that the leasing cost is time independent, therefore, there is no need to hire surplus empty containers in preparation for the demand of the next time period, which will incur more inventory cost.

(ii) Evaluation of neighbour solutions
To decide which move should be performed, we need to evaluate the objective value of each neighbour solution [50]. In the distribution phase, the handling cost, inventory cost, leasing cost, and CO2 emission-related cost are straightforward and easy to calculate. However, since the routes for transporting empty containers are derived from the routing phase, which is implemented after the distribution phase, the transportation cost is hard to estimate. In order to get a relatively exact

Distribution Phase
In the distribution phase, a tabu search algorithm is developed to determine how many empty containers are required to reposition from surplus nodes to deficit nodes, as well as where and how many empty containers to be leased at deficit nodes, i.e., distribution plan and leasing plan. The main elements are presented, next, in detail.

(i) Neighbourhood structure
In the procedure of tabu search, neighbourhood structure is a key factor for the performance and efficiency of this algorithm, which transforms the current solution into one of its neighbour solutions in the search space [49]. The neighbourhood structure of distribution phase is obtained by considering add or drop moves based on the current solution. Add moves are defined as the moves that increase the amount of repositioned empty containers between one surplus node and one deficit node, while drop moves are described as the moves that decrease the number of repositioned empty containers between two nodes.
For each decision variable of distribution plan x t ij , x t pq , or x t iq , a random number is generated from uniform distribution on the interval (0, 1) to compare with the predetermined selection possibility, which indicates whether an add or drop move is performed for the variable. Then, at each iteration of distribution phase, an add/drop move sample is introduced to represent the move status of all distribution variables. The move status is sorted into three categories: drop move, no move, and add move, corresponding to values −1, 0, and 1, respectively, which the add/drop move sample can take, and the variability range of add or drop moves follows a uniform distribution between 1 TEU and 5 TEUs. Accordingly, all neighbour solutions are obtained, based on the add/drop move sample. For the leasing plan, the corresponding decision variables, z t i and z t p (i.e., the amount of leasing containers at each deficit railway station or port), can be calculated by the difference between the number of repositioned and deficit empty containers. Considering that the leasing cost is time independent, therefore, there is no need to hire surplus empty containers in preparation for the demand of the next time period, which will incur more inventory cost.

(ii) Evaluation of neighbour solutions
To decide which move should be performed, we need to evaluate the objective value of each neighbour solution [50]. In the distribution phase, the handling cost, inventory cost, leasing cost, and CO 2 emission-related cost are straightforward and easy to calculate. However, since the routes for transporting empty containers are derived from the routing phase, which is implemented after the distribution phase, the transportation cost is hard to estimate. In order to get a relatively exact evaluation, we assume that each empty container is shipped on the route with the minimal cost. Thus, the evaluation of neighbour solutions is the sum of transportation cost estimate and other component costs.
(iii) Tabu attributes, tabu duration, and aspiration criterion For add or drop moves in distribution phase, the tabu attribute is the increase or decrease of repositioned empty containers from surplus nodes to deficit nodes, respectively. To avoid the cycling search, a dynamic tabu duration is introduced, which varies with the number of decision variables of the distribution plan. In addition, an aspiration criterion is used to ensure that a tabu add or drop move can be applied when a better solution is found [51].

Routing Phase
The routing phase proposes a tabu search algorithm to specify the routes for transporting the empty containers from surplus nodes to deficit nodes, i.e., a routing plan. The details are stated as follows.
(i) Neighbourhood structure Considering the characteristics of the binary decision variables of routing phase, swap moves are introduced to obtain the neighbour structure, which transit the current solution to its neighbour solutions in the search space. Swap moves are defined as the moves that specify one candidate route to replace the current route for transporting empty containers. This type of move can ensure that all neighbour solutions satisfy the route selection constraints. Similar to the distribution phase, a swap move sample is employed in this phase to indicate whether a swap move is performed for each decision variable of routing phase, i.e., y k,t ij , y r,t pq , and y k,t iq .

(ii) Evaluation of neighbour solutions
Compared with the distribution phase, the evaluation of neighbour solutions in routing phase is more straightforward. For the given distribution plan and leasing plan, routing phase yields the corresponding optimal routing plan, which only causes the change of transportation cost. Therefore, in order to simplify the calculation and select the best swap move in a highly efficient manner, we use the transportation cost as the evaluation value of each neighbour solution.
(iii) Tabu attributes, tabu duration, and aspiration criterion For swap moves, the tabu attribute is the changes of routes for transporting empty containers. Since the number of decision variables of distribution phase is relevant with that of routing phase, tabu duration needs to change dynamically to maintain algorithm performance. In the process of searching, if a better solution is encountered, resulting from one tabu swap move, the tabu status of this move will be updated as nontabu, which is defined as an aspiration criterion.

Computational Study
In order to demonstrate the proposed chance-constrained nonlinear mixed integer programming model, and evaluate the two-phase tabu search algorithm integrating SAA method, a numerical example is studied in Section 4, considering railway and maritime transportation. In addition, we compare the results of deterministic and stochastic cases, and conduct the sensitivity analyses of model parameters and random variables.

Deterministic Case
In this section, a sea-rail intermodal network is considered, which consists of seven railway stations and three ports, as shown in Figure 2. The planning horizon T is three time periods, i.e., t = 1, 2, 3. At the beginning of the planning horizon, the numbers of empty containers at all nodes are assumed to be zero. Additionally, the unit handling cost, the unit inventory cost, the unit leasing cost, and the unit CO 2 emission-related cost are set to be 15 (US$/TEU), 5.6 (US$/TEU), 200 (US$/TEU), and 2 (US$/kg), respectively. and the unit CO2 emission-related cost are set to be 15 (US$/TEU), 5.6 (US$/TEU), 200 (US$/TEU), and 2 (US$/kg), respectively.
The railway network information is listed in Table 3, including the distance, cost, and CO2 emission when empty containers are transported by train, while the maritime network information is displayed in Table 4, including the port calling sequence of each ship route and the cost and CO2 emission between two ports when empty containers are transported by ship. Figure 2. A sea-rail intermodal network.  The railway network information is listed in Table 3, including the distance, cost, and CO 2 emission when empty containers are transported by train, while the maritime network information is displayed in Table 4, including the port calling sequence of each ship route and the cost and CO 2 emission between two ports when empty containers are transported by ship. It is assumed that three stations (i.e., S1, S2, S3) and three ports (i.e., P1, P2, P3), in this case, are treated as surplus railway station/port or deficit railway station/port, which can generate the demand and supply of empty containers. In other words, the other four stations S4, S5, S6, and S7, are treated as a balanced station/port without repositioning and leasing containers. Tables 5 and 6 give the demand and supply of empty containers at these six nodes during the planning horizon.  The key parameters of two-phase tabu search algorithm are described as follows. For the distribution phase, the maximum number of iterations max_d and tabu length are set as 100 and 4. For the routing phase, the maximum number of iterations max_r is set as 50, while the tabu length is dependent on the number of routing decision variables at each time period. The solution methodology, incorporating SAA method and two-phase tabu search algorithm, is coded in MATLAB R2012a. The program is run on a 2.5GHz dual core PC with 4GB RAM.
Given the weights of repositioning cost and CO 2 emission-related cost ω 1 = 1, ω 2 = 1, we can obtain the repositioning plan with a total cost of $65,991, including the distribution plan, routing plan, and leasing plan at each time period, as shown in Table 7. Table 8 lists the costs of the deterministic case, consisting of transportation cost, handling cost, inventory cost, leasing cost, CO 2 emission-related cost, and total cost.  To assess the influence of the weights of repositioning cost and CO 2 emission-related cost, we test the deterministic case under different weight combinations. As shown in Figure 3a,b,e, for a given weight ω 1 , transportation cost, handling cost and CO 2 emission-related cost decline with the increasing value of weight ω 2 . As CO 2 emission is only incurred by repositioning empty containers, less empty containers are transported when CO 2 emission-related cost dominates the total weighted cost. In other words, when the value of weight ω 2 goes up, more empty containers are leased to satisfy the demand because it makes no contribution to CO 2 emission. As a consequence, the three component costs drop to a low level with the decreasing number of repositioned empty containers, and vice versa.
On the contrary, inventory cost and leasing cost present an increasing tendency, as illustrated in Figure 3c,d. As discussed above, less empty containers are repositioned resulting from the growth of weight ω 2 . Consequently, it leads to more empty containers accumulating at surplus railway stations or ports, which incur high inventory cost. On the other hand, more leased empty containers undoubtedly result in the rise of leasing cost. As depicted in Figure 3f, compared with other component costs, total cost ascends significantly with the weights ω 1 and ω 2 both declining from 0.6 to 0. Nevertheless, it grows slowly when the two weights are larger than 0.6. We conjecture that this can be explained by the interaction of multiple component costs and relevant model parameters. Hence, it is can be concluded that total cost and other component costs are strongly dependent on the weight combinations.
Next, we examine the effect of unit leasing cost on the number of leased empty containers under three different weight combinations for each time period and the whole planning horizon. As expected, the amounts of leased empty containers for the three weight settings all descend with the growth of unit leasing cost in Figure 4. However, for the weight setting ω 1 = 0, ω 2 = 1, the number of leased empty containers is less sensitive to the unit leasing cost, and fluctuates fiercely because the change of unit leasing cost cannot affect the total weighted cost when ω 1 = 0. Furthermore, we note that the number of leased empty containers for the weight setting ω 1 = 0, ω 2 = 1 is generally more than the other two weight settings. Since the value of ω 1 for the other two weight settings is 1, that means the leasing cost accounts for a larger percentage of total weighted cost, and freight companies are preferable to repositioned empty containers rather than leased empty containers to reduce high leasing cost.  To assess the influence of the weights of repositioning cost and CO2 emission-related cost, we test the deterministic case under different weight combinations. As shown in Figure 3a, b, e, for a given weight 1 ω , transportation cost, handling cost and CO2 emission-related cost decline with the increasing value of weight 2 ω . As CO2 emission is only incurred by repositioning empty containers, less empty containers are transported when CO2 emission-related cost dominates the total weighted cost. In other words, when the value of weight 2 ω goes up, more empty containers are leased to satisfy the demand because it makes no contribution to CO2 emission. As a consequence, the three component costs drop to a low level with the decreasing number of repositioned empty containers, and vice versa.
On the contrary, inventory cost and leasing cost present an increasing tendency, as illustrated in Figure 3c,d. As discussed above, less empty containers are repositioned resulting from the growth of weight 2 ω . Consequently, it leads to more empty containers accumulating at surplus railway stations or ports, which incur high inventory cost. On the other hand, more leased empty containers undoubtedly result in the rise of leasing cost. As depicted in Figure 3f, compared with other component costs, total cost ascends significantly with the weights 1 ω and 2 ω both declining from 0.6 to 0. Nevertheless, it grows slowly when the two weights are larger than 0.6. We conjecture that this can be explained by the interaction of multiple component costs and relevant model parameters.
Hence, it is can be concluded that total cost and other component costs are strongly dependent on the weight combinations.
(a) (b) Next, we examine the effect of unit leasing cost on the number of leased empty containers under three different weight combinations for each time period and the whole planning horizon. As expected, the amounts of leased empty containers for the three weight settings all descend with the growth of unit leasing cost in Figure 4. However, for the weight setting 1 0 ω = , 2 1 ω = , the number of leased empty containers is less sensitive to the unit leasing cost, and fluctuates fiercely because the change of unit leasing cost cannot affect the total weighted cost when 1 0 ω = . Furthermore, we note that the number of leased empty containers for the weight setting 1 0 ω = , 2 1 ω = is generally more than the other two weight settings. Since the value of 1 ω for the other two weight settings is 1, that means the leasing cost accounts for a larger percentage of total weighted cost, and freight companies are preferable to repositioned empty containers rather than leased empty containers to reduce high leasing cost.  Figure 5 illustrates the influence of unit leasing cost on CO 2 emission-related cost under three different weight combinations. As indicated in Figure 5, CO 2 emission-related cost grows with the increase of unit leasing cost. This can be explained by less leased empty containers, which forces more empty containers to be repositioned as a complementary way to satisfy the demand. At the same time, more CO 2 is emitted during the process of repositioning, resulting in more CO 2 emission-related cost.    Figure 5 illustrates the influence of unit leasing cost on CO2 emission-related cost under three different weight combinations. As indicated in Figure 5, CO2 emission-related cost grows with the increase of unit leasing cost. This can be explained by less leased empty containers, which forces more empty containers to be repositioned as a complementary way to satisfy the demand. At the same time, more CO2 is emitted during the process of repositioning, resulting in more CO2 emission-related cost.

Stochastic Case
As mentioned in Section 2.1, the uncertain demand of empty containers  N , and replication number M, to be 500, 10,000, and 10, respectively. Confidence levels 1 ε and 2 ε are both assumed to be 0.5. In addition, other input parameters are the same with those of deterministic case. Table 9 reports the distribution plan, routing plan, and leasing plan of stochastic case at each time period, while Table  10 shows the associated costs.

Stochastic Case
As mentioned in Section 2.1, the uncertain demand of empty containers d t i and d t p both follow uniform distributions, i.e., d t i ∼ U(m i , m i + λ 1 ), with m i the demand of deterministic case at railway station i, and d t p ∼ U m p , m p + λ 1 , with m p the demand of deterministic case at port p, where coefficient λ 1 is defined as the variability of uncertain demand. Similarly, the uncertain supply of empty containers s t i and s t p also follow uniform distributions, i.e., s t i ∼ U(n i , n i + λ 2 ) and s t p ∼ U n p , n p + λ 2 , where n i and n p , the supplies of deterministic case at railway station i and port p, respectively, as well as λ 2 , the variability of uncertain supply. Here, λ 1 and λ 2 both take value 6. For SAA method parameters, we set sample sizes N, N , and replication number M, to be 500, 10,000, and 10, respectively. Confidence levels ε 1 and ε 2 are both assumed to be 0.5. In addition, other input parameters are the same with those of deterministic case. Table 9 reports the distribution plan, routing plan, and leasing plan of stochastic case at each time period, while Table 10 shows the associated costs.  As compared with the results of deterministic case, stochastic demand and supply result in a different repositioning plan. With respect to the distribution plan, repositioned empty containers rise dramatically from 418 TEUs of deterministic case to 516 TEUs of stochastic case. Similarly, for the leasing plan, leased empty containers grow relatively steadily, from 88 TEUs to 106 TEUs. By contrast, there is no change in the routing plan. Due to the rise of repositioned and leased empty containers, the total cost and other component costs go up considerably in comparison to those of the deterministic case.
To further investigate the effect of stochasticity, we test the stochastic case with different λ 1 and λ 2 , i.e., different variabilities of stochastic demand and supply. As shown in Figure 6, the leasing cost and total cost ascend with coefficient λ 1 . In Figure 6a, we can also observe that the leasing cost with λ 2 = 1 is always more than those with λ 2 = 5 and λ 2 = 9, when λ 1 varies from 1 to 10. This implies that lower supply variability provides less available empty containers than the other two situations, so more leased containers are required to satisfy the demand, which incurs more leasing cost.
To further investigate the effect of stochasticity, we test the stochastic case with different 1 λ and 2 λ , i.e., different variabilities of stochastic demand and supply. As shown in Figure 6, the leasing cost and total cost ascend with coefficient 1 λ . In Figure 6a, we can also observe that the leasing cost with 2 1 λ = is always more than those with 2 5 λ = and 2 9 λ = , when 1 λ varies from 1 to 10. This implies that lower supply variability provides less available empty containers than the other two situations, so more leased containers are required to satisfy the demand, which incurs more leasing cost.
(a) (b) Likewise, the effect of coefficient 2 λ is also analysed. It can be seen from Figure 7 that the inventory cost rises noticeably, while the total cost fluctuates within a certain range and goes up very slightly. In addition, we notice that the total cost with 1 9 λ = is more than other two situations in Figure 7b, which indicates that more operation cost is requested to cope with higher variability in empty container demand. Likewise, the effect of coefficient λ 2 is also analysed. It can be seen from Figure 7 that the inventory cost rises noticeably, while the total cost fluctuates within a certain range and goes up very slightly. In addition, we notice that the total cost with λ 1 = 9 is more than other two situations in Figure 7b, which indicates that more operation cost is requested to cope with higher variability in empty container demand.

Conclusions
This paper has presented an empty container repositioning problem (ECRP) allowing for CO2 emissions, and stochastic demand and supply. To formulate this problem, a chance-constrained mixed integer programming model is proposed to minimise the expected value of total weighted cost in a sea-rail intermodal network. In order to solve it, sample average approximation (SAA) method is used to convert it into its deterministic equivalents. Two-phase tabu search algorithm is then developed to obtain the optimal solution of the converted SAA problem. Specifically, the first phase yields the distribution plan and leasing plan, whilst the second phase produces the routing plan. Finally, a computational study has been carried out to validate the proposed model and developed solution methodology. The results indicate that the repositioning plan and total cost are affected not only by stochastic factors, but also by weight combinations. To explore the effect of them, a series of sensitivity analyses is conducted. We conclude that stochastic demand and supply can both increase the transportation cost, handling cost, inventory cost, leasing cost, CO2 emission-related cost, and total cost compared with the deterministic case. At the same time, the amounts of repositioned and

Conclusions
This paper has presented an empty container repositioning problem (ECRP) allowing for CO 2 emissions, and stochastic demand and supply. To formulate this problem, a chance-constrained mixed integer programming model is proposed to minimise the expected value of total weighted cost in a sea-rail intermodal network. In order to solve it, sample average approximation (SAA) method is used to convert it into its deterministic equivalents. Two-phase tabu search algorithm is then developed to obtain the optimal solution of the converted SAA problem. Specifically, the first phase yields the distribution plan and leasing plan, whilst the second phase produces the routing plan. Finally, a computational study has been carried out to validate the proposed model and developed solution methodology. The results indicate that the repositioning plan and total cost are affected not only by stochastic factors, but also by weight combinations. To explore the effect of them, a series of sensitivity analyses is conducted. We conclude that stochastic demand and supply can both increase the transportation cost, handling cost, inventory cost, leasing cost, CO 2 emission-related cost, and total cost compared with the deterministic case. At the same time, the amounts of repositioned and leased empty containers also present an increasing tendency. Moreover, the impact of variabilities of stochastic factors on costs is examined. Higher variability of stochastic demand results in more leasing cost and total cost, while higher variability of stochastic supply leads to the increased inventory cost and fluctuating total cost. With respect to the weights of repositioning cost and CO 2 emission-related cost, it is found that both of them have a positive or negative relationship with the costs. Transportation cost, handling cost, and CO 2 emission-related cost decline with the increasing value of weight ω 2 and the decrease of weight ω 1 . By contrast, inventory cost and leasing cost present an opposite trend. Additionally, the results also show that the rise of unit leasing cost can lead to the decrease of the number of leased empty containers, and the increase of CO 2 emission-related cost. These findings in this paper can provide implications for decision-makers when an operation plan for repositioning and leasing empty containers needs to be determined in a stochastic environment. Furthermore, the decision-maker can obtain the different repositioning plans under different weight combinations, according to their preferences.
There are two directions which are worth investigating for future research. The first direction is to integrate the routing problem of loaded containers in the decision process of ECRP. The second direction is to consider uncertain factors from the perspective of fuzzy theory.