Fuzzy Programming Approaches for Modeling a Customer-Centred Freight Routing Problem in the Road-Rail Intermodal Hub-and-Spoke Network with Fuzzy Soft Time Windows and Multiple Sources of Time Uncertainty

Abstract

In this study, we systematically investigate a road-rail intermodal routing problem the optimization of which is oriented on the customer demands on transportation economy, timeliness and reliability. The road-rail intermodal transportation system is modelled as a hub-and-spoke network that contains time-flexible container truck services and scheduled container train services. The transportation timeliness is optimized by using fuzzy soft time windows associated with the service level of the transportation. Reliability is enhanced by considering multiple sources of time uncertainty, including road travel time and loading/unloading time. Such uncertainty is modelled by using fuzzy set theory. Triangular fuzzy numbers are adopted to represent the uncertain time. Under the above consideration, we first establish a fuzzy mixed integer nonlinear programming model with a weighted objective that includes minimizing the costs and maximizing the service level for accomplishing transportation orders. Then we use the fuzzy expected value model and fuzzy chance-constrained programming separately to realize the defuzzification of the fuzzy objective and use fuzzy chance-constrained programming to deal with the fuzzy constraint. After defuzzification and linearization, an equivalent mixed integer linear programming (MILP) model is generated to enable the problem to be solved by mathematical programming software. Finally, a numerical case modified from our previous study is presented to demonstrate the feasibility of the proposed fuzzy programming approaches. Sensitivity analysis and fuzzy simulation are comprehensively utilized to discuss the effects of the fuzzy soft time windows and time uncertainty on the routing optimization and help decision makers to better design a crisp transportation plan that can effectively make tradeoffs among economy, timeliness and reliability.

1. Introduction

Globalization and accompanying international trade enable companies to take advantage of the rich market resources of the entire world, including outsourcing businesses to professional partners to reduce production cost and extending markets to seek more customers to make more profit [1,2,3]. As a result, international freight transportation has a positive trend and its volume will increase fourfold by 2050 [4]. Although the international trade and global commodity circulation bring great opportunities for the growth of companies, challenges also exist. One of the biggest of these is from the transportation industry [5,6]. With the rapid expansion of the companies’ businesses, the distribution channels for their raw materials and products are extended significantly [7]. The long-distance distribution channels enhance the difficulty of transportation organizations, and thereby increase both the cost of logistics and time of accomplishing the transportation orders of companies. As a result, improving the logistics performance is widely acknowledged to be a crucial approach for companies in order to maintain competitiveness in the worldwide market [8,9].

In response to the increasing demand of companies for reducing the logistics cost and time, an advanced transportation mode, namely intermodal transportation, has been widely adopted by large numbers of companies to transport their goods in international trade [10]. Intermodal transportation can be defined as the transportation of containerized cargoes from their origins to associated destinations by using more than one transportation mode, including air, rail, road and water [11,12]. The combination of various transportation modes can form a seamless door-to-door chain that can fully make use of the respective advantages of different modes, which can help enterprises to reduce the cost created in the transportation process [13]. Furthermore, by using the ISO standard containers to carry goods, mechanized operations can be promoted in intermodal transportation. Therefore, the timeliness can be enhanced to improve the service level of the transportation. Currently, intermodal transportation has been widely used in North America [14] and Asia [15]. In Europe, for example in Italy [16], although the road industry is still the main means of freight transportation, freight volume accomplished by intermodal transportation is sustainably growing.

Intermodal transportation is considered as a promising means of efficiently improving the logistics performance [17]. Among the diverse forms of transportation, road-rail intermodal transportation integrates the time-flexible road transportation implemented by container trucks and scheduled rail transportation. It therefore enjoys both the good mobility of trucks on short/medium-distance collection and delivery and the cost efficiency as well as large capacity of rail transportation on long-distance distribution [18,19]. Thus, the road-rail intermodal transportation is the most representative form of the diverse intermodal transportation modes and gets more and more popular in inland transportation. Therefore, in this study, we focus on road-rail intermodal transportation planning. With the road-rail intermodal transportation network becoming mature, how to optimally utilize the existing transportation facilities and equipment in the network to accomplish the transportation orders draws considerable attention from both transportation demanders (e.g., companies with transportation demands) and transportation managers (e.g., intermodal transportation operators) [20]. Consequently, the road-rail intermodal routing problem becomes the forefront in the intermodal transportation planning field [18].

The road-rail intermodal routing problem aims at designing the best origin-to-destination routes that combine container trucks and container trains to enable customers to accomplish their transportation orders. It is much more complex than the famous vehicle routing problem, since two different transportation modes, i.e., road transportation (container trucks) and rail transportation (container trains), should be optimized in a combinatorial way in the same transportation network [21]. Satisfying customer demands is the foundation of the routing optimization, especially when the traditional transportation industry is trying its best to develop into the modern service industry [22]. Therefore, the road-rail intermodal routing investigated by this study is a customer-centred optimization. Generally, the customer demands can be summarized as in Figure 1 [18].

First of all, since the cost created in the logistics activities research is up to nearly 30–50% of the companies’ total production cost [5,8], reducing logistics cost is considered as an effective way for companies to make profits. This motivates the minimization of the costs paid for accomplishing the transportation orders as the optimization objective of the road-rail intermodal routing modelling. Such an objective is established in the modelling by all the relative literature.

Secondly, improving the transportation timeliness to realize on-time transportation is important for companies that need to distribute their materials or products through the extensive intermodal transportation network. In the era when large numbers of companies resort to just-in-time (JIT) strategy to minimize inventory, minimizing time does not always lead to the minimization of costs [23]. Therefore, besides reducing costs paid for accomplishing their transportation orders, they also expect goods delivery at the right time instead of traditional delivery using the least time, i.e., avoiding both early and late delivery [9]. Solving the question of how to formulate and further improve the customer demand on timeliness is thus an important goal related to enhancing the service level of the intermodal transportation and its routing optimization.

Last but not least, customers attach great importance to transportation reliability so that their transportation orders can be accomplished without disruptions, in order that they can reduce opportunity costs. Transportation reliability is significantly influenced by the uncertainty of the operations of the intermodal transportation network [24]. The operation uncertainty leads to time uncertainty. In the road-rail intermodal transportation network, the container trains are operated strictly by fixed schedules and usually get less disruptions [1,25,26]. Consequently, the time uncertainty of rail transportation can be neglected. As a result, in this study, road travel time and loading/unloading time are considered as the sources of time uncertainty. Road travel time uncertainty emerges due to traffic congestion, bad weather and accidents [1,27,28], while loading/unloading time uncertainty results from the unstable proficiency and state of staff that conduct loading/unloading operations, technical issues of the loading/unloading equipment and unpredictable tasks that occupy staff and equipment. These two sources of time uncertainty will not only influence the goods delivery but also disrupt the transshipment between road and rail. They should therefore be modeled in the road-rail intermodal routing optimization to improve the routing reliability.

Above all, in this study, we explore the road-rail intermodal routing problem that is directly oriented on satisfying the customer demands on reducing costs, improving timeliness and enhancing reliability. The contributions made by this study are fivefold.

(1) Fuzzy soft time windows are employed to model the due dates for accomplishing transportation orders. Maximizing the service level associated with the fuzzy soft time windows is set as part of the weighted objective of the road-rail intermodal routing model.

(2) Multiple sources of time uncertainty, i.e., road travel time uncertainty and loading/unloading time uncertainty, are comprehensively modeled by using fuzzy set theory.

(3) A hub-and-spoke network is utilized to model the road-rail intermodal transportation system with time-flexible container truck services and scheduled container train services.

(4) A fuzzy mixed integer nonlinear programming model is constructed to formulate the road-rail intermodal routing problem with fuzzy soft time windows and multiple sources of time uncertainty, and an exact solution approach combining defuzzification and linearization is developed.

(5) Sensitivity analysis and fuzzy simulation are adopted to quantify the effects of the fuzzy soft time windows and the uncertainty of road travel time and loading/unloading time on the road-rail intermodal routing optimization.

The remaining sections of this study are organized as follows. In Section 2, the existing literature on the intermodal routing problem is reviewed to demonstrate our improvements. In Section 3, we present the methods used to model the multiple sources of time uncertainty, the fuzzy soft time windows with respect to the due dates of accomplishing transportation orders and the road-rail intermodal hub-and-spoke transportation system. Based on the modelling foundation proposed in Section 3, we establish a fuzzy mixed integer nonlinear programming model in Section 4 for the road-rail intermodal routing problem that fully considers to satisfy the customer demands on costs, timeliness and reliability. In Section 5, considering the fuzziness of the model, defuzzification is first of all conducted to get a crisp model by using the fuzzy expect value model and fuzzy chance-constrained programming, so that decision makers can obtain crisp road-rail intermodal route planning. Then using the linearization technique developed in our previous study [1,25], linear reformulation of the nonlinear model is realized, so that global optimal solutions to the road-rail intermodal routing problem can be effectively obtained by using an exact solution algorithm that can be implemented by standard mathematical programming software. In Section 6, computational experiment is designed to verify the feasibility of the proposed methods. The effects of the fuzzy soft time windows and the uncertainty of road travel time and loading/unloading time on the road-rail intermodal routing optimization are discussed by using sensitivity analysis and fuzzy simulation. Finally, the conclusions of this study are drawn in Section 7.

2. Literature Review

There are large numbers of research articles that focus on the intermodal routing optimization. These studies have been very successful and so provide a solid foundation for us to continue this study. In order to fully understand the current progress of the existing literature on the customer-centred intermodal routing optimization, a systematic literature review on how existing literature deals with the customer demands is presented. Additionally, how to model the intermodal transportation system is a fundamental issue to optimize the associated routing problem. Consequently, the approaches for such issue will also be reviewed in this section.

2.1. Review on Modeling Customer Demand on Economy in the Intermodal Routing

Lowering costs reflects the customer demand on improving transportation economy and is always the most important objective of intermodal routing optimization. Even in the multi-objective optimization on the intermodal routing problem, for example, the hazardous materials routing problem [29,30] and fresh food routing problem [31], the economic objective plays an important role in the optimization. A few studies on the intermodal routing problem, e.g., Chang [8] and Moccia et al. [32], adopt a piecewise linear cost function that shows the effect of economies of scale to calculate the economic objective. However, the piecewise linear cost function is only dependent on the transported weight [8]. Actually, the costs are not only associated with transported weight, but also related to the time and distance. Moreover, it is difficult to determine the various parameters of the piecewise linear cost function in transportation practice. Therefore, it is not completely applicable to use the piecewise linear cost function to represent the economic objective function. Currently, the majority of the existing studies employ generalized costs to formulate the charges for accomplishing the transportation orders, e.g., Sun et al. [1,18,25,30], Ayar and Yaman [31], Hrušovský et al. [2] and Demir et al. [3]. Generalized costs have a structure shown as Figure 2 [30,33,34] and thus cover all the payment created in the activities of the intermodal transportation. In this study, we continue to use generalized costs to formulate the economic objective of road-rail intermodal routing problem.

2.2. Review on Modeling Customer Demand on Timeliness in the Intermodal Routing

As for the customer demand on improving transportation timeliness, there are various studies on intermodal routing problem that consider the “less time, the better” concept, and hence set minimizing the transportation time as their optimization objective, e.g., Chang [8], Vale and Ribeiro [21], Sun and Chen [35], Cai et al. [36], Xiong and Wang [37] and Xiong and Dong [38]. In these studies, a multi-objective optimization is used to deal with the intermodal routing problem, and minimizing the transportation time is of one of the multiple objectives of the modelling. Tradeoff between minimizing the transportation time and other objectives (e.g., minimizing transportation costs [8,35,36,37,38] or minimizing carbon emissions [21]) is made by using Pareto analysis. However, as emphasized by Dua and Sinha [9], minimizing transportation time is not always the best option to improve transportation timeliness, since nowadays the customers highly expect that their goods can be delivered at the right time to avoid both early and late delivery. In this case, the majority of the existing literature consider improving transportation timeliness by formulating due dates of accomplishing transportation orders. Currently, the methods to formulate the due dates can be summarized as in Figure 3.

As for the formulation methods shown in Figure 3, at the initial stage, using hard time points to describe the due dates is very popular in the intermodal routing optimization. Sun et al. [19,39], Verma and Verter [40] and Wang and Han [41] formulate the due date constraint based on the hard time point formulation that requires that the instants of accomplishing transportation orders by intermodal transportation should not be later than the prescribed due dates. Although using hard time points to formulate the due date constraint can avoid lateness of the delivery, early delivery might occur. Consequently, some studies on the intermodal routing problem propose the concept of soft due dates, so that penalty can be used to reduce the degree of both earliness and lateness of the delivery. Hrušovský et al. [2] and Demir et al. [3] use soft time points to represent the due dates. In their studies, the penalty cost is linear with the degree of the earliness or lateness. By lowering the penalty cost caused by earliness and lateness of the delivery through minimizing the economic objective, the on-time transportation can be enhanced.

As claimed by Sun et al. [42], time windows are more suitable and flexible to model the due dates by using lower and upper bounds to describe the customers’ opinions on the delivery that is neither too early nor too late. Currently, Zhang et al. [43] develop hard time windows to model the due dates and require that the transportation time of each order must fall into the time window. In fact, customers also accept the violation of the due date time windows to a certain degree [42,44,45]. Under this situation, soft time windows receive a lot of attention from the intermodal routing literature when discussing the formulation of due dates. Sun et al. [1,18,30] and Fazayeli et al. [46] explore the intermodal location-routing problem with due dates denoted by soft time windows. In Fazayeli et al.’s study, the penalty cost strategy that is the same as that of Hrušovský et al. [2] and Demir et al. [3] is used to optimize the timeliness of the intermodal transportation. In many cases, however, violation of the soft due date time windows does not lead to a penalty cost [44]. Moreover, it is difficult to determine the value of the unit penalty cost in transportation practice.

The due dates claimed by customers usually involve their subjective opinions [44]. The service level of the intermodal routing on timeliness significantly influences the customer satisfaction. When formulating soft time windows and the corresponding penalty cost strategy is infeasible, it is worthwhile to try to model the service level of the intermodal routing in order to further improve the timeliness by constructing an associated constraint or objective. The most popular method of building the function of customer satisfaction associated with subjective opinions is fuzzy set theory [45]. From that viewpoint of fuzzy set theory, the due dates can be modeled as fuzzy soft time windows that are expressed by trapezoidal fuzzy numbers [44,45]. The customer satisfaction on the instants when the transportation orders are accomplished can be measured by the fuzzy membership of the trapezoidal fuzzy numbers. Such measurement quantifies the service level of the intermodal routing. Currently, to the best of our knowledge, there is no existing literature that investigates the intermodal routing problem with fuzzy soft time windows. But such an idea gets widely discussed in the vehicle routing problem, e.g., Tang et al. [44], Lopez-Castro and Montoya-Torres [45], Ghannadpour et al. [47] and Xu et al. [48]. In this study, we will introduce the fuzzy soft time windows into the intermodal routing problem to improve the service level of the routing optimization. The optimization of the service level will be realized by setting service objective and constructing service level constraint.

2.3. Review on Modeling Customer Demand on Reliability in the Intermodal Routing

Transportation reliability is a measure of the containers being successfully delivered to the receivers by using the planned intermodal routes in actual transportation. It is highly influenced by uncertainty [49]. As explained in Section 1, operational time is the main source of the uncertainty. In the majority of the existing studies on the intermodal routing problem, travel time and loading/unloading time are treated as deterministic parameters. Sun et al. [50] pointed out that both overestimation and underestimation might exist if uncertain parameters are valued in a deterministic way, which will reduce the reliability of the intermodal routing.

However, the intermodal routing problem with operational time uncertainty (specifically, road travel time and loading/unloading time uncertainty in this study) does not receive enough of the highlights it deserves. A few relative studies can be found so far. The current research progress on the intermodal routing problem with operational time uncertainty is illustrated by Figure 4. Hrušovský et al. [2] and Demir et al. [3] explore the green intermodal routing problem with travel time uncertainty. Stochastic programming is employed by their study to model the uncertain travel time. Similar study on the sea-rail intermodal routing problem with stochastic travel time is then conducted by Zhao et al. [51], in which a chance-constrained stochastic approach is provided under the assumption that the travel time follows a normal distribution. Uddin and Huynh [24,52] also underline the importance of formulating road travel time uncertainty in improving the reliability of the road-rail intermodal routing, and develop a multicommodity model based on stochastic programming. Besides stochastic programming, Sun et al. [1] adopt the time-dependent travel time represented by the continuous piecewise linear function to describe the road travel time uncertainty caused by traffic congestion. Further they integrate the time-dependent road travel time with the intermodal routing problem. Resat and Turkay [17], Liu et al. [26] and Tang et al. [53] use the classical BPR (Bureau of Public Roads) equation to estimate the road travel time in road-rail intermodal transportation.

The above studies only consider the uncertainty of travel time, while neglecting the loading/unloading time uncertainty that also exists in the routing decision-making stage. Considering multiple sources of uncertainty however might help decision makers find improved intermodal routes [1]. As a result, in this study, we will consider both road travel time uncertainty and loading/unloading time uncertainty when modeling the road-rail intermodal routing, so that road-rail intermodal routes with improved reliability can be provided to decision makers.

Moreover, as indicated by Figure 4, stochastic programming, time-dependent programming and the BPR equation all have their own disadvantages in dealing with the uncertainty [27,54], which considerably reduces their feasibility in dealing with the optimization problems with uncertainty. To the best of our knowledge, fuzzy programming as an alternative method to address uncertainty has not been considered in the intermodal routing literature that involves uncertain issues. Based on the fuzzy set theory introduced by Zadeh [55], fuzzy programming has been widely considered as a useful tool to deal with uncertainty [56]. In the fuzzy programming, the uncertain parameters or variables can be effectively valued by the pessimistic, optimistic and most likely estimations from experts and decision makers based on their knowledge and experience [50]. Consequently, compared with the rest three methods, fuzzy programming does not rely on the size of the existing data or the real-time information of a system and is hence more flexible in practical decision making. Although less attention has been paid on the fuzzy programming in dealing with the intermodal routing problem with time uncertainty, it has already been employed by Sun et al. [1] to model capacity uncertainty and by Sun et al. [18] to formulate demand uncertainty in the intermodal routing problem, and shows good feasibility. Additionally, it has been widely adopted by the vehicle routing problem with travel time uncertainty, e.g., Zhang and Liu [27], Zarandi et al. [54] and Djadane et al. [57]. Since fuzzy programming can combine limited historical data with the expert knowledge and experience to effectively estimate the uncertainty, in this study, we will use fuzzy programming to model the multiple sources of time uncertainty that exist in the road-rail intermodal routing problem.

2.4. Review on Modeling Intermodal Transportation System

How to model the intermodal transportation system is the foundation of the routing optimization. First of all, an efficient intermodal transportation network where the routes are planned should be constructed. The consolidation of the intermodal transportation network is the most important factor that influences the efficiency of the transportation [6]. The four kinds of consolidation networks, i.e., point-to-point, line, collection-and-distribution and hub-and-spoke networks, can be seen from Macharis and Bontekoning’s work [6]. The hub-and-spoke network is widely acknowledged to be the most suitable consolidation network to construct the road-rail intermodal transportation system [18]. The superiority of the road-rail intermodal hub-and-spoke network is its systematic integration of the different advantages of the time-flexible road transportation and scheduled rail transportation in the pre haul-long haul-end haul transportation chain, which is stressed by Wang et al. [19] and Sun et al. [1,18,30].

However, modeling the road-rail transportation system as a hub-and-spoke network does not receive enough attention from researchers in the intermodal routing field. A few studies, e.g., Sun et al. [1,18] and Uddin Huynh [24], consider the hub-and-spoke structure as part of the transportation network when modeling the intermodal transportation system. Sun et al. [30] adopted the hub-and-spoke network to represent the road-rail intermodal transportation network when exploring the bi-objective optimization on the hazardous materials routing problem. Although such consolidation network is not very popular in the road-rail intermodal routing problem, it has had great importance attached to the intermodal transportation design problem, e.g., Wang et al. [19], Meng and Wang [58], Yang et al. [59], Lin and Chen [60,61] and Konings et al. [62]. Consequently, in this study, we will model the road-rail intermodal transportation system as a hub-and-spoke network and further optimize the associated road-rail intermodal routing problem.

Furthermore, the routing optimization is used to support the practical decision making on the transportation organization. As a result, the modeling of the road-rail intermodal transportation system should match the real world. As stated in Section 1, the intermodal routing problem is more complex than the vehicle routing problem, since it involves more than one transportation mode and different transportation modes are operated in different ways. In the road-rail transportation system, road transportation is time flexible, while rail transportation should be operated with fixed schedules [6,8]. Therefore, when modeling the road-rail intermodal routing problem, a combination of the two different transportation modes by transshipping operations at rail terminal (hubs) in the road-rail intermodal hub-and-spoke network should therefore follow the restrictions of the schedules of the selected container trains. Moreover, the arrival and departure of the containers carried by container trains should also observe the schedules.

At the beginning of the intermodal routing optimization, schedules are not formulated in the intermodal routing problem. The majority of the existing studies, e.g., Yiping et al. [36], Xiong and Wang [37], Xiong and Dong [38], Zhang et al. [43], Çakır [63], Chang et al. [64], consider that after arriving at the transshipping node by one transportation mode, the containers are immediately unloaded, then loaded on to another transportation mode once the unloading operation is finished, and finally depart from the transshipping node when the loading operation is finished. Since the schedules of some transportation modes really exist in the transportation practice and need to be modeled in the intermodal routing problem, modeling schedules become part of the routing optimization. A few studies, e.g., Hrušovský et al. [2], Demir et al. [3], Liu et al. [65] and Lin [66], formulate the scheduled departure instant constraint that the arrival instant of containers at the transshipping node should not be later than the scheduled departure instant of the successive transportation mode with schedules. Limited studies can be found that consider the scheduled service time window of transportation modes with schedules when optimizing the intermodal routing problem, e.g., Chang [8], Ayar and Yaman [31] and Moccia et al. [32]. In transportation practice, schedules regulate more than the departure instant and service/operation time window. As summarized by our previous study [1,25,30], in the road-rail intermodal transportation scenario, the schedules of rail transportation contain operation time windows, arrival instants and departure instants of freight trains at and from the rail terminals covered on their running routes as well as their operational periods. Therefore, in order to establish a road-rail intermodal transportation system that represents the real-world transportation organization, we will model the above-mentioned contents of schedules by referring to our previous study.

2.5. Review Summaries

As we can see from Section 2.1, Section 2.2, Section 2.3 and Section 2.4, the existing literature on the intermodal routing problem has shown that sustainable accomplishments have been achieved, especially in optimizing the customer demand on transportation economy. But as we can see from the literature review, research gaps still exist in the following aspects:

(1) Optimizing the customer demand on transportation efficiency by formulating due dates of accomplishing transportation orders in a more applicable way.

(2) Optimizing the customer demand on transportation reliability by modeling multiple sources of time uncertainty, i.e., road travel time uncertainty and loading/unloading time uncertainty, in a more feasible programming method.

(3) Modeling an efficient road-rail intermodal transportation system that also matches the real-world transportation organization to provide a solid foundation for the routing modelling.

In this study, we focus on the customer-centred freight routing optimization in the road-rail intermodal transportation network and aim at bridging these research gaps through the way that is presented in Section 1 (see the fivefold contributions made by this study).

3. Methodology

In this study, we extend the road-rail intermodal routing problem by modeling multiple sources of time uncertainty, i.e., road travel time and loading/unloading time uncertainty, to improve reliability and considering fuzzy soft time windows to improve service level associated with timeliness. The routing problem is further oriented on a road-rail intermodal hub-and-spoke transportation system. In this section, we present the methods for modeling time uncertainty, fuzzy soft time windows and the road-rail intermodal transportation system.

3.1. Modeling Time Uncertainty by Fuzzy Set Theory

As mentioned in Section 1, considering the limitation of using stochastic programming to deal with uncertainty, in this study, we adopt fuzzy set theory to model the two kinds of time uncertainty by using fuzzy numbers. In practical decision-making, different decision makers might hold different viewpoints on estimating the values of fuzzy parameters, including pessimistic, optimistic and most likely estimations [30,67].

Interval, triangular and trapezoidal fuzzy numbers can be used to represent fuzzy parameters [18,68]. In this study, we select triangular fuzzy numbers to model the above characteristics of the estimation of the fuzzy parameters (specifically, road travel time and loading/unloading time in this study) due to the following two reasons.

(1) Triangular fuzzy numbers are simpler and more flexible in the fuzzy arithmetic operations than the other two kinds of fuzzy numbers [69,70].

(2) Triangular fuzzy numbers can match the estimations held by different decision makers to fully reflect the practical decision-making scenario under fuzzy environment, which can be seen in Figure 5.

For a triangular fuzzy number $\tilde{t}=\left(t_1, t_2, t_3\right)$ representing fuzzy road travel time or fuzzy loading/unloading time, its three prominent points are defined as follows.

(1) $t_1$ is the most optimistic estimation, and corresponds to the best case that the traffic conditions of the road transportation are extremely good, or the loading/loading operations are conducted quite smoothly.

(2) $t_2$ is the most likely estimation, and corresponds to the most likely case that shows what the traffic conditions of the road transportation or the loading/unloading operations usually are.

(3) $t_3$ is the most pessimistic estimation, and corresponds to the worst case that the traffic conditions of the road transportation are extremely bad (for example, severe congestion occurs), or there are severe technical or operational issues that happen to disrupt the loading/unloading operations (for example, equipment breakdown happens).

3.2. Modeling the Due Dates of Transportation Orders by Fuzzy Soft Time Windows

As stated in Section 1, in this study, we employ fuzzy soft time windows to model the due dates of accomplishing transportation orders, so that the timeliness of road-rail intermodal routing can be effectively improved. First of all, there need time intervals to represent the instants of accomplishing transportation orders that are neither too early nor too late for customers. When the instant of accomplishing a transportation order falls into such a time interval, the satisfaction of the customer reaches the highest. Otherwise the satisfaction will be reduced. Moreover, since the customers cannot accept that the transportation orders are accomplished too early or too lately. Therefore, the fuzzy soft time windows have lower bounds and upper bounds separately representing the earliest and latest instants that the customers can accept. Above all, the fuzzy soft time window has four prominent points, and shows the same representation as the trapezoidal fuzzy numbers. The fuzzy membership function of the trapezoidal fuzzy number can therefore be used to measure the customers’ satisfaction level [44]. A fuzzy soft time window can be denoted by $\tilde{T}=\left(T_1, T_2, T_3, T_4\right)$ the fuzzy membership of which is illustrated by Figure 6. The four prominent points of a fuzzy soft time window are given as follows.

(1) $T_1$ is the endurable earliest instant of accomplishing a transportation order. The customer cannot accept the instant of accomplishing a transportation order that is earlier than $T_1$.

(2) $\left[T_2, T_3\right]$ is the preferred instant range of accomplishing a transportation order. The customer considers the transportation timeliness reaches its maximum, in other words, his or her satisfaction level reaches maximum (100%), if the instant of accomplishing a transportation order falls into $\left[T_2, T_3\right]$.

(3) $T_4$ is the endurable latest instant of accomplishing a transportation order. The customer cannot accept the instant of accomplishing a transportation order that is later than $T_4$.

Given an instant of accomplishing a transportation order denoted by $T$, the corresponding satisfaction level $\mu \left(T\right)$ is calculated by Equation (1) [44].

$$\mu \left(T\right)=\{\begin{array}{c}\begin{array}{c}\begin{array}{cc}\frac{T-T_1}{T_2-T_1}& T_1\le T<T_2\end{array}\\ \begin{array}{cc}1& T_2\le T\le T_3\end{array}\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{T_4-T}{T_4-T_3}& T_3<T\le T_4\end{array}\\ \begin{array}{cc}0& \mathrm{otherwise}\end{array}\end{array}\end{array}$$

(1)

Under the above setting, it is acceptable that the instant of accomplishing a transportation order falls into $\left[T_1, T_2\right]$ or $\left[T_3, T_4\right]$, which however reduces the customer’s satisfaction level. If the customer requires a satisfaction level that is no less than $\eta$, the endurable instant range of accomplishing the transportation order is $\left[T_{2\eta }, T_{3\eta }\right]$ (see Figure 6), where $T_{2\eta }=\eta ·\left(T_2-T_1\right)+T_1$ and $T_{3\eta }=T_4-\eta ·\left(T_4-T_3\right)$ according to Equation (1).

In this study, the service level associated with fuzzy soft time windows in the road-rail intermodal routing is optimized by the following two aspects.

(1) Maximizing the total service levels of all the transportation orders is set as part of the weighted objective of the routing model.

(2) A service level constraint is established to ensure the service level of each transportation order is not lower than a satisfaction degree requested by corresponding customer.

3.3. Modeling the Road-Rail Intermodal Transportartion System

As claimed in Section 2, a hub-and-spoke network shown as Figure 7 is used to model the road-rail intermodal transportation system. Rail terminals installed with loading/unloading equipment are the hub where transshipment between container trucks and container trains is realized. Origins and destinations are the spokes, and are connected with rail terminals by container truck transportation [19,30]. If containers from various transportation orders at a rail terminal are gathered together, the operations of the rail terminal can be centralized to lead to economies of scale [19], which is the main advantage of the road-rail intermodal hub-and-spoke network.

Another issue is how to coordinate various transportation modes with different operational characteristics by routing to generate origin-to-destination intermodal routes that are feasible in both space and time. In the road-rail intermodal transportation system, road transportation is time flexible and thus offers good mobility, while rail transportation is operated by fixed schedules that regulate the operation time windows, arrival and departure instants of a container train at the rail terminals covered on its running route as well as its operational period [1,25,30].

Therefore, in order to design feasible road-rail intermodal routes, the operations of container trucks should coordinate with the schedules of the container trains, especially during the “Arriving → Transshipping → Departing” operation process at rail terminals. Figure 8 indicates the road-rail intermodal transportation process in a hub-and-spoke network that captures the operational characteristics of road and rail transportation [30]. Such process will be formulated in the modelling.

In the transportation process illustrated by Figure 8, there are three kinds of costs that are created in the road-rail intermodal transportation, i.e., travel costs, loading and unloading costs and storage costs. The economic objective of the road-rail intermodal routing optimization is also to minimize the sum of these costs paid for accomplishing all the transportation orders.

As we can see from Figure 8, the arrival instants of containers by trucks are related to the road travel time, and the storage time of containers at rail terminals are determined by the arrival instants of containers and the loading/unloading time. Obviously, the uncertainty of road travel time and loading/unloading time lead to uncertainty of the arrival instant and storage time. Therefore, there are two parameters (i.e., road travel time and loading/unloading time) and two decision variables (i.e., arrival instants of containers by trucks and storage time of containers at rail terminals) in the optimization model that are uncertain.

4. A Fuzzy Mixed Integer Nonlinear Programming Model

4.1. Modelling Characteristics

The optimization on the routing problem is based on a given framework that contains a series of modelling characteristics [10,30,71]. According to Sun et al. [30,71] and Kumar et al. [10], the modelling of a freight routing problem should comprehensively consider the characteristics listed in Figure 9. It also systematically indicates how we select the modelling characteristics to establish our optimization model.

4.2. Symbols Used to Establish the Model

The symbols used to establish the optimization model are similar to the ones presented in our previous study that provides a solid modelling foundation for us to carry out this study [30] and are introduced in

Table 1. Notations.

Categories	Symbols	Representations
TRANSPORTATION ORDERS	$P$	Transportation order sets served by the road-rail intermodal routing.
	$p\in P$	Index of transportation orders.
	$o_p$ & $d_p$	Indexes of the origin and destination of the container transportation, respectively.
	$q_p$	Volume in TEU of transportation order $p$.
	$t_p^0$	Release instant of the containers of transportation order $p$ at origin $o_p$ (also the instant when the containers of transportation order $p$ is prepared to get loaded to start the pre-haul process).
	$\left[T{W}_p^1, T{W}_p^2, T{W}_p^3, T{W}_p^4\right]$	Fuzzy soft time window of transportation order $p$ claimed by the customer that requests this transportation order.
	${\mu }_p(\ast )$	Service level of transportation order $p$ whose input is the loading completion instant of containers of transportation order $p$ at its destination.
	${\eta }_p\in \left[0, 1\right]$	Endurable lowest service level of transportation order $p$.
ROAD-RAIL INTERMODAL TRANSPORTATION NETWORK	$V=\left(N, A, S\right)$	Node-arc-based network, where $N$, $A$ and $S$ are the node set, arc set and transportation service set of the network.
	$i$, $j$ & $k\in N$	Indexes of the nodes.
	$m$ & $n\in S$	Indexes of transportation services.
	$\left(i, j\right)\in A$	A direct arc from $i$ and $j$, on which there are transportation services.
	$S_{ij}$, $Truc{k}_{ij}$ & $Rai{l}_{ij}$	Transportation service set, road (container truck) service set and rail (container train) service set on arc $\left(i, j\right)$, respectively, and $S_{ij}=Truc{k}_{ij}{\displaystyle \cup }Rai{l}_{ij}\subseteq S$.
	$N_i^{-}$ & $N_i^{+}\subseteq N$	Predecessor node set and successor node set to node $i$, respectively.
PARAMETERS	$d_{ijm}$	Travel distance in km in TEU of transportation service $m$ on arc $\left(i, j\right)$
	$ca{p}_{ijm}$	Transportation capacity in TEU of transportation service $m$ on arc $\left(i, j\right)$.
	${\tilde{t}}_i^m=\left(t_1{}_i^m,t_2{}_i^m, t_3{}_i^m\right)$	Fuzzy loading/unloading time in hour per TEU of transportation $m$ at node $i$.
	$\left[t{w}_i^{m-}, t{w}_i^{m+}\right]$	Operation time window of rail service $m$ at rail terminal $i$, where $t{w}_i^{m-}$ is the operation start instant of the loading/unloading operation and $t{w}_i^{m+}$ is the corresponding operation cutoff instant (as we can see in Figure 8, although schedules of container trains regulate many parameters that control their operations, only operation time windows of container trains at rail terminals influence the routing).
	${\tilde{t}}_{ijm}=\left(t_{ijm}^1,t_{ijm}^2, t_{ijm}^3\right)$	Fuzzy travel time of road service $m$ in hour on arc $\left(i, j\right)$.
	$c_m$	Travel costs per TEU per km of transportation service $m$,
	$e_m$	Separate loading/unloading costs per TEU of transportation service $m$.
	$f_{store}$	Storage costs per TEU per hour at rail terminal.
DECISION VARIABLES	$x_{ijm}^p$	If transportation service $m$ on arc $\left(i, j\right)$ is used to move containers of transportation order $p$, $x_{ijm}^p=1$, otherwise $x_{ijm}^p=0.$
	${\tilde{y}}_i^p=\left(y_1{}_i^p, y_2{}_i^p, y_3{}_i^p\right)$	Non-negative fuzzy variable: fuzzy arrival instant of the containers of transportation order $p$ at node $i$ by container trucks.
	${\tilde{z}}_{ijm}^p=\left(z_1{}_{ijm}^p, z_2{}_{ijm}^p, z_3{}_{ijm}^p\right)$	Non-negative fuzzy variable: fuzzy storage time of container of transportation order $p$ at node $i$ before being moved by container train $m$ on arc $\left(i, j\right)$.

. It should be noted that there are two fuzzy decision variables in the model, and the their fuzziness has been explained in Section 3.3.

4.3. Objective Functions of the Road-Rail Intermodal Routing Problem

• Economic Objective

  $$\begin{array}{c}\mathrm{minimize} {\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{c}_m·q_p·d_{ijm}·x_{ijm}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{i\in N_i^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{e}_m·q_p·x_{ijm}^p+{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{k\in N_i^{+}}}{\displaystyle {\displaystyle \sum }_{n\in S_{jk}}}{e}_n·q_p·x_{jkn}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·{\tilde{z}}_{ijm}^p\end{array}$$

  (2)
  The economic objective is to minimize the total generalized costs of accomplishing all the transportation orders that contains travel costs, loading and unloading costs and storage cost.
• Service Objective

  $$\mathrm{maximize} {\displaystyle \sum }_{p\in P}{\mu }_p\left(E\left[{\tilde{y}}_{d_p}^p+{\displaystyle {\displaystyle \sum }_{i\in N_{d_p}^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Truc{k}_{i{d}_p}}}{\tilde{t}}_{d_p}^m·q_p·x_{i{d}_{p}m}^p\right]\right)$$

  (3)
  The service objective is to maximize the service levels of all the transportation orders to improve the customers’ satisfaction degrees. Since the road travel time and unloading time are fuzzy, the instants of accomplishing transportation orders are also fuzzy. In this study, we use the fuzzy expected value of the instant of accomplishing a transportation order as the input to calculate the service level.
• Weighted Objective

  $$\begin{array}{c}\mathrm{minimize} {\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{c}_m·q_p·d_{ijm}·x_{ijm}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{i\in N_i^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{e}_m·q_p·x_{ijm}^p+{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{k\in N_i^{+}}}{\displaystyle {\displaystyle \sum }_{n\in S_{jk}}}{e}_n·q_p·x_{jkn}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·{\tilde{z}}_{ijm}^p\\ -W·{\displaystyle {\displaystyle \sum }_{p\in P}}{\mu }_p\left(E\left[{\tilde{y}}_{d_p}^p+{\displaystyle {\displaystyle \sum }_{i\in N_{d_p}^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Truc{k}_{i{d}_p}}}{\tilde{t}}_{d_p}^m·q_p·x_{i{d}_{p}m}^p\right]\right)\end{array}$$

  (4)
  The economic objective and service objective can be linearly combined together by distributing a weight to the service objective. The combination is shown as Equation (4). The value of the weight is determined by the decision makers according to their preference to the service objective.

4.4. Constraints of the Road-Rail Intermodal Routing Problem

$$\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{i\in N_i^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{x}_{ijm}^p-{\displaystyle {\displaystyle \sum }_{k\in N_i^{+}}}{\displaystyle {\displaystyle \sum }_{n\in S_{jk}}}{x}_{jkn}^p=\{\begin{array}{cc}-1& \forall j=o_p\\ 0& \forall j\in N\backslash \left\{o_p, d_p\right\}\\ 1& \forall j=d_p\end{array}& \forall p\in P\end{array}$$

(5)

$$\begin{array}{ccc}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{x}_{ijm}^p\le 1& \forall p\in P& \forall \left(i, j\right)\in A\end{array}$$

(6)

Equation (5) is the flow conservation constraint to ensure the integrity of an origin-to-destination route for each transportation order. The combination of Equation (5) and Equation (6) ensures the containers in each transportation order are unsplittable.

$$\begin{array}{cc}\begin{array}{cc}\left(t_p^0+t_1{}_{o_p}^n·q_p+t_{o_{p}kn}^1-y_1{}_k^p\right)·x_{o_{p}kn}^p=0& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(7)

$$\begin{array}{cc}\begin{array}{cc}\left(t_p^0+t_2{}_{o_p}^n·q_p+t_{o_{p}kn}^2-y_2{}_k^p\right)·x_{o_{p}kn}^p=0& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(8)

$$\begin{array}{cc}\begin{array}{cc}\left(t_p^0+t_3{}_{o_p}^n·q_p+t_{o_{p}kn}^3-y_3{}_k^p\right)·x_{o_{p}kn}^p=0& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(9)

$$\begin{array}{c}\left[{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_1{}_j^m·q_p\right)·x_{ijm}^p+t_1{}_j^n·q_p+t_{j{d}_{p}n}^1-y_1{}_{d_p}^p\right]·x_{j{d}_{p}n}^p=0\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(10)

$$\begin{array}{c}\left[{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_2{}_j^m·q_p\right)·x_{ijm}^p+t_2{}_j^n·q_p+t_{j{d}_{p}n}^2-y_2{}_{d_p}^p\right]·x_{j{d}_{p}n}^p=0\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(11)

$$\begin{array}{c}\left[{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_3{}_j^m·q_p\right)·x_{ijm}^p+t_3{}_j^n·q_p+t_{j{d}_{p}n}^3-y_3{}_{d_p}^p\right]·x_{j{d}_{p}n}^p=0\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(12)

$$\begin{array}{cc}\begin{array}{cc}{\tilde{y}}_k^p=\left(y_1{}_k^p, y_2{}_k^p, y_3{}_k^p\right)& \forall p\in P\end{array}& \forall k\in N/\left(\left\{o_p\right\}{\displaystyle \cup }{N}_{d_p}^{-}\right)\end{array}$$

(13)

Equations (7)–(9) are used to calculate the three prominent points of the fuzzy arrival instant of containers at the rail terminal by container trucks after departing from the origin (pre-haul). Equations (10)–(13) calculate the three prominent points of the fuzzy arrival instant of containers at the destination by container trucks after departing from the rail terminal (end-haul). After obtaining the prominent points, the fuzzy arrival instant represented by a triangular fuzzy number can be determined by Equation (13).

$$\begin{array}{c}\left[\max \left\{t{w}_j^{n-}-\left(y_1{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_1{}_j^m·q_p·x_{ijm}^p\right), 0\right\}-z_3{}_{jkn}^p\right]·x_{jkn}^p=0\\ \begin{array}{ccc}\forall p\in P& \forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(14)

$$\begin{array}{c}\left[\max \left\{t{w}_j^{n-}-\left(y_2{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_2{}_j^m·q_p·x_{ijm}^p\right), 0\right\}-z_2{}_{jkn}^p\right]·x_{jkn}^p=0\\ \begin{array}{ccc}\forall p\in P& \forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(15)

$$\begin{array}{c}\left[\max \left\{t{w}_j^{n-}-\left(y_3{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_3{}_j^m·q_p·x_{ijm}^p\right), 0\right\}-z_1{}_{jkn}^p\right]·x_{jkn}^p=0\\ \begin{array}{ccc}\forall p\in P& \forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(16)

$$\begin{array}{cc}\begin{array}{cc}{\tilde{z}}_{jkn}^p=\left(z_1{}_{jkn}^p, z_2{}_{jkn}^p, z_3{}_{jkn}^p\right)& \forall p\in P\end{array}& \begin{array}{cc}\forall \left(j, k\right)\in A& \forall n\in R_{ij}\end{array}\end{array}$$

(17)

The three prominent points of the fuzzy storage time at the rail terminal can be further calculated by Equations (14)–(16) to obtain the triangular fuzzy number representation of the fuzzy storage time by Equation (17).

$$\begin{array}{c}{\tilde{y}}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{\tilde{t}}_j^m·q_p·x_{ijm}^p+{\tilde{z}}_{jkn}^p+{\tilde{t}}_j^n·q_p\le t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)\\ \begin{array}{ccc}\forall p\in P& \forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(18)

In Equation (18), $ϵ$ is a large enough positive constant. Equation (18) ensures that the loading completion instant of the containers at the rail terminal should not be later than the operation cutoff instant (upper bound of the operation time window) of the selected container train. Obviously, Equation (18) is a fuzzy constraint.

$$\begin{array}{ccc}{\displaystyle {\displaystyle \sum }_{p\in P}}{q}_p·x_{ijm}^p\le ca{p}_{ijm}& \forall \left(i, j\right)\in A& \forall m\in S_{ij}\end{array}$$

(19)

Equation (19) is the transportation capacity constraint that ensures the total volume loaded on a common container train or a common group of truck fleets does not exceed the allowable capacity.

$$\begin{array}{cc}{\mu }_p\left(E\left[{\tilde{y}}_{d_p}^p+{\displaystyle {\displaystyle \sum }_{i\in N_{d_p}^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Truc{k}_{i{d}_p}}}{\tilde{t}}_{d_p}^m·q_p·x_{i{d}_{p}m}^p\right]\right)\ge {\eta }_p& \forall p\in P\end{array}$$

(20)

Equation (20) is the service level constraint that ensures the service level associated with the expected instant of accomplishing each transportation order should not be lower than the customer’s satisfaction degree.

$$\begin{array}{cc}\begin{array}{cc}{x}_{ijm}^p\in \left\{0, 1\right\}& \forall p\in P\end{array}& \begin{array}{cc} \forall \left(i, j\right)\in A& \forall m\in S_{ij}\end{array}\end{array}$$

(21)

$$\begin{array}{ccc}{y}_3{}_i^p\ge y_2{}_i^p\ge y_1{}_i^p\ge 0& \forall p\in P& \forall i\in N\end{array}$$

(22)

$$\begin{array}{cc}\begin{array}{cc}{z}_3{}_{jkn}^p\ge z_2{}_{jkn}^p\ge z_1{}_{jkn}^p\ge 0& \forall p\in P\end{array}& \begin{array}{cc}\forall \left(i, j\right)\in A& \forall m\in Rai{l}_{ij}\end{array}\end{array}$$

(23)

Equations (21)–(23) are the domain constraints of the variables.

5. Solution Approaches

The aim of designing the approaches is to generate an equivalent mixed integer linear programming model, so that the global optimal solutions to the road-rail intermodal routing problem can be obtained by using an exact solution algorithm run by mathematical programming software. Consequently, the solution approaches contain the following two parts:

(1) Defuzzification that is used to generate a crisp optimization model that can provide decision makers with crisp route planning.

(2) Linearization that is used to generate an equivalent linear programming model that can be effectively solved by using an exact solution algorithm to get the global optimal solutions.

5.1. Defuzzification of the Fuzzy Constraint

Equation (18) is a fuzzy constraint, since it contains fuzzy decision variables ${\tilde{y}}_j^p$ and ${\tilde{z}}_{jkn}^p$ as well as fuzzy parameters ${\tilde{t}}_j^m$ and ${\tilde{t}}_j^n$. To realize the defuzzification of Equation (18), fuzzy chance-constrained programming is applied in this study to reformulate the fuzzy constraint. Compared with the fuzzy possibility and fuzzy necessity measures, the fuzzy credibility measure is self-dual and can ensure that a fuzzy event must hold if its credibility equals 1 while it must fail if its credibility equals 0 [27,50]. Therefore, fuzzy credibility is more suitable than the other two measures to be used to construct a fuzzy chance constraint.

Based on the fuzzy credibility measure, the fuzzy chance constraint of Equation (18) is as Equation (24). In this equation, $\mathsf{\alpha }$ is a confidence level determined by decision makers subjectively. Equation (24) ensures the credibility of the fuzzy event in $\mathrm{Cr}\{\ast \}$ should not be lower than a given credibility level.

$$\begin{array}{c}\mathrm{Cr}\left\{{\tilde{y}}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{\tilde{t}}_j^m·q_p·x_{ijm}^p+{\tilde{z}}_{jkn}^p+{\tilde{t}}_j^n·q_p\le t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)\right\}\ge \mathsf{\alpha }\\ \begin{array}{ccc}\forall p\in P& \forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(24)

The left-hand function of Equation (24) can be rewritten as Equation (25).

$$\mathrm{Cr}\left\{\right(\begin{array}{c}t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)-\left(y_3{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_3{}_j^m·q_p·x_{ijm}^p+z_3{}_{jkn}^p+t_3{}_j^n·q_p\right),\\ t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)-\left(y_2{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_2{}_j^m·q_p·x_{ijm}^p+z_2{}_{jkn}^p+t_2{}_j^n·q_p\right),\\ t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)-\left(y_1{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_1{}_j^m·q_p·x_{ijm}^p+z_1{}_{jkn}^p+t_1{}_j^n·q_p\right)\end{array})\ge 0\}$$

(25)

As indicated by Cao and Lai [72] and Zheng and Liu [27], given a deterministic number $r$ and a triangular fuzzy number $\tilde{d}=\left(d_1, d_2, d_3\right)$, there exists Equation (26).

$$\mathrm{Cr}\left\{\tilde{d}\ge r\right\}=\{\begin{array}{c}\begin{array}{cc}1& \mathrm{if} r\le d_1\\ \frac{2d_2-d_1-r}{2\left(d_2-d_1\right)}& \mathrm{if} d_1\le r\le d_2\end{array}\\ \begin{array}{cc}\frac{d_3-r}{2\left(d_3-d_2\right)}& \mathrm{if} d_2\le r\le d_3\\ 0& \mathrm{if} r\ge d_3\end{array}\end{array}$$

(26)

Furthermore, based on Equation (26), $\mathrm{Cr}\left\{\tilde{d}\ge r\right\}\ge \alpha$ equals Equations (27) and (28) according to proof proposed by Sun et al. [1] and Wang et al. [19].

$$\begin{array}{cc}2\mathsf{\alpha }·d_2-\left(2\mathsf{\alpha }-1\right)·d_3\ge a& \mathrm{if} 0\le \mathsf{\alpha }\le \frac{1}{2}\end{array}$$

(27)

$$\begin{array}{cc}2\left(1-\mathsf{\alpha }\right)·d_2+\left(2\mathsf{\alpha }-1\right)·d_1\ge a& \mathrm{if} \frac{1}{2}<\mathsf{\alpha }\le 1\end{array}$$

(28)

Accordingly, the fuzzy chance constraint Equation (24) can be reformulated as Equations (29) and (30) that are crisp and linear.

$$\begin{array}{c}2\mathsf{\alpha }·\left[t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)-\left(y_2{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_2{}_j^m·q_p·x_{ijm}^p+z_2{}_{jkn}^p+t_2{}_j^n·q_p\right)\right]\\ \ge \left(2\mathsf{\alpha }-1\right)·\left[t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)-\left(y_1{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_1{}_j^m·q_p·x_{ijm}^p+z_1{}_{jkn}^p+t_1{}_j^n·q_p\right)\right]\\ \begin{array}{ccc}\begin{array}{cc}\mathrm{if} 0\le \mathsf{\alpha }\le \frac{1}{2}& \forall p\in P\end{array}& \forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(29)

$$\begin{array}{c}2\left(1-\mathsf{\alpha }\right)·\left[t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)-\left(y_2{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_2{}_j^m·q_p·x_{ijm}^p+z_2{}_{jkn}^p+t_2{}_j^n·q_p\right)\right]\\ \ge \left(1-2\mathsf{\alpha }\right)·\left[t{w}_j^{n+}·x_{jkn}^p+ϵ·\left(1-x_{jkn}^p\right)-\left(y_3{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_3{}_j^m·q_p·x_{ijm}^p+z_3{}_{jkn}^p+t_3{}_j^n·q_p\right)\right]\\ \begin{array}{ccc}\begin{array}{cc}\mathrm{if} \frac{1}{2}<\mathsf{\alpha }\le 1& \forall p\in P\end{array}& \forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(30)

5.2. Defuzzification of the Fuzzy Objective

The optimization objective is fuzzy, because it contains fuzzy decision variables ${\tilde{y}}_j^p$ and ${\tilde{z}}_{jkn}^p$. Therefore, it is also necessary to realize the defuzzification of the fuzzy objective in order to obtain crisp solutions to the road-rail intermodal routing problem. Currently, there are two defuzzification methods that can be used to deal with the fuzzy objective, including fuzzy expected value model and the fuzzy chance-constrained programming. However, there is no solid evidence demonstrating that one method is better than the other or the two methods have same performance. As a result, in this study, we will use both defuzzification methods to address the fuzzy objective and also test their performance in dealing with the road-rail intermodal routing problem with time uncertainty in the case study section.

5.2.1. Using the Fuzzy Expected Value Model

The expected value of a triangular fuzzy number $\tilde{d}=\left(d_1, d_2, d_3\right)$ is expressed by Equation (31) [73].

$$E\left[\tilde{d}\right]=\frac{d_1+2d_2+d_3}{4}$$

(31)

The fuzzy expected value model is widely acknowledged to be an effective approach to deal with a fuzzy objective [18]. It aims at minimizing or maximizing the expected value of the fuzzy objective. According to Equation (30), the fuzzy expected value of the total storage costs for all transportation orders can be obtained by Equation (32).

$$\begin{array}{c}E\left[{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·{\tilde{z}}_{ijm}^p\right]\\ ={\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·E\left[{\tilde{z}}_{ijm}^p\right]\\ ={\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·\frac{z_1{}_{jkn}^p+2z_2{}_{jkn}^p+z_3{}_{jkn}^p}{4}\end{array}$$

(32)

The fuzzy expected value of the weighted objective is therefore given by Equations (33) and (34) by replacing the fuzzy storage costs with the corresponding fuzzy expected value.

$$\begin{array}{c}\mathrm{minimize} {\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{c}_m·q_p·d_{ijm}·x_{ijm}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{i\in N_i^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{e}_m·q_p·x_{ijm}^p+{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{k\in N_i^{+}}}{\displaystyle {\displaystyle \sum }_{n\in S_{jk}}}{e}_n·q_p·x_{jkn}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·\frac{z_1{}_{jkn}^p+2z_2{}_{jkn}^p+z_3{}_{jkn}^p}{4}-W·{\displaystyle {\displaystyle \sum }_{p\in P}}{\mu }_p\left(\frac{v_p}{4}\right)\end{array}$$

(33)

$$\begin{array}{c}{v}_p=y_1{}_{d_p}^p+{\displaystyle {\displaystyle \sum }_{i\in N_{d_p}^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Truc{k}_{i{d}_p}}}{t}_1{}_{d_p}^m·q_p·x_{i{d}_{p}m}^p+2\left(y_2{}_{d_p}^p+{\displaystyle {\displaystyle \sum }_{i\in N_{d_p}^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Truc{k}_{i{d}_p}}}{t}_2{}_{d_p}^m·q_p·x_{i{d}_{p}m}^p\right)\\ \begin{array}{cc}+y_3{}_{d_p}^p+{\displaystyle {\displaystyle \sum }_{i\in N_{d_p}^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Truc{k}_{i{d}_p}}}{t}_3{}_{d_p}^m·q_p·x_{i{d}_{p}m}^p& \forall p\in P\end{array}\end{array}$$

(34)

5.2.2. Using Fuzzy Chance-Constrained Programming

Fuzzy chance-constrained programming can be also used to deal with the fuzzy objective [74]. In this case, Equation (4) can be reformulated as Equations (35) and (36) where $\psi$ is a non-negative auxiliary variable.

$$\begin{array}{c}\mathrm{minimize} {\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{c}_m·q_p·d_{ijm}·x_{ijm}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{i\in N_i^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{e}_m·q_p·x_{ijm}^p+{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{k\in N_i^{+}}}{\displaystyle {\displaystyle \sum }_{n\in S_{jk}}}{e}_n·q_p·x_{jkn}^p\\ +\psi -W·{\displaystyle {\displaystyle \sum }_{p\in P}}{\mu }_p\left(\frac{v_p}{4}\right)\end{array}$$

(35)

$$\mathrm{Cr}\left\{{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·{\tilde{z}}_{ijm}^p\le \psi \right\}\ge \mathsf{\alpha }$$

(36)

Similar to the crisp reformulation of Equation (24), Equation (36) is equivalent to Equations (37) and (38).

$$\begin{array}{cc}2\mathsf{\alpha }·{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·z_2{}_{ijm}^p+\left(1-2\mathsf{\alpha }\right)·{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·z_1{}_{ijm}^p\le \psi & \mathrm{if} 0\le \mathsf{\alpha }\le \frac{1}{2}\end{array}$$

(37)

$$\begin{array}{cc}2\left(1-\mathsf{\alpha }\right)·{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·z_2{}_{ijm}^p+\left(2\mathsf{\alpha }-1\right)·{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·z_3{}_{ijm}^p\le \psi & \mathrm{if} \frac{1}{2}<\mathsf{\alpha }\le 1\end{array}$$

(38)

5.3. Linear Reformulation of the Service Objective

The service level function ${\mu }_p\left(\frac{v_p}{4}\right)$ is a continuous piecewise linear function. Therefore, the service objective function ${{\displaystyle \sum }}_{p\in P}{\mu }_p\left(\frac{v_p}{4}\right)$ is nonlinear.

Constrained by Equation (20), for $\forall p\in P$, $\frac{v_p}{4}$ must fall into range $\left[T{W}_p^1, T{W}_p^4\right]$, and can only fall into only one of its three sub ranges shown as Figure 10. First of all, we define a 0–1 auxiliary variable $w_g^p$: if $\frac{v_p}{4}$ falls into range $g$ ($g\in G=\left\{1, 2, 3\right\}$), $w_g^p=1$, otherwise, $w_g^p=0$. Hence we have Equations (39) and (40).

$$\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{g\in G}}{w}_g^p=1& \forall p\in P\end{array}$$

(39)

$$\begin{array}{ccc}{w}_g^p\in \left\{0, 1\right\}& \forall p\in P& \forall g\in \left\{1, 2, 3\right\}\end{array}$$

(40)

Then we distribute non-negative auxiliary variables ${\xi }_g^{p-}$ and ${\xi }_g^{p+}$ to the lower bound and upper bound of range $g$, and so we have Equations (41)–(43).

$$\begin{array}{ccc}{\xi }_g^{p-}+{\xi }_g^{p+}=w_g^p& \forall p\in P& \forall g\in \left\{1, 2, 3\right\}\end{array}$$

(41)

$$\begin{array}{ccc}{\xi }_g^{p-}\ge 0& \forall p\in P& \forall g\in \left\{1, 2, 3\right\}\end{array}$$

(42)

$$\begin{array}{ccc}{\xi }_g^{p+}\ge 0& \forall p\in P& \forall g\in \left\{1, 2, 3\right\}\end{array}$$

(43)

Furthermore, $\frac{v_p}{4}$ can be expressed by Equation (44).

$$\begin{array}{cc}\frac{v_p}{4}={\displaystyle {\displaystyle \sum }_{g\in G}}\left({\xi }_g^{p-}·T{W}_p^g+{\xi }_g^{p+}·T{W}_p^{g+1}\right)& \forall p\in P\end{array}$$

(44)

By using Equation (44) to replace $\frac{v_p}{4}$ in the service objective function ${{\displaystyle \sum }}_{p\in P}{\mu }_p\left(\frac{v_p}{4}\right)$, we have the following formula.

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{p\in P}}{\mu }_p\left(\frac{v_p}{4}\right)={\displaystyle {\displaystyle \sum }_{p\in P}}{\mu }_p\left({\displaystyle {\displaystyle \sum }_{g\in G}}\left({\xi }_g^{p-}·T{W}_p^g+{\xi }_g^{p+}·T{W}_p^{g+1}\right)\right)\\ ={\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{g\in G}}{\mu }_p\left({\xi }_g^{p-}·T{W}_p^g+{\xi }_g^{p+}·T{W}_p^{g+1}\right)={\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{g\in G}}\left[{\mu }_p\left({\xi }_g^{p-}·T{W}_p^g\right)+{\mu }_p\left({\xi }_g^{p+}·T{W}_p^{g+1}\right)\right]\\ ={\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{g\in G}}\left[{\xi }_g^{p-}·{\mu }_p\left(T{W}_p^g\right)+{\xi }_g^{p+}·{\mu }_p\left(T{W}_p^{g+1}\right)\right]={\displaystyle {\displaystyle \sum }_{p\in P}}\left({\xi }_1^{p+}+{\xi }_2^{p-}+{\xi }_2^{p+}+{\xi }_3^{p-}\right)\end{array}$$

(45)

As a result, the linearized crisp weighted objective function under the fuzzy expected value model is as in Equation (46).

$$\begin{array}{c}\mathrm{minimize} {\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{c}_m·q_p·d_{ijm}·x_{ijm}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{i\in N_i^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{e}_m·q_p·x_{ijm}^p+{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{k\in N_i^{+}}}{\displaystyle {\displaystyle \sum }_{n\in S_{jk}}}{e}_n·q_p·x_{jkn}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}{f}_{store}·q_p·\frac{z_1{}_{jkn}^p+2z_2{}_{jkn}^p+z_3{}_{jkn}^p}{4}-W·{\displaystyle {\displaystyle \sum }_{p\in P}}\left({\xi }_1^{p+}+{\xi }_2^{p-}+{\xi }_2^{p+}+{\xi }_3^{p-}\right)\end{array}$$

(46)

Under fuzzy chance-constrained programming, the linearized crisp weighted objective function is as in Equation (47).

$$\begin{array}{c}\mathrm{minimize} {\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{\left(i, j\right)\in A}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{c}_m·q_p·d_{ijm}·x_{ijm}^p\\ +{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{i\in N_i^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{e}_m·q_p·x_{ijm}^p+{\displaystyle {\displaystyle \sum }_{p\in P}}{\displaystyle {\displaystyle \sum }_{j\in N}}{\displaystyle {\displaystyle \sum }_{k\in N_i^{+}}}{\displaystyle {\displaystyle \sum }_{n\in S_{jk}}}{e}_n·q_p·x_{jkn}^p\\ +\psi -W·{\displaystyle {\displaystyle \sum }_{p\in P}}\left({\xi }_1^{p+}+{\xi }_2^{p-}+{\xi }_2^{p+}+{\xi }_3^{p-}\right)\end{array}$$

(47)

5.4. Linear Reformulation of the Nonlinear Constraints

According to Sun et al. [25], nonlinear constraints Equations (7)–(12) can be linearized as Equations (48)–(59).

$$\begin{array}{cc}\begin{array}{cc}{t}_p^0+t_1{}_{o_p}^n·q_p+t_{o_{p}kn}^1-y_1{}_k^p\ge ϵ·\left(x_{o_{p}kn}^p-1\right)& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(48)

$$\begin{array}{cc}\begin{array}{cc}{t}_p^0+t_1{}_{o_p}^n·q_p+t_{o_{p}kn}^1-y_1{}_k^p\le ϵ·\left(1-x_{o_{p}kn}^p\right)& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(49)

$$\begin{array}{cc}\begin{array}{cc}{t}_p^0+t_2{}_{o_p}^n·q_p+t_{o_{p}kn}^2-y_2{}_k^p\ge ϵ·\left(x_{o_{p}kn}^p-1\right)& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(50)

$$\begin{array}{cc}\begin{array}{cc}{t}_p^0+t_2{}_{o_p}^n·q_p+t_{o_{p}kn}^2-y_2{}_k^p\le ϵ·\left(1-x_{o_{p}kn}^p\right)& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(51)

$$\begin{array}{cc}\begin{array}{cc}{t}_p^0+t_3{}_{o_p}^n·q_p+t_{o_{p}kn}^3-y_3{}_k^p\ge ϵ·\left(x_{o_{p}kn}^p-1\right)& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(52)

$$\begin{array}{cc}\begin{array}{cc}{t}_p^0+t_3{}_{o_p}^n·q_p+t_{o_{p}kn}^3-y_3{}_k^p\le ϵ·\left(1-x_{o_{p}kn}^p\right)& \forall p\in P\end{array}& \begin{array}{cc}\forall k\in N_{o_p}^{+}& \forall n\in Truc{k}_{o_{p}k}\end{array}\end{array}$$

(53)

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_1{}_j^m·q_p\right)·x_{ijm}^p+t_1{}_j^n·q_p+t_{j{d}_{p}n}^1-y_1{}_{d_p}^p\ge ϵ·\left(x_{j{d}_{p}n}^p-1\right)\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(54)

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_1{}_j^m·q_p\right)·x_{ijm}^p+t_1{}_j^n·q_p+t_{j{d}_{p}n}^1-y_1{}_{d_p}^p\le ϵ·\left(1-x_{j{d}_{p}n}^p\right)\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(55)

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_2{}_j^m·q_p\right)·x_{ijm}^p+t_2{}_j^n·q_p+t_{j{d}_{p}n}^2-y_2{}_{d_p}^p\ge ϵ·\left(x_{j{d}_{p}n}^p-1\right)\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(56)

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_2{}_j^m·q_p\right)·x_{ijm}^p+t_2{}_j^n·q_p+t_{j{d}_{p}n}^2-y_2{}_{d_p}^p\le ϵ·\left(1-x_{j{d}_{p}n}^p\right)\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(57)

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_3{}_j^m·q_p\right)·x_{ijm}^p+t_3{}_j^n·q_p+t_{j{d}_{p}n}^3-y_3{}_{d_p}^p\ge ϵ·\left(x_{j{d}_{p}n}^p-1\right)\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(58)

$$\begin{array}{c}{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in Rai{l}_{ij}}}\left(t{w}_j^{m-}+t_3{}_j^m·q_p\right)·x_{ijm}^p+t_3{}_j^n·q_p+t_{j{d}_{p}n}^3-y_3{}_{d_p}^p\le ϵ·\left(1-x_{j{d}_{p}n}^p\right)\\ \begin{array}{cc}\begin{array}{cc}\forall p\in P& \forall j\in N_{d_p}^{-}\end{array}& \forall n\in Truc{k}_{j{d}_p}\end{array}\end{array}$$

(59)

According to Sun et al. [25], nonlinear constraints Equations (14)–(16) are equivalent to the following linear equations.

$$\begin{array}{ccc}{z}_3{}_{jkn}^p\ge t{w}_j^{n-}-\left(y_1{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_1{}_j^m·q_p·x_{ijm}^p\right)+ϵ·\left(x_{jkn}^p-1\right)& \forall p\in P& \begin{array}{cc}\forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(60)

$$\begin{array}{cc}\begin{array}{cc}{z}_3{}_{jkn}^p\le ϵ·x_{jkn}^p& \forall p\in P\end{array}& \begin{array}{cc}\forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(61)

$$\begin{array}{ccc}{z}_2{}_{jkn}^p\ge t{w}_j^{n-}-\left(y_2{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_2{}_j^m·q_p·x_{ijm}^p\right)+ϵ·\left(x_{jkn}^p-1\right)& \forall p\in P& \begin{array}{cc}\forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(62)

$$\begin{array}{cc}\begin{array}{cc}{z}_2{}_{jkn}^p\le ϵ·x_{jkn}^p& \forall p\in P\end{array}& \begin{array}{cc}\forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(63)

$$\begin{array}{ccc}{z}_1{}_{jkn}^p\ge t{w}_j^{n-}-\left(y_3{}_j^p+{\displaystyle {\displaystyle \sum }_{i\in N_j^{-}}}{\displaystyle {\displaystyle \sum }_{m\in S_{ij}}}{t}_3{}_j^m·q_p·x_{ijm}^p\right)+ϵ·\left(x_{jkn}^p-1\right)& \forall p\in P& \begin{array}{cc}\forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(64)

$$\begin{array}{cc}\begin{array}{cc}{z}_1{}_{jkn}^p\le ϵ·x_{jkn}^p& \forall p\in P\end{array}& \begin{array}{cc}\forall \left(j, k\right)\in A& \forall n\in Rai{l}_{ij}\end{array}\end{array}$$

(65)

The service level constraint Equation (20) is also nonlinear, since it contains continuous piecewise linear function. Based on the description in Section 3.2, it can be reformulated as the following two equations that are both linear.

$$\begin{array}{cc}\frac{v_p}{4}\ge {\eta }_p·\left(T{W}_p^2-T{W}_p^1\right)+T{W}_p^1& \forall p\in P\end{array}$$

(66)

$$\begin{array}{cc}\frac{v_p}{4}\le T{W}_p^4-{\eta }_p·\left(T{W}_p^4-T{W}_p^3\right)& \forall p\in P\end{array}$$

(67)

5.5. Equivalent Mixed Integer Linear Programming Model

After defuzzification and linearization, the fuzzy mixed integer nonlinear model can be modified into two kinds of MILP models. The process that generates the two models is illustrated by Figure 11.

As we can see from Figure 11, two different crisp linear programming models can be obtained after using two different defuzzification methods to deal with the fuzzy objective and employing same linearization techniques to remove the nonlinearity. We use the fuzzy expected value model to defuzzy the initial fuzzy objective to obtain the MILP model I, and the fuzzy chance-constrained programming to generate the MILP model II. The differences of the two models are that they have different crisp objective functions and the MILP model II has two more auxiliary constraints (i.e., Equations (37) and (38)) deriving from the process of using fuzzy chance-constrained programming to deal with the fuzzy objective. Apart from Equations (37) and (38), the two MILP models share the same constraints.

6. Computational Experiment

We modify the numerical case designed in our previous study [30] to make it match the specific routing problem discussed by this study, and use the modified case to demonstrate the feasibility of proposed methods in dealing with the road-rail intermodal routing problem with fuzzy soft time windows and multiple sources of time uncertainty. The detailed description of the modified numerical case is presented in Appendix A (see Figure A1 and

Table A1. The schedules of the container trains in the network.

Container Trains	Train Routes	Operation Start Instants at Origins	Operation Cutoff Instants at Origins	Operation Start Instants at Destinations	Capacities in TEU	Periods in Train per Day	Travel Distances in km
1	(4, 7)	15	30	40	300	1	184
2	(4, 7)	6	23	32	350	1	184
3	(4, 8)	5	20	31	380	1	210
4	(4, 8)	3	21	34	400	1	210
5	(4, 9)	9	26	42	320	1	280
6	(4, 9)	14	31	46	300	1	280
7	(5, 7)	8	24	35	335	1	193
8	(5, 7)	13	29	41	355	1	193
9	(5, 8)	7	25	35	400	1	185
10	(5, 8)	10	26	37	325	1	185
11	(5, 9)	0	18	32	292	1	178
12	(5, 9)	6	25	35	270	1	178
13	(6, 7)	11	30	47	290	1	315
14	(6, 7)	14	29	50	315	1	315
15	(6, 8)	4	21	32	362	1	236
16	(6, 8)	8	26	40	345	1	236
17	(6, 9)	12	27	39	262	1	291
18	(6, 9)	15	33	45	285	1	291

,

Table A2. The capacities, travel times, and travel costs of the container truck fleet groups in the network.

Container Truck Fleet Groups	Arcs	Capacities in TEU	Fuzzy Travel Times in Hour	Travel Distances in km
19	(1, 4)	150	1.0, 1.5, 2.8	68
20	(1, 5)	245	1.4, 2.0, 4.2	85
21	(1, 6)	200	3.6, 5.0, 8.3	120
22	(2, 4)	190	2.4, 3.0, 4.8	90
23	(2, 5)	230	1.7, 2.8, 5.2	105
24	(2, 6)	145	1.2, 1.8, 3.8	75
25	(3, 4)	160	3.3, 5.6, 6.5	114
26	(3, 5)	205	1.6, 2.2, 3.6	94
27	(3, 6)	170	0.8, 2.0, 3.8	100
28	(7, 10)	165	1.7, 2.3, 3.5	95
29	(7, 11)	205	2.5, 3.2, 5.5	118
30	(7, 12)	150	4.6, 6.0, 9.2	130
31	(8, 10)	175	3.5, 4.3, 5.5	106
32	(8, 11)	210	1.0, 1.7, 3.1	64
33	(8, 12)	221	1.6, 2.4, 4.0	93
34	(9, 10)	144	4.6, 6.2, 8.5	122
35	(9, 11)	195	2.1, 3.7, 5.8	102
36	(9, 12)	132	0.9, 1.4, 3.1	70

and

Table A3. The values of the cost and operation parameters defined in Table 3.

Transportation Mode	Transportation Costs in CNY per TEU per km	Separate Loading and Unloading Costs in CNY per TEU	Storage Costs in CNY per TEU per Hour	Separate Fuzzy Loading and Unloading Time in CNY per TEU
Rail	2.025	195	3.125	0.05, 0.1, 0.15
Road	6	25	-	0.1, 0.2, 0.25

for details). The following works will also be untaken in this section in addition to demonstrating the feasibility:

(1) Exploring whether and how fuzzy soft time windows and fuzziness of both road travel time and loading/unloading time influence the road-rail intermodal routing optimization.

(2) Comparing the fuzzy expected value model and fuzzy chance-constrained programming in dealing with the fuzzy objective.

(3) Helping decision makers to identify the optimal value of confidence level $\mathsf{\alpha }$, so that a crisp road-rail intermodal transportation scheme can be provided to them.

6.1. Computational Environment

In this study, we use the standard Branch-and-Bound algorithm that is a famous exact solution algorithm in operations research to solve the two equivalent MILP models that formulate the specific routing problem discussed by this study. The Branch-and-Bound algorithm is run by the mathematical programming software LINGO version 12.0 [75]. All the computation is performed on a ThinkPad Laptop with Intel Core i5-5200U, 2.20 GHz CPU, and 8 GB RAM.

In this study, we provide two MILP models to solve problem. MILP model I uses fuzzy expected value model to deal with the fuzzy objective, while MILP model II employs fuzzy chance-constrained programming. Apart from the differences deriving from the defuzzification of the fuzzy objective, the two MILP models yield the same formulations. When using two MILP models to optimize the road-rail intermodal problem, the computational scale is shown in

Table 2. Computational scale of the routing problem formulated by two MILP models.

Models	Variables	Integer Variables	Constraints
MILP model I	2547	696	5999
MILP model II	2549	696	6001

. As we can see from Table 2, MILP model II has two more variables than MILP model I, because non-negative auxiliary variable $\psi$ emerges twice (one in the crisp objective Equation (47), and one in constraint Equation (37) or (38)) during the defuzzification of the fuzzy objective to obtain MILP model II, while there is no auxiliary variable generated when obtaining the crisp objective of MILP model I. Moreover, MILP model II has two more constraints than MILP model I, since the defuzzification for generating the crisp objective of MILP model II lead to two auxiliary constraints, including the constraint regarding the non-negativity of auxiliary variable $\psi$ and the auxiliary constraint Equation (37) (or (38)). However, MILP model I does not generate any auxiliary constraints during the defuzzification of fuzzy objective.

When the confidence level $\mathsf{\alpha }$ is set to 0.9, the service levels ${\eta }_p$ ($\forall p\in P$) are set to 0.7 and the weight distributed to the service objective is set to 1000, the computational times of the models in 10 times computation are presented as follows. MILP model I has computational times that vary from 33 s to 34 s with an average value of 33.4 s in the 10 times computation. As for MILP model II, the maximum and minimum computational times are 26 s and 25 s, respectively. The average computational time of MILP model II is 25.6 s, which saves ~23.4% of the computational time compared with MILP model I. Therefore, using the Branch-and-Bound algorithm to solve the two MILP models to generate the best solutions to the routing problem can be accomplished within quite a short time. MILP model II that uses fuzzy chance-constrained programming to deal with the fuzzy objective shows better efficiency than MILP model I that uses the fuzzy expect value model to address the fuzzy objective.

6.2. Sensitivity Analysis of the Routing Optimization with Respect to the Weight of the Service Objective

The weight distributed to the service objective reflects the decision makers’ preference for the improvement on the transportation timeliness. In this study, we analyze the sensitivity of the routing optimization with respect to the weight associated with the service objective to reveal the relationship between the two optimization objectives. The following analysis is under the setting that confidence level $\mathsf{\alpha }$ is set to 0.9 and service levels ${\eta }_p$ ($\forall p\in P$) are set to 0.7. The changes of the two objectives with respect to the weight, i.e., the Pareto frontiers to the bi-objective optimization on the road-rail intermodal routing problem, are indicated by Figure 12.

It can be observed from Figure 12 that by enhancing the weight distributed to the service objective, the value of the economic objective increases, which means that the transportation economy is reduced. Contrary to the transportation economy, the increase of the weight leads to the improved transportation timeliness, since the value of the service objective increases. Consequently, the economic objective and service objective of the road-rail intermodal routing problem cannot reach their respective optimum simultaneously, and enhancing one objective will worsen the other. A solution with minimum costs does not yield maximum service level. As a result, the road-rail intermodal routing problem yields non dominated solutions (i.e., Pareto solutions) instead of dominated solutions.

As shown in Figure 12, it is impossible to simultaneously satisfy the customer demands on transportation economy and timeliness by planning the road-rail intermodal routes. A tradeoff should therefore be made between the economic and service objectives based on the decision makers’ preference for the service improvement when planning the road-rail intermodal routes. The Pareto frontiers shown in Figure 12 can help decision makers to select the preferred solution. In practical decision making, decision makers can utilize the multiple-criteria decision-making methods, e.g., AHP method [76] and TOPSIS method [77], to preciously select the Pareto solution that best matches the decision-making situation.

Moreover, it can be noticed that MILP model II yields larger values of the economic objective than MILP model I, while the two models have the same values of the service objective. Although the two models have different values of the economic objective, the planned best road-rail intermodal routes given by them are exactly the same, which is demonstrated in the following Section 6.5.2. This leads to the same values of the service objective of the two models, since they share the same service objective function. However, the economic objective functions of the two models are different, which results in the difference of the values of the economic objective between the two models.

6.3. Sensitivity Analysis of the Routing Optimization with Respect to the Service Level

Besides the weight distributed to the service objective, the service level is also set by the decision makers manually. Its settings might influence the optimization results. In this study, we also use sensitivity analysis to explore its effect on the routing. We first of all set the weight as 1000 and the confidence level as 0.9. Then we modify service levels ${\eta }_p$ ($\forall p\in P$) from 0.1 to 1.0 with a step of 0.1, and calculate the best solutions corresponding to each value of the service level by solving MILP model I and MILP model II. The results of the sensitivity analysis are shown as Figure 13.

Figure 13 indicates that the road-rail intermodal routing optimization is sensitive to the service level. As we can see from Figure 13, by increasing the service level ${\eta }_p$, both the values of the economic objective and service objectives show a growth tendency, which means that the transportation timeliness of the routing is getting improved, while its transportation economy is getting worse. Enhancing the service level shows an effect on the routing optimization that is the same as improving the weight distributed to the service objective.

Therefore, when customers prefer or propose a higher service level for accomplishing their transportation orders, enhancing the service level ${\eta }_p$ will help decision makers to find the best road-rail intermodal routes with improved transportation timeliness, on condition that customers could accept more costs paid for accomplishing transportation orders. With the help of Figure 13, decision makers can determine the value of the service level ${\eta }_p$ in a more objective way in order to avoid that the service level is overestimated and the economy of the routing is considerably reduced.

As well as the tendency shown in Figure 12, Figure 13 also shows that MILP model II provides larger values of the economic objective than MILP model I and the two models have same values of service objective. The reasons that lead to the tendency indicated by Figure 13 are same to the ones that result in Figure 12. Figure 13 points out that although the variation of the service level ${\eta }_p$ leads to the modification of the planned best routes, the two models always generate the same solutions under each value of the service level ${\eta }_p$.

6.4. Sensitivity Analysis of the Routing Optimization with Respect to the Confidence Level

The confidence level $\mathsf{\alpha }$ is another important parameter that is determined by decision makers according to their subjective preference. In this study, we set the weight distributed to the service objective as 1000 and the service level ${\eta }_p$ ($\forall p\in P$) as 0.5. Under the above situation, we analyze the sensitivity of the routing optimization with respect to the confidence level. We vary the confidence level $\mathsf{\alpha }$ from 0.3 to 0.9 with a step of 0.1, and generate the best solutions to the routing problem by solving the two MILP models. The results of the sensitivity analysis are illustrated by Figure 14. It should be noted that there is no feasible solution that can be found when the confidence level $\mathsf{\alpha }$ is set to 1.0.

As we can see from Figure 14, the road-rail intermodal routing optimization is sensitive to the confidence level, and increasing the confidence level leads to the increase of both economic objective and service objective. Since the larger confidence level means higher transportation reliability, Figure 14 also reveals that improving the transportation reliability will help to improve its transportation timeliness.

However, the value of the economic objective also increases with respect to the confidence level, which means that the transportation economy of the routing is in conflict with its transportation reliability. If customers prefer a road-rail intermodal route plan with higher transportation reliability, they should spend more on accomplishing their transportation orders. Figure 14 provides a solid foundation for decision makers to make effective tradeoffs among the economic objective, service objective and transportation reliability in order to avoid that the overestimation of one objective significantly reduces other objectives in the practical decision making.

Moreover, similar to Figure 12 and Figure 13, in the sensitivity analysis summarized in Figure 14, the two MILP models yield same service objective. However, their economic objectives do not provide the same value, and the value of the economic objective of MILP model II is not always larger than that of MILP model I. Figure 14 shows such tendency due to the same reasons that motivate the tendencies shown in Figure 12 and Figure 13. Moreover, as we can see from Figure 14, the two models yield the same route plan under each value of the confidence level, which is similar to the tendency shown in Figure 13.

6.5. Fuzzy Simulation to Gain the Crisp Road-Rail Intermodal Route Plan

In practical routing decision making, the decision makers who include transportation experts, transportation companies and customers will first of all determine the weight distributed to the service objective and the lowest service level that customers can accept. In this case, we set ${\eta }_p$ as 0.5 and $W$ as 1000. Then if they want to obtain a crisp road-rail intermodal route plan, the confidence level that influences whether the routes are feasible in practice should be determined. Furthermore, the decision makers also need to determine which MILP model is more suitable. Considering the above decision making requirements, we design a fuzzy simulation method to help decision makers overcome the above issues in order to generate the best crisp road-rail intermodal route plan.

In this study, we randomly generate 10 sets of the deterministic road travel times and loading/unloading times according to the fuzzy membership function of the triangular fuzzy number as shown in Figure 5. Then we can obtain 10 deterministic cases that can be used to test the performances of different confidence levels and the two MILP models. The fuzzy simulation is run according to Figure 15 [1,50].

6.5.1. Testing the Feasibility of the Planned Best Road-Rail Intermodal Routes in the Deterministic Cases

In this study, we first of all test the feasibility of the road-rail intermodal routes provided by the fuzzy programming models under different confidence levels in the 10 deterministic cases. When moving the containers along the planned road-rail intermodal routes in the deterministic cases, the route plan is feasible if the constraints of the upper bounds of the operation time windows of covered container trains (i.e., the deterministic formulation of Equation (18)) and the capacity constraints are satisfied. Otherwise, the route plan is infeasible and failed, which means that the transportation orders cannot be accomplished by using the planned routes. The feasible ratios of the route plans provided by the two MILP models under difference confidence levels in the 10 deterministic cases are shown in Figure 16.

As we can see from Figure 16, the route plans generated by solving the two MILP models show the same feasibility ratios in the 10 deterministic cases, regardless of the values of the confidence level. Furthermore, Figure 16 shows that with the increase of the confidence level, the feasibility ratios of the route plans, i.e., the transportation reliability of the routing, tend to be improved. However, in the cases set by this study, the improvement is not very sensitive.

When the confidence level is set to 0.5, 0.6, 0.7, 0.8 or 0.9, the route plans are feasible in all the deterministic cases. In the practical routing decision making, customers always demand that their transportation orders can be accomplished successfully by the planned routes to avoid rerouting that might cause much higher costs [1]. Therefore, as for the case study in this article, considering the customer demand on transportation reliability, the confidence level can be 0.5, 0.6, 0.7, 0.8 or 0.9.

6.5.2. Analyzing the Gaps Between the Planned Best Road-Rail Intermodal Routes and the Actual Best Routes in the Deterministic Cases

In the 10 deterministic cases, we replace the fuzzy parameters and fuzzy variables in the model presented in Section 4 with their deterministic forms, and thus get a deterministic road-rail intermodal routing model that can be linearized by using the linearization techniques given in Section 5. Consequently, the actual best road-rail intermodal routes in the 10 deterministic cases can be generated by solving the deterministic model.

We can also calculate the values of the economic objective and the service objective of the planned best road-rail intermodal routes when they are used to move containers in the 10 deterministic cases. Then we can compare the gaps between planned best road-rail intermodal routes with the actual best routes in the 10 deterministic cases. The computational results are indicated by Figure 17 and Figure 18.

Shown as Figure 17 and Figure 18, we notice that the two MILP models yield the same computational results. Consequently, we can draw the conclusion that the performances of the two MILP models in dealing with the specific routing problem are the same. Either of them can be used to optimize the problem.

Based on the computation results shown in Figure 17 and Figure 18, we can calculate the root mean square (RMS) of the objectives of the planned best routes provided by MILP model I with respect to the objectives of the actual best routes in the 10 deterministic cases. The calculation results which are based on the associated RMS are given in

Table 3. Objective gaps between the planned best routes provided by MILP model I and the actual best routes.

Confidence Levels	RMS in 10 Thousand CNY Regarding Economic Objective	RMS Regarding Service Objective
0.5 and 0.6	1.112	0.981
0.7, 0.8 and 0.9	0.384	1.047

.

Compared with the confidence level of 0.5 and 0.6, the MILP model I with the confidence levels of 0.7, 0.8 or 0.9 can significantly bridge the economic objective gap by ~65.5% by slightly extending the service objective gap by ~6.7%. As a result, when using the MILP model I to optimize the road-rail intermodal routing problem, the confidence level is recommended to be 0.7, 0.8 or 0.9. Above all, as for the numerical case introduced in the Appendix A, the best road-rail intermodal routes that can be used in the real-world transportation under ${\eta }_p$ = 0.5 and $W$ = 1000 can be planned by solving any one of the two MILP models illustrated by Figure 11 with confidence level of 0.7, 0.8 or 0.9.

Furthermore, when confidence level is set to 0.5, MILP model II is converted into the deterministic optimization model used by the existing literature to deal with the road-rail intermodal routing with time certainty [25]. In the studies on the deterministic intermodal routing problem, the road travel time and loading/unloading time are valued in a deterministic way, i.e., by most likely estimations [18,23,25,30,37,43]. Therefore, indicated by Table 3, compared with the solutions generated by the deterministic optimization model, MILP models (i.e., models for the road-rail intermodal routing problem with time fuzziness) can obtain the crisp best road-rail intermodal routes that better match the actual best situation in the transportation practice. And considering the multiple sources of time uncertainty helps to improve the overall performance of the road-rail intermodal routing.

7. Conclusions

In this study, we focus on modeling and optimizing a customer-centred freight routing problem in the road-rail intermodal hub-and-spoke network with fuzzy soft time windows and multiple sources of time uncertainty. The following contributions are made by this study in order to improve the problem optimization:

(1) We employ fuzzy soft time windows to represent the due dates of accomplishing transportation orders. Maximizing the service level associated with the fuzzy soft time windows is set as the optimization objective and a service level constraint is also established, so that the transportation timeliness can be improved.

(2) We simultaneously consider the road travel time uncertainty and loading/unloading time uncertainty in the road-rail intermodal routing problem. The combination of the multiple sources of time uncertainty helps to enhance the transportation reliability of the routing.

(3) We model the road-rail intermodal transportation system as a hub-and-spoke network that contains time-flexible road transportation and scheduled rail transportation.

(4) We propose a bi-objective mixed integer nonlinear programming model the objectives of which are addressed by a weighting method to formulate the specific routing problem discussed in this study. By using the fuzzy expected value model, fuzzy chance-constrained programming and linearization techniques, two equivalent mixed integer linear programming models are generated and can be effectively solved by mathematical programming software.

In the case study, we make full use of sensitivity analysis and fuzzy simulation to quantify the effects of the fuzzy soft time windows and time uncertainty on the routing optimization. The following managerial implications can be deduced.

(1) The economic objective and service objective are in conflict with each other, i.e., the routing optimization cannot satisfy the customer demands on economy and timeliness simultaneously. An effective tradeoff between the two objectives can be made by using the Pareto frontier to the bi-objective optimization problem.

(2) When customers prefer or propose higher service levels for accomplishing their transportation orders, enhancing service level ${\eta }_p$ will help decision makers to find the best road-rail intermodal routes with improved transportation timeliness, on condition that customers could accept more costs paid for accomplishing transportation orders.

(3) Time uncertainty (fuzziness in this study) has significant effect on the routing optimization. A larger confidence level that means higher transportation reliability leads to improved service level (i.e., transportation timeliness), however, worsens the transportation economy of the routing.

(4) The fuzzy expected value model has the same performance as the fuzzy chance-constrained programming in dealing with the fuzzy objectives.

(5) By using the sensitivity analysis and fuzzy simulation designed in this study, decision makers can identify the best value(s) of the confidence level to provide a crisp road-rail intermodal route plan to the customers.

In the future work, we will first devote to exploring the green road-rail intermodal routing problem that considers to optimize the carbon dioxide emissions. We will discuss the green routing problem under different carbon emission regulations, including carbon tax regulation and carbon cap-and-trade regulation. The performances of the two regulations will be compared. Additionally, we will try to add other sources of uncertainty, e.g., demand uncertainty and capacity uncertainty, into the routing problem to further enhance the transportation reliability. Last but not least, developing a heuristic algorithm to efficiently solve the large-scale road-rail intermodal routing problem that is an NP-hard problem is also a direction of our future work.