Bi-Objective Modelling for Hazardous Materials Road–Rail Multimodal Routing Problem with Railway Schedule-Based Space–Time Constraints

The transportation of hazardous materials is always accompanied by considerable risk that will impact public and environment security. As an efficient and reliable transportation organization, a multimodal service should participate in the transportation of hazardous materials. In this study, we focus on transporting hazardous materials through the multimodal service network and explore the hazardous materials multimodal routing problem from the operational level of network planning. To formulate this problem more practicably, minimizing the total generalized costs of transporting the hazardous materials and the social risk along the planned routes are set as the optimization objectives. Meanwhile, the following formulation characteristics will be comprehensively modelled: (1) specific customer demands; (2) multiple hazardous material flows; (3) capacitated schedule-based rail service and uncapacitated time-flexible road service; and (4) environmental risk constraint. A bi-objective mixed integer nonlinear programming model is first built to formulate the routing problem that combines the formulation characteristics above. Then linear reformations are developed to linearize and improve the initial model so that it can be effectively solved by exact solution algorithms on standard mathematical programming software. By utilizing the normalized weighted sum method, we can generate the Pareto solutions to the bi-objective optimization problem for a specific case. Finally, a large-scale empirical case study from the Beijing–Tianjin–Hebei Region in China is presented to demonstrate the feasibility of the proposed methods in dealing with the practical problem. Various scenarios are also discussed in the case study.


Introduction
Transportation of hazardous materials has always been a highlighted problem in both the public security and transportation planning fields. Hazardous materials have five characteristics: corrosivity, toxicity, ignitability, reactivity, and infectivity [1]. Therefore, unlike regular goods or materials transportation, the transportation of hazardous materials often involves many risk factors [2]. Once accidents happen in the transportation of hazardous materials, severe consequences, such as loss of life, public damage (e.g., environment pollution, infrastructure unavailability, and economic losses) [3], and social instability, will emerge. In China, hazardous materials are plentiful. The freight volume of hazardous materials reaches more than 4 billion tons per year [4]. However, accidents happen frequently, e.g., transportation accidents involving chemicals account for 30%-40% of total the problem exactly and how to test the performance of the heuristic algorithm are not indicated in their studies.
In this study, in order to further develop the road-rail multimodal routing problem of hazardous materials, we discuss it from the operational level of network planning, and comprehensively consider the following formulation characteristics.
In particular, the restrictions of railway schedules on the routing decision are formulated as the railway schedule-based space-time constraints. (4) Environmental risk constraint to lower and balance the environmental risk under a given threshold. (5) Bi-objective optimization, including minimizing the total generalized costs of transporting the multiple hazardous materials and the social risk along the planned routes.
The remaining sections of this study are organized as follows. In Section 2, the specific transportation scenario for the hazardous materials and some background information are introduced and described, and the railway schedule-based space-time constraints are presented and explained. In Section 3, we evaluate and formulate the transportation risk of hazardous materials from the viewpoints of society and the environment. In Section 4, we build a bi-objective mixed integer nonlinear programming to model the road-rail multimodal routing problem of hazardous materials that comprehensively considers the above five characteristics. In Section 5, linear reformulations are developed to gain the equivalent linear programming of the initial model so that it can be effectively solved by exact solution algorithms on standard mathematical programming software. The normalized weighted sum method is also introduced in this section to address the bi-objective optimization. In Section 6, an empirical case study from the Beijing-Tianjin-Hebei Region in China is designed to demonstrate the feasibility of the proposed methods in dealing with the practical problem. Some discussions on various scenarios are also presented in this section. Finally, the conclusions of this study are drawn in Section 7.

Transportation Scenario Description
The road-rail multimodal service network is composed of three kinds of nodes, including origins, destinations, and hazardous materials stations. The hazardous materials stations are also the transshipment nodes, where transshipments between rail services and road services and between different rail services are conducted. We focus on combining the rail services and road services to generate the optimal routes for transporting the multiple hazardous materials flows. The two kinds of services adopt different operation modes in the transportation practice. The operation of a rail service is based on its prescribed schedule, while a road service is a kind of flexible service [26].
For a certain rail service, its schedule formulates its operation from both time and space viewpoints. These include its route, loading/unloading operation time windows at the goods yards, and classification/disassembly operation windows at the marshalling yards of the stations on the scheduled route. Arrival times and departure times at and from these stations are also included. The two kinds of operation time windows are intervals from operation start times to corresponding cutoff times. In addition, due to the limitations of the effective length of the rail tracks and of the locomotives' tractions, rail services are a kind of capacitated service. On the contrary, the organization of the road service is relatively flexible and is not restricted in terms of time and space. It is also easy to assign enough road vehicles to carry the hazardous materials. Hence, we formulate road services as a kind of uncapacitated service. Additionally, in this study, we neglect the operation times of loading/unloading hazardous materials on/from the trains and road vehicles as well as the operation times of classifying/disassembling rail wagons carrying the hazardous materials into/from the trains. That is, the hazardous materials will complete the unloading operation once they arrive at a node by road service, or complete unloading/disassembly operation at the corresponding unloading/disassembly start time at the node when arriving by the rail service.
As described above, the prescribed schedules of rail services restrict their operations from the space-time viewpoint. Cases 1-3 describe the railway schedule-based space-time constraints, and should be comprehensively modelled in hazardous materials routing. Case 1. When adopting a certain rail service to transport hazardous materials from the current node to the successor node, the arrival time of the hazardous materials at the current node should not be later than the operation (loading/classification) cutoff time of the rail service at the same node. If the arrival time is earlier than the operation (loading/classification) start time, it should wait until that time at the current node, and then get processed further. This case has two sub-cases as follows: Case 1.1. If the hazardous materials arrive at the current node by road service, they should be loaded on the rail service at the goods yard current node. Consequently, in Case 1, the operation start time and operation cutoff time separately refer to the loading start time and loading cutoff time of the rail service. In this case, the selection of a rail service is constrained by its loading operation time window at the current node, and the inventory costs and loading/unloading costs will be charged. Case 1.2. If the hazardous materials arrive at the current node by another rail service, carried by rail wagons, they will get classified into the rail service at the marshalling yards of the current node, and there is no loading operation. Consequently, in Case 1, the operation start time and operation cutoff time separately refer to the classification start time and classification cutoff time of the rail service. In this case, the selection of a rail service is constrained by its classification operation time window at the current node. Contrary to Case 1.1, because there are no loading/unloading operations and inventory services in this case, corresponding costs do not exist. Case 2.
The hazardous materials will not depart from the current node once the loading/classification operation is completed. They will wait until the scheduled departure time of the rail service at the same node, and then depart from the current node along with the train. Case 3. The hazardous materials will arrive at the successor node covered on the scheduled route along with the rail service at its scheduled arrival time at the same node. After arriving at the node, they cannot get unloaded/disassembled immediately and should wait until the unloading/disassembled operation start time of the rail service at the same node.
Contrary to the rail service, the organization of the road service is time-flexible. Transshipment is unnecessary if the hazardous materials both arrive at and then depart from the node by road services. A road service route can be entirely or partly covered in the hazardous materials multimodal routes. For the convenience of modelling, we can divide the road service routes into several segments, e.g., a road service route (1, 2, 3, 4) is divided into (1,2), (1,3), (1,4), (2,3), (2,4), and (3,4). The adjacent segments should not be covered in a hazardous materials route at the same time to avoid extra loading/unloading costs and transportation delay. Furthermore, there are no such constraints (Case 1 to Case 3) above when utilizing the road service to transport hazardous materials from the current node to the successor node. Consequently, there exist the following Cases 4 and 5 for road service in the multimodal service network. Case 4. The hazardous materials can immediately be loaded onto the road vehicles once they arrive at the current node. Once the loading operation is finished, the hazardous materials can immediately depart from the current node.
Case 5. The hazardous materials can immediately be unloaded from the road vehicles once they arrive at the successor node.

Risk Evaluation and Modelling
In order to manage and control the potential risk, risk evaluation on the planned route for the transportation of hazardous materials is quite important and necessary in decision-making. People and the environment are the key factors and indexes to evaluate the risk.

Social Risk Evaluation
The purpose of social risk evaluation is to determine how many people along the transportation routes of hazardous materials will be potentially impacted and how much the degree of potential risk will be. The degree can be reflected by the volume of the hazardous materials: the larger the volume, the higher the risk degree will be. We can multiply the population exposure and its corresponding volume of the hazardous materials in order comprehensively to express the two aspects of the social risk quantitatively, i.e., Social risk " population exposureˆvolume of the hazardous materials The population exposure can be calculated from the population density and exposure area (area of the potential dispersion zone of the hazardous materials once an accident occurs). The exposure area is determined by the dispersion distance threshold in the accidental release of hazardous materials. The dispersion distance is related to the spatial distribution of the toxic concentration level. In addition, such spatial distribution is influenced by various factors, including the types of the hazardous materials, atmospheric and geographical conditions of the accident, and so on [27]. Specifically, in this study, we consider that hazardous materials such as chlorine and methane become airborne in the accidental release.
There are many models for estimating such spatial distribution, among which the Gaussian plume model is the one commonly used [18]. Also, many standards formulated by the U.S. Environmental Protection Agency refer to this model [28,29]. So, in this study, we adopt the Gaussian plume model to estimate the spatial distribution of the toxic concentration level. The Gaussian plume model comprehensively considers the distance to the release source, the release rate, and the wind speed and direction in order to evaluate the toxic concentration level of the hazardous materials in the accidental release [24]. It assumes that the release rate of the hazardous materials and the surrounding atmospheric conditions remain constant during the dispersion process. When the release source and impact zone are at zero elevation, the concentration C pxq (unit: mg/m 3 ) at downwind distance x to the release source is formulated by Equation (1) [17]: In Equation (1), Q is the source release rate of the hazardous materials (unit: mg/s); u is the average downwind speed (unit: m/s); σ y and σ z are separately horizontal and vertical dispersion coefficients, and are determined by the atmospheric stability category and downwind distance x. There exist σ y " a¨x b and σ z " c¨x d , where a, b, c, and d are the dispersion parameters. In China, their values are formulated by the Technical Methods for Making Local Emission standards for Air Pollutants (GB/T 3840-91) according to the atmospheric stability category and downwind distance. Q is determined by various factors, including the type of the hazardous materials, the size of the damaged container, the pressure and temperature inside the container, and so on. Let VOL and T denote the maximal volume that a container can carry and the release duration, respectively. The estimation of Q can be simplified by Equation (2) [17]: For each hazardous material, the concentration threshold C that is immediately dangerous to life and health (IDLH concentration threshold) is known. By using this threshold value and modifying Equation (1), we can get the impact distance threshold x from Equation (3): In this study, we consider that the potential dispersion direction of the hazardous materials in the accidental release is 360˝. When evaluating the population exposure, we should consider the worst cases at the nodes and on the arcs. For node i, its population exposure is given by Equation (4), where ρ i is the population density around node i: For arc pi, jq, when using service s to transport the hazardous materials, its population exposure is given by Equation 5, where ρ ijs is the population density along path pi, j, sq and d ijs is the transportation distance of service s on arc pi, jq:

Environmental Risk Evaluation
Besides people, the environment will also be potentially impacted by the transportation of hazardous materials. For airborne hazardous materials, once accidental release occurs, surrounding ambient air and some sensitive zones (e.g., nature reserves, forests, rivers, lakes, and farmland) will be polluted. If the release volume of the hazardous materials in the ambient air exceeds the environmental capacity, the environment will be permanently damaged and so is difficult to clean up. Therefore, the environmental capacity should be considered in the hazardous materials road-rail multimodal routing problem. When evaluating the environmental risk, we first calculate the environmental capacities of the nodes and of the arcs.
We simplify by assuming that the potential impact spaces at nodes and at arcs are as in the box model, i.e., the potential impact space at a node is hemispherical while that at an arc is semi-cylindrical. Under this assumption, we can easily estimate the environmental capacity as follows. For node i, environmental capacity is given by Equation (6): For arc pi, jq, when using service s to transport the hazardous materials, its environmental capacity is given by Equation (7): where C˚is the allowed emission concentration threshold. It is formulated by the Integrated Emission Standard of Air Pollutants (GB 16297-1996). R˚is the potential impact radius of the hazardous materials when accidental release occurs. After determining the environmental capacity, we can then evaluate the environmental risk according to the definition proposed by Zhao in her doctoral dissertation [30]. For node i, environment risk is given by Equation (8), where q i is the cumulative hazardous materials volume that node i takes: For arc pi, jq, when using service s to transport the hazardous materials, environmental risk is given by Equation (9), where q ijs is cumulative hazardous materials volume that path pi, j, sq takes: The values of ER i and of ER ijs reflect the aggregation degrees of the risk at a certain node and on a certain path, respectively, and can also reflect the environmental risk distribution in the multimodal service network.

Notations
In this study, G " pN, A, Sq represents a road-rail multimodal service network, where N, A, and S are the node set, the directed arc set, and transportation service set in the network, respectively. In addition, for hazardous materials flow k (k P K, where K is the flow set), its transportation demands, including origin o k , destination d k , release time at origin t k release , due date T k , and the volume q k (unit: ton), are all determined and known in the location-routing problem. The remaining symbols in the model and their representations are listed as follows. Some of the following symbols are similar to our previous study [26].
‚ Indices h, i, j´index of the nodes, and h, i, j P N. s, r´index of the transportation services, and s, r P S. pi, jq´directed arc from node i to node j, and pi, jq P A. pi, j, sq´path from node i to node j served by service s.

‚ Sets
S ij´s et of transportation services on arc pi, jq, S ij Ď S and S ij " Γ ij Y Ω ij , where Γ ij and Ω ij are separately the rail service set and road service set on arc pi, jq.
δ´piq´predecessor node set to terminal i, and δ´piq Ď N. δ`piq´successor node set to terminal i, and δ`piq Ď N. l s i , u s i ‰´l oading/unloading operation time window of rail service s at node i, where l s i is the operation start time and u s i is the operation cutoff time, and i P N 2 . SA s i´s cheduled arrival time of rail service s at node i. SD s i´s cheduled departure time of rail service s from node i. Q ijs´c apacity of rail service s on arc pi, jq, unit: ton. d ijs´t ransportation distance of path pi, j, sq, unit: km. t ijs´t ransportation time of service s on arc pi, j, q, unit: h. POP ijs´p opulation exposure along path pi, j, sq. POP i´p opulation exposure around node i. CAP ijs´e nvironmental capacity of path pi, j, sq, unit: ton. CAP i´e nvironmental capacity of node i, unit: ton. ER max´e nvironmental risk threshold for the multimodal service network. c ijs´u nit costs of transporting hazardous materials flow on arc pi, jq by service s, unit: ¥/ton. c s´u nit costs of loading/unloading operation of service s at the node, unit: ¥/ton. Specifically: c rail´u nit costs of loading/unloading operation of rail service at nodes, unit: ¥/ton. τ´inventory period free of charge, unit: h. c store´u nit inventory costs of rail service at node, unit: ¥/ton-h. M´a significant large positive number.

Decision Variables
W k i´0 -1 decision variable. If node i is included in the route for hazardous materials flow k, W k i = 1; otherwise W k i = 0. X k ijs´0 -1 decision variable. If hazardous materials flow k is transported along path pi, j, sq, X k ijs = 1; otherwise X k ijs = 0. Y k i´a rrival time of hazardous materials flow k at node i. Z k ijs´c harged inventory time of hazardous materials flow k at node i before transported by rail service s on arc pi, jq, unit: h.

A Node-Arc-Based Model
Using the symbols above, the specific routing problem is first formulated as a bi-objective mixed integer nonlinear programming model in Equations/Components (10)- (31).
Objective 1 is to minimize the total generalized costs of routing the multiple hazardous material flows. Objective 1 consists of three components including the transportation costs (Component (10)), loading/unloading operation costs (Component (11)), and inventory costs (Component (12)) at the nodes. Objective 1 reflects the private goal of transportation demanders to pay the least capital to realize their transportation demands. Note that loading/unloading costs are only created (1) at the goods yards of the hazardous materials stations, where there are transshipments between road services and rail services; and (2) at the origins and destinations of the flows of hazardous materials. Therefore, the extra loading and unloading costs calculated in the transshipments between different rail services should be deducted. However, the inventory costs are only created in the first case. Consequently, we formulate the total loading/unloading costs and inventory costs as Components (11) and (12), respectively.
Objective 2 is to minimize the cumulative social risk along the planned routes for the transportation of hazardous materials. Objective 2 consists of two components, including the cumulative social risk at the nodes (Component (13)) and the cumulative social risk on the paths (Component (14)). Objective 2 reflects the public goal of lowering the potential impacted population and the potential risk degree in the multimodal service network.
Equation (15) is the general flow conservation constraint. The combination of Equations (15) and (16) ensures that each hazardous materials flow is non-bifurcated [31].
Equation (17) ensures that after dividing the road services route into several segments according to the method in Section 3.1, transshipments between road services are extra and unnecessary and should be avoided in the routing [26].
Equation (18) ensures the compatibility requirements between decision variables W k i and X k ijs . Because the origins and destinations are always included in the routes, for i P to k , d k u and @k P K, there exists W k i = 1. For the other nodes, if and only if a path linked with node i is selected for hazardous materials flow k, i.e., ř jδ`piq ř sPS ij X k ijs = 1, there exists W k i = 1, otherwise W k i = 0. Consequently, we formulate Equation (18) as above.
Equation (19) is the environmental risk constraint. It ensures that the total environmental risk in the multimodal service network does not exceed a given threshold, so that the environmental risk distribution in the multimodal service network will be balanced and not excessively aggregated.
Equation (20) is the rail service capacity constraint. It ensures that the loaded and classified volume of the hazardous materials should not exceed its available carrying capacity.
Equation (21) is the operation time window constraint of the rail service that is adopted to move the hazardous materials from the current node to the successor node. This constraint contains three independent cases: , which means the arrival time of the hazardous materials at the node is not restricted by the unselected rail services.
If X k ijs = 1 and ř hPδ´piq ř sPΩ ij X k hir = 1, according to Equation (16), ř hPδ´piq ř sPΓ ij X k hir = 0, i.e., the hazardous materials k arrive at node i by road service, then according to Equation (21), We assume that the arrival times of the flows of hazardous materials at their origins are their release times. Therefore, we have Equation (22) as above.
Equations (23) and (24) ensure the compatibility requirements between decision variables X k ijs and Y k i . They calculate the arrival times of the flow of each hazardous material at the selected nodes. Especially for Equation (24), when adopting rail service r to transport hazardous materials k from predecessor node h to current node i (X k hir = 1): (i) If the following transportation is by road service ( ř jPδ`piq ř sPΩ ij X k ijs " 1, and according to Equation (15), there exists ř jPδ`piq ř sPΓ ij X k ijs = 0); we are far more concerned about the unloading start time l s i instead of SA s i , because the following operations cannot be conducted until l s i . Therefore, in such a case, If the following transportation is by rail service ( ř jPδ`piq ř sPΓ ij X k ijs " 1, and according to Equation (16), there exists ř Equation (25) is the due date constraint. It ensures that the arrival times of the flows of hazardous materials at their respective destinations should be later than the due dates claimed by customers.
Equation (26) ensures the compatibility requirements among decision variables X k ijs , Y k i and Z k ijs . It calculates the charged inventory times of the flow of each hazardous material at the nodes. Note that, when X k ijs " 0, @Z k ijs ě 0 all satisfy Equation (26), but the minimization of Component (12) in Objective 1 will finally restrict Z k ijs " 0.

Improved Linear Reformulation
The initial model built in Section 4.2 is easy to understand. However, due to the nonlinear Components (11) and (12) in Objective 1 and nonlinear Equations (21), (23), (24), and (26), which contain both multiplications of decision variables and the maximum function, this model is a typical nonlinear programming. These nonlinear functions make the model quite difficult to be solved by the exact solution algorithms. So, in this section, we design equivalent linear forms for the nonlinear components in order to reformulate the initial model. Then the model can be easily solved with the help of standard mathematical programming software, in which many exact solution algorithms can be effectively implemented. Although the equivalent linear forms are not very straightforward, the effectiveness of solving the reformulated mixed integer linear programming model is significantly improved. For detailed linear reformulation proofs, readers can refer to Appendix A: Model Linear Reformulation.
By replacing these nonlinear objective functions and constraints with the equivalent linear formulations given in Appendix A, the initial model is transformed to a mixed integer linear programming model shown as follows. We can then straightforwardly utilize the standard mathematical programming software (e.g., Lingo) to implement the exact solution algorithm (e.g., Branch-and-Bound Algorithm) to solve the model effectively.

Normalized Weighted Sum Method for the Bi-Objective Optimization
Multi-objective optimization problems have two kinds of solutions, dominated solutions and Pareto solutions (non-dominated). For the multi-objective optimization with dominated solutions, its multiple objectives will reach their respective optimum simultaneously, i.e., the optimal solution for a certain objective is also the optimal one for the remaining objectives. Consequently, there only exists one optimal solution for such a multi-objective optimization. However, in most cases, the objectives in the multi-objective optimization are in mutual conflict, and the objectives cannot reach their optimum simultaneously, i.e., the respective optimal solutions for the objectives are different. There usually exists a group of solutions to such multi-objective optimization, namely Pareto solutions or non-dominated solutions. A group of Pareto solutions forms the Pareto frontier of the multi-objective optimization problem. An example of the Pareto solutions to the multi-objective optimization can be seen in Figure 1 [32]. In multi-objective decision making, with the help of the Pareto frontier of the problem, decision-makers can make a trade-off between different objectives.
Multi-objective optimization problems have two kinds of solutions, dominated solutions and Pareto solutions (non-dominated). For the multi-objective optimization with dominated solutions, its multiple objectives will reach their respective optimum simultaneously, i.e., the optimal solution for a certain objective is also the optimal one for the remaining objectives. Consequently, there only exists one optimal solution for such a multi-objective optimization. However, in most cases, the objectives in the multi-objective optimization are in mutual conflict, and the objectives cannot reach their optimum simultaneously, i.e., the respective optimal solutions for the objectives are different. There usually exists a group of solutions to such multi-objective optimization, namely Pareto solutions or non-dominated solutions. A group of Pareto solutions forms the Pareto frontier of the multi-objective optimization problem. An example of the Pareto solutions to the multi-objective optimization can be seen in Figure 1 [32]. In multi-objective decision making, with the help of the Pareto frontier of the problem, decision-makers can make a trade-off between different objectives. There are many methods that can be utilized to generate Pareto solutions to a multi-objective optimization problem. One of the most classical and extensively used methods is the weighted sum method. In this method, different weights are distributed to the objectives, and then the weighted objectives are combined linearly to transform the multi-objective optimization problem into a singleobjective one. It is a simple method that is easy to understand and utilize. Hence it has already received wide application, e.g., Xie et al. [22], Sheu [33], Samanlioglu [34], and Rakas et al. [35]. Specifically for the bi-objective optimization model constructed in this study, after weighted summation, there are: There are many methods that can be utilized to generate Pareto solutions to a multi-objective optimization problem. One of the most classical and extensively used methods is the weighted sum method. In this method, different weights are distributed to the objectives, and then the weighted objectives are combined linearly to transform the multi-objective optimization problem into a single-objective one. It is a simple method that is easy to understand and utilize. Hence it has already received wide application, e.g., Xie et al. [22], Sheu [33], Samanlioglu [34], and Rakas et al. [35]. Specifically for the bi-objective optimization model constructed in this study, after weighted summation, there are: minimize λ 1¨f1`λ2¨f2 (31) where f 1 and f 2 represent the objective functions of minimizing the generalized costs and of minimizing the social risk along the planned routes, λ 1 and λ 2 are the distributed weights to the two objectives, and their values are determined by decision makers before model simulation. We can gain different solutions when different values are assigned to λ 1 and λ 2 , and finally generate the Pareto frontier for the problem.
In the classical weighted sum method, the objectives usually have different units and magnitudes, which will affect the optimization result when directly utilizing this method to address the multi-objective optimization problem. Consequently, in this study, we adopt an improved weighted sum method, namely the normalized weighted sum method, to solve our bi-objective optimization problem.
The normalized weighted sum method first solves the multi-objective optimization problem as a series of single-objective ones, and obtains their respective optimal solution. For the specific model in this study, let f1 and f2 separately denote the optimal values of Objective 1 and of Objective 2. After dividing by their respective optimal values, the objectives become normalized and there is no difference regarding the unit and magnitude among the different objectives. Therefore, in the normalized weighted sum method, the objective function is as below instead of as in Equation (31): By using the normalized weighted sum method introduced in this section, we can generate the Pareto frontier for the case that will be discussed in the next section.

Case Description and Parameter Setting
In this section, we design an empirical case regarding liquefied chlorine transportation to demonstrate the feasibility of the proposed methods in dealing with a practical problem. The empirical case is based on a Beijing-Tianjin-Hebei Region scenario in China, in which four chlorine-alkali chemical factories in Hebei Province supply liquefied chlorine to several waterworks factories in the Beijing-Tianjin-Hebei Region. The locations of the nodes in the region are shown in Figure 2. The population density distribution and average wind speed distribution of the region are also illustrated in Figure 2. They can be used to estimate the population exposure in the empirical case. The detailed topological structure of the multimodal service network for this empirical case is shown as Figure 3. It is composed of 45 nodes (including four origin nodes, 13 destination nodes, and 28 hazardous materials station nodes) and 132 directed arcs (including 31 rail service arcs and 101 road service arcs). The corresponding information on the road services in the multimodal service network is presented in Table B2. The schedules of the rail services in the multimodal service network as well as their risk parameters are given in Table B3. In Table B3, all the time data are discretized into real numbers (e.g., 10:30 on the first day is transformed into 10.5; on the second day, it is transformed into 34.5), and the same rail services operated in different periods are treated as different ones. In addition, the operation periods of all the rail services are one day/train.
The cost parameters in the road-rail multimodal routing problem of hazardous materials are In the large-scale empirical case, we assume that the liquefied chlorine is carried in 40-ft tank containers (in Equation (2), VOL = 30,480 kg). The most common consequence in the liquefied chlorine transportation accidents is that the liquefied chlorine constantly releases from the damaged tank container and disperses into the surroundings as gas. This kind of accident accounts for 95% of total liquefied chlorine transportation accidents [36]. In such a transportation accident, the release of chlorine usually lasts about one hour (in Equation (2), T = 3600 s) [37]. Therefore, according to Equation (2), when evaluating the population exposure, the release rate Q = 8.5 kg/s (8.5ˆ10 6 mg/s). The IDLH concentration threshold C in Equations (4) and (5) is 2500 mg/m 3 [37]. Considering the conservative evaluation of the social risk, the atmospheric stability category is set to Level D. Consequently, in Equations (4) and (5), a = 0.15, b = 0.89, c = 0.40, and d = 0.63. Based on the information provided above, we can obtain the population exposure and environmental capacity of each node shown in Appendix B, Table B1. When evaluating the environmental capacity, in Equations (6) and (7), the allowed emissions concentration threshold C˚is 80 mg/m 3 according to GB 16297-1996, and the potential impact radius R˚= 8 km [38]. Therefore, the environmental capacity of each node is 85,743 tons.
The detailed topological structure of the multimodal service network for this empirical case is shown as Figure 3. It is composed of 45 nodes (including four origin nodes, 13 destination nodes, and 28 hazardous materials station nodes) and 132 directed arcs (including 31 rail service arcs and 101 road service arcs). The corresponding information on the road services in the multimodal service network is presented in Table B2. The schedules of the rail services in the multimodal service network as well as their risk parameters are given in Table B3. In Table B3, all the time data are discretized into real numbers (e.g., 10:30 on the first day is transformed into 10.5; on the second day, it is transformed into 34.5), and the same rail services operated in different periods are treated as different ones. In addition, the operation periods of all the rail services are one day/train.   The cost parameters in the road-rail multimodal routing problem of hazardous materials are stated as follows. For the rail service, its transportation unit costs can be determined by c ijs " c 1 rail`c 2 rail¨d ijs , where c 1 rail = 27.9 ¥/ton and c 2 rail = 0.206 ¥/ton-km. Its unit loading/unloading costs c rail = 5.8 ¥/ton. The unit inventory costs c store = 0.1 ¥/ton-h. The free-of-charge inventory period is 48 h. For the road service, its transportation unit costs can be determined by c ijs " c 2 road¨d ijs , where c 2 road = 0.76 ¥/ton-km. Its unit loading/unloading costs c road = 5.5 ¥/ton. In this empirical case, there are 25 hazardous material flows that need to be routed. Their origins, destinations, volumes, release times, and due dates are all listed in Table 1. In Table 1, the release times and due dates are all real numbers.

Simulation Environment
We set the environmental risk threshold ER max = 0.6. Then this empirical case, which is formulated as a linear programming problem, can be solved by the exact solution algorithms. In this study, we utilize the Branch-and-Bound Algorithm to solve the specific routing problem. This algorithm is implemented by the mathematical programming software Lingo 12 (LINDO Systems Inc., Chicago, IL, USA) on a ThinkPad Laptop with Intel Core i5-5200U 2.20 GHz CPU 8 GB RAM. The scale of the empirical case problem is indicated in Table 2.

Multimodal Routes Illustration
First of all, we use the optimal routes with minimal generalized costs as examples to indicate how to transport hazardous materials in the road-rail multimodal service network. In this empirical case, the least-costs road-rail multimodal routes are presented in Table 3, where items such as #40103 and #32123 represent the train numbers of the rail services. As we can see from Table 3, there are two kinds of routes for the transportation of hazardous materials in the multimodal service network, including single road service routes and road-rail multimodal routes. In the latter, rail services are dominant, and road services are helpful supplements to realize the through transportation from shippers to receivers. Such routes are quite suitable for transporting hazardous materials with longer lead times. The minimal generalized costs are 850,192 ¥ in this empirical case. The structure of these costs is shown in Figure 4. As we can see from Figure 4, costs en route account for the majority of the total generalized costs (in this case, approximately 90%). The planning horizon of the empirical case is limited; the inventory times of the hazardous materials are hence shorter than the free-of-charge inventory period. Consequently, as for this empirical case, there are no inventory costs in the minimal generalized costs.
The minimal generalized costs are 850,192 ¥ in this empirical case. The structure of these costs is shown in Figure 4. As we can see from Figure 4, costs en route account for the majority of the total generalized costs (in this case, approximately 90%). The planning horizon of the empirical case is limited; the inventory times of the hazardous materials are hence shorter than the free-of-charge inventory period. Consequently, as for this empirical case, there are no inventory costs in the minimal generalized costs.

Single-Objective Scenarios
First we address the two single-objective optimization problems that separately aim at minimizing the total generalized costs (unit: ¥) and the social risk (unit: 10 4 people-ton). The respective optimal solutions to the three objectives are given in Table 4. Table 4 clearly indicates that the two objectives cannot reach their respective optimal values simultaneously, i.e., the optimal solution to a certain objective is not the optimal one for the remaining objective. Therefore, there exist Pareto solutions for the bi-objective optimization problem.

Bi-Objective Scenario
In bi-objective decision making, because the two objectives are conflicting, decision-makers should make a trade-off between lowering the generalized costs and reducing the social risk. By using the normalized weighted sum method introduced in Section 5, we can generate Pareto solutions (Pareto frontier) to the bi-objective optimization problem as in Figure 5. Figure 5 shows that improving the economy of the multimodal routing will lead to an increase in the social risk along the planned routes, and reducing the social risk along the planned routes will lower the economy of the multimodal routing. Figure 5 provides various candidate solutions to the bi-objective optimization problem that will be helpful for decision-makers trying to balance the conflict between cost objectives and risk objectives. For example, in the trade-off between generalized costs and social risk, if decision-makers

Single-Objective Scenarios
First we address the two single-objective optimization problems that separately aim at minimizing the total generalized costs (unit: ¥) and the social risk (unit: 10 4 people-ton). The respective optimal solutions to the three objectives are given in Table 4. Table 4 clearly indicates that the two objectives cannot reach their respective optimal values simultaneously, i.e., the optimal solution to a certain objective is not the optimal one for the remaining objective. Therefore, there exist Pareto solutions for the bi-objective optimization problem.

Bi-Objective Scenario
In bi-objective decision making, because the two objectives are conflicting, decision-makers should make a trade-off between lowering the generalized costs and reducing the social risk. By using the normalized weighted sum method introduced in Section 5, we can generate Pareto solutions (Pareto frontier) to the bi-objective optimization problem as in Figure 5. Figure 5 shows that improving the economy of the multimodal routing will lead to an increase in the social risk along the planned routes, and reducing the social risk along the planned routes will lower the economy of the multimodal routing. Figure 5 provides various candidate solutions to the bi-objective optimization problem that will be helpful for decision-makers trying to balance the conflict between cost objectives and risk objectives. For example, in the trade-off between generalized costs and social risk, if decision-makers consider that it is acceptable as long as the normalized social risk is lower than 1.05, then the last four Pareto solutions in Figure 5 can be their candidate multimodal route schemes. Considering minimization of the generalized costs, decision-makers will finally select the multimodal route corresponding to the 7th Pareto solution to transport the hazardous materials.
consider that it is acceptable as long as the normalized social risk is lower than 1.05, then the last four Pareto solutions in Figure 5 can be their candidate multimodal route schemes. Considering minimization of the generalized costs, decision-makers will finally select the multimodal route corresponding to the 7th Pareto solution to transport the hazardous materials.

Sensitivity of the Pareto Solutions with Respect to the Environmental Risk Threshold
In this section, we conduct sensitivity analysis of the Pareto solutions with respect to the environmental risk threshold. Because the minimal environmental risk in the multimodal service network is 0.135, we generate three Pareto frontiers corresponding to the thresholds 0.2, 0.4, and 0.6. Figure 6 shows the variation of the Pareto frontiers when the environmental risk threshold increases. In Figure 6, the generalized costs and social risk are separately normalized by 850,192 and 506,362, which are the optimal values of generalized costs and social risk, respectively, when the environmental risk threshold is 0.6.

Sensitivity of the Pareto Solutions with Respect to the Environmental Risk Threshold
In this section, we conduct sensitivity analysis of the Pareto solutions with respect to the environmental risk threshold. Because the minimal environmental risk in the multimodal service network is 0.135, we generate three Pareto frontiers corresponding to the thresholds 0.2, 0.4, and 0.6. Figure 6 shows the variation of the Pareto frontiers when the environmental risk threshold increases. In Figure 6, the generalized costs and social risk are separately normalized by 850,192 and 506,362, which are the optimal values of generalized costs and social risk, respectively, when the environmental risk threshold is 0.6.
consider that it is acceptable as long as the normalized social risk is lower than 1.05, then the last four Pareto solutions in Figure 5 can be their candidate multimodal route schemes. Considering minimization of the generalized costs, decision-makers will finally select the multimodal route corresponding to the 7th Pareto solution to transport the hazardous materials.

Sensitivity of the Pareto Solutions with Respect to the Environmental Risk Threshold
In this section, we conduct sensitivity analysis of the Pareto solutions with respect to the environmental risk threshold. Because the minimal environmental risk in the multimodal service network is 0.135, we generate three Pareto frontiers corresponding to the thresholds 0.2, 0.4, and 0.6. Figure 6 shows the variation of the Pareto frontiers when the environmental risk threshold increases. In Figure 6, the generalized costs and social risk are separately normalized by 850,192 and 506,362, which are the optimal values of generalized costs and social risk, respectively, when the environmental risk threshold is 0.6. As we can see from Figure 6, the Pareto frontier of the bi-objective optimization problem will move to the left with the increase of the environmental risk threshold. That is to say, we can lower the generalized costs and the social risk simultaneously if we relax the environmental risk constraint by enlarging the environmental risk threshold. The extrapolated Figure 6 can also provide a reference for decision-makers when they are not sure of the setting of the environmental risk threshold for the multimodal service network.

Conclusions
In this study, we explore the road-rail multimodal routing problem of hazardous materials from the operational level of network planning. The main contribution made in this study is its comprehensive formulation of the following characteristics, so that the specific routing problem can be addressed from an improved practical viewpoint: (1) specific customer demands from transportation economy and efficiency; (2) multiple hazardous materials flows; (3) railway schedule-based space-time constraints; (4) environmental risk constraint; and (5) bi-objective optimization from the viewpoint of generalized costs vs. social risk management. Consideration of all these formulation characteristics significantly enhances the feasibility of the proposed model in dealing with practical problems, which is also a good development of current studies on the routing problem of hazardous materials.
Besides the improvement in model formulation, this study also develops concise linear reformulations that transform the initial nonlinear model into its equivalent linear programming. On the one hand, the linear reformulations enable the problem to be solved effectively by an exact solution algorithm that can be easily implemented by standard mathematical programming software. On the other hand, the linearization method can also provide an exact benchmark in order to systematically test the performance of various heuristic algorithms in dealing with the specific routing problem. In addition, the exact solution algorithm, e.g., the Branch-and-Bound Algorithm, can also test whether the model itself is mathematically logical or not.
In this large-scale empirical case study, by using the Branch-and-Bound Algorithm in Lingo 12, the optimal solution for a certain scenario can be generated within 3 min. Although the computational time for a single solution is acceptable, a large amount of computational time will be consumed when several Pareto solutions are required and a sensitivity analysis needs to be conducted. Therefore, it is still necessary to design a heuristic algorithm to address the problem more efficiently. Consequently, our further study will focus on the design of the heuristic algorithm. As claimed above, this study has provided an exact benchmark for testing the heuristic algorithm, which has provided a solid foundation for us to carry on such a study.

Appendix A. Model Linear Reformulation
. Consequently, the equivalency above is verified. ‚ Reformulation 2. By using a non-negative decision variable V k ijs to replace the nonlinear function Z k ijs¨řhPδ´piq ř rPΩ ij X k hir and adding two linear Equations (A5) and (A6), the nonlinear Component (12) in Objective 1 can be reformulated as a linear function (A4).
In this reformulation, the following two cases obviously exist, which proves the equivalency above.
If ř hPδ´piq ř rPΩ ij X k hir = 0, i.e., hazardous materials k arrive at node i by the rail service, there are V k ijs ě´M and V k ijs ě 0 according to Equations (A5) and (A6), and the minimization of Component (A4) will restrict V k ijs = 0, which matches the setting that there are no inventory costs at the node if the hazardous materials are transshipped between different rail services. (ii) If ř hPδ´piq ř rPΩ ij X k hir = 1, i.e., hazardous materials k arrive at node i by road service, there are V k ijs ě Z k ijs and V k ijs ě 0, and minimization of Component (A4) will restrict V k ijs " Z k ijs , which matches the setting that inventory costs will be charged at the node according to the charged inventory time if the hazardous materials are transshipped from road service to rail service. (21) is equivalent to linear Equation (A7) as follows.

Reformulation 3. Nonlinear Equation
It is obvious that we can gain the same results as Equation (21) after we separately input the three given cases, Case (i)-Case (iii), from Equation (21) into Equation (A7). Therefore, the two constraints are equivalent to each other. (23), (24) and (26) are separately equivalent to linear Equations (A8)-(A13). For detailed proofs of the above three equivalencies, readers can refer to our previous study [26].