A Fuzzy Multi-Objective Routing Model for Managing Hazardous Materials Door-to-Door Transportation in the Road-Rail Multimodal Network With Uncertain Demand and Improved Service Level

A routing optimization with risk control and reduction provides a quantitative support for an efficient organization of the hazardous materials transportation where road-rail multimodal transportation plays an important role. Oriented on the viewpoint of the multimodal operator, this study investigates a hazardous materials routing problem for their door-to-door transportation through a road-rail multimodal hub-and-spoke network. It is difficult to determine the exact demand of hazardous materials for each transportation order during the advanced routing decision making. Therefore, demand uncertainty is introduced into the problem, in which the uncertain demands are modeled by triangular fuzzy numbers. To enhance the quality of service of the transportation, the hazardous materials routing optimization in this study uses fuzzy soft time windows to optimize the last-mile delivery of hazardous materials that determines the timeliness of transportation. Accordingly, this study proposes a fuzzy multi-objective routing model that is linear and considers the constraints of the customers’ preferred service levels. To solve the proposed model, its crisp reformulation is first achieved by using the fuzzy expected value method and Jimenez’s fuzzy ranking method. Then, the $\varepsilon $ -constraint method is adopted by this study to obtain the Pareto frontier for the multi-objective optimization problem. Finally, an experimental case study is carried out to verify the proposed modeling of this study. By using sensitivity analysis and fuzzy simulation, the conflicting relationship among the economic, risk, reliability and timeliness objectives that improving one objective worsens the others is clearly indicated, and the insights that helps decision makers to make effective tradeoffs among these objectives to get a most suitable hazardous materials road-rail multimodal transportation plan are also summarized.


I. INTRODUCTION
The production and transportation of hazardous materials is rapidly growing in many industrialized countries, e.g., U.S. [1] and China [2]. The dangerous nature of hazardous materials indicated in Zhao et al.'s work [3]

results in
The associate editor coordinating the review of this manuscript and approving it for publication was Wen-Sheng Zhao . extensive risks in the hazardous materials transportation. Once hazardous materials transportation accidents occur, negative affects on human health, environment sustainability and social stability are severe. However, it is impossible to avoid hazardous materials transportation, since it is closely related to the industrial and agricultural production and people's daily life [1]. Nowadays, road freight plays a principal role in the hazardous materials transportation. However, the frequent road traffic accidents increase the risks of hazardous materials transportation and thereby considerably affect its sustainable development [4].
As an alternative to the road freight that is a traditional unimodal transportation mode, a multimodal transportation that combines road and rail transportation has been extensively promoted in hazardous materials transportation [1], [5]. Compared with road transportation, rail transportation is more cost-efficient and yields better safety, especially in the long-distance bulk transportation settings. Transferring hazardous materials between road and rail can take advantage of the above characteristics of rail transportation and at the same time keep the good mobility of road transportation in goods pickup and delivery to provide door-to-door services [6]. Consequently, hazardous materials road-rail multimodal transportation is being promoted by both governments and transportation industries globally.
In the transportation planning field, researchers consider that routing optimization provides a quantitative support to effectively manage the transportation activities in unimodal and multimodal transportation systems [7]. It is also found to be significantly effective in reducing risks from hazardous materials transportation [8]. Generally, hazardous materials road-rail multimodal routing problem aims at accomplishing hazardous materials transportation orders according to customers' requirements and refers to combining road and rail services based on their organization characteristics to generate the best routes to transport hazardous materials from origins to destinations. In the routing optimization, with the help of the multi-objective optimization techniques, a most suitable route plan can be identified by balancing the conflicting objectives of improving economy and reducing risk of hazardous materials transportation. In recent years, the study on the hazardous materials multimodal routing becomes a highlight in the transportation planning field.
Currently, the hazardous materials transportation planning can be classified into six categories, which can be seen in Figure 1. As stated by Erkut et al. [9], compared with hazardous materials transportation planning in road transportation system, multimodal transportation is not the dominant area that is studied. However, there are many research articles on the hazardous materials multimodal (location-) routing problem that can be found in recent decades.
In the existing literature, Xie et al. [10] model a locationrouting problem to plan and manage the hazardous materials multimodal transportation, in which location refers to siting the suitable nodes to transfer hazardous materials between road and rail. In this study, the optimization objective weighted combination of costs and risks. By linearizing the model to generate an equivalent linear formulation, the authors use CPLEX optimizer to solve a real-world case. This study provides a solid foundation for future works on this topic. Verma and Verter [11] systematically introduce a Gaussian plume model based method to determine the population exposure along the planned routes, and formulate it as the risk objective for routing hazardous materials shipments in the road-rail multimodal network. Then, population exposure becomes the commonly used indicator measuring the risks of hazardous materials multimodal transportation. Taking the minimization of the population exposure as the risk objective, Verma et al. [12] explore the optimization framework for routing hazardous materials road-rail multimodal shipments with lead time constraint and develop a solution approach based on tabu search to solve the large-scale real-world cases. To extend the study by Verma et al. [12], Assadipour et al. [13] discusses the hazardous materials multimodal routing problem that considers congestion at yards and equipment capacity selection. The same group of authors also presents a bi-level programming approach for the problem that government's toll-setting policy on using terminals [14]. Both the two works design multi-objective heuristic algorithms to solve the real-world cases. Recently, Ke [15] investigates a hazardous materials road-rail multimodal routing problem with random yard disruptions and establishes a scenario-based robust optimization model.
Chinese scholars also contribute a lot to this research topic, especially under the background that the annual hazardous materials transportation volume has been up to about 400 million tons in China [2]. Based on the shortest path problem, Xin et al. [16] concerns a hazardous materials multimodal routing problem under a time-varying environment and design a modified Dijkstra algorithm. Zhao et al. [17] VOLUME 8,2020 formulate the hazardous materials multimodal locationrouting problem. In their work, different ranks are granted to different types of hazardous materials, and the transshipment nodes will first handle these with higher ranks. Cao et al. [18] use the theory of conditional value at risk to estimate the risk of the hazardous materials routing and build a multi-objective routes selection model based on the multimodal network transformation. Huang and Shuai [19] focus on modeling the bi-objective multi-commodity hazardous materials multimodal routing problem, and present multi-objective algorithm integrating path finding algorithm and non-dominated sorting principle to obtain the Pareto solutions. Li et al. [20] carry out the research on the multi-period hazardous materials multimodal path selection and demonstrate that considering the multiple periods reduces the costs of hazardous materials transportation compared with single-period modeling. Sun et al. [21] construct a bi-objective nonlinear model and draw the equivalent linear reformulations for the hazardous materials routing problem in the multimodal network, and use normalized weighing method to get the Pareto solutions on the platform of LINGO optimizer. Using a modeling method similar to [21], a hazardous materials road-rail multimodal routing problem under risk uncertainty and sustainability consideration is addressed by a recent article by Sun et al. [22], in which the transportation system is formulated by a hub-and-spoke network.
Above studies provide helpful thoughts and tools to plan and manage hazardous materials multimodal transportation by routing optimization. However, these researches are interested in the routing decision under deterministic environment, in which the demands of hazardous materials between the origins and the destinations are certain. In fact, as stressed by Sun et al.'s work [23], the hazardous materials routing decision making should be carried out after the customers (i.e., shippers and receivers) propose transportation orders, during which the actual transportation associated with these orders does not start. Because of the rapidly changing markets and less effective communication among shippers, receivers, transportation providers and the multimodal operators [24], the exact demands for all the hazardous materials transportation orders are difficult to be determined, which results in the demand uncertainty during the decision making phase. In the deterministic modeling, the routes are probably infeasible due to capacity constraint if there are underestimations on the deterministic values of hazardous materials demands, while are not the actual optimum if overestimations happen [23]. Since demand uncertainty influences the reliability of routing [25], it is necessary to pay attention to this issue when modeling the hazardous materials routing problem. Therefore, improving the reliability of hazardous materials road-rail multimodal routing optimization through formulating demand uncertainty is a target of this study.
Moreover, the current relevant literature pays less attention to the timeliness of hazardous materials transportation. In some studies [18], [20], the timeliness of transportation is not considered, so the routing is not restricted by the due dates (i.e., deadlines) and it is acceptable that the hazardous materials can arrive at the destinations at any time. Other studies [11], [21] consider that the hazardous materials transportation ought to follow prescribed due dates, in which a hard time constraint is proposed to ensure that deliveries of hazardous materials to their destinations should be accomplished before deadlines. However, this setting can only avoid delayed delivery while may lead to early delivery that possibly causes extra costs for processing and storing hazardous materials. Sun et al. [22] consider the due dates as soft time windows, in which penalty costs will create if the arrivals of hazardous materials at destinations fall out of the windows. However, as claimed by Tang et al. [26], violation of the time windows is not directly related to penalty. Additionally, it is challenging to objectively set the applicable penalty costs per ton per hour, since the value is oriented on decision makers' subjective opinions. Besides the traditional economic and risk objectives, improving the timeliness of hazardous materials transportation is associated with the quality of service of transportation. Under the background that more and more companies demand on-time delivery to implement just-intime production strategy [27], the improvement on the service level of hazardous materials road-rail multimodal routing that shows customers' subjective satisfaction is targeted by this study.
Above all, besides improving the economy and reducing the risks that are the traditional objectives, this study continues to explore the hazardous materials routing problem in with the emerging objectives of optimizing its reliability and timeliness, which is illustrated by Figure 2. In order to achieve these objectives, following works has been carried out by this study.
(1) Demand uncertainty is introduced into the hazardous materials routing problem in the road-rail multimodal network. According to the fuzzy set theory, uncertain demands are represented by triangular fuzzy numbers to describe their impreciseness.
(2) Service level, i.e., the customers' satisfaction on the delivery services, of the hazardous materials transportation is quantified based on the fuzzy soft time windows. Accordingly, the customers' endurable least service level is ensured by building the service level constraint in the hazardous materials routing modeling.
(3) A fuzzy multi-objective linear model is proposed to address the hazardous materials routing problem explored in this study. The fuzzy expected value method and Jimenez's fuzzy ranking method are used to deal with the fuzzy objectives and fuzzy constraint of the proposed model. Then, the Pareto solutions to the problem are generated by using the ε-constraint method.
(4) By using sensitivity analysis and fuzzy simulation, the effects of improving the reliability and service level on the hazardous materials routing optimization are summarized. Some insights that are helpful to deal with the objectives of the hazardous materials routing are drawn. It is the multimodal operator's work that addresses the hazardous materials road-rail multimodal routing problem and organize the practical hazardous materials transportation based on the customer demands, which will be clarified in Section II. For the customers, they would like to use least expenditure to accomplish the hazardous materials transportation, while also seek for improved reliability and timeliness. The public and governments are concerned with the safety of the hazardous materials transportation and risk reduction is their interest. The application of this study is thereby to provide a quantitative tool to help the multimodal operator plan the hazardous materials road-rail multimodal transportation by making effective tradeoffs among the economic, risk, reliability and timeliness objectives concerned by different actors and stakeholders.
The rest sections are presented as follows. Introduction to the modeling of the hazardous materials road-rail multimodal transportation system is systematically presented in Section 2, in which the consolidation network, the transportation process combining two transportation modes, and the risk estimation on the network are given based on the author's previous studies. In Section 3, after presenting the formulation of demand fuzziness by triangular fuzzy numbers and service level based on fuzzy soft time windows, a fuzzy multi-objective routing model is established. Then, Section 4 designs a solution method that combines the fuzzy expected value method, Jimenez's fuzzy ranking method and ε-constraint method to obtain the Pareto solutions to the hazardous materials routing problem. In Section 5, an experimental case study is given to first verify the proposed modeling of this study and then reveal the effects of the demand fuzziness and service level improvement on the multi-objective routing optimization for hazardous materials transportation. Finally, Section 6 summarizes the conclusions of this study.

II. MODELING THE HAZARDOUS MATERIALS TRANSPORTATION SYSTEM
This section systematically introduces the modeling of the hazardous materials road-rail multimodal transportation system, which provides a solid foundation for the problem optimization carried out in this study.

A. MODELING THE ORGANIZATION AND CONSOLIDATION NETWORK OF THE TRANSPORTATION SYSTEM
Modeling the hazardous materials transportation system is the foundation of carrying out the study on the associated routing problem. Various actors and stakeholders in the system take part in the hazardous materials routing that is an operational-level planning, which are listed as follows [28], [29]. Moreover, since the hazardous materials involve with extensive risks, governments should regulate the hazardous materials transportation for the sake of risk control.
(1) Shippers and receivers (i.e., customers), who propose hazardous materials transportation orders that include the origin-to-destination pairs, demands, release times and due dates.
(2) Drayage operators, who undertake the short/mediumdistance pickups (from shippers to multimodal terminals) and deliveries (from multimodal terminals to receivers) by using road services (i.e., trucks).
(3) Terminal operators, who carry out the transshipment of hazardous materials between road and rail according to the schedules of freight trains.
(4) Network operators, who provide freight trains with fixed schedules to move hazardous materials between multimodal terminals. The fixed schedules include the freight trains' routes, operation time windows at multimodal terminals, departure times from the terminals, arrival times at the terminals and operation periods [30]. (5) Multimodal operators, who represent the customers to coordinate drayage operators, terminal operators and network operators to accomplish the hazardous materials transportation.
It is the multimodal operators' concern to select the routes by integrating the services provided by other operators [28]. When planning the routes, the multimodal operator should first of all determine the consolidation network of the transportation system, which is the basis for managing the hazardous materials road-rail multimodal transportation. In both theory and practice, a hub-and-spoke network is most suitable for a multimodal transportation system that consists of road and rail transportation modes [31]. As indicated by Wang et al. [32], it has various advantages compared with the point-to-point network, line network, and collection-anddistribution network. It can improve the customer service level and the multimodal transportation frequency, centralize the operations at the multimodal terminals to achieve economies of scales, and dominate more market share in a particular region.
However, the existing literature does not attach importance to a hub-and-spoke network for hazardous materials transportation. Among the articles reviewed in Section 1, only Sun et al. [22] and Verma et al. [12] consider such a structure in their studies. Therefore, this study focuses on the hazardous materials routing problem in a transportation network illustrated by Figure 3.

B. MODELING THE HAZARDOUS MATERIALS DOOR-TO-DOOR TRANSPORTATION PROCESS
When using the transportation system indicated by Figure 3 to organize the hazardous materials transportation, the drayage operators and network operators provide information on available road and rail services to the multimodal operator, respectively. Based on the given information, the multimodal operator selects suitable transportation services to carry out pickups, long-distance transportation and deliveries of hazardous materials. Since the freight trains follow fixed schedules, the transshipments of hazardous materials conducted by terminal operators should consider the restrictions of freight trains' fixed operation time windows at the multimodal terminals.
When a freight train is selected, the upper bound of its fixed operation time window at the multimodal terminal regulates the loading operation cutoff time, and restricts that the arrival of the hazardous materials by trucks should not be later than such a bound. Therefore, road pickup services whose arrival times at the multimodal terminal are later than the upper bound (e.g., road pickup service 1 in Figure 4) should not be selected by the multimodal operator.
Road pickup services whose arrival times at the multimodal terminal falls into operation time window (e.g., road pickup service 2 in Figure 4) are optional. Under this situation, the transshipment of the hazardous materials from road to rail can be immediately carried out by the terminal operator.
The corresponding lower bound regulates the loading operation start time of the freight train at the multimodal terminal. Road pickup services whose arrival times at the multimodal terminal are earlier than such a bound (e.g., road pickup service 3 in Figure 4) are also applicable. Under this situation, the transshipment of the hazardous materials from road to rail should be wait until the lower bound, and then undertaken by the terminal operator. During the waiting, storage of hazardous materials will be charged by the terminal operator. After the freight train arrives at the multimodal terminal, the hazardous materials cannot be unloaded immediately until the freight train's unloading operation start time. Once such a time is reached, the hazardous materials will get transshipped from rail to road to start deliveries.
Road transportation usually implements a full-truck-load strategy [23]. Therefore, the road pickups and deliveries are not same to their operations in the vehicle routing problem that relates to the less-than-truck-load transportation. In this case, one truck only loads the hazardous materials from one transportation order. There usually needs a fleet of trucks to load all the hazardous materials of one transportation order. During the transportation, the multimodal operator can integer all the trucks from different drayage operators who provide road services on the same line into a group of trucks [22]. Then, based on the assignment of hazardous materials transportation orders according to the routing, the multimodal operator will flexibly organize truck fleets to conduct the road transportation.
Considering the flexibility and good mobility of road transportation, in the road pickup services, the truck fleets carrying the hazardous materials depart from their origins at the release time indicated in the transportation orders. When the hazardous materials arrive at the goods yard of the multimodal terminals by tucks, their unloading can be immediately carried out. As for the road delivery services, transshipments of hazardous materials from rail to road will be undertaken at the lower bounds of the freight trains' operation time windows. After transshipments, the hazardous materials will be distributed to their destinations by truck fleets.

C. MODELING THE RISKS OF THE HAZARDOUS MATERIALS TRANSPORTATION SYSTEM
Shown as Figure 3, the hazardous materials road-rail multimodal transportation mainly includes in-transit process where hazardous materials get moved from one node to anther by one transportation mode and transshipment process where transportation modes of hazardous materials are changed. The two processes are operated on the arcs and at the nodes of the transportation network, respectively. Both the two processes may suffer transportation accidents that cause the accidental release of hazardous materials. The accidental release will affect the people living along the arcs and around the nodes. Consequently, according to Verma and Verter [11], the risks are represented by population exposure that defines the number of people exposed to the possibilities of hazardous materials transportation accidents happened on the planned road-rail multimodal routes.
As stated above, the risks exist on the arcs and at the nodes. In large numbers of the existing literature, minimizing population exposure is used as the risk objective of the hazardous materials transportation planning problems. There are also some studies that consider the multiplication of the population exposure on the arcs and around the nodes with corresponding possibilities of transportation accidents as the risks [33]. However, in order to avoid the situation that slight possibilities of transportation accidents lead to severe consequences, this study directly uses population exposure to represent the risks.
Verma and Verter [11] have proposed an effective method to determine the population exposure of the multimodal routes based on Gaussian plume model. Therefore, this study VOLUME 8, 2020 considers the population exposure due to moving or transshipping one ton of hazardous materials as known parameter.

III. OPTIMIZATION MODEL
This section presents the modeling of the demand fuzziness and service level in the hazardous materials road-rail multimodal transportation setting and proposes a fuzzy multi-objective routing model.

A. MODELING DEMAND FUZZINESS
When modeling the optimization problems under uncertainty, Zheng and Liu [34], Zarandi et al. [35] and many other studies acknowledge that it is difficult to attain enough data to represent the stochasticity of uncertain parameters in most cases. Furthermore, as clarified by Fazayeli et al. [36], stochastic programming uses scenario-based method to model optimization problems under uncertainty, in which uncertainty is formulated by a lot of scenarios. Such a characteristic of stochastic programming makes its computation time consuming. Consequently, as an alternative, this study employs fuzzy programming to model the hazardous materials routing problem under demand uncertainty.
First of all, according to the fuzzy set theory, fuzzy programming adopts fuzzy numbers to represent uncertain parameters. As the most widely used form [36], triangular fuzzy numbers are employed by this study to represent the fuzzy demands, in which there are three prominent points that describe a fuzzy demand. This study takesq = (q 1 , q 2 , q 3 ) shown as Figure 5 to demonstrate a fuzzy demand, in which there is q 3 ≥ q 2 ≥ q 1 . In the fuzzy demand, q 1 is the minimum possible values of the fuzzy demand, while q 3 is the maximum. They are related to the optimistic and pessimistic estimations, respectively. q 2 is the most likely value that appears in most cases. The fuzzy membership degree indicated in Figure 5 shows the possibilities of these values appearing in the practical transportation. Although demand fuzziness has not received enough attention from hazardous materials multimodal routing, limited articles on the similar topic for regular goods that consider such a characteristic can be found [6], [23], [36].

B. ESTABLISHING THE SERVICE LEVEL OF THE HAZARDOUS MATERIALS ROAD-RAIL MULTIMODAL ROUTING
Oriented on the customers' following subjective opinions, this study uses a fuzzy soft time window to represent a due date. Furthermore, the service level of the hazardous materials routing can be established.
(1) The customers would like on-time transportation, which means that the deliveries of hazardous materials to the destinations should be accomplished at a right time that is neither too early nor too late. Consequently, such a preference needs a time window to reflect. The customers are satisfied when the arrivals of the hazardous materials at the destinations fall into the time windows. Under this situation, the service level of the hazardous materials transportation reaches the highest.
(2) The customers accept early or late deliveries to a certain degree, for example, if the costs can be reduced. These situations will worsen the service level. But they cannot tolerate that the deliveries are too early or too late. As a result, they have the endurable earliest and latest times for receiving their orders. It is unacceptable that deliveries are accomplished out of the endurable time ranges.
Considering above descriptions, this study uses 0 to 1.0 to measure the service level. The service level in the first situation is up 1.0. For the second, the service level will decrease from 1.0 to 0 with the violation of the preferred time window aggravates. Therefore, the service level of hazardous materials routing can be expressed by Figure 6.
The general function for the service level of hazardous materials routing is illustrated by functions (1a) to (1c) shown in Figure 6. Its expression thereby is a piecewise function shown as Equation (1) [37] where [ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 ] is a fuzzy soft time window and y is the time when the hazardous materials arrive at the destination by using the planned route, i.e., the time when the delivery is accomplished.
In Equation (1), α and β are the equation constants. In order to avoid the nonlinearity of the modeling, both α and β are set to 1. Consequently, the service level is as represented by functions (2a) to (2c) in Figure 6, which is same to the fuzzy membership function of trapezoidal fuzzy numbers.
During the routing decision making, the customers should first of all propose the least service level that they can accept. The service level is determined by the time when the delivery of hazardous materials is accomplished. According to the planned routes, the service level should not be lower than the prescribed lower bound. Otherwise, the planned routes are infeasible. According to Figure 6, functions (2a) and (2c) decide the lower and upper bounds of the acceptable arrival times of hazardous materials at destinations. Supposing the least service level is represented by γ , there are y−ϕ 1

C. NOTATIONS
This study classifies the symbols used for the modeling into four categories. The symbols representing the sets, indices, parameters and decision variables are defined and presented in Table 1.

D. A FUZZY MULTI-OBJECTIVE ROUTING MODEL
In order to present a rigorous problem modeling, some general assumptions that are associated with the multimodal routing should be followed and are listed as below [23]: (1) The hazardous materials of each transportation order depart from the corresponding origin at the prescribed release time.
(2) Each hazardous materials transportation order is considered to be accomplished when the associated hazardous materials arrive at the corresponding destination.
(3) The hazardous materials of each transportation order cannot be splitted during the entire transportation process, so that each hazardous materials transportation order is integrated.
(4) The hazardous materials road-rail multimodal transportation system is under deterministic environment, which means that all the parameters regarding the road-rail multimodal network are known and deterministic.
Based on above assumptions and using the symbols defined in Section 3.3, this section presents a fuzzy multi-objective routing model whose framework is oriented on the author's previous works [21], [30]. The formulation of the model presented below is linear, which means that it is suitable be solved by exact solution algorithms after its defuzzification.
Economic objective: Equation (2), as shown at the bottom of the next page, is the economic objective that minimize the total costs charged by the multimodal operator for all the hazardous materials transportation orders. The total costs include travel costs, handling costs and storage costs that successively correspond to the three formulas in Equation (2).
Risk objective: Equation (3), as shown at the bottom of the next page, is the risk objective for reducing the total population exposure due to the hazardous materials door-to-door transportation by using the planned routes. In Equation (3), the first formula is the cumulative population exposure at the nodes, while the second one refers to the cumulative population exposure along the arcs. Equations (2) and (3) are fuzzy objectives, since they contain fuzzy parameterq k in their objective functions.
Constraint Set: Equation (4) is the hazardous materials flow conservation constraint. For each hazardous materials transportation order, ''-1'' means that there is only outbound flow at the origin, while ''1'' shows that there is only inbound flow at the destination. ''0'' explains that both inbound and outbound flows whose volumes are same exist at the multimodal terminals that serve as transshipment nodes. By the above setting, the routes continuity can be ensured.
Equation (5) ensures that the hazardous materials of one transportation order should not be splitted during the transportation process.
The combination of Equations (4) and (5) keep the space continuity of the planned hazardous materials road-rail multimodal routes consisting of nodes and arcs, as well as the integrity of hazardous materials transportation orders.
Equation (6) ensures that the hazardous materials of each transportation order depart from the origin at the release time during the process of road pickup service.
Equations (7) and (8) calculate the arrival times of hazardous materials of each transportation order at the nodes on the route when using road services. When a transportation service on a directed arc is selected for a hazardous materials transportation order, the above two equations equal to y k i + t ijs − y k j = 0. The arrival times at the nodes can be then determined.
Equations (9) and (10) have similar structures to Equations (7) and (8). They calculate the arrival times of hazardous materials of each transportation order at the nodes on the routes when using rail services. After arriving at the nodes at the scheduled arrive times of the freight train, the hazardous materials cannot get unloaded until the scheduled unloading operation start time. As a result, it is more suitable to formulate such lower bounds as the arrival times of the hazardous materials at the nodes by freight trains.
Equations (11) and (12) calculate the waiting periods of hazardous materials of each transportation order at the multimodal terminals before loaded on freight trains. According to the above equations, when a freight train is selected for moving a certain transportation order, there is z k ijs ≥ l s i − y k i . If y k i is smaller than l s i (e.g., road pickup service 3 in Figure 4), there is z k ijs ≥ l s i − y k i > 0. Minimization of Equation (2) results in that z k ijs = l s i − y k i . However, if y k i is larger than l s i (e.g., road pickup service 2 in Figure 4), z k ijs ≥ l s i − y k i < 0. Considering the domain constraint that z k ijs ≥ 0, minimization of Equation (2) results in that z k ijs = 0. (13) Equation (13) is the operation time window constraint to ensure the feasible transshipment of hazardous materials from road to rail, i.e., the arrival time of hazardous materials at the multimodal terminal by road should be later than the scheduled loading operation cutoff time of the selected freight train.
Equations (6) to (13) ensure the time validity of the planned hazardous materials road-rail multimodal routes, so that the transshipments of hazardous materials between road and rail at the multimodal terminals on the routes can be realized smoothly according to the schedules of rail services. As a result, the coordination of the two transportation modes can be enhanced.
Equations (14) and (15) are the service level constraints. They ensure that the service level of hazardous materials transportation satisfy the customers' endurable least level.
Equation (16) is the capacity constraint requiring that the total demands of hazardous materials loaded on a transportation service should not exceed its available capacity. As for a rail service, its capacity is the available payload of a freight train, while the capacity of a road service is associated with the group of trucks on the road transportation line. Equation (16)'s left formula contains fuzzy parameterq k . Therefore, it is a fuzzy inequality constraint.
Equations (17) to (19) are the domain constrains of decision variables that regulate their values.

IV. SOLUTION APPROACH
As already proved by Ayar and Yaman [38], the multimodal routing problem is a NP-hard problem. It is widely acknowledged that heuristic algorithms are generally more suitable to solve such a kind of problems than exact solution algorithms, especially in the large-scale case setting. However, as stated by Sun et al. [32], the design and application of an exact solution algorithm to solve the multimodal routing problem has two advantages.
First of all, the exact solution algorithms strictly follow the structure of the mathematical model to obtain solutions and can verify if the model itself is logically sound or not, while heuristic algorithm depends less on the model and more on the natural language describing the model. Second, it is difficult to make sure that the solutions generated by heuristic algorithms are global or local optimum. However, the exact solution algorithm can provide an exact benchmark for the problem by giving their global optimal solutions.
In the existing literature, there are various studies that employ exact solutions algorithms to solve the multimodal routing problem and demonstrate their feasibility in dealing with empirical cases, e.g., Xie et al. [10] and Sun [6], Sun et al. [21], Sun and Lang [30], [32]. As a result, this study continues to seek for the application of exact solution algorithms in solving the hazardous materials road-rail multimodal routing problem, in which mathematical programming software can be adopted to run the algorithms efficiently. VOLUME 8, 2020 Under this situation, this study needs first to convert the initial model proposed in Section III into a crisp single-objective linear programming model.

A. DEFUZZIFICATION OF THE INITIAL FUZZY MODEL
The optimization model proposed in Section 3 contain imprecise information. In order to obtain crisp optimization results to support the practical hazardous materials transportation organization, defuzzification should be undertaken to make the model solvable. The widely used fuzzy expected value method is adopted by this study to address the fuzzy economic and risk objectives. As for fuzzy demandq k = q 1 k , q 2 k , q 3 k that is a triangular fuzzy number, the expected value is as Equation (20).
By using this method, the economic objective, i.e., Equation (2) can be converted into Equation (21), as shown at the bottom of the next page, that minimizes the expected costs for accomplishing all the hazardous materials transportation orders.
The linearity of fuzzy expected value operator ensures that is a fuzzy parameter represented by triangular fuzzy number and x i is a variable [39]. As a result, the crisp economic objective is as Equation (22), as shown at the bottom of the next page. Similarly, the fuzzy risk objective, i.e., Equation (2), can be rewritten as Equation (23), as shown at the bottom of the next page, that minimizes the fuzzy expected value of all the population exposure.
By using this method, defuzzification of fuzzy inequality constraints can be also realized. As demonstrated by Sun [6], using this method reduces the flexibility of the fuzzy optimization. Therefore, Jimenez's fuzzy ranking method [40] is employed by this study as an alternative to carry out the defuzzification of Equation (16).
In the computing, this method preserves the linearity of the constraints and does not increase the number of constraints [40]. It is thereby simpler and yields better efficiency in addressing the fuzzy inequality constraints in comparison to the widely used fuzzy chance-constrained programming. Such a method has been adopted in the transportation, logistics, and supply chain planning problems under a fuzzy environment [4], [36], [41].
Consider a fuzzy constraintã · x ≥b wherẽ b = (b 1 , b 2 , b 3 ) is another triangular fuzzy parameter. According to Jimenez [40], such a fuzzy constraint can be converted into where λ ∈ [0, 1] is the confidence level that defines the least degree thatã · x is not smaller thanb. Accordingly, Equation (16) can be rewritten as Equation (24).

B. THE E -CONSTRAINT METHOD TO GENERATE THE PARETO FRONTIER
Hazardous materials routing problem is a kind of multi-objective optimization, since it has two objectives that cannot be simultaneously optimized by a single optimal solution [42]. In a multi-objective optimization, an objective can be only improved with worsening at least one other objective [43], [44]. Consequently, there are non-dominated Pareto solutions in multi-objective optimizations, instead of dominated solution. All the Pareto solutions form the Pareto frontier. A suitable Pareto solution can be determined as the transportation plan based on the decision makers' preferences [45].
The weighting method (also known as weighed sum method) and the ε-constraint method are both adopted by the existing literature on the multi-objective hazardous materials routing problems, in which the utilization of the first method is wider [46]. The weighting method is a priori method, while the ε-constraint method is posterior. Although the weighting method is easy to be conducted in the computing, it has many disadvantages. For example: (1) In the decision making practice, as a kind of priori method, the weighting method requires the predetermination of the weights to different objectives before solving the multi-objective optimization problems. However, during practical decision making, the relative importance of different objectives may be unclear, which reduces the effectiveness of this method. Furthermore, this method is not effective in the posteriori decision making [47]. However, the ε-constraint method does not need pre-determined weights.
(2) In the computation, the weighting method often generates extreme Pareto solutions, especially when the feasible solution region is non-convex [42], [47]. However, the ε-constraint method can result in evenly distributed Pareto solutions. Furthermore, the weighting method may need to try a lot of weight combinations to obtain the Pareto solutions, and sometime different weights lead to same solutions, which makes the weighting method time-consuming. However, the ε-constraint method yields better efficiency.
A detailed introduction to the disadvantages can be seen in Mavrotas's research [42]. As a result, this study uses the ε-constraint method, a posteriori method, to get the evenly distributed Pareto solutions to the specific hazardous materials routing problem in this study. The hazardous materials routing model contains minimization objectives. Therefore, the ε-constraint method optimizes one objective with converting the other objectives into inequality constraints that their objective function values should exceed respective upper bounds. By relaxing the upper bounds with fixed steps in acceptable ranges and solving the single-objective optimization model iteratively, the Pareto solutions can be obtained [43], [48].
This study considers the economic objective denoted whose function is denoted by G 1 as the priori objective, and let the risk objective function denoted by G 2 be the constraint. By solving the single-objective model whose objectives are Equation (22) and constraint set is Equations (4) to (15), (17) to (19), and (24), G * 1 , the optimal value of G 1 , and G 2 , the corresponding value of G 2 , can be obtained. Similarly, by solving the model sharing the same constraint set and taking Equation (23) as the objective, the optimal value G * 2 can be generated.
The decision makers will first set the number of Pareto solutions that they want to get. Let M denote the number of solutions. To generate Pareto solution m where m = 1, 2, . . . , M , the constraint converted from the risk objective is as Equation (25).
By solving the single-objective liner programming model whose objectives are Equation (22) and constraint set is Equations (4) to (15), (17) to (19), (24), and (25) for m = 1, 2, . . . , M , the Pareto solutions can be obtained. Mathematical programming software is a useful tool to automatically run exact solution algorithm to solve such a kind of linear model by outputting its global optimal solutions.

V. EXPERIMENTAL CASE STUDY
This section gives an experimental case study to first demonstrate the feasibility of the proposed method, and then draw some insights based on sensitivity analysis and fuzzy simulation to help to better manage the hazardous materials road-rail multimodal transportation.

A. CASE DESCRIPTION
An experimental case designed in Sun et al.'s study [22] is introduced to demonstrate the feasibility of this study. The hazardous materials road-rail multimodal network is illustrated by Figure 7, in which the numbers beside the nodes and arcs are the corresponding population exposure in thousand people.
The information of the road and rail services in the network can be seen in the appendix in Sun et al.'s work [22], and is also presented in Tables S1 and S2 in the supplementary materials of this study. In China, the cost rates on the hazardous materials transportation that are the main controlling parameters of hazardous materials routing problem are regulated by the China State Railway Group Company (the former Ministry of Railways of China) and the Ministry of Transportation of China [6]. Although this study focuses on an experimental case, the values of these controlling parameters are still determined based on the practical transportation situation in order to make the optimization results on the experimental case and resulting conclusions feasible in practice. This study    sets the cost rates as Table 2 by referring to an empirical case presented in Sun et al.'s study [21]. The hazardous materials transportation orders are given in Table 3. In order to improve the computational efficiency, the release times and fuzzy soft time windows are transformed into real numbers for the convenience of computation, e.g., 12:00 in the second day of the planning horizon is equal to 36 (12+24=36).

B. OPTIMIZATION RESULTS
This study uses LINGO optimizer, a standard mathematical programming software, to run the branch-and-bound algorithm to solve the problem. First of all, this study assumes that: (1) The customers attach considerable importance to the hazardous materials routing and would not like to bear any chance that the capacity constraint is violated during the actual transportation that uses the planned routes. As a result, λ in Equation (24) is set to 1.0.
(2) The customers require on-time transportation and accept to pay more charges to have their orders accomplished within their preferred time windows. In this case, γ k in Equations (14) and (15) is also set to 1.0.
The scale of the experimental case shown in Figure 7 is shown in Table 4.
This study runs the branch-and-bound algorithm in the LINGO optimizer on a ThinkPad Laptop with Intel Core i5-5200U, 2.20 GHz CPU, and 8 GB RAM. The computational time consumed to generate each Pareto solution is given in Table 5.
As we can see from Table 5, the proposed solution method can solve the problem by generating its each Pareto solutions within 3.0 seconds and hence shows good computational efficiency in solving the experimental case discussed in this study. According to customers' above requirements, the multimodal operator can use the proposed methods presented in Section IV to generate the Pareto frontier shown in Figure 8 for the experimental case.  According to Figure 8, the economic and risk objectives of the hazardous materials road-rail multimodal routing problem are conflicting. As a result, it is impossible to find an optimal solution that can make the two objectives simultaneously reach their optimum. Improving the economic objective will worsen the risk objective, and vice versa. Therefore, during the hazardous materials routing decision making, the multimodal operator should make a tradeoff between the two objectives.
The Pareto frontier in Figure 8 provides the multimodal operator with diverse candidate plans. In practical decision making, when the multimodal operator considers the risk VOLUME 8, 2020 objective lower than a certain value is acceptable, the Pareto solution yielding the minimum costs and satisfying the preferred condition is the best choice. For example, if it is accepted by the multimodal operator that the risk objective should be lower than 341 (10 million people), the 7 th Pareto solution (41.2, 340.3) is the most suitable plan for organizing the hazardous materials transportation.
In this section, this study also employs the widely used weighting method to generate the Pareto solutions to the hazardous materials road-rail multimodal routing problem under the setting that λ = 1.0 and γ k = 1.0. After normalizing the economic objective (i.e., Equations (22)) and risk objective (i.e., Equation (23)) by using their respective optimal values, this method distributes different weights to normalized objectives and then combines them into Equation (26), where θ ∈ [0, 1.0] is the weight distributed to the normalized economic objective and its value is determined by decision makers before optimization. When the value of θ is set, the corresponding Pareto solution to the problem can be obtained by solving the single-objective linear programming model whose objective is Equation (26) and constraints include Equations (4) to (15), (17) to (19) and (24). This study changes θ from 0.01 to 1.0 with a step of 0.01. By solving the above model under each value of θ, the Pareto solutions to the problem can be obtained and are also shown in Figure 8 where the weights distributed to the two objectives are indicated in the brackets.
As we can see from Figure 8, although the weighting method is widely used by the literature, it misses some Pareto solutions that can be identified by the ε-constraint method. Although these missing solutions may be found by further decreasing the step of changing θ, it will consume much more time to finish the entire computation. Moreover, indicated by LINGO optimizer, the computational time consumed by solving the above model to generate a Pareto solution falls into [5.0, 13.0] seconds. Therefore, the weighting method yields lower computational efficiency compared with the ε-constraint method. As a result, the ε-constraint method is more suitable to be applied to solve the multi-objective optimization for the hazardous materials road-rail, since it can efficiently output evenly distributed Pareto solutions, which matches its advantages explained in Section IV.

C. SENSITIVITY ANALYSIS ON THE HAZARDOUS MATERIALS ROUTING IN REGARD TO THE SERVICE LEVEL
Generally, service level γ i influences the hazardous materials routing optimization results. By using sensitivity analysis, this study discusses the effects of improving the service level on the Pareto frontier for the problem. In this section, setting λ as 1.0, this study presents the Pareto frontiers under service levels of 0.6, 0.8 and 1.0. The optimization results are illustrated by Figure 9.
As we can see from Figure 9, the diversity of the Pareto solutions is sensitive to the service levels. By setting a lower service level, the difference between the two extreme solutions are more significant, and the solution space covered by the Pareto frontier are larger. Consequently, the multimodal operator can obtain more diverse plans. Moreover, from the horizontal viewpoint, with the service level increasing, the Pareto frontier for the problem tends to move to the upper right. However, the tendency is not significant.
The reason that leads to such a sensitivity is as follows. With the improvement of the service level, Equation (14) will restrict a higher lower bound for the arrival times of hazardous materials at the destinations that are decision variables, while Equation (15) will reduce the upper bound of these variables. As a result, the feasible solution space of the hazardous materials routing problem becomes smaller, which is clearly illustrated by Figure 9. With the service level improving, some feasible solutions with smaller economic or risk objectives will become infeasible due to the fact that the corresponding arrival times at destinations cannot satisfy the service level constraints. As a result, these solutions will be excluded during the optimization process. Indicated by Figure 9, it is not wise for customers to look for extremely high service level. Take the Pareto frontiers under service levels of 0.8 and 1.0 for demonstration. Compared with the service level of 1.0, when the service level is changed to 0.8, Pareto solutions with fewer costs and risks can be provided to the multimodal operator and customers. Meanwhile, a relative high service level, i.e., 0.8, can be still guaranteed. Consequently, the multimodal operator should reveal this insight to customers and assist them to determine a suitable service level through sensitivity analysis.

D. SENSITIVITY ANALYSIS ON THE HAZARDOUS MATERIALS ROUTING IN REGARD TO THE CONFIDENCE LEVEL
The demand fuzziness also influences the routing optimization results through confidence level λ associated with the capacity constraint. λ is also determined by the multimodal operator and customers subjectively. In this section, assuming γ k = 0.8, this study gives the Pareto frontier for the problem under confidence levels of 0.5, 0.6, 0.8 and 1.0. The sensitivity in this regard is shown in Figure 10.
As shown in Figure 10, the sensitivity of the hazardous materials road-rail multimodal routing optimization in regard to the confidence level is significant, especially when the confidence level changes from 0.6 to 1.0. The Pareto frontier for the problem moves to the upper right with the confidence level increasing. In theory, higher value of the confidence level represents better reliability. Therefore, the economic, risk and reliability of the hazardous materials routing problem explored in this study are conflicting with each other. Improving the reliability worsens both the economic and risk objectives of the routing. Therefore, the multimodal operator and customers should work together to determine a suitable confidence level.
The reason that leads to the sensitivity shown as Figure 10 comes from the restriction of Equation (24). As we can see From Equation (24), with the confidence level improves from 0 to 1.0, (1 − λ) · q 1 k +q 2 k 2 + λ · q 2 k +q 3 k 2 , the parameter multiplied with x k ijs that is a decision variable, increases from Considering a fixed upper bound, i.e., cap ijs , increasing the value of such a parameter will reduce the feasible solution space of decision variable x k ijs . Consequently, some solutions that used to be feasible under smaller values of confidence level become infeasible. Therefore, the Pareto frontier for the hazardous materials routing problem moves to the upper right when the confidence level.
In existing literature on hazardous materials multimodal routing problem, the demands are modelled as deterministic parameters. In the deterministic modeling, demands are determined by their most likely values [49]. Therefore, the deterministic version of the hazardous materials routing model yields objectives as Equations (27) and (28), as shown at the bottom of the next page.
The corresponding deterministic capacity constraint is as Equation (29).
Consequently, the deterministic version of the hazardous materials routing model contains Equations (27) and (28) as its objectives and Equations (4) to (15), (17) to (19) and (29) as the constraint. By using the ε-constraint method and keeping γ k as 0.8, the Pareto solutions to the hazardous materials road-rail multimodal routing problem under demand certainty can be also generated and are also shown in Figure 10. As we can see from Figure 10, the Pareto frontier of the deterministic modeling is nearest to the one given by the fuzzy routing optimization under confidence level of 0.5 and is on its left. According to the sensitivity analysis in this section, the deterministic routing optimization has the reliability that is smaller than or equal to that of the fuzzy routing optimization under confidence level of 0.5.
This study uses the upper left and lower right Pareto solutions to measure the gaps between two routing optimization results. Compared with the routing optimization under confidence level of 1.0, the deterministic routing optimization only slightly reduces the economic objective values of the two Pareto solutions by 1.28% and 0.92%, and cuts down on the risk objective values by 1.69% and 1.17%. However, the reliability of the routing optimization is significantly reduced when employing the deterministic modeling. Consequently, from a comprehensive tradeoff among the economy, risk and reliability objectives, the fuzzy routing optimization is more feasible than the deterministic version, since it can remarkably improve the reliability with a slight sacrifice on the economy and risk objectives.

E. TESTING THE RELIABILITY OF THE HAZARDOUS MATERIALS ROUTING OPTIMIZATION RESULTS IN THE SIMULATED ACTUAL TRANSPORTATION
In the last section, the conclusion that fuzzy routing optimization for the hazardous materials road-rail multimodal transportation has better reliability than the deterministic optimization is only based on the theory. The performances of these optimization results in the actual transportation is difficult to determine, since the optimization should be carried out before the actual transportation where the demands of hazardous materials are known exactly.
In this case, using the simulation technique is the only way to represent the actual transportation. In the studies given by the author [6] and Cao and Lai [50], a fuzzy simulation based method is developed to simulate the actual transportation by randomly generating the deterministic values of fuzzy parameters based on the significance of their fuzzy membership degrees. The fuzzy simulation for the actual transportation can be implemented several times, and the optimization results given by the fuzzy programming models can be tested (28) in large numbers of simulated actual transportation cases, so that their performances can be evaluated more objectively and convincingly. As for the hazardous materials road-rail multimodal routing problem under demand uncertainty, the actual demands in the actual transportation can be simulated according to the above method. The fuzzy membership function of a fuzzy demandq k = q 1 k , q 2 k ,q 3 k is as Equation (30).
In the w th simulation (where w = 1, 2, . . . , W ), the simulated actual demand of hazardous materials transportation order k is q w k . The fuzzy membership degree of q w k in regard toq k is µq k q w k . Then, the values of an indicator that falls into [0, 1] and measures the significance of q w k is randomly generated. If µq k q w k is bigger than the indicator's value, q w k can be considered as the actual demand of hazardous materials transportation order k. Otherwise, the random generation of q w k should be implemented again until its significance is satisfied. After undertaking above steps for k ∈ K in all W simulations, W actual transportation cases can be obtained. Then, this study can test if the optimization results given by the fuzzy routing optimization satisfy Equation (31).
A routing optimization result can be considered as a successful plan in actual transportation case w if it satisfies Equation (31). Otherwise, it fails in this case. This study tests if a routing optimization result is successful in all the simulated actual transportation cases, and obtains the success ratio that can reflect the reliability of transportation in a more specific way.
In this study, the fuzzy simulation is carried out 20 times to generate 20 actual transportation cases. The simulated actual demands are given in Table S3 in the supplementary material. For demonstration, this study sets the service level (γ k ) as 0.8, and uses the optimization results given by the deterministic modeling indicated at the end of last section and the fuzzy routing model under confidence levels of 0.5, 0.6, 0.8 and 1.0 for comparisons. Since the optimization results contain various Pareto solutions, this study only considers the one that yields minimum Euclidean distance, i.e., Equation (32), to the ideal point G * 1 , G * 2 after objective normalization. Such a solution makes a balanced tradeoff between the economic and risk objectives of the hazardous materials road-rail multimodal routing problem.
Accordingly, the selected Pareto solutions of deterministic modeling and the fuzzy routing model under confidence level of 0.5, 0.6, 0.8 and 1.0 are presented in Table 6 where the units of the economic and risk objective values are 10 thousand CNY and 10 million people. As we can see from Table 6, the values of both economic and risk objectives of the balanced solution increase with the confidence level improves, and the deterministic model yields a balanced solution that has minimum objective values. Such a variation matches the sensitivity analysis illustrated by Figure 10. Then, this study tests these solutions in the 20 simulated actual transportation cases, and gets their success ratios regarding the deterministic capacity constraint shown as Equation (32). The test result is shown in Figure 11. As we can see from Figure 11, the success ratio of the deterministic modeling is 55%, and is same to that of the fuzzy routing model under confidence level of 0.5. However, when the confidence level changes from 0.5 to 0.6, the success ratio significantly improves from 55% to 90% by 63.6% with considerably slight increases of the economic and risk objectives. This improvement indicates that the consideration of fuzzy demands enhances the reliability of the hazardous materials routing optimization, which further verifies the conclusion drawn at the end of last section.
In the practical decision making, assume that the multimodal operator and customers pay equal attention to the economic and risk objectives and accept a plan whose success ratio should be larger than 80%. Under this situation, the balanced solution given by the fuzzy routing optimization under confidence level of 0.6 will be selected as the hazardous materials road-rail multimodal transportation plan, since it can satisfy the decision makers' requirement on transportation reliability and meanwhile yield the minimum economic and risk objective values compared with other optimization results.
Above all, based on the insights drawn from the sensitivity analysis and fuzzy simulation, there need a negotiation between the multimodal operator and the customers to determine their preferred service level and confidence level in a reasonable way. Then, Pareto optimization can be carried out to generate the Pareto solutions to serve as hazardous materials transportation plans. In practical decision making, multiple group decision making [51], [52] is an effective method to help the multimodal operator and the customers to accurately select a Pareto solution that most fit a given decision making situation.

VI. CONCLUSION
This study investigates a freight routing problem for hazardous materials transportation in a road-rail multimodal huband-spoke network. Besides paying attention to the economic and risk objectives, this study also considers to improve the reliability and timeliness (service level) of the hazardous materials routing. Following contributions are made by this study to enhance the systematics of the goals of the routing.
(1) Demand uncertainty that influences the optimization reliability is introduced into the hazardous materials routing, in which triangular fuzzy numbers are used to represent the uncertain demands of hazardous materials and fuzzy programming is employed to formulate hazardous materials routing model.
(2) Service level of hazardous materials routing from the perspective of on-time transportation is modeled based on the fuzzy soft time windows, and the service level constraint that ensures the customers' endurable least level to be satisfied is accordingly constructed.
(3) A hazardous materials road-rail multimodal transportation system is systematically established, in which the associated actors and stakeholders, consolidation network, transportation modes and process, and risks are presented in detail.
(4) With the help of a fuzzy multi-objective routing model, the effects of enhancing the reliability (by considering the demand uncertainty) and improving timeliness (by formulating the service level improvement) on the hazardous materials routing optimization are summarized based on the sensitivity analysis and fuzzy simulation.
The insights based on the experimental case study are summarized as follows.
(1) The economic and risk objectives of the hazardous materials road-rail multimodal routing problem are conflicting with each other and cannot reach their own optimum simultaneously. The Pareto solutions to the problem can provide various candidates, and decision makers can select a suitable one according to the specific decision making situation and their preferences.
(2) Improving the service level will increase both the costs and risk of the hazardous materials road-rail multimodal transportation. Therefore, it is not very wise for the customers to seek for excessive service levels. The multimodal operator and the customers should work together to determine a reasonable service level, in which the sensitivity analysis shown as Figure 10 is very helpful.
(3) It is necessary to formulate the demand uncertainty in the hazardous materials road-rail multimodal routing problem, since such a formulation can remarkably enhance the reliability of hazardous materials transportation. Fuzzy programming provides a useful tool to deal with the hazardous materials routing optimization with uncertain demands.
(4) Similar to the effect of service level improvement on the hazardous materials routing optimization, enhancing the reliability also leads to the increase of both the costs and risk of the hazardous materials road-rail multimodal transportation. Meanwhile, through fuzzy simulation to simulate the actual transportation and test the performance of the routing optimization under different confidence levels, this study finds that improving the confidence level does not lead to a constant improvement on the transportation reliability, which is illustrated by Figure 11.
(5) With the help of the performance analysis in the fuzzy simulation illustrated in Figure 11, a suitable confidence level however can be determined by making tradeoffs among the economy, risk and reliability objectives according to the preferences of decision makers.
However, there are still limitations of this study. First of all, this study only considers the uncertainty from the side of customers and does not take the uncertainty of the transportation system into account. Consequently, the future work will model an uncertain environment that covers transportation system and its users for the hazardous materials road-rail multimodal routing problem. Furthermore, the hazardous materials road-rail multimodal routing problem is NP hard. Although the proposed solution method based on exact solution algorithm can provide an exact benchmark for solving the problem, its performance in dealing with the large-scale real-world case is doubtful and should be tested in the future work. If the computational efficiency is not satisfying, design on intelligent algorithms (e.g., heuristic algorithm) should be a direction of the future work. Moreover, it is also interesting to use other forms of fuzzy numbers (e.g., Type-2 fuzzy numbers) to represent the fuzziness and discuss the effects of such a setting on the routing optimization results.