# A Time-Dependent Fuzzy Programming Approach for the Green Multimodal Routing Problem with Rail Service Capacity Uncertainty and Road Traffic Congestion

Authors: Yan Sun; Martin Hrušovský; Chen Zhang; Maoxiang Lang
Journal: Complexity
 (2018)

---

## Abstract

This study explores an operational-level container routing problem in the road-rail multimodal service network. In response to the demand for an environmentally friendly transportation, we extend the problem into a green version by using both emission charging method and bi-objective optimization to optimize the CO2 emissions in the routing. Two uncertain factors, including capacity uncertainty of rail services and travel time uncertainty of road services, are formulated in order to improve the reliability of the routes. By using the triangular fuzzy numbers and time-dependent travel time to separately model the capacity uncertainty and travel time uncertainty, we establish a fuzzy chance-constrained mixed integer nonlinear programming model. A linearization-based exact solution strategy is designed, so that the problem can be effectively solved by any exact solution algorithm on any mathematical programming software. An empirical case is presented to demonstrate the feasibility of the proposed methods. In the case discussion, sensitivity analysis and bi-objective optimization analysis are used to find that the bi-objective optimization method is more effective than the emission charging method in lowering the CO2 emissions for the given case. Then, we combine sensitivity analysis and fuzzy simulation to identify the best confidence value in the fuzzy chance constraint. All the discussion will help decision makers to better organize the green multimodal transportation.


---

## Body

## 1. Introduction

With the rapid development of globalization, companies are seeking for specialized partners all over the world to outsource their businesses and also for new markets to make more profit. International trade and global commodity circulation thus get significantly motivated during the last decades [1]. Supported by the transportation industry, this prosperity results in a remarkable expansion of the transportation network from a limited region to the entire world. The expanded transportation network extends the distribution channels and enhances the difficulty of transportation organization from the viewpoint of economy, efficiency, reliability, and environmental protection. The traditional unimodal transportation is no longer suitable for this situation. Meanwhile, an advanced transportation mode, namely, the multimodal transportation, emerges in the transportation industry. It combines different transportation modes to generate origin-to-destination routes and integrates their respective advantages in the transportation process. Furthermore, by using standardized containers, multimodal transportation can be highly mechanized and the efficiency of the operations at terminals is greatly improved [2]. Many empirical case studies, for example, Bookbinder and Fox [3], Janic [4], and Liao et al. [5], have demonstrated the superiority of the multimodal transportation in the long-haul transportation setting. Nowadays, multimodal transportation is getting more and more popular and plays a crucial role in international trade. Since the transportation chains in this area are usually long combining different actors, effective solutions for their coordination are necessary in order to ensure that the goods are delivered according to the customer’s requirements. In this context, the multimodal routing problem has gained interest of both researchers and practitioners [6, 7]. This operational problem aims at selecting the best routes to move containers from their origins to destinations through the multimodal service network according to customer demands [2]. In order to achieve this objective, a number of different factors have to be considered.First of all, multimodal transportation combines different modes with different operational constraints that have to be represented in the mathematical model. In general, it can be distinguished between time-flexible modes (e.g., road) that can be used whenever they are needed and schedule-based modes (e.g., rail) that operate according to fixed schedules planned in advance. These schedules regulate the route as well as the (un-)loading, departure, and arrival times of the transportation service. Since these schedules cannot be easily changed, they should be considered when planning the routes of the containers and their transshipment between different transportation services; otherwise, the multimodal plans might be infeasible in reality.Secondly, the actors and stakeholders involved in multimodal transportation systems have usually different interests which are often in conflict. Whereas transportation operators usually tend to minimize only the transportation costs, customers are also interested in the reliability of the routes in order to avoid situations where a plan becomes infeasible due to late arrival of a service or insufficient capacity. Since the multimodal routing is done before the start of transportation, when the information about the transportation services (e.g., traffic situation, travel time, and capacity) is not completely known, possible sources of uncertainty should be considered already during the planning [8]. As an example, the scheduled rail services have usually limited capacity due to the limited length of the tracks in the railway stations and limited locomotive power [9]. However, the available capacities of rail services are uncertain since they are also dependent on other transportation plans that are difficult to predict [10]. If only deterministic capacities are considered in the model, two consequences may happen in practice: the planned routes are either infeasible if the capacities are estimated to be too large, or the planned routes are not the actual best ones if the capacity estimation is too small. Therefore, capacity uncertainty is a factor that should be considered in the decision-making process. In addition to that, road services frequently suffer from various disruptions, such as natural disasters (e.g., bad weather) and manmade disruptions (e.g., traffic accidents and congestion) [11]. Among these disruptions, traffic congestion is very common due to the rapidly increasing traffic volume in recent years [12]. Traffic congestion results in the travel time uncertainty and delayed arrivals of the containers at the nodes, which disrupts the transshipment if the containers have to take another scheduled transportation service at the same node and also possibly violates the due date constraint at the destination. Thus, congestion is another factor that should be also considered in the decision-making process.Thirdly, the environmental aspects of transportation have drawn considerable attention from the public and the government [13]. Among the various environmental issues, global warming caused by greenhouse gas emissions has been considered as a tremendous threat to human sustainable development. CO2 accounts for approximately 80% of the entire greenhouse gas emissions [5]. During recent decades, large numbers of countries have proposed regulations to ease the global warming by reducing CO2 emissions. For example, in 2012, the National Development and Reform Commission of China proposed a pilot project on trading carbon emissions [2]. Since the transportation sector is one of the biggest contributors to the CO2 emissions [14, 15], the transportation industry should accept this responsibility and consider controlling and optimizing the CO2 emissions during the transportation planning stage to promote the development of an environmentally friendly transportation. In order to respond to these environmental requirements, in this study, we extend the multimodal routing problem into a green version by considering the CO2 emissions.The inclusion of uncertain factors and the combination of multiple objectives (i.e., economic and environmental) in multimodal transport planning is an emerging field which has been studied in recent years [13–15]. However, the authors usually consider only one source of uncertainty or model the CO2 emissions in form of charges which limits their influence on the resulting routes. Therefore, we present a multimodal routing planning approach which deals with the following question:How can the sustainability and reliability of multimodal routing be improved by considering environmental factors and uncertainty?The contributions of the paper are fivefold:(1)
A mixed integer nonlinear optimization model for multimodal routing including multiple objectives (i.e., economic and environmental) and multiple uncertainty factors (i.e., road traffic congestion and rail service capacity) is defined.
(2)
A fuzzy logic is proposed to model the capacity uncertainty, and a piecewise linear function is used to represent the road traffic congestion.
(3)
The influence of CO2 emissions on the routes is analyzed comparing the CO2 emissions charging method and a bi-objective approach in which minimizing CO2 emissions is set as an independent objective.
(4)
Since the proposed optimization model is nonlinear, a linearization approach is proposed in order to reduce the complexity of the model.
(5)
The planning approach is applied to a real-world case study from China with specific network design and parameters. In this case, the network is limited to road and rail services; however, other modes with fixed schedules (e.g., inland waterway and deep sea shipping) could be also treated in a similar way as it is the case of rail.The paper is organized as follows: In Section2, we review the current literature on multimodal transportation routing problem. In Section 3, we present a detailed modelling methodology. The improvements of our settings compared to the existing multimodal routing literature are indicated by comparing the two sections. In Section 4, based on the proposed modelling framework, we establish a fuzzy chance-constrained mixed-integer nonlinear programming model to mathematically describe the multimodal routing problem explored in this study. Since the proposed model is nonlinear, in Section 5 we design a linearization-based exact solution strategy to solve the problem. In Section 6, an empirical case study based on the Chinese scenario is given to demonstrate the feasibility of the methods in dealing with the practical problem. Also in the case study section, sensitivity analysis, bi-objective optimization scenario analysis on CO2 emissions, and fuzzy simulation are adopted to draw conclusions which can further help the decision makers to better organize the multimodal transportation. Finally, conclusions of this study are drawn in Section 7.


## 2. Literature Review

The multimodal routing problem is a combinatorial optimization problem which deals with finding the optimal routes for customer demands respecting the constraints of the network [16]. Customer demands are represented by transportation orders characterized by their origin, destination, volume of containers, release time at origin, and due date at destination. Whereas the release time is usually considered as fixed, the due date can be either considered as a hard constraint (see, e.g., Verma et al. [17], Sun and Lang [2, 18], and Sun et al. [19]) or delayed deliveries can be accepted causing additional penalty costs (see, e.g., Demir et al. [20], Hrušovský et al. [14], and Zhang et al. [21].However, considering only late arrivals might not be sufficient since early deliveries can also cause problems with the storage of the goods until they are needed. Therefore, various studies consider due date time windows instead of time points [8, 22], which is also the case in this study. Moreover, modelling a soft time window is more suitable than a hard one, because customers allow the violation to the time window (i.e., early and delayed deliveries) within a certain degree.The containers from transportation orders need to be transported in the multimodal service network by various transportation services with their specific characteristics. One of these characteristics is the usually fixed schedule of the high-capacity and environmentally friendly transportation modes (e.g., railway or inland waterway) in contrast to the time-flexible (nonscheduled) road service. However, the majority of the related articles neglects the schedules and considers rail services to be also time-flexible. This facilitates the modelling since only transshipment times between services in terminals need to be considered (Sun et al. [23], Jiang and Lu [24], Cai et al. [25], Zhang et al. [21], and Xiong and Wang [26]). However, this does not correspond to reality where railway services are usually planned well in advance [2, 19]. Therefore, the schedules are usually fixed in the tactical network design phase (see, e.g., Wieberneit [27], Andersen et al. [28], and SteadieSeifi et al. [29] for this area) and have to be accepted in the operative multimodal routing considered in this paper (Liu et al. [30], Lin [31], Demir et al. [20], and Hrušovský et al. [14]).The routes can be optimized according to different criteria. Traditionally, the optimization models concentrate on costs or time [32]. However, various surveys show that multiple criteria, among others also sustainability and reliability, are important for transportation planners [33, 34]. In case of sustainability, vast majority of available literature only considers economic and environmental criteria whereas social aspects of transportation are neglected [35, 36]. As a result, the comparison of transportation costs and CO2 emissions is decisive for finding the optimal routes.The CO2 emissions can be either integrated into the cost objective function as emission charges or modelled as an independent objective in multiobjective optimization. The majority of the articles on optimizing CO2 emissions in transportation adopts the first approach (e.g., Chang et al. [37], Zhang et al. [21], Wen et al. [38], Sun and Lang [2], Chen et al. [39], Demir et al. [20], and Hrušovský et al. [14]), whereas only a few studies used the second approach (e.g., Sun and Chen [40] and Qu et al. [41]). However, it is difficult to decide which method is better to use since the charging method is easier to implement but might not always be applicable, especially in cases where the emission charges are too low in comparison to transportation costs and therefore do not have any influence on the resulting routes [22].Another important decision is the choice of the emission calculation method since most of the models require a lot of parameters (e.g., fuel type, vehicle and route characteristics, and traffic condition) that cannot be always estimated in advance [20, 42–45]. As an alternative, the activity-based method that multiplies activity intensity (unit: TEU-km) with the CO2 emission factor (unit: g/TEU-km) to calculate the total CO2 emissions [5] yields a better feasibility in transportation practice and has been already applied in various studies [2, 5, 37, 39, 41, 46].The reliability of transportation is influenced by various sources of uncertainty (e.g., weather, accidents, congestion, and technical issues) that can cause disruptions in transportation processes. Despite this fact, most of the existing models consider all factors as deterministic and do not take into account uncertainty [47]. In order to improve the reliability of transportation plans, the uncertainty should be considered during the planning phase. Despite the variety of possible reasons for disruption, their effects influence either the demand, the travel time, or the supply (capacity) of transportation. Whereas extensive literature about demand uncertainty in multimodal routing is available (see, e.g., Lium et al. [48] and Ho et al. [49]), the consideration of travel time and capacity uncertainty is rather limited [50, 51]. Moreover, usually the models deal with only one source of uncertainty. However, the consideration of multiple uncertainties at the same time might lead to improved routes and therefore in our study the travel time uncertainty for road services and capacity uncertainty for rail services should be combined.Road traffic congestion leads to travel time uncertainty of road services. Only a few studies, for example, Demir et al. [20], Sun and Chen [40], and Uddin and Huynh [11], adopted stochastic programming to formulate the multimodal service network problem under disruptions. However, in practical planning, due to the lack of reliable historical data [52], it is quite difficult to fit the probability distributions to the travel times of all road services in the multimodal service network. Such disadvantage of stochastic programming in addressing the uncertainty is emphasized by Zarandi et al. [52] in a location-routing problem. Furthermore, the fuzzy programming that uses fuzzy variables (usually triangular or trapezoidal numbers) to evaluate the uncertainty fails to represent the variation of the travel time of a road service en route with respect to the time of the day [53]. Therefore, it is worthwhile to find an applicable approach (e.g., time-dependent travel time that is widely discussed in the vehicle routing problems [53, 54]) to formulate the road traffic congestion and introduce it to the multimodal routing modelling.The capacity uncertainty has only been considered by a few studies in different transportation planning problems, for example, multiobjective solid minimal cost flow problem by Wang [55], multiobjective multi-item solid transportation problem by Kaur and Kumar [56], and resource allocation problem for containerized cargo transportation networks by Kundu et al. [57]. Whereas Wang [55] used stochastic programming to address this issue, Kaur and Kumar [56] and Kundu et al. [57] formulated the capacity uncertainty from a fuzzy programming viewpoint. The fuzzy programming approach has been also used in the context of vehicle routing problem [52, 58] and shortest path problem [59], since it can avoid the infeasibility of the stochastic programming in case of lack of historical data [58]. However, to the best of our knowledge, there exists no study which would consider capacity uncertainty in multimodal routing problem.The considered aspects of the multimodal routing approach increase the complexity of the problem. The modelling of constraints for schedule-based services leads to the nonlinearity of the model, which can be seen in the models proposed by Liu et al. [30], Lin [31], and Chang [60]. Moreover, the modelling of uncertainty sources might also lead to nonlinearity of the model. As for the problem formulated by a nonlinear model, Xie et al. [61] pointed out that the nonlinearity makes the problem very difficult to solve, while using linearization technique to get the equivalent linear model can enable the problem to be relatively easily solved. The linearization technique should be designed according to the specific nonlinearity of the nonlinear model. Its main idea is to introduce new linear forms to replace the nonlinear components in the model equivalently.As the literature review shows, the interest of researchers in multimodal routing has been increasing during the last years, considering different aspects of the studied planning problem. However, there still exists potential for improvement especially in the area of sustainability and reliability of transportation. Therefore, our paper contributes to these aspects of the multimodal routing since it introduces and compares two different methods for modelling CO2 emissions. Moreover, two different sources of uncertainty (travel time and capacity) are considered within the planning process which should reduce the influence of potential disruptions on created transportation plans and improve the reliability. In case of rail service capacity, the fuzzy programming approach is used in order to better represent the possible capacity variations. Although the application of the model to the Chinese case study requires some specific assumptions, the general model can be used for any multimodal service network.In order to give an overview of the current state of the literature, the most important articles with the covered aspects are shown in Table1.Table 1
Literature review of the articles in the multimodal routing problem.
Articles
Due date consideration
Rail schedule constraints
CO2 emission optimization
Road traffic congestion
Rail service capacity
Time point
Time window
Yes
No
Charging method
Bi-objective optimization
Yes
No
Uncapacitated
Capacitated
Chang (2008) [60]
American scenario
★
Deterministic
Sun et al. (2008) [23]
★
★
Deterministic
Jiang and Lu (2009) [24]
Hard
★
★
★
Çakır (2009) [62]
★
★
Deterministic
Cai et al. (2010) [25]
★
★
Deterministic
Chang et al. (2010) [37]
★
★
★
Deterministic
Liu et al. (2011) [30]
Hard
Scheduled departure time
★
Deterministic
Zhang et al. (2011) [21]
Soft
★
★
★
★
Moccia et al. (2011) [63]
Italian scenario
★
Deterministic
Ayar and Yaman (2012) [64]
Hard
★
★
Deterministic
Verma et al. (2012) [17]
Hard
★
★
Deterministic
Sun and Chen (2013) [40]
★
★
★
★
Wen et al. (2013) [38]
Hard
★
★
★
Deterministic
Xiong and Wang (2014) [26]
★
★
★
Demir et al. (2015) [20]
Soft
Scheduled departure time
★
Stochastic
Deterministic
Lin (2015) [31]
Scheduled departure time
★
Deterministic
Sun and Lang (2015) [2]
Hard
Chinese scenario
★
★
Deterministic
Uddin and Huynh (2016) [11]
Hard
★
Stochastic
Deterministic
Hrušovský et al. [14]
Soft
Scheduled departure time
★
Stochastic
Deterministic
Qu et al. (2016) [41]
★
★
★
Deterministic
Sun et al. (2016) [19]
Hard
Chinese scenario
★
Deterministic
Zhang et al. (2017) [22]
Hard
★
★
★
★


## 3. Modelling Methodology

As already mentioned, the objective of the multimodal routing problem in this study is to select the best green routes to move containers from multiple transportation orders through the multimodal service network with road traffic congestion and rail service capacity uncertainty. The multimodal service network is composed of nodes and arcs, where nodes represent the terminals serving as origins and destinations of transportation orders as well as transshipment points between rail and road services. The arcs represent the services connecting the nodes. Related information on road and rail services, such as schedules, transportation distance, and capacity, is regarded as the parameters of the nodes and arcs.In the studied network, rail services are the backbone for multimodal container transportation, while road services play an important role in origin pickup and destination delivery. In some cases, road services also take direct origin-to-destination transportation. Figure1 shows the basic topological structure of the multimodal service network. Overall, the network structure is similar to a hub-and-spoke network with multiple hubs where origins and destinations represent spokes while railway container stations represent hubs.Figure 1
Basic topological structure of the multimodal service network.


### 3.1. Modelling Transportation Services

Service characteristics are dependent on the chosen transportation mode. In case of road service, the trucks are assumed to be fully flexible which means that they are available every time they are needed and it is possible to assign as many trucks as necessary to carry the containers. As a result, the road services are uncapacitated and their departure time can be adjusted according to the customer requirements and the current state of the network. This means that the departure time can be either immediately after the containers are loaded on the vehicle, but can be also postponed in order to avoid traffic congestion. However, it has to be ensured that the truck arrives to its destination before the loading cutoff time of the next planned service in order to avoid infeasibility of the transportation plan. Besides that, too early arrivals should be also avoided as inventory costs can be charged for containers waiting for the next service at a terminal. Therefore, it is important to determine the planned departure time of a truck from the current node.In contrast to that, rail services are modelled as block container trains driving between two nodes according to fixed schedules with the following characteristics:(1)
Operation time window at the loading/unloading organization station that regulates the operation start time and cutoff time regarding loading/unloading containers on/off the train
(2)
Scheduled departure/arrival time from/at the loading/unloading organization station
(3)
Operation period. For the convenience of modelling, same rail services in different operation periods are considered as different services.Since the loading/unloading operations of both road and rail services are highly mechanized, the time needed for loading/unloading is relatively low. Consequently, we assume that these times can be neglected in the multimodal routing model. Besides that, the capacities of rail services are limited. Taking into account these assumptions, the origin-to-destination transportation process combining road and rail services is depicted in Figure2 and can be described as follows:
(1)
The containers will first depart from the origin (node 1) at the planned departure time that is not earlier than the earliest release time defined by the customer, and then arrive at the loading organization station (node 2).
(2)
In order to be able to load the containers on the train operated from node 2 to node 3, their arrival time at node 2 should be no later than the operation start time. If the arrival time is earlier than the operation start time, there will be an inventory period from the arrival time to the operation start time at node 2, and the containers should wait until the operation start time and then get loaded. Moreover, there exists a free-of-charge period for inventory costs [64]. Only when the inventory period is longer than this period will inventory costs be created.
(3)
When loaded on the train, the containers will wait until the scheduled departure time and then leave node 2 along with the train.
(4)
The containers will arrive at the unloading organization station (node 3) at the scheduled arrival time of the train. However, they should wait until the operation start time and then get unloaded. Therefore, the effective arrival time of the containers at node 3 is the operation start time instead of the scheduled arrival time.
(5)
After unloading from the train, the transportation organizer can select a planned departure time that is no earlier than the operation start time to move containers to their destination (node 4). In case of an early arrival to the destination, the containers have to be stored again until the planned delivery date in a warehouse provided by the third-party inventory companies that can have different conditions in comparison to warehouses used in intermodal terminals [65].Figure 2
Diagram of a simple road-rail multimodal route.


### 3.2. Modelling Uncertainty

Considering the uncertainty factors, traffic congestion is an important factor in road transportation that influences the travel times. Congestion can occur either regularly due to insufficient road capacity or unexpectedly due to, for example, an accident. Since the latter one is hard to predict, only regularly occurring congestion is considered in this study. However, this congestion can vary with the time of the day or the section considered, therefore piecewise linear functions can be used to model such situations [53, 54]. By combining limited historical data and experiences from the experts, such travel time functions can be established. Figure 3 shows an example of the time-dependent travel time function. There are two peak hours from 9 am to 10 am and from 6 pm to 8 pm on the arc served by road services.Figure 3
Time-dependent travel time formulated by piecewise linear function.Besides the road congestion, the available capacity of a rail service can be also uncertain due to various reasons (e.g., order cancellation, vehicle breakdown, and lack of available vehicles). Since the reasons are manifold and the historical data about capacity variation might not be available, using the stochastic programming for capacity uncertainty might not be feasible. As a consequence, capacity uncertainty can be effectively expressed in a fuzzy way by using the triangular fuzzy variable represented by the minimum of the possible capacities, most possible capacity and maximum of the possible capacities. Similar to the traffic congestion issue, fuzzy capacity can be attained by limited historical data and expert experiences.


### 3.3. Modelling and Optimizing the CO2 Emissions

As stated in Section2, there are two approaches to optimize the CO2 emissions, including the CO2 emission charging method and the bi-objective optimization method. The optimization model is firstly constructed in Section 4 using the charging method, and the bi-objective method is then introduced in Section 6 where both methods are compared. For calculating the CO2 emissions, the activity-based method is used. This method calculates the total CO2 emissions for a transportation service using the formula (EM·Q·L) where EM is the CO2 emission factor of the service (unit: g/TEU-km), Q is the number of the containers carried by the transportation service (unit: TEU), and L is the corresponding transportation distance (unit: km).


### 3.4. Modelling the Customer Demands

Transportation orders are considered to be unsplittable which means that all containers of an order have to use the same route. The release times of orders are fixed whereas the delivery can be done within a due date time window limited by a lower and upper bound. For each transportation order, it is best if the arrival time of the containers falls within the time window. In this case, there will be no extra costs at the destination. However, if the containers arrive at the destination earlier than the lower bound, inventory costs will be created, because the goods need to be stored before being processed at the destination. If the arrival time is later than the upper bound, penalty costs should be charged for the delayed delivery.


## 3.1. Modelling Transportation Services

Service characteristics are dependent on the chosen transportation mode. In case of road service, the trucks are assumed to be fully flexible which means that they are available every time they are needed and it is possible to assign as many trucks as necessary to carry the containers. As a result, the road services are uncapacitated and their departure time can be adjusted according to the customer requirements and the current state of the network. This means that the departure time can be either immediately after the containers are loaded on the vehicle, but can be also postponed in order to avoid traffic congestion. However, it has to be ensured that the truck arrives to its destination before the loading cutoff time of the next planned service in order to avoid infeasibility of the transportation plan. Besides that, too early arrivals should be also avoided as inventory costs can be charged for containers waiting for the next service at a terminal. Therefore, it is important to determine the planned departure time of a truck from the current node.In contrast to that, rail services are modelled as block container trains driving between two nodes according to fixed schedules with the following characteristics:(1)
Operation time window at the loading/unloading organization station that regulates the operation start time and cutoff time regarding loading/unloading containers on/off the train
(2)
Scheduled departure/arrival time from/at the loading/unloading organization station
(3)
Operation period. For the convenience of modelling, same rail services in different operation periods are considered as different services.Since the loading/unloading operations of both road and rail services are highly mechanized, the time needed for loading/unloading is relatively low. Consequently, we assume that these times can be neglected in the multimodal routing model. Besides that, the capacities of rail services are limited. Taking into account these assumptions, the origin-to-destination transportation process combining road and rail services is depicted in Figure2 and can be described as follows:
(1)
The containers will first depart from the origin (node 1) at the planned departure time that is not earlier than the earliest release time defined by the customer, and then arrive at the loading organization station (node 2).
(2)
In order to be able to load the containers on the train operated from node 2 to node 3, their arrival time at node 2 should be no later than the operation start time. If the arrival time is earlier than the operation start time, there will be an inventory period from the arrival time to the operation start time at node 2, and the containers should wait until the operation start time and then get loaded. Moreover, there exists a free-of-charge period for inventory costs [64]. Only when the inventory period is longer than this period will inventory costs be created.
(3)
When loaded on the train, the containers will wait until the scheduled departure time and then leave node 2 along with the train.
(4)
The containers will arrive at the unloading organization station (node 3) at the scheduled arrival time of the train. However, they should wait until the operation start time and then get unloaded. Therefore, the effective arrival time of the containers at node 3 is the operation start time instead of the scheduled arrival time.
(5)
After unloading from the train, the transportation organizer can select a planned departure time that is no earlier than the operation start time to move containers to their destination (node 4). In case of an early arrival to the destination, the containers have to be stored again until the planned delivery date in a warehouse provided by the third-party inventory companies that can have different conditions in comparison to warehouses used in intermodal terminals [65].Figure 2
Diagram of a simple road-rail multimodal route.


## 3.2. Modelling Uncertainty

Considering the uncertainty factors, traffic congestion is an important factor in road transportation that influences the travel times. Congestion can occur either regularly due to insufficient road capacity or unexpectedly due to, for example, an accident. Since the latter one is hard to predict, only regularly occurring congestion is considered in this study. However, this congestion can vary with the time of the day or the section considered, therefore piecewise linear functions can be used to model such situations [53, 54]. By combining limited historical data and experiences from the experts, such travel time functions can be established. Figure 3 shows an example of the time-dependent travel time function. There are two peak hours from 9 am to 10 am and from 6 pm to 8 pm on the arc served by road services.Figure 3
Time-dependent travel time formulated by piecewise linear function.Besides the road congestion, the available capacity of a rail service can be also uncertain due to various reasons (e.g., order cancellation, vehicle breakdown, and lack of available vehicles). Since the reasons are manifold and the historical data about capacity variation might not be available, using the stochastic programming for capacity uncertainty might not be feasible. As a consequence, capacity uncertainty can be effectively expressed in a fuzzy way by using the triangular fuzzy variable represented by the minimum of the possible capacities, most possible capacity and maximum of the possible capacities. Similar to the traffic congestion issue, fuzzy capacity can be attained by limited historical data and expert experiences.


## 3.3. Modelling and Optimizing the CO2 Emissions

As stated in Section2, there are two approaches to optimize the CO2 emissions, including the CO2 emission charging method and the bi-objective optimization method. The optimization model is firstly constructed in Section 4 using the charging method, and the bi-objective method is then introduced in Section 6 where both methods are compared. For calculating the CO2 emissions, the activity-based method is used. This method calculates the total CO2 emissions for a transportation service using the formula (EM·Q·L) where EM is the CO2 emission factor of the service (unit: g/TEU-km), Q is the number of the containers carried by the transportation service (unit: TEU), and L is the corresponding transportation distance (unit: km).


## 3.4. Modelling the Customer Demands

Transportation orders are considered to be unsplittable which means that all containers of an order have to use the same route. The release times of orders are fixed whereas the delivery can be done within a due date time window limited by a lower and upper bound. For each transportation order, it is best if the arrival time of the containers falls within the time window. In this case, there will be no extra costs at the destination. However, if the containers arrive at the destination earlier than the lower bound, inventory costs will be created, because the goods need to be stored before being processed at the destination. If the arrival time is later than the upper bound, penalty costs should be charged for the delayed delivery.


## 4. Mathematical Model for the Multimodal Routing Problem

### 4.1. Notations

In this study,G=N,A,S represents a road-rail multimodal service network in form of a directed graph, where N, A, and S separately denote the node set, arc set, and transportation service set in the network. For each transportation order k (∀k∈K where K is the transportation order set), its five attributes are represented as origin ok, destination dk, volume qk (unit: TEU), earliest release time treleasek at origin, and due date tduek=tkL,tkU where tkL and tkU are separately the lower bound and upper bound of the due date time window, respectively. The rest of the symbols in the mathematical model and their respective representations are summarized as follows in Table 2.Table 2
Notations.
Indices
Representations
i, j, h
Indices of the nodes,i, j, and h∈N.
s, r
Indices of the transportation services,s and r∈S.
h,i
Directed arc from nodeh to node i, and h,i∈A.
i,j
Directed arc from nodei to node j, and i,j∈A.
Sets
Representations
Γij
Set of the rail services on arci,j.
Ωij
Set of the road services on arci,j.
Sij
Set of the transportation services on arci,j, Sij=Γij∪Ωij, and Sij⊆S.
δ−i
Set of the predecessor nodes to nodei, and δ−i⊆N.
δ+i
Set of the successor nodes to nodei, and δ+i⊆N.
Variables
Representations
wijsk
Charged inventory time for transportation orderk at node i before its containers are moved on arc i,j by service s, unit: h.
xijsk
A binary variable.xijsk=1 if service s on arc i,j is used for transportation order k; xijsk=0 otherwise.
yik
Arrival time of the containers of transportation orderk at node i.
zik
Planned departure time of the containers of transportation orderk from node i.
nik
A nonnegative integer variable that indicates how many periods (24 hours)zik exceed.
tijsk
Travel time of road services on arc i,j when used for transportation order k, unit: h. tijsk=fijszik where fijs is the piecewise linear function of the travel time with respect to the corresponding planned departure time. Note that in this function, the input departure time zik should fall within range [0, 24]. Otherwise, it should be first normalized into such range before being input into the function.
Parameters
Representations
lis
Loading/unloading operation start time of rail services at node i.
uis
Loading/unloading operation cutoff time of rail services at node i.
v~ijs
Fuzzy available capacity of rail services on arc i,j and v~ijs=vijsmin,vijsM,vijsmax where vijsmin<vijsM<vijsmax, unit: TEU.
dijs
Travel distance of services on arc i,j, unit: km.
cijs
Unit transportation cost when services is used to move containers on arc i,j, unit: ¥/TEU.
cs
Unit loading/unloading operation costs of services, unit: ¥/TEU.
cis
Unit inventory costs of transportation services at node i, unit: ¥/TEU-h.
cdestk
Unit destination inventory costs for transportation orderk, unit: ¥/TEU-h.
cpenk
Unit penalty costs for transportation orderk caused by delayed delivery at its destination, unit: ¥/h.
cco2
Unit CO2 emission costs, unit: ¥/g.
ems
CO2 emission factor for service s, unit: g/TEU-km.
τis
Free-of-charge inventory period of services at node i, unit: h.
M
A significant large positive number.
α
Confidence in the fuzzy chance constraint.


### 4.2. A Fuzzy Chance-Constrained MINLP Model

The framework of the mathematical model is derived from the classical arc-node-based formulation that has been widely used in the routing problems of various research fields, for example, telecommunication systems, manufacturing systems, and internet service systems [1]. The mathematical model for the multimodal routing problem explored in this study is presented as follows. The proposed model is a combination of the fuzzy chance-constrained programming and the mixed integer nonlinear programming (MINLP).Objective function:(1)minimize∑k∈K∑i,j∈A∑s∈Sijcijsk·qk·xijsk,(2)+∑k∈K∑i∈N∑h∈δ−i∑r∈Shicr·qk·xhirk+∑j∈δ+i∑s∈Sijcs·qk·xijsk,(3)+∑k∈K∑i,j∈A∑s∈Sijcis·qk·wijsk+∑k∈Kcdestk·qk·maxtkL−ydk,0,(4)+∑k∈Kcpenk·maxydk−tkU,0,(5)+∑k∈K∑i,j∈A∑s∈Sijcco2·ems·qk·dijs·xijsk,Eqs. (1), (2), (3), (4), and (5) are successively the transportation costs en route, loading and unloading costs at the nodes, inventory costs, penalty costs caused by delayed deliveries at destinations, and CO2 emission costs. Their linear summation represents the generalized costs. The objective function minimizes the generalized costs for multiple transportation orders using multimodal transportation. This setting reflects the objective of satisfying customer demands and minimizing monetary expenditures. Therefore, the multimodal routing in this study is a kind of service-oriented planning.Subject to(6)∑j∈δ+i∑s∈Sijxijsk−∑h∈δ−i∑r∈Shixhirk=1i=ok0∀i∈N\ok,dk−1i=dk∀k∈K,∀i∈N,(7)∑s∈Sijxijsk≤1∀k∈K,∀i,j∈A.Eq. (6) is the flow conservation constraint which ensures the continuity of the nodes and arcs on the planned routes, which is a common constraint in the node-arc-based formulation [1]. Eq. (7) means that no more than one transportation service on one arc can be utilized to serve a transportation order. The combination of the two equations above ensures that containers in a transportation order are unsplittable, which is claimed in Section 3.4.
(8)Cr∑k∈Kqk·xijsk≤v~ijs≥α∀i,j∈A,∀s∈Γij.Eq. (8) is the fuzzy chance constraint of the rail service capacity based on the fuzzy credibility measure [66]. In this equation, Cr· represents the credibility of the event in ·. Eq. (8) ensures that the credibility of the event that the entire loaded container volume on a rail service does not exceed its available capacity should not be lower than a predetermined confidence α∈0,1. The value of α indicates the decision makers’ preference for addressing the corresponding decision-making.
(9)yokk=treleasek∀k∈K.Eq. (9) assumes that the arrival time of the containers of transportation order k at the origin equals its earliest release time.
(10)yik≤zik∀k∈K,∀i∈N\dk,(11)nik=zik24∀k∈K,∀i∈N\dk,(12)tijsk=fijszik−24·nik∀k∈K,i,j∈A,∀s∈Ωij,(13)zik+tijsk−yjk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij.Eqs. (10), (11), (12), and (13) compute the arrival times of the containers of transportation order k at the nodes on the route by road services. First, we determine the planned departure time of the containers at the current node i by (10). Then, we normalize such departure time into range [0, 24] by using the floor integral function (11) and normalization formula zik−24·nik, so that we can further generate the corresponding travel time tijsk of the containers on arc i,j by road service s according to the piecewise linear function (12). Finally, we can gain the arrival time of container flow k at the following node j by (13).
(14)lir−yik·xhirk=0∀k∈K,∀h,i∈A,∀r∈Γhi.Eq. (14) computes the arrival times of the containers of transportation order k at the nodes on the route by rail services. As claimed in Section 3.1, lir is the effective arrival time of containers at node i. Because we assume that loading/unloading time can be neglected, yik is also the earliest departure time of containers at transshipping node.
(15)yik≤uis·xijsk+M·1−xijsk∀k∈K,∀i,j∈A,∀s∈Γij.Eq. (15) is the rail service operation cutoff time constraint. It ensures that the arrival time of the containers of transportation order k at node i should not exceed the operation cutoff time of rail service s at the same node, if this rail service has to be used in order to move these containers on arc i,j.
(16)maxlis−yik−τis,0−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Γij,(17)maxzik−yik−τis,0−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij.Eqs. (16) and (17) compute the charged inventory time of the containers of transportation order k at node i before being moved on arc i,j by a rail service and a road service, respectively. In (17), especially if i=ok and s∈Ωij, τis=0.
(18)nik∈0,1,2,…∀k∈K,∀i∈N,(19)tijsk≥0∀k∈K,i,j∈A,∀s∈Ωij,(20)wijsk≥0∀k∈K,∀i,j∈A,∀s∈Sij,(21)xijsk∈0,1∀k∈K,∀i,j∈A,∀s∈Sij,(22)yik≥0∀k∈K,∀i∈N,(23)zik≥0∀k∈K,∀i∈N\dk.Eqs. (18), (19), (20), (21), (22), and (23) are the variable domain constraints.The above mathematical model for the multimodal routing problem is summarized in Figure4. Figure 4 classifies all equations (excluding variable domain constraints) into three categories, including “Modelling the network” which realizes the methodology proposed in Section 3, “Setting the objective” which reflects that satisfying customer demands is the core of the routing problem, and “Basic assumptions” which are the basic conditions of the model. This should help the readers to better understand the functions and roles of various equations in the mathematical model.Figure 4
Framework of the mathematical model.By using the solution strategy designed in the next section, the proposed model can be used to optimize any green multimodal routing problems that are composed of scheduled rail services with uncertain capacities and time-flexible road services with uncertain travel times. It can plan the best multimodal routes to accomplish multiple transportation orders with soft due date time windows.


## 4.1. Notations

In this study,G=N,A,S represents a road-rail multimodal service network in form of a directed graph, where N, A, and S separately denote the node set, arc set, and transportation service set in the network. For each transportation order k (∀k∈K where K is the transportation order set), its five attributes are represented as origin ok, destination dk, volume qk (unit: TEU), earliest release time treleasek at origin, and due date tduek=tkL,tkU where tkL and tkU are separately the lower bound and upper bound of the due date time window, respectively. The rest of the symbols in the mathematical model and their respective representations are summarized as follows in Table 2.Table 2
Notations.
Indices
Representations
i, j, h
Indices of the nodes,i, j, and h∈N.
s, r
Indices of the transportation services,s and r∈S.
h,i
Directed arc from nodeh to node i, and h,i∈A.
i,j
Directed arc from nodei to node j, and i,j∈A.
Sets
Representations
Γij
Set of the rail services on arci,j.
Ωij
Set of the road services on arci,j.
Sij
Set of the transportation services on arci,j, Sij=Γij∪Ωij, and Sij⊆S.
δ−i
Set of the predecessor nodes to nodei, and δ−i⊆N.
δ+i
Set of the successor nodes to nodei, and δ+i⊆N.
Variables
Representations
wijsk
Charged inventory time for transportation orderk at node i before its containers are moved on arc i,j by service s, unit: h.
xijsk
A binary variable.xijsk=1 if service s on arc i,j is used for transportation order k; xijsk=0 otherwise.
yik
Arrival time of the containers of transportation orderk at node i.
zik
Planned departure time of the containers of transportation orderk from node i.
nik
A nonnegative integer variable that indicates how many periods (24 hours)zik exceed.
tijsk
Travel time of road services on arc i,j when used for transportation order k, unit: h. tijsk=fijszik where fijs is the piecewise linear function of the travel time with respect to the corresponding planned departure time. Note that in this function, the input departure time zik should fall within range [0, 24]. Otherwise, it should be first normalized into such range before being input into the function.
Parameters
Representations
lis
Loading/unloading operation start time of rail services at node i.
uis
Loading/unloading operation cutoff time of rail services at node i.
v~ijs
Fuzzy available capacity of rail services on arc i,j and v~ijs=vijsmin,vijsM,vijsmax where vijsmin<vijsM<vijsmax, unit: TEU.
dijs
Travel distance of services on arc i,j, unit: km.
cijs
Unit transportation cost when services is used to move containers on arc i,j, unit: ¥/TEU.
cs
Unit loading/unloading operation costs of services, unit: ¥/TEU.
cis
Unit inventory costs of transportation services at node i, unit: ¥/TEU-h.
cdestk
Unit destination inventory costs for transportation orderk, unit: ¥/TEU-h.
cpenk
Unit penalty costs for transportation orderk caused by delayed delivery at its destination, unit: ¥/h.
cco2
Unit CO2 emission costs, unit: ¥/g.
ems
CO2 emission factor for service s, unit: g/TEU-km.
τis
Free-of-charge inventory period of services at node i, unit: h.
M
A significant large positive number.
α
Confidence in the fuzzy chance constraint.


## 4.2. A Fuzzy Chance-Constrained MINLP Model

The framework of the mathematical model is derived from the classical arc-node-based formulation that has been widely used in the routing problems of various research fields, for example, telecommunication systems, manufacturing systems, and internet service systems [1]. The mathematical model for the multimodal routing problem explored in this study is presented as follows. The proposed model is a combination of the fuzzy chance-constrained programming and the mixed integer nonlinear programming (MINLP).Objective function:(1)minimize∑k∈K∑i,j∈A∑s∈Sijcijsk·qk·xijsk,(2)+∑k∈K∑i∈N∑h∈δ−i∑r∈Shicr·qk·xhirk+∑j∈δ+i∑s∈Sijcs·qk·xijsk,(3)+∑k∈K∑i,j∈A∑s∈Sijcis·qk·wijsk+∑k∈Kcdestk·qk·maxtkL−ydk,0,(4)+∑k∈Kcpenk·maxydk−tkU,0,(5)+∑k∈K∑i,j∈A∑s∈Sijcco2·ems·qk·dijs·xijsk,Eqs. (1), (2), (3), (4), and (5) are successively the transportation costs en route, loading and unloading costs at the nodes, inventory costs, penalty costs caused by delayed deliveries at destinations, and CO2 emission costs. Their linear summation represents the generalized costs. The objective function minimizes the generalized costs for multiple transportation orders using multimodal transportation. This setting reflects the objective of satisfying customer demands and minimizing monetary expenditures. Therefore, the multimodal routing in this study is a kind of service-oriented planning.Subject to(6)∑j∈δ+i∑s∈Sijxijsk−∑h∈δ−i∑r∈Shixhirk=1i=ok0∀i∈N\ok,dk−1i=dk∀k∈K,∀i∈N,(7)∑s∈Sijxijsk≤1∀k∈K,∀i,j∈A.Eq. (6) is the flow conservation constraint which ensures the continuity of the nodes and arcs on the planned routes, which is a common constraint in the node-arc-based formulation [1]. Eq. (7) means that no more than one transportation service on one arc can be utilized to serve a transportation order. The combination of the two equations above ensures that containers in a transportation order are unsplittable, which is claimed in Section 3.4.
(8)Cr∑k∈Kqk·xijsk≤v~ijs≥α∀i,j∈A,∀s∈Γij.Eq. (8) is the fuzzy chance constraint of the rail service capacity based on the fuzzy credibility measure [66]. In this equation, Cr· represents the credibility of the event in ·. Eq. (8) ensures that the credibility of the event that the entire loaded container volume on a rail service does not exceed its available capacity should not be lower than a predetermined confidence α∈0,1. The value of α indicates the decision makers’ preference for addressing the corresponding decision-making.
(9)yokk=treleasek∀k∈K.Eq. (9) assumes that the arrival time of the containers of transportation order k at the origin equals its earliest release time.
(10)yik≤zik∀k∈K,∀i∈N\dk,(11)nik=zik24∀k∈K,∀i∈N\dk,(12)tijsk=fijszik−24·nik∀k∈K,i,j∈A,∀s∈Ωij,(13)zik+tijsk−yjk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij.Eqs. (10), (11), (12), and (13) compute the arrival times of the containers of transportation order k at the nodes on the route by road services. First, we determine the planned departure time of the containers at the current node i by (10). Then, we normalize such departure time into range [0, 24] by using the floor integral function (11) and normalization formula zik−24·nik, so that we can further generate the corresponding travel time tijsk of the containers on arc i,j by road service s according to the piecewise linear function (12). Finally, we can gain the arrival time of container flow k at the following node j by (13).
(14)lir−yik·xhirk=0∀k∈K,∀h,i∈A,∀r∈Γhi.Eq. (14) computes the arrival times of the containers of transportation order k at the nodes on the route by rail services. As claimed in Section 3.1, lir is the effective arrival time of containers at node i. Because we assume that loading/unloading time can be neglected, yik is also the earliest departure time of containers at transshipping node.
(15)yik≤uis·xijsk+M·1−xijsk∀k∈K,∀i,j∈A,∀s∈Γij.Eq. (15) is the rail service operation cutoff time constraint. It ensures that the arrival time of the containers of transportation order k at node i should not exceed the operation cutoff time of rail service s at the same node, if this rail service has to be used in order to move these containers on arc i,j.
(16)maxlis−yik−τis,0−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Γij,(17)maxzik−yik−τis,0−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij.Eqs. (16) and (17) compute the charged inventory time of the containers of transportation order k at node i before being moved on arc i,j by a rail service and a road service, respectively. In (17), especially if i=ok and s∈Ωij, τis=0.
(18)nik∈0,1,2,…∀k∈K,∀i∈N,(19)tijsk≥0∀k∈K,i,j∈A,∀s∈Ωij,(20)wijsk≥0∀k∈K,∀i,j∈A,∀s∈Sij,(21)xijsk∈0,1∀k∈K,∀i,j∈A,∀s∈Sij,(22)yik≥0∀k∈K,∀i∈N,(23)zik≥0∀k∈K,∀i∈N\dk.Eqs. (18), (19), (20), (21), (22), and (23) are the variable domain constraints.The above mathematical model for the multimodal routing problem is summarized in Figure4. Figure 4 classifies all equations (excluding variable domain constraints) into three categories, including “Modelling the network” which realizes the methodology proposed in Section 3, “Setting the objective” which reflects that satisfying customer demands is the core of the routing problem, and “Basic assumptions” which are the basic conditions of the model. This should help the readers to better understand the functions and roles of various equations in the mathematical model.Figure 4
Framework of the mathematical model.By using the solution strategy designed in the next section, the proposed model can be used to optimize any green multimodal routing problems that are composed of scheduled rail services with uncertain capacities and time-flexible road services with uncertain travel times. It can plan the best multimodal routes to accomplish multiple transportation orders with soft due date time windows.


## 5. An Exact Solution Strategy

### 5.1. Determining the Solution Strategy

The proposed mathematical model in Section4.2 is easily understood, because it is established straightforwardly according to the settings of the specific multimodal routing problem. However, the mathematical model contains three categories of nonlinearity, including (1) fuzzy chance constraint as (8), (2) piecewise linear constraint as (12), and (3) general nonlinear components as (3), (4), (13), (14), (16), and (17). Such nonlinearity increases the difficulty of the problem to be effectively solved by the exact solution algorithms (e.g., branch-and-bound algorithm and simplex algorithm). Especially for the large-scale real-world problem, due to the nonlinearity, the solution to the problem will usually get into the local optimum, and massive computational time will be consumed to generate the solution. Therefore, if we can eliminate the nonlinearity of the model by designing the equivalent linear reformulations to the nonlinear equations in the model, then the problem can be solved with the help of the mathematical programming software (e.g., LINGO, CPLEX, and GAMS) where exact solution algorithms can be implemented efficiently.We devote this study to developing a linearization-based exact solution strategy to address the multimodal routing problem. Although it is widely acknowledged that heuristic algorithms have better performance than the exact solution algorithms in dealing with the multimodal routing, which is known as an NP-hard problem, exact solution algorithms have two significant advantages. The first one is that they can test if the model observes the mathematical logic. The second is that they can provide an exact benchmark to verify the efficiency and accuracy of the heuristic algorithms [67]. Fazayeli et al. [68] also highlighted that the development of an exact algorithm and the comparison between the algorithms will be their future work for the location-routing problem on a multimodal transportation network.


### 5.2. Model Linearization Technique

#### 5.2.1. Crisp Equivalent and Linearization of the Fuzzy Chance Constraint

From the definition of the fuzzy creditability, (24) can be obtained where a is a deterministic number while b~ is a triangular fuzzy number described by b1,b2,b3 and b1<b2<b3. [69]:
(24)Cra≤b~=1,ifa≤b12b2−b1−a2b2−b1,ifb1≤a≤b2b3−a2b3−b2,ifb2≤a≤b30,ifa≥b3.According to (24), we can gain the crisp equivalency to the fuzzy chance constraint (8) by replacing a with ∑k∈Kqk·xijsk and b~ with v~ijs. However, the crisp equivalency is a piecewise linear function which is still essentially a nonlinear formula. Thus, we conduct following modifications shown as (25) and (26) to transform it into a pure linear function. After such modification, during the simulation process, we can use (25) to substitute (8) if α∈0.5,1, (26) to substitute the same equation otherwise.
(25)21−α·vijsM+2α−1·vijsmin≥∑k∈Kqk·xijskif0.5≤α≤1,∀i,j∈A,∀s∈Γij,(26)2α·vijsM−2α−1·vijsmax≥∑k∈Kqk·xijskif0≤α≤0.5,∀i,j∈A,∀s∈Γij.


#### 5.2.2. Linearization of the Piecewise Linear Function

Similarly to (24), the piecewise linear function fijs· is also a nonlinear expression where · represents zik−24·nik for the sake of conciseness. Let Pijs represent the set of the linear segments of fijs·, p the index of the segments (p=1,…,Pijs), and ttijskp−,ttijskp+ the departure time interval corresponding to the linear segment fijsp·. Then, fijs· has a universal formulation shown as (27).
(27)fijs·=fijs1·,·∈ttijs1−,ttijs1+⋮fijsp·,·∈ttijsp−,ttijsp+∀i,j∈A,∀s∈Ωij.⋮fijsPijs·,·∈ttijsPijs−,ttijsPijs+To linearize the piecewise linear function (27), the following procedures are conducted.Step 1.
Define a 0-1 variableφijskp. φijskp=1, if zik−24·nik falls within ttijsp−,ttijsp+. φijsp=0 otherwise. And φijsp should observe the following constraint.
(28)∑p∈Pijsφijskp=1∀k∈K,∀i,j∈A,∀s∈Ωij.Step 2.
Define two nonnegative variablesλijskp− and λijskp+, and distribute them to the lower bound ttijsp− and upper bound ttijsp+ of the corresponding departure time interval ttijsp−,ttijsp+ as their weights. Hence, we have the following equations.
(29)zik−24·nik=∑p∈Pijsλijsp−·ttijsp−+λijsp+·ttijsp+∀k∈K,∀i,j∈A,∀s∈Ωij,λijskp−+λijskp+=φijskp∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs,λijskp+≥0∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs,λijskp−≥0∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs.Step 3.
Modify the piecewise linear function as follows.(30)fijszik−24·nik=fijs∑p∈Pijsλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsfijsλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsfijspλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsλijskp−·fijspttijsp−+λijskp−·fijspttijsp+,where fijspttijsp− and fijspttijsp+ for ∀i,j∈A,∀s∈Ωij,∀p∈Pijs are all known. Finally, the pure linear equivalency to (12) is shown as (31).
(31)tijsk=∑p∈Pijsλijsp−·fijspttijsp−+λijsp−·fijspttijsp+∀k∈K,∀i,j∈A,∀s∈Ωij.


#### 5.2.3. Linearization of Other Nonlinear Components

(1) Linear Reformulation 1. Eq. (3) with nonlinear component ∑k∈Kcdestk·qk·maxtkL−ydk,0 can be linearized by (32), (33), and (34) where gk is a newly introduced variable.
(32)∑k∈K∑i,j∈A∑s∈Sijcis·qk·wijsk+∑k∈Kcdestk·qk·gk,(33)gk≥tkL−ydk∀k∈K,(34)gk≥0∀k∈K.If the containers arrive at the destination earlier thantkL, tkL−ydk>0. Consequently, maxtkL−ydk,0=tkL−ydk. According to (32) and (33), gk≥tkL−ydk. Minimization of (32) will then lead to gk=tkL−ydk. In this case, ∑k∈Kcdestk·qk·maxtkL−ydk,0 is equivalent to (32) and (33).If the containers arrive at the destination later thantkL, tkL−ydk<0. Consequently, maxtkL−ydk,0=0. According to (32) and (33), gk≥0. Minimization of (32) will then lead to gk=0. In this case, ∑k∈Kcdestk·qk·maxtkL−ydk,0 is also equivalent to (32) and (33).Above all, the equivalency above is verified. Similarly, we have the following linear reformulation.(2) Linear Reformulation 2. Nonlinear component ∑k∈Kcpenk·maxydk−tkU,0 in (4) can be linearized by (35), (36), and (37) where σk is also a newly introduced variable.
(35)∑k∈Kcpenk·σk,(36)σk≥ydk−tkU∀k∈K,(37)σk≥0∀k∈K.(3) Linear Reformulation 3. Nonlinear constraint nik=zik/24∀k∈K,∀i∈N\dk (12) can be reformulated by linear constraints (38) and (39).
(38)nik≤zik24∀k∈K,∀i∈N\dk,(39)nik>zik24−1∀k∈K,∀i∈N\dk.We takezik=60 as example to verify the equivalency above. According to (12), nik=60/24=2.5=2. Meanwhile, according to (45) and (46), there exists 1.5<nik≤2.5. Because nik is a nonnegative integer, finally nik=2. Consequently, the equivalency above is tenable.(4) Linear Reformulation 4. By using the same linearization technique, nonlinear constraints zik+tijsk−yjk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij (13) and lir−yik·xhirk=0∀k∈K,∀h,i∈A,∀r∈Γhi (14) can be reformulated by linear constraints (40) and (41) and (42) and (43), respectively.
(40)zik+tijsk−yjk≥M·xijsk−1∀k∈K,∀i,j∈A,∀s∈Ωij,(41)zik+tijsk−yjk≤M·1−xijsk∀k∈K,∀i,j∈A,∀s∈Ωij,(42)lir−yik≥M·xhirk−1∀k∈K,∀h,i∈A,∀r∈Γhi,(43)lir−yik≤M·1−xhirk∀k∈K,∀h,i∈A,∀r∈Γhi.(5) Linear Reformulation 5. By using the same linearization form, nonlinear constraints max0,lis−yik−τ−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Γij (16) and max0,zik−yik−τ−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij (17) can be reformulated by linear constraints (44) and (45) and (46) and (47).
(44)wijsk≥M·xijsk−1+lis−yik−τis∀k∈K,∀i,j∈A,∀s∈Γij,(45)wijsk≤M·xijsk∀k∈K,∀i,j∈A,∀s∈Γij,(46)wijsk≥M·xijsk−1+zik−yik−τis∀k∈K,∀i,j∈A,∀s∈Ωij,(47)wijsk≤M·xijsk∀k∈K,∀i,j∈A,∀s∈Ωij.For detailed proofs of Linear Reformations 4 and 5, readers can refer to our previous study [2].


### 5.3. Improved Linear Reformulation

For the sake of readability, the improved linear reformulation derived from the mathematical model in Section4.2 is summarized as follows.Objective function:Minimize (1)+(2)+(32)+(35)+(5)Subject to:Constraints (6), (7), (9), (10), (15), (18–23), (25), (26), (28), (29), (31), (33), (34), (36–47)


## 5.1. Determining the Solution Strategy

The proposed mathematical model in Section4.2 is easily understood, because it is established straightforwardly according to the settings of the specific multimodal routing problem. However, the mathematical model contains three categories of nonlinearity, including (1) fuzzy chance constraint as (8), (2) piecewise linear constraint as (12), and (3) general nonlinear components as (3), (4), (13), (14), (16), and (17). Such nonlinearity increases the difficulty of the problem to be effectively solved by the exact solution algorithms (e.g., branch-and-bound algorithm and simplex algorithm). Especially for the large-scale real-world problem, due to the nonlinearity, the solution to the problem will usually get into the local optimum, and massive computational time will be consumed to generate the solution. Therefore, if we can eliminate the nonlinearity of the model by designing the equivalent linear reformulations to the nonlinear equations in the model, then the problem can be solved with the help of the mathematical programming software (e.g., LINGO, CPLEX, and GAMS) where exact solution algorithms can be implemented efficiently.We devote this study to developing a linearization-based exact solution strategy to address the multimodal routing problem. Although it is widely acknowledged that heuristic algorithms have better performance than the exact solution algorithms in dealing with the multimodal routing, which is known as an NP-hard problem, exact solution algorithms have two significant advantages. The first one is that they can test if the model observes the mathematical logic. The second is that they can provide an exact benchmark to verify the efficiency and accuracy of the heuristic algorithms [67]. Fazayeli et al. [68] also highlighted that the development of an exact algorithm and the comparison between the algorithms will be their future work for the location-routing problem on a multimodal transportation network.


## 5.2. Model Linearization Technique

### 5.2.1. Crisp Equivalent and Linearization of the Fuzzy Chance Constraint

From the definition of the fuzzy creditability, (24) can be obtained where a is a deterministic number while b~ is a triangular fuzzy number described by b1,b2,b3 and b1<b2<b3. [69]:
(24)Cra≤b~=1,ifa≤b12b2−b1−a2b2−b1,ifb1≤a≤b2b3−a2b3−b2,ifb2≤a≤b30,ifa≥b3.According to (24), we can gain the crisp equivalency to the fuzzy chance constraint (8) by replacing a with ∑k∈Kqk·xijsk and b~ with v~ijs. However, the crisp equivalency is a piecewise linear function which is still essentially a nonlinear formula. Thus, we conduct following modifications shown as (25) and (26) to transform it into a pure linear function. After such modification, during the simulation process, we can use (25) to substitute (8) if α∈0.5,1, (26) to substitute the same equation otherwise.
(25)21−α·vijsM+2α−1·vijsmin≥∑k∈Kqk·xijskif0.5≤α≤1,∀i,j∈A,∀s∈Γij,(26)2α·vijsM−2α−1·vijsmax≥∑k∈Kqk·xijskif0≤α≤0.5,∀i,j∈A,∀s∈Γij.


### 5.2.2. Linearization of the Piecewise Linear Function

Similarly to (24), the piecewise linear function fijs· is also a nonlinear expression where · represents zik−24·nik for the sake of conciseness. Let Pijs represent the set of the linear segments of fijs·, p the index of the segments (p=1,…,Pijs), and ttijskp−,ttijskp+ the departure time interval corresponding to the linear segment fijsp·. Then, fijs· has a universal formulation shown as (27).
(27)fijs·=fijs1·,·∈ttijs1−,ttijs1+⋮fijsp·,·∈ttijsp−,ttijsp+∀i,j∈A,∀s∈Ωij.⋮fijsPijs·,·∈ttijsPijs−,ttijsPijs+To linearize the piecewise linear function (27), the following procedures are conducted.Step 1.
Define a 0-1 variableφijskp. φijskp=1, if zik−24·nik falls within ttijsp−,ttijsp+. φijsp=0 otherwise. And φijsp should observe the following constraint.
(28)∑p∈Pijsφijskp=1∀k∈K,∀i,j∈A,∀s∈Ωij.Step 2.
Define two nonnegative variablesλijskp− and λijskp+, and distribute them to the lower bound ttijsp− and upper bound ttijsp+ of the corresponding departure time interval ttijsp−,ttijsp+ as their weights. Hence, we have the following equations.
(29)zik−24·nik=∑p∈Pijsλijsp−·ttijsp−+λijsp+·ttijsp+∀k∈K,∀i,j∈A,∀s∈Ωij,λijskp−+λijskp+=φijskp∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs,λijskp+≥0∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs,λijskp−≥0∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs.Step 3.
Modify the piecewise linear function as follows.(30)fijszik−24·nik=fijs∑p∈Pijsλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsfijsλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsfijspλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsλijskp−·fijspttijsp−+λijskp−·fijspttijsp+,where fijspttijsp− and fijspttijsp+ for ∀i,j∈A,∀s∈Ωij,∀p∈Pijs are all known. Finally, the pure linear equivalency to (12) is shown as (31).
(31)tijsk=∑p∈Pijsλijsp−·fijspttijsp−+λijsp−·fijspttijsp+∀k∈K,∀i,j∈A,∀s∈Ωij.


### 5.2.3. Linearization of Other Nonlinear Components

(1) Linear Reformulation 1. Eq. (3) with nonlinear component ∑k∈Kcdestk·qk·maxtkL−ydk,0 can be linearized by (32), (33), and (34) where gk is a newly introduced variable.
(32)∑k∈K∑i,j∈A∑s∈Sijcis·qk·wijsk+∑k∈Kcdestk·qk·gk,(33)gk≥tkL−ydk∀k∈K,(34)gk≥0∀k∈K.If the containers arrive at the destination earlier thantkL, tkL−ydk>0. Consequently, maxtkL−ydk,0=tkL−ydk. According to (32) and (33), gk≥tkL−ydk. Minimization of (32) will then lead to gk=tkL−ydk. In this case, ∑k∈Kcdestk·qk·maxtkL−ydk,0 is equivalent to (32) and (33).If the containers arrive at the destination later thantkL, tkL−ydk<0. Consequently, maxtkL−ydk,0=0. According to (32) and (33), gk≥0. Minimization of (32) will then lead to gk=0. In this case, ∑k∈Kcdestk·qk·maxtkL−ydk,0 is also equivalent to (32) and (33).Above all, the equivalency above is verified. Similarly, we have the following linear reformulation.(2) Linear Reformulation 2. Nonlinear component ∑k∈Kcpenk·maxydk−tkU,0 in (4) can be linearized by (35), (36), and (37) where σk is also a newly introduced variable.
(35)∑k∈Kcpenk·σk,(36)σk≥ydk−tkU∀k∈K,(37)σk≥0∀k∈K.(3) Linear Reformulation 3. Nonlinear constraint nik=zik/24∀k∈K,∀i∈N\dk (12) can be reformulated by linear constraints (38) and (39).
(38)nik≤zik24∀k∈K,∀i∈N\dk,(39)nik>zik24−1∀k∈K,∀i∈N\dk.We takezik=60 as example to verify the equivalency above. According to (12), nik=60/24=2.5=2. Meanwhile, according to (45) and (46), there exists 1.5<nik≤2.5. Because nik is a nonnegative integer, finally nik=2. Consequently, the equivalency above is tenable.(4) Linear Reformulation 4. By using the same linearization technique, nonlinear constraints zik+tijsk−yjk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij (13) and lir−yik·xhirk=0∀k∈K,∀h,i∈A,∀r∈Γhi (14) can be reformulated by linear constraints (40) and (41) and (42) and (43), respectively.
(40)zik+tijsk−yjk≥M·xijsk−1∀k∈K,∀i,j∈A,∀s∈Ωij,(41)zik+tijsk−yjk≤M·1−xijsk∀k∈K,∀i,j∈A,∀s∈Ωij,(42)lir−yik≥M·xhirk−1∀k∈K,∀h,i∈A,∀r∈Γhi,(43)lir−yik≤M·1−xhirk∀k∈K,∀h,i∈A,∀r∈Γhi.(5) Linear Reformulation 5. By using the same linearization form, nonlinear constraints max0,lis−yik−τ−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Γij (16) and max0,zik−yik−τ−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij (17) can be reformulated by linear constraints (44) and (45) and (46) and (47).
(44)wijsk≥M·xijsk−1+lis−yik−τis∀k∈K,∀i,j∈A,∀s∈Γij,(45)wijsk≤M·xijsk∀k∈K,∀i,j∈A,∀s∈Γij,(46)wijsk≥M·xijsk−1+zik−yik−τis∀k∈K,∀i,j∈A,∀s∈Ωij,(47)wijsk≤M·xijsk∀k∈K,∀i,j∈A,∀s∈Ωij.For detailed proofs of Linear Reformations 4 and 5, readers can refer to our previous study [2].


## 5.2.1. Crisp Equivalent and Linearization of the Fuzzy Chance Constraint

From the definition of the fuzzy creditability, (24) can be obtained where a is a deterministic number while b~ is a triangular fuzzy number described by b1,b2,b3 and b1<b2<b3. [69]:
(24)Cra≤b~=1,ifa≤b12b2−b1−a2b2−b1,ifb1≤a≤b2b3−a2b3−b2,ifb2≤a≤b30,ifa≥b3.According to (24), we can gain the crisp equivalency to the fuzzy chance constraint (8) by replacing a with ∑k∈Kqk·xijsk and b~ with v~ijs. However, the crisp equivalency is a piecewise linear function which is still essentially a nonlinear formula. Thus, we conduct following modifications shown as (25) and (26) to transform it into a pure linear function. After such modification, during the simulation process, we can use (25) to substitute (8) if α∈0.5,1, (26) to substitute the same equation otherwise.
(25)21−α·vijsM+2α−1·vijsmin≥∑k∈Kqk·xijskif0.5≤α≤1,∀i,j∈A,∀s∈Γij,(26)2α·vijsM−2α−1·vijsmax≥∑k∈Kqk·xijskif0≤α≤0.5,∀i,j∈A,∀s∈Γij.


## 5.2.2. Linearization of the Piecewise Linear Function

Similarly to (24), the piecewise linear function fijs· is also a nonlinear expression where · represents zik−24·nik for the sake of conciseness. Let Pijs represent the set of the linear segments of fijs·, p the index of the segments (p=1,…,Pijs), and ttijskp−,ttijskp+ the departure time interval corresponding to the linear segment fijsp·. Then, fijs· has a universal formulation shown as (27).
(27)fijs·=fijs1·,·∈ttijs1−,ttijs1+⋮fijsp·,·∈ttijsp−,ttijsp+∀i,j∈A,∀s∈Ωij.⋮fijsPijs·,·∈ttijsPijs−,ttijsPijs+To linearize the piecewise linear function (27), the following procedures are conducted.Step 1.
Define a 0-1 variableφijskp. φijskp=1, if zik−24·nik falls within ttijsp−,ttijsp+. φijsp=0 otherwise. And φijsp should observe the following constraint.
(28)∑p∈Pijsφijskp=1∀k∈K,∀i,j∈A,∀s∈Ωij.Step 2.
Define two nonnegative variablesλijskp− and λijskp+, and distribute them to the lower bound ttijsp− and upper bound ttijsp+ of the corresponding departure time interval ttijsp−,ttijsp+ as their weights. Hence, we have the following equations.
(29)zik−24·nik=∑p∈Pijsλijsp−·ttijsp−+λijsp+·ttijsp+∀k∈K,∀i,j∈A,∀s∈Ωij,λijskp−+λijskp+=φijskp∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs,λijskp+≥0∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs,λijskp−≥0∀k∈K,∀i,j∈A,∀s∈Ωij,∀p∈Pijs.Step 3.
Modify the piecewise linear function as follows.(30)fijszik−24·nik=fijs∑p∈Pijsλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsfijsλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsfijspλijskp−·ttijsp−+λijskp+·ttijsp+=∑p∈Pijsλijskp−·fijspttijsp−+λijskp−·fijspttijsp+,where fijspttijsp− and fijspttijsp+ for ∀i,j∈A,∀s∈Ωij,∀p∈Pijs are all known. Finally, the pure linear equivalency to (12) is shown as (31).
(31)tijsk=∑p∈Pijsλijsp−·fijspttijsp−+λijsp−·fijspttijsp+∀k∈K,∀i,j∈A,∀s∈Ωij.


## 5.2.3. Linearization of Other Nonlinear Components

(1) Linear Reformulation 1. Eq. (3) with nonlinear component ∑k∈Kcdestk·qk·maxtkL−ydk,0 can be linearized by (32), (33), and (34) where gk is a newly introduced variable.
(32)∑k∈K∑i,j∈A∑s∈Sijcis·qk·wijsk+∑k∈Kcdestk·qk·gk,(33)gk≥tkL−ydk∀k∈K,(34)gk≥0∀k∈K.If the containers arrive at the destination earlier thantkL, tkL−ydk>0. Consequently, maxtkL−ydk,0=tkL−ydk. According to (32) and (33), gk≥tkL−ydk. Minimization of (32) will then lead to gk=tkL−ydk. In this case, ∑k∈Kcdestk·qk·maxtkL−ydk,0 is equivalent to (32) and (33).If the containers arrive at the destination later thantkL, tkL−ydk<0. Consequently, maxtkL−ydk,0=0. According to (32) and (33), gk≥0. Minimization of (32) will then lead to gk=0. In this case, ∑k∈Kcdestk·qk·maxtkL−ydk,0 is also equivalent to (32) and (33).Above all, the equivalency above is verified. Similarly, we have the following linear reformulation.(2) Linear Reformulation 2. Nonlinear component ∑k∈Kcpenk·maxydk−tkU,0 in (4) can be linearized by (35), (36), and (37) where σk is also a newly introduced variable.
(35)∑k∈Kcpenk·σk,(36)σk≥ydk−tkU∀k∈K,(37)σk≥0∀k∈K.(3) Linear Reformulation 3. Nonlinear constraint nik=zik/24∀k∈K,∀i∈N\dk (12) can be reformulated by linear constraints (38) and (39).
(38)nik≤zik24∀k∈K,∀i∈N\dk,(39)nik>zik24−1∀k∈K,∀i∈N\dk.We takezik=60 as example to verify the equivalency above. According to (12), nik=60/24=2.5=2. Meanwhile, according to (45) and (46), there exists 1.5<nik≤2.5. Because nik is a nonnegative integer, finally nik=2. Consequently, the equivalency above is tenable.(4) Linear Reformulation 4. By using the same linearization technique, nonlinear constraints zik+tijsk−yjk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij (13) and lir−yik·xhirk=0∀k∈K,∀h,i∈A,∀r∈Γhi (14) can be reformulated by linear constraints (40) and (41) and (42) and (43), respectively.
(40)zik+tijsk−yjk≥M·xijsk−1∀k∈K,∀i,j∈A,∀s∈Ωij,(41)zik+tijsk−yjk≤M·1−xijsk∀k∈K,∀i,j∈A,∀s∈Ωij,(42)lir−yik≥M·xhirk−1∀k∈K,∀h,i∈A,∀r∈Γhi,(43)lir−yik≤M·1−xhirk∀k∈K,∀h,i∈A,∀r∈Γhi.(5) Linear Reformulation 5. By using the same linearization form, nonlinear constraints max0,lis−yik−τ−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Γij (16) and max0,zik−yik−τ−wijsk·xijsk=0∀k∈K,∀i,j∈A,∀s∈Ωij (17) can be reformulated by linear constraints (44) and (45) and (46) and (47).
(44)wijsk≥M·xijsk−1+lis−yik−τis∀k∈K,∀i,j∈A,∀s∈Γij,(45)wijsk≤M·xijsk∀k∈K,∀i,j∈A,∀s∈Γij,(46)wijsk≥M·xijsk−1+zik−yik−τis∀k∈K,∀i,j∈A,∀s∈Ωij,(47)wijsk≤M·xijsk∀k∈K,∀i,j∈A,∀s∈Ωij.For detailed proofs of Linear Reformations 4 and 5, readers can refer to our previous study [2].


## 5.3. Improved Linear Reformulation

For the sake of readability, the improved linear reformulation derived from the mathematical model in Section4.2 is summarized as follows.Objective function:Minimize (1)+(2)+(32)+(35)+(5)Subject to:Constraints (6), (7), (9), (10), (15), (18–23), (25), (26), (28), (29), (31), (33), (34), (36–47)


## 6. Empirical Case Study Based on the Chinese Scenario

### 6.1. Case Description

In this section, we present a real-world case to demonstrate the feasibility of the proposed model and exact solution strategy in optimizing the practical problem. This empirical case is oriented on a Chinese scenario. In this scenario, several consignments of containers need to be transported from the western inland cities (e.g., Lanzhou and Hohhot) to the eastern sea ports (e.g., Qingdao) through the road-rail multimodal service network.In the multimodal service network, “road services” refer to the container trucks and “rail services” refer to the block container trains. The topological structure of the multimodal service network is shown in Figure5. In the multimodal service network, there are 10 periodically operated rail services (block container trains) and 15 road services (container trucks). The schedules of the rail services and the time-dependent travel times of the road services are separately presented in Appendix A (available here) and Appendix B (available here) in the supplementary material. Note that all the values in Appendix A (available here) and Appendix B (available here) are expressed as real numbers; for example, 10 : 30 is rewritten as 10.5 and 5 : 00 at the second day of the planning horizon is 29 (24+5). Other times in different days can be expressed in the same way.Figure 5
Multimodal service network for the empirical case.There are 10 transportation orders within a short-term planning horizon. The containers in these orders need to be transported from Lanzhou and Hohhot (inland cities in Northwest China) to Jiaozhou and Huangdaogang (sea ports in Qingdao, Shandong Province) and Lianyungang (a sea port in Jiangsu Province) through the multimodal service network shown in Figure5. The detailed information on the multiple transportation orders is given in Appendix C (available here) in the supplementary material. Note that in Appendix C (available here), if the destination inventory costs equal zero, it means that the receivers have their own inventory facilities and there is no need to rent them from other warehouse companies.


### 6.2. Parameter Setting

The unit transportation costs en route (unit: ¥/TEU) regarding the two transportation services are determined by (48).
(48)cijsk=crail1+crail2·dijss∈Γijcroad2·dijss∈Ωij∀k∈K,∀i,j∈A,∀s∈Sij.As can be seen in (48), the unit transportation costs for rail services en route include two parameters: crail1 which is the costs related to the container volume (unit: ¥/TEU) and crail2 which is the costs related to the turnover of the containers (the multiplication of the container volume and corresponding transportation distance, unit: ¥/TEU-km). For road services, such costs only relate to one parameter croad2 which is similar to crail2 (unit: ¥/TEU-km). The values of all these parameters above are regulated by the Ministry of Transport of China and China Railway Corporation [70].20 ft ISO containers are utilized to contain goods in the empirical case. According to 2015 China Railway Statistical Bulletin [71], the national railways consumed 15.69 million tons of the standard coal and accomplished 33,503.67 hundred million ton-kilometers of freight transportation. The unit CO2 emission rate of the standard coal is 0.69, and approximately 92% carbon is conversed into CO2. Assuming that the 20 ISO containers are all fully loaded (each container weights 24 tons), we can then calculate the CO2 emission factor of rail service as the following formula [72] where 44 is the molecular weight of one CO2 molecule and 12 is the atomic weight of one C atom.
(49)15.69∗2433503.67∗102∗0.69∗92%∗4412∗106=262gTEU−km.Similarly, by using the data from 2015 China Transportation Industry Development Statistical Bulletin [73], we can also gain the CO2 emission factor of road services which is 1064 g/TEU-km. The values of all the parameters in the empirical case are summarized as shown in Table 3. In addition, the values of the unit inventory costs of transportation services at the nodes (cis) and the free-of-charge inventory periods of services at the nodes (τis) are all presented in Figure 5. When the nodes are specifically the railway stations, all the corresponding cis=3.125 ¥/TEU and τis=48h. When the nodes are the origins of the transportation orders, cis varies based on the locations of the nodes and τis=0h if s∈Ωij.Table 3
Values of the parameters in the empirical case.
Parameter
Description
Value
Unit
crail1
Railway cost parameter regarding the volume of the containers
500
¥/TEU
crail2
Railway cost parameter regarding the turnover of the containers
2.025
¥/TEU-km
croad2
Road cost parameter regarding the turnover of the containers
6
¥/TEU-km
cs
Unit loading/unloading operation costs of the transportation services
25∗
¥/TEU
195∗∗
cco2
Unit CO2 emission costs
50∗∗∗
¥/ton
α
Confidence in the fuzzy chance constraint (8)
0.9
—
ems
CO2 emission factors for the transportation services
1064∗
g/TEU-km
262∗∗
∗Road service; ∗∗rail service. ∗∗∗Suggested by the research group at the Chinese Academy of Environmental Planning [74].


### 6.3. Simulation Environment

Described by the linearized programming model shown in Section5.3, the empirical case can be solved by any exact solution algorithms, for example, branch-and-bound algorithm and simplex algorithm. In this study, we adopt the branch-and-bound algorithm to solve the problem. The algorithm is conducted by the mathematical programming software LINGO 12 designed by LINDO Systems Inc., Chicago, IL, USA [75], on a ThinkPad Laptop with Intel Core i5-5200U 2.20 GHz CPU 8 GB RAM. The scale of the empirical case is indicated by Table 4.Table 4
Scale of the empirical case problem.
Total variables
Integer variables
Constraints
5829
1980
8417


### 6.4. Optimization Result and Multimodal Route Illustration

The optimization result of the empirical case and the performance of the exact solution strategy are shown in Table5.Table 5
Optimization result and strategy performance.
Solution algorithm
Solver state
Best solution
Computational time
Standard branch-and-bound algorithm
Global optimum
¥ 1,575,559
1 min 06 sec∗
∗Average value of 10 times simulation.The structure of the minimal generalized costs for the multiple transportation orders is illustrated by Figure6 where C1 to C5 are successively the transportation costs en route, loading/unloading operation costs, inventory costs, penalty costs, and CO2 emission costs, which also correspond to (1) to (5) in the objective function.Figure 6
Structure of the minimal generalized costs.The best routes for the containers in the multiple transportation orders are given in Table6. According to the planned routes, the accomplishments of 4 transportation orders satisfy their due date time windows exactly.Table 6
Best routes in the empirical case.
Number
Best multimodal routes
Arrival time
1
Lanzhou {28} — (14) — Xinzhu — (4) — Jiaozhou
130.4
2
Lanzhou {37} — (14) — Xinzhu — (4) — Jiaozhou
130.4∗∗
3
Lanzhou {111} — (14) — Xinzhu — (5) — Huangdaogang
214.6
4
Lanzhou {15} — (17) — Putian — (10) — Lianyungang
79∗∗
5
Lanzhou {90} — (17) — Putian — (10)— Lianyungang
175∗∗
6
Hohhot — (7) — Xingang {96} — (23) — Jiaozhou
106.9∗∗
7
Hohhot — (7) — Xingang {101.5} — (23) — Jiaozhou
106.9
8
Hohhot — (7) — Xingang {148} — (23) — Jiaozhou
72∗
9
Hohhot — (7) — Xingang {192} — (23) — Jiaozhou
72∗
10
Hohhot {72} — (19) — Baotou — (6) — Xingang {171.5} — (25) — Lianyungang
186
∗Early delivery; ∗∗delayed delivery; numbers in {}: planned departure time.


### 6.5. Case Discussion

#### 6.5.1. Sensitivity of the Multimodal Routing with respect to the CO2 Emissions

First of all, we analyze the effect of the CO2 emissions on the multimodal routing. Let the unit emission costs vary from 50 to 100 ¥/ton, the sensitivity of the multimodal routing result with respect to this variation can be seen in Figure 7.Figure 7
Sensitivity of the multimodal routing with respect to the CO2 emissions.Figure7 shows that the CO2 emission costs increase linearly due to the fact that the emission volume is a constant (87.1 tons), while the rest of the generalized costs stay constant with the unit emission costs increasing. Even though we continue the simulation until the unit emission costs reach up to 1000 ¥/ton, which is obviously infeasible in practice (according to the formulation of the Chinese Academy of Environmental Planning [73], the feasible unit emission costs should reach up to 100 ¥/ton at 2030), the sensitivity still keeps the same tendency as shown in Figure 6. It means that the best routes in the empirical case do not modify when the emission costs vary and are hence insensitive to this variation. Therefore, charging for CO2 emissions within a reasonable range is not particularly helpful to promote the green multimodal transportation, at least in this empirical case.


#### 6.5.2. Bi-Objective Optimization regarding CO2 Emissions

Furthermore, we modify the proposed model into a bi-objective optimization model with the objectives of minimizing the generalized costs (Z1=Eq.1+Eq.2+Eq.3+Eq.4) and of minimizing the CO2 emissions (Z2=Eq.51).
(50)minimize∑k∈K∑i,j∈A∑s∈Sijems·qk·dijs·xijsk.Using the widely utilized weighted sum method [43] defined in (51) and (52), we can solve the bi-objective optimization problem and generate its Pareto solutions shown in Figure 7.
(51)minimizeω·Z1+1−ω·Z2,(52)0<ω≤1.Especially, in Figure8, compared with the last Pareto solution, the second last one can lower the CO2 emissions from 87.1 to 83.7 by 3.90% and cause only a slight increase of the generalized costs from 1571.2 to 1586.0 by 0.94%. By comparing Figure 8 with Figure 7, we can conclude that a bi-objective optimization approach is much more effective than charging for CO2 emissions in helping decision makers to control and optimize the CO2 emissions in the multimodal service network.Figure 8
Pareto solutions to the bi-objective optimization problem.


#### 6.5.3. Sensitivity of the Multimodal Routing with respect to the Confidence

Confidenceα in the fuzzy chance constraint (8) is set by the decision makers and reflects their preference for the reliability of the multimodal routing such that all the loaded volume of a rail service does not exceed its carrying capacity. Its value has a remarkable effect on the optimization result of the multimodal routing, which can be investigated by using sensitivity analysis. Let α vary from 0.1 to 1.0 with a step size of 0.1, we can obtain the sensitivity of the multimodal routing result (indexed by the minimal generalized costs) with respect to such variation, which can be seen in Figure 9.Figure 9
Sensitivity of the multimodal routing with respect to the confidence.We can also draw some helpful insights from Figure9 summarized as follows.
(1)
Overall, larger values of confidenceα will lead to larger generalized costs for the best routes.
(2)
The variation of the generalized costs with respect to the confidence is stepwise. In this case study, the multimodal routing is sensitive to the confidence when its value exceeds 0.3.
(3)
The economy and reliability of the multimodal routing cannot reach their respective optimum simultaneously. The economy will be sacrificed if the decision makers decide to improve the reliability of the planned best routes, and vice versa.


#### 6.5.4. Fuzzy Simulation

During the decision-making process regarding the multimodal routing, it is challenging for decision makers to determine an objective confidence value. In order to help them select reasonable confidence values, fuzzy simulation is adopted in this study.The fuzzy simulation simulates the actual deterministic transportation scenario by randomly generating the carrying capacities of rail services according to the membership degree of the triangular fuzzy numbers. The simulation process is shown in Figure10 [8].Figure 10
Fuzzy simulation process.After the simulation, we can first gain deterministic carrying capacitiesvijs for i,j∈A, s∈Γij. Eq. (8) can be consequently transformed into a deterministic carrying capacity constraint as (53).
(53)∑k∈Kqk·xijsk≤vijs∀i,j∈A,∀s∈Γij.Then, we can check if the planned routes given by the fuzzy programming model satisfy (15) or not. If the planned routes satisfy the constraint, we define them as a successful plan; otherwise, a failed plan which means that the planned routes are infeasible due to the violation of the capacity constraint when adopted to move containers in the practical transportation. If the transportation starts and a plan is found to be infeasible due to the insufficient actual capacities of some rail services, an alternative plan should be generated according to the following principle in order to reduce the additional costs and CO2 emissions caused by the alternative: the rail service with insufficient capacity should be preferentially used by the transportation orders with larger container volumes initially assigned to it as many as possible on the condition that (53) must be satisfied, that is, the transportation orders with larger container volumes should be accomplished by the initial plan as far as possible, and the rest transportation orders usually with smaller container volumes assigned to the rail service should be accomplished by the origin-to-destination road services. The fuzzy simulation should be conducted several times. In this study, we ran the simulation 50 times, and the fuzzy simulation result is shown in Figure 11.Figure 11
Fuzzy simulation result regarding confidenceα.As we can see from Figure11, the successful times of the planned routes increase remarkably when the confidence value changes from 0.7 to 0.8. Although the generalized costs increase approximately by 1.29%, the success ratio of the planned routes considerably increases from 26% to 100% by more than 3.8 times when the confidence value is set to 0.8 instead of 0.7.Note that whenα=0.5, according to (25) and (26), the fuzzy chance constraint can be rewritten as ∑k∈Kqk·xijsk≤vijsM which is the deterministic capacity constraint widely employed by the current literature (e.g., Sun and Lang [2] and Sun et al. [19]). So we know that the successful ratio of the multimodal routing with deterministic capacity consideration is only 12% in the 50 simulations for the specific case. As we can see from Figure 11, the successful ratio of the multimodal routing can be significantly improved by more than 8.3 times by considering the capacity uncertainty and setting the confidence to 0.8, 0.9, or 1.0. Therefore, it is of great value to formulate the capacity uncertainty when planning the multimodal routes, so that higher reliability of the multimodal routing decision-making can be guaranteed.We can further achieve the optimization results of the 50 deterministic problems formulated by the linearized model with (25) and (26) replaced by (53). The optimization results are presented in Appendix D (available here) in the supplementary material. Such results can be considered as the actual best routes. By comparing the actual best routes with the planned ones, we can find the confidence value that best matches the actual transportation.In the 50 fuzzy simulations, the best actual routes whose generalized costs are equal to ¥ 1,575,559 emerge 38 times that account for 76% of the entire simulation times (see in Figure12). That is to say, such routes are the most likely best ones in the practical transportation. The planned best routes given by the fuzzy programming model with confidence values of 0.8, 0.9, and 1.0 match the most likely best actual routes better than others. Above all, the best values of the confidence in this empirical case are recommended to be 0.8, 0.9, and 1.0.Figure 12
Ratios of the optimization results in the 50 fuzzy simulations.


### 6.6. Managerial Implications

It should be noted that the optimization results and conclusions in this section are drawn based on a specific case whose setting is shown in Sections6.1 and 6.2. The results and conclusions are not universal and are sensitive to the case itself and might change if the case changes. Therefore, the above analysis and discussion in this section only provide a procedure and a demonstration on how to adopt the proposed methods to deal with the practical problems. Although the cases might be different, we can always use such procedure to address them to draw corresponding optimization results and conclusions.Despite this fact, the results show important implications for the practical implementation of the model:(1)
Although the multimodal service network is complex, the linearized optimization model can deliver the optimal solution in a very short time so that new plans for incoming orders can be found very fast.
(2)
The objective of the optimization is to minimize the total costs consisting of different cost categories (i.e., transportation, transshipment, inventory, penalty, and emission costs). However, there is a clear dominance of transportation costs over all other cost categories which makes it difficult to influence the optimal solution by changing other cost factors than transportation costs.
(3)
Based on the previous implication, it is clear that the charging method for CO2 emissions does not improve the environmental sustainability of the multimodal routes if realistic values are used for emission costs. Therefore, rather a bi-objective approach should be used in order to be able to compare the cost-optimal and emission-optimal solution. In this case, often a solution can be found where one objective is significantly improved whereas there is only a small increase for the other one.
(4)
The minimal values of total costs can be usually achieved under the deterministic setting since no buffer times or capacities need to be included. However, uncertainty plays an important role in practice and therefore should be considered also in planning. Although the total costs are increasing with the increasing confidence level, this also reduces the risk of disruptions during transportation which can be seen in a lower ratio of failed plans. As a result, more transportations can be performed as originally planned and the need for re-planning for failed plans is minimized. Since such re-planning usually leads to additional costs, reducing the share of failed plans is of high interest.


## 6.1. Case Description

In this section, we present a real-world case to demonstrate the feasibility of the proposed model and exact solution strategy in optimizing the practical problem. This empirical case is oriented on a Chinese scenario. In this scenario, several consignments of containers need to be transported from the western inland cities (e.g., Lanzhou and Hohhot) to the eastern sea ports (e.g., Qingdao) through the road-rail multimodal service network.In the multimodal service network, “road services” refer to the container trucks and “rail services” refer to the block container trains. The topological structure of the multimodal service network is shown in Figure5. In the multimodal service network, there are 10 periodically operated rail services (block container trains) and 15 road services (container trucks). The schedules of the rail services and the time-dependent travel times of the road services are separately presented in Appendix A (available here) and Appendix B (available here) in the supplementary material. Note that all the values in Appendix A (available here) and Appendix B (available here) are expressed as real numbers; for example, 10 : 30 is rewritten as 10.5 and 5 : 00 at the second day of the planning horizon is 29 (24+5). Other times in different days can be expressed in the same way.Figure 5
Multimodal service network for the empirical case.There are 10 transportation orders within a short-term planning horizon. The containers in these orders need to be transported from Lanzhou and Hohhot (inland cities in Northwest China) to Jiaozhou and Huangdaogang (sea ports in Qingdao, Shandong Province) and Lianyungang (a sea port in Jiangsu Province) through the multimodal service network shown in Figure5. The detailed information on the multiple transportation orders is given in Appendix C (available here) in the supplementary material. Note that in Appendix C (available here), if the destination inventory costs equal zero, it means that the receivers have their own inventory facilities and there is no need to rent them from other warehouse companies.


## 6.2. Parameter Setting

The unit transportation costs en route (unit: ¥/TEU) regarding the two transportation services are determined by (48).
(48)cijsk=crail1+crail2·dijss∈Γijcroad2·dijss∈Ωij∀k∈K,∀i,j∈A,∀s∈Sij.As can be seen in (48), the unit transportation costs for rail services en route include two parameters: crail1 which is the costs related to the container volume (unit: ¥/TEU) and crail2 which is the costs related to the turnover of the containers (the multiplication of the container volume and corresponding transportation distance, unit: ¥/TEU-km). For road services, such costs only relate to one parameter croad2 which is similar to crail2 (unit: ¥/TEU-km). The values of all these parameters above are regulated by the Ministry of Transport of China and China Railway Corporation [70].20 ft ISO containers are utilized to contain goods in the empirical case. According to 2015 China Railway Statistical Bulletin [71], the national railways consumed 15.69 million tons of the standard coal and accomplished 33,503.67 hundred million ton-kilometers of freight transportation. The unit CO2 emission rate of the standard coal is 0.69, and approximately 92% carbon is conversed into CO2. Assuming that the 20 ISO containers are all fully loaded (each container weights 24 tons), we can then calculate the CO2 emission factor of rail service as the following formula [72] where 44 is the molecular weight of one CO2 molecule and 12 is the atomic weight of one C atom.
(49)15.69∗2433503.67∗102∗0.69∗92%∗4412∗106=262gTEU−km.Similarly, by using the data from 2015 China Transportation Industry Development Statistical Bulletin [73], we can also gain the CO2 emission factor of road services which is 1064 g/TEU-km. The values of all the parameters in the empirical case are summarized as shown in Table 3. In addition, the values of the unit inventory costs of transportation services at the nodes (cis) and the free-of-charge inventory periods of services at the nodes (τis) are all presented in Figure 5. When the nodes are specifically the railway stations, all the corresponding cis=3.125 ¥/TEU and τis=48h. When the nodes are the origins of the transportation orders, cis varies based on the locations of the nodes and τis=0h if s∈Ωij.Table 3
Values of the parameters in the empirical case.
Parameter
Description
Value
Unit
crail1
Railway cost parameter regarding the volume of the containers
500
¥/TEU
crail2
Railway cost parameter regarding the turnover of the containers
2.025
¥/TEU-km
croad2
Road cost parameter regarding the turnover of the containers
6
¥/TEU-km
cs
Unit loading/unloading operation costs of the transportation services
25∗
¥/TEU
195∗∗
cco2
Unit CO2 emission costs
50∗∗∗
¥/ton
α
Confidence in the fuzzy chance constraint (8)
0.9
—
ems
CO2 emission factors for the transportation services
1064∗
g/TEU-km
262∗∗
∗Road service; ∗∗rail service. ∗∗∗Suggested by the research group at the Chinese Academy of Environmental Planning [74].


## 6.3. Simulation Environment

Described by the linearized programming model shown in Section5.3, the empirical case can be solved by any exact solution algorithms, for example, branch-and-bound algorithm and simplex algorithm. In this study, we adopt the branch-and-bound algorithm to solve the problem. The algorithm is conducted by the mathematical programming software LINGO 12 designed by LINDO Systems Inc., Chicago, IL, USA [75], on a ThinkPad Laptop with Intel Core i5-5200U 2.20 GHz CPU 8 GB RAM. The scale of the empirical case is indicated by Table 4.Table 4
Scale of the empirical case problem.
Total variables
Integer variables
Constraints
5829
1980
8417


## 6.4. Optimization Result and Multimodal Route Illustration

The optimization result of the empirical case and the performance of the exact solution strategy are shown in Table5.Table 5
Optimization result and strategy performance.
Solution algorithm
Solver state
Best solution
Computational time
Standard branch-and-bound algorithm
Global optimum
¥ 1,575,559
1 min 06 sec∗
∗Average value of 10 times simulation.The structure of the minimal generalized costs for the multiple transportation orders is illustrated by Figure6 where C1 to C5 are successively the transportation costs en route, loading/unloading operation costs, inventory costs, penalty costs, and CO2 emission costs, which also correspond to (1) to (5) in the objective function.Figure 6
Structure of the minimal generalized costs.The best routes for the containers in the multiple transportation orders are given in Table6. According to the planned routes, the accomplishments of 4 transportation orders satisfy their due date time windows exactly.Table 6
Best routes in the empirical case.
Number
Best multimodal routes
Arrival time
1
Lanzhou {28} — (14) — Xinzhu — (4) — Jiaozhou
130.4
2
Lanzhou {37} — (14) — Xinzhu — (4) — Jiaozhou
130.4∗∗
3
Lanzhou {111} — (14) — Xinzhu — (5) — Huangdaogang
214.6
4
Lanzhou {15} — (17) — Putian — (10) — Lianyungang
79∗∗
5
Lanzhou {90} — (17) — Putian — (10)— Lianyungang
175∗∗
6
Hohhot — (7) — Xingang {96} — (23) — Jiaozhou
106.9∗∗
7
Hohhot — (7) — Xingang {101.5} — (23) — Jiaozhou
106.9
8
Hohhot — (7) — Xingang {148} — (23) — Jiaozhou
72∗
9
Hohhot — (7) — Xingang {192} — (23) — Jiaozhou
72∗
10
Hohhot {72} — (19) — Baotou — (6) — Xingang {171.5} — (25) — Lianyungang
186
∗Early delivery; ∗∗delayed delivery; numbers in {}: planned departure time.


## 6.5. Case Discussion

### 6.5.1. Sensitivity of the Multimodal Routing with respect to the CO2 Emissions

First of all, we analyze the effect of the CO2 emissions on the multimodal routing. Let the unit emission costs vary from 50 to 100 ¥/ton, the sensitivity of the multimodal routing result with respect to this variation can be seen in Figure 7.Figure 7
Sensitivity of the multimodal routing with respect to the CO2 emissions.Figure7 shows that the CO2 emission costs increase linearly due to the fact that the emission volume is a constant (87.1 tons), while the rest of the generalized costs stay constant with the unit emission costs increasing. Even though we continue the simulation until the unit emission costs reach up to 1000 ¥/ton, which is obviously infeasible in practice (according to the formulation of the Chinese Academy of Environmental Planning [73], the feasible unit emission costs should reach up to 100 ¥/ton at 2030), the sensitivity still keeps the same tendency as shown in Figure 6. It means that the best routes in the empirical case do not modify when the emission costs vary and are hence insensitive to this variation. Therefore, charging for CO2 emissions within a reasonable range is not particularly helpful to promote the green multimodal transportation, at least in this empirical case.


### 6.5.2. Bi-Objective Optimization regarding CO2 Emissions

Furthermore, we modify the proposed model into a bi-objective optimization model with the objectives of minimizing the generalized costs (Z1=Eq.1+Eq.2+Eq.3+Eq.4) and of minimizing the CO2 emissions (Z2=Eq.51).
(50)minimize∑k∈K∑i,j∈A∑s∈Sijems·qk·dijs·xijsk.Using the widely utilized weighted sum method [43] defined in (51) and (52), we can solve the bi-objective optimization problem and generate its Pareto solutions shown in Figure 7.
(51)minimizeω·Z1+1−ω·Z2,(52)0<ω≤1.Especially, in Figure8, compared with the last Pareto solution, the second last one can lower the CO2 emissions from 87.1 to 83.7 by 3.90% and cause only a slight increase of the generalized costs from 1571.2 to 1586.0 by 0.94%. By comparing Figure 8 with Figure 7, we can conclude that a bi-objective optimization approach is much more effective than charging for CO2 emissions in helping decision makers to control and optimize the CO2 emissions in the multimodal service network.Figure 8
Pareto solutions to the bi-objective optimization problem.


### 6.5.3. Sensitivity of the Multimodal Routing with respect to the Confidence

Confidenceα in the fuzzy chance constraint (8) is set by the decision makers and reflects their preference for the reliability of the multimodal routing such that all the loaded volume of a rail service does not exceed its carrying capacity. Its value has a remarkable effect on the optimization result of the multimodal routing, which can be investigated by using sensitivity analysis. Let α vary from 0.1 to 1.0 with a step size of 0.1, we can obtain the sensitivity of the multimodal routing result (indexed by the minimal generalized costs) with respect to such variation, which can be seen in Figure 9.Figure 9
Sensitivity of the multimodal routing with respect to the confidence.We can also draw some helpful insights from Figure9 summarized as follows.
(1)
Overall, larger values of confidenceα will lead to larger generalized costs for the best routes.
(2)
The variation of the generalized costs with respect to the confidence is stepwise. In this case study, the multimodal routing is sensitive to the confidence when its value exceeds 0.3.
(3)
The economy and reliability of the multimodal routing cannot reach their respective optimum simultaneously. The economy will be sacrificed if the decision makers decide to improve the reliability of the planned best routes, and vice versa.


### 6.5.4. Fuzzy Simulation

During the decision-making process regarding the multimodal routing, it is challenging for decision makers to determine an objective confidence value. In order to help them select reasonable confidence values, fuzzy simulation is adopted in this study.The fuzzy simulation simulates the actual deterministic transportation scenario by randomly generating the carrying capacities of rail services according to the membership degree of the triangular fuzzy numbers. The simulation process is shown in Figure10 [8].Figure 10
Fuzzy simulation process.After the simulation, we can first gain deterministic carrying capacitiesvijs for i,j∈A, s∈Γij. Eq. (8) can be consequently transformed into a deterministic carrying capacity constraint as (53).
(53)∑k∈Kqk·xijsk≤vijs∀i,j∈A,∀s∈Γij.Then, we can check if the planned routes given by the fuzzy programming model satisfy (15) or not. If the planned routes satisfy the constraint, we define them as a successful plan; otherwise, a failed plan which means that the planned routes are infeasible due to the violation of the capacity constraint when adopted to move containers in the practical transportation. If the transportation starts and a plan is found to be infeasible due to the insufficient actual capacities of some rail services, an alternative plan should be generated according to the following principle in order to reduce the additional costs and CO2 emissions caused by the alternative: the rail service with insufficient capacity should be preferentially used by the transportation orders with larger container volumes initially assigned to it as many as possible on the condition that (53) must be satisfied, that is, the transportation orders with larger container volumes should be accomplished by the initial plan as far as possible, and the rest transportation orders usually with smaller container volumes assigned to the rail service should be accomplished by the origin-to-destination road services. The fuzzy simulation should be conducted several times. In this study, we ran the simulation 50 times, and the fuzzy simulation result is shown in Figure 11.Figure 11
Fuzzy simulation result regarding confidenceα.As we can see from Figure11, the successful times of the planned routes increase remarkably when the confidence value changes from 0.7 to 0.8. Although the generalized costs increase approximately by 1.29%, the success ratio of the planned routes considerably increases from 26% to 100% by more than 3.8 times when the confidence value is set to 0.8 instead of 0.7.Note that whenα=0.5, according to (25) and (26), the fuzzy chance constraint can be rewritten as ∑k∈Kqk·xijsk≤vijsM which is the deterministic capacity constraint widely employed by the current literature (e.g., Sun and Lang [2] and Sun et al. [19]). So we know that the successful ratio of the multimodal routing with deterministic capacity consideration is only 12% in the 50 simulations for the specific case. As we can see from Figure 11, the successful ratio of the multimodal routing can be significantly improved by more than 8.3 times by considering the capacity uncertainty and setting the confidence to 0.8, 0.9, or 1.0. Therefore, it is of great value to formulate the capacity uncertainty when planning the multimodal routes, so that higher reliability of the multimodal routing decision-making can be guaranteed.We can further achieve the optimization results of the 50 deterministic problems formulated by the linearized model with (25) and (26) replaced by (53). The optimization results are presented in Appendix D (available here) in the supplementary material. Such results can be considered as the actual best routes. By comparing the actual best routes with the planned ones, we can find the confidence value that best matches the actual transportation.In the 50 fuzzy simulations, the best actual routes whose generalized costs are equal to ¥ 1,575,559 emerge 38 times that account for 76% of the entire simulation times (see in Figure12). That is to say, such routes are the most likely best ones in the practical transportation. The planned best routes given by the fuzzy programming model with confidence values of 0.8, 0.9, and 1.0 match the most likely best actual routes better than others. Above all, the best values of the confidence in this empirical case are recommended to be 0.8, 0.9, and 1.0.Figure 12
Ratios of the optimization results in the 50 fuzzy simulations.


## 6.5.1. Sensitivity of the Multimodal Routing with respect to the CO2 Emissions

First of all, we analyze the effect of the CO2 emissions on the multimodal routing. Let the unit emission costs vary from 50 to 100 ¥/ton, the sensitivity of the multimodal routing result with respect to this variation can be seen in Figure 7.Figure 7
Sensitivity of the multimodal routing with respect to the CO2 emissions.Figure7 shows that the CO2 emission costs increase linearly due to the fact that the emission volume is a constant (87.1 tons), while the rest of the generalized costs stay constant with the unit emission costs increasing. Even though we continue the simulation until the unit emission costs reach up to 1000 ¥/ton, which is obviously infeasible in practice (according to the formulation of the Chinese Academy of Environmental Planning [73], the feasible unit emission costs should reach up to 100 ¥/ton at 2030), the sensitivity still keeps the same tendency as shown in Figure 6. It means that the best routes in the empirical case do not modify when the emission costs vary and are hence insensitive to this variation. Therefore, charging for CO2 emissions within a reasonable range is not particularly helpful to promote the green multimodal transportation, at least in this empirical case.


## 6.5.2. Bi-Objective Optimization regarding CO2 Emissions

Furthermore, we modify the proposed model into a bi-objective optimization model with the objectives of minimizing the generalized costs (Z1=Eq.1+Eq.2+Eq.3+Eq.4) and of minimizing the CO2 emissions (Z2=Eq.51).
(50)minimize∑k∈K∑i,j∈A∑s∈Sijems·qk·dijs·xijsk.Using the widely utilized weighted sum method [43] defined in (51) and (52), we can solve the bi-objective optimization problem and generate its Pareto solutions shown in Figure 7.
(51)minimizeω·Z1+1−ω·Z2,(52)0<ω≤1.Especially, in Figure8, compared with the last Pareto solution, the second last one can lower the CO2 emissions from 87.1 to 83.7 by 3.90% and cause only a slight increase of the generalized costs from 1571.2 to 1586.0 by 0.94%. By comparing Figure 8 with Figure 7, we can conclude that a bi-objective optimization approach is much more effective than charging for CO2 emissions in helping decision makers to control and optimize the CO2 emissions in the multimodal service network.Figure 8
Pareto solutions to the bi-objective optimization problem.


## 6.5.3. Sensitivity of the Multimodal Routing with respect to the Confidence

Confidenceα in the fuzzy chance constraint (8) is set by the decision makers and reflects their preference for the reliability of the multimodal routing such that all the loaded volume of a rail service does not exceed its carrying capacity. Its value has a remarkable effect on the optimization result of the multimodal routing, which can be investigated by using sensitivity analysis. Let α vary from 0.1 to 1.0 with a step size of 0.1, we can obtain the sensitivity of the multimodal routing result (indexed by the minimal generalized costs) with respect to such variation, which can be seen in Figure 9.Figure 9
Sensitivity of the multimodal routing with respect to the confidence.We can also draw some helpful insights from Figure9 summarized as follows.
(1)
Overall, larger values of confidenceα will lead to larger generalized costs for the best routes.
(2)
The variation of the generalized costs with respect to the confidence is stepwise. In this case study, the multimodal routing is sensitive to the confidence when its value exceeds 0.3.
(3)
The economy and reliability of the multimodal routing cannot reach their respective optimum simultaneously. The economy will be sacrificed if the decision makers decide to improve the reliability of the planned best routes, and vice versa.


## 6.5.4. Fuzzy Simulation

During the decision-making process regarding the multimodal routing, it is challenging for decision makers to determine an objective confidence value. In order to help them select reasonable confidence values, fuzzy simulation is adopted in this study.The fuzzy simulation simulates the actual deterministic transportation scenario by randomly generating the carrying capacities of rail services according to the membership degree of the triangular fuzzy numbers. The simulation process is shown in Figure10 [8].Figure 10
Fuzzy simulation process.After the simulation, we can first gain deterministic carrying capacitiesvijs for i,j∈A, s∈Γij. Eq. (8) can be consequently transformed into a deterministic carrying capacity constraint as (53).
(53)∑k∈Kqk·xijsk≤vijs∀i,j∈A,∀s∈Γij.Then, we can check if the planned routes given by the fuzzy programming model satisfy (15) or not. If the planned routes satisfy the constraint, we define them as a successful plan; otherwise, a failed plan which means that the planned routes are infeasible due to the violation of the capacity constraint when adopted to move containers in the practical transportation. If the transportation starts and a plan is found to be infeasible due to the insufficient actual capacities of some rail services, an alternative plan should be generated according to the following principle in order to reduce the additional costs and CO2 emissions caused by the alternative: the rail service with insufficient capacity should be preferentially used by the transportation orders with larger container volumes initially assigned to it as many as possible on the condition that (53) must be satisfied, that is, the transportation orders with larger container volumes should be accomplished by the initial plan as far as possible, and the rest transportation orders usually with smaller container volumes assigned to the rail service should be accomplished by the origin-to-destination road services. The fuzzy simulation should be conducted several times. In this study, we ran the simulation 50 times, and the fuzzy simulation result is shown in Figure 11.Figure 11
Fuzzy simulation result regarding confidenceα.As we can see from Figure11, the successful times of the planned routes increase remarkably when the confidence value changes from 0.7 to 0.8. Although the generalized costs increase approximately by 1.29%, the success ratio of the planned routes considerably increases from 26% to 100% by more than 3.8 times when the confidence value is set to 0.8 instead of 0.7.Note that whenα=0.5, according to (25) and (26), the fuzzy chance constraint can be rewritten as ∑k∈Kqk·xijsk≤vijsM which is the deterministic capacity constraint widely employed by the current literature (e.g., Sun and Lang [2] and Sun et al. [19]). So we know that the successful ratio of the multimodal routing with deterministic capacity consideration is only 12% in the 50 simulations for the specific case. As we can see from Figure 11, the successful ratio of the multimodal routing can be significantly improved by more than 8.3 times by considering the capacity uncertainty and setting the confidence to 0.8, 0.9, or 1.0. Therefore, it is of great value to formulate the capacity uncertainty when planning the multimodal routes, so that higher reliability of the multimodal routing decision-making can be guaranteed.We can further achieve the optimization results of the 50 deterministic problems formulated by the linearized model with (25) and (26) replaced by (53). The optimization results are presented in Appendix D (available here) in the supplementary material. Such results can be considered as the actual best routes. By comparing the actual best routes with the planned ones, we can find the confidence value that best matches the actual transportation.In the 50 fuzzy simulations, the best actual routes whose generalized costs are equal to ¥ 1,575,559 emerge 38 times that account for 76% of the entire simulation times (see in Figure12). That is to say, such routes are the most likely best ones in the practical transportation. The planned best routes given by the fuzzy programming model with confidence values of 0.8, 0.9, and 1.0 match the most likely best actual routes better than others. Above all, the best values of the confidence in this empirical case are recommended to be 0.8, 0.9, and 1.0.Figure 12
Ratios of the optimization results in the 50 fuzzy simulations.


## 6.6. Managerial Implications

It should be noted that the optimization results and conclusions in this section are drawn based on a specific case whose setting is shown in Sections6.1 and 6.2. The results and conclusions are not universal and are sensitive to the case itself and might change if the case changes. Therefore, the above analysis and discussion in this section only provide a procedure and a demonstration on how to adopt the proposed methods to deal with the practical problems. Although the cases might be different, we can always use such procedure to address them to draw corresponding optimization results and conclusions.Despite this fact, the results show important implications for the practical implementation of the model:(1)
Although the multimodal service network is complex, the linearized optimization model can deliver the optimal solution in a very short time so that new plans for incoming orders can be found very fast.
(2)
The objective of the optimization is to minimize the total costs consisting of different cost categories (i.e., transportation, transshipment, inventory, penalty, and emission costs). However, there is a clear dominance of transportation costs over all other cost categories which makes it difficult to influence the optimal solution by changing other cost factors than transportation costs.
(3)
Based on the previous implication, it is clear that the charging method for CO2 emissions does not improve the environmental sustainability of the multimodal routes if realistic values are used for emission costs. Therefore, rather a bi-objective approach should be used in order to be able to compare the cost-optimal and emission-optimal solution. In this case, often a solution can be found where one objective is significantly improved whereas there is only a small increase for the other one.
(4)
The minimal values of total costs can be usually achieved under the deterministic setting since no buffer times or capacities need to be included. However, uncertainty plays an important role in practice and therefore should be considered also in planning. Although the total costs are increasing with the increasing confidence level, this also reduces the risk of disruptions during transportation which can be seen in a lower ratio of failed plans. As a result, more transportations can be performed as originally planned and the need for re-planning for failed plans is minimized. Since such re-planning usually leads to additional costs, reducing the share of failed plans is of high interest.


## 7. Conclusions

Nowadays, the multimodal transportation has been widely promoted all over the world, and its routing optimization has attached great importance by both researchers and practitioners. Meanwhile, the public concern on environmental protection raises, and environmentally friendly development gets tremendous highlights in almost every industry. Therefore, it is valuable to combine multimodal routing with the environmental views. In addition to that, different uncertainty factors can influence the performance of multimodal transportation and therefore should be considered already in the planning phase. However, all these factors increase the complexity of multimodal routing modelling and therefore make the transportation planning very challenging in practice.In order to respond to these challenges, we developed a nonlinear mixed-integer optimization model covering multiple objectives (i.e., economic and environmental) and multiple uncertainty factors (i.e., traffic congestion and rail service capacity). Moreover, in order to facilitate the solution of this problem, a linearization approach was proposed. The application of the proposed model to a real-world case study showed that it can deliver optimal results in a relatively short time which is very helpful for the practical transportation planning. Besides that, it was shown that the inclusion of environmental criteria is meaningful only when the bi-objective method is used, since the weight of the emission costs is very low in comparison to transportation costs if the charging method is used. Regarding the reliability of transportation, including the uncertainty factors into planning and increasing confidentiality levels significantly reduce the risk of disruptions during the execution phase but also lead to increased total costs. However, the reduced risk of disruption also reduces the need for re-planning and therefore helps to avoid additional costs connected with this process.Compared with the models proposed by other studies listed in the literature review section, our model can comprehensively formulate the issues of road traffic congestion and rail service capacity uncertainty that actually exist in the routing decision-making process and thus have better feasibility in dealing with the practical problem. Another improvement is that by modifying the objective function, our model can easily use two different approaches, including charging for CO2 emissions and bi-objective optimization, to separately optimize the CO2 emissions and identify which one is more effective in planning green routes. Although the model is nonlinear, it can be linearized, which enables the problem to be effectively solved by any exact solution algorithm on any mathematical programming software.This study mainly focuses on the modelling and the utilization of a classical exact solution algorithm based on a linearized model and does not involve any intelligent algorithms. The exact solution strategy can efficiently optimize the presented case, but its performance is doubtful when the scale of the case is increasingly expanded. As a result, for further improvement to this study, we will test the feasibility of the exact solution strategy in dealing with a larger-scale empirical case. If the test result is not satisfactory, we will try to design some advanced intelligent solution algorithms (e.g., heuristic algorithms) with higher solution efficiency, so that the problem can be effectively optimized and routing suggestions can be provided to the decision makers in time when the problem scale gets larger.



---

*Source: 8645793-2018-06-24.xml*