Covering models are the most widespread location models for locating emergency facilities. The objective of the covering models is to provide “coverage” to the demand points. The “coverage” is deemed as the service provided by facilities to demand points, and the distance between a demand point and the corresponding facility is less than a pre-set distance limit. In order to provide full coverage of all demand points, Toregas et al. first proposed the location set covering problem (LSCP) to minimize number of facilities needed [8C. Toregas, R. Swain, C. ReVelle, and L. Bergman, "The location of emergency service facilities", Oper. Res., vol. 19, no. 6, pp. 1363-1373, 1971.
[http://dx.doi.org/10.1287/opre.19.6.1363] ]. However, since full coverage is required in LSCP, the required facilities number may be larger than the actual available amount, resulting in an unreasonable situation. Aiming at this problem, Church and ReVelle [9R. Church, and C. Revelle, "The maximal covering location problem", Pap. Reg. Sci., vol. 32, no. 1, pp. 101-118, 1974.
[http://dx.doi.org/10.1007/BF01942293] ] and White and Case [10J.A. White, and K.E. Case, "On covering problems and the central facilities location problem", Geogr. Anal., vol. 6, no. 3, pp. 281-294, 1974.
[http://dx.doi.org/10.1111/j.1538-4632.1974.tb00513.x] ] develop the maximum covering location problem (MCLP), which does not require full coverage to all demand points. Instead, the model seeks to provide as much coverage as possible based on a given number of facilities.

Since LSCP aims to minimize the facility number and MCLP to maximize the coverage under limited facilities, they both take no consideration of service efficiency of transportation network. Service efficiency could be evaluated by the total distance between the demand points and selected facilities. A larger value of total distance indicates that the accessibility and effectiveness of facilities is worse. The P-median problem takes efficiency as an objective and seek to find a solution with the minimum total distance between demand points and facilities. So, the models are often applied in occasions where proximity is desirable, such as supermarkets, post offices, as well as emergency service facilities location. Dawson et al. used a hybrid method integrating P-median and P-center models to find locations with minimized travelling distance and minimized maximum distance between missile site and required security forces for the Minuteman Intercontinental Ballistic Missile (ICBM) nuclear weapon system [11M.C. Dawson, J.E. Bell, and J.D. Weir, "A hybrid location method for missile security team positioning", J. Business Manage., vol. 13, pp. 5-17, 2007.]. Serra et al. introduced classic scenario approach to solve uncertain demand and travel time in P-median problem when locating public facilities [12D. Serra, and V. Marianov, "The p-median problem in a changing network: The case of Barcelona", Location Sci., vol. 6, pp. 383-394, 1998.
[http://dx.doi.org/10.1016/S0966-8349(98)00049-7] ]. Afshartous et al. proposed a simulation-optimization methodology to address uncertain distress call in P-median model facing with US Coast Guard air station location problem [13D. Afshartous, Y. Guan, and A. Mehrotra, "US Coast Guard air station location with respect to distress calls: A spatial statistics and optimization based methodology", Eur. J. Oper. Res., vol. 196, no. 3, pp. 1086-1096, 2009.
[http://dx.doi.org/10.1016/j.ejor.2008.04.010] ]. Dantrakul et al. compared greedy algorithm based method, p-median method and p-center method when solving facility location problem where the objective function is to minimize total transportation and setup costs [14S. Dantrakul, C. Likasiri, and R. Pongvuthithum, "Applied p-median and p-center algorithms for facility location problems", Expert Syst. Appl., vol. 41, no. 8, pp. 3596-3604, 2014.
[http://dx.doi.org/10.1016/j.eswa.2013.11.046] ].

In contrast to P-median problem models, which concentrate on optimizing the overall performance of the whole system, the P-center model attempts to minimize the worst performance of the system, thus it is also known as minimax model. The P-center model considers that a demand point is served by its nearest facility and therefore full coverage to all demand points is always achieved.

2.1.3. Decision Variables

The formulation for the original P-median model is as follows:

(1)

(2)

(3)

(4)

(5)

(6)

where Objective (1) is to minimize the total distance between DPs and DCs; Constraint (2) ensures that only the selected DCs could serve the DPs; Constraint (3) ensures that each DP is served by one and only one DC; Constraint (4) states that there are N DCs to be located at the candidate locations; finally, constraints (5) and (6) enforce the integrality of variables x_ij and y_j.

After earthquake disasters, the causes of road damage are multiple. According to the statistics and analysis results of historical seismic road damage data, 7 factors have great influence on road damage degree. They are seismic intensity (α₁), roadbed soil type (α₂), site classification (α₃), subgrade failure (α₄), roadbed type (α₅), roadbed height difference (α₆), seismic fortification intensity (α₇). The categories of the 7 factors and corresponding quantized values are given in Table 1.

Table 1
Seismic road damage factors and corresponding quantized values.

The 7 factors impact road damage at the same time, so road damage index (RDI) is proposed to reflect the integrated influence [15Y. Chen, and X. Chen, "An introduction to computer aid system for damage prediction of highway transportation system", Build. Sci., vol. 04, pp. 71-75, 1994. [In Chinese]., 16S. Jiang, and C. Liu, "Vulnerability analysis of transportation system in Sanya city", World Earthquake Eng., vol. 03, pp. 23-27, 2005. [In Chinese].]. The relationship between RDI and 7factors is:

(7)

Different values of RDI indicate different road damage degrees. The higher value of RDI represents severer seismic damage, and lower RDI declares slighter damage. The ranges of RDI and corresponding damage degrees are obtained according to analysis of historical data and given in Table 2.

Table 2
The ranges of RDI and corresponding damage degrees.

Therefore, to obtain the damage degree of roads in seismic affected areas, the RDIs should be computed by the product of the 7 factors. However, the accurate value of these factors could not be acquired exactly in time after earthquake. Immediately after an earthquake disaster, α₁ could be approximately obtained by spatial analysis of road net layer and estimated seismic intensity layer. α₂, α₃, α₅, α₆, α₇ might be available through queries to traffic information database in some developed areas with perfect traffic informatization level, but they could not be easily got in every road at the range of a nation or all over the world. However α₄ could only be accessed by field investigation which is a time-consuming work.

As a consequence, assuming α₁ is deterministic and α₂, α₃, α₅, α₆, α₇ are stochastic values obeying certain distributions, the values in Table 1 could be analyzed according to different values of seismic intensity (α₁). The statistical results are shown in Table 3. The max_in, min_in and mean_in are the maximum, minimum and mean values of RDI_in when α₁=in.

Table 3
Statistical results of RDI_in under different seismic intensity degrees.

The historical experience in earthquake disaster shows that RDIs under different seismic intensity degrees obey normal distribution RDI_in~N(μ_in, σ_in²). According to the unique characteristics of normal distribution, the probability that the value of RDI_in locates in the interval (-3σ_in+μ_in, μ_in+3σ_in) is 0.9973, and the probability beyond the interval is 0.0027, which is considered as small probability event (Fig. 2). So it could be assumed that the values of RDI_in completely locate in (-3σ_in+μ_in, μ_in+3σ_in) and the relationships between max_in, min_in, mean_in and μ_in, σ_in are:

(8)

(9)

Fig. (2) The ±3σ+μ interval of normal distribution.

The computational results of normal distributions parameters for RDIs in different seismic intensity degree areas are given in Table 4 and the corresponding probability density function graphs are shown in Fig. (3).

Table 4
The parameters of normal distributions for RDIs in different seismic intensity degree areas.

Fig. (3) The probability density function graphs for RDIs in different seismic intensity degree areas.

Since different RDIs reflect different road damage situations and higher value of RDI indicate severer damage, the traffic speed after seismic disaster would vary depending on road damage degrees, which is quantified as RDI. As a low value of RDI stands for slight road damage, the real traffic speed is close to designed traffic speed. However, when RDI is high, the road is severely damaged and the real traffic speed is much lower than the designed traffic speed.

In P-median model, the distance between facility and demand point is the only measurement in the objective. The value of distance is obtained under the condition that road is intact and real traffic speed is equal to design traffic speed. However, after earthquake disasters, roads are usually damaged in different degrees, thus leading to irrationality of original P-median model. So, improvements should be made to adjust P-median model to be suitable for facility location problems after earthquake disasters.

Because road damage would result in the decline of vehicle traffic speed, and lengthen traffic time as a result. So, an artificial increase of distance could offset the influence of speed decline being equivalent to the effect of road damage.

In this study, a distance penalty coefficient ξ is introduced to increase road distance, as an equivalent method to model the influence of road damage. Assuming the undamaged road distance is d, the following equations would be satisfied.

(10)

(11)

dⁱⁿ is the road distance exposed i,so the sum of dⁱⁿ is the total travel distance between a DP and a DC as illustrated in Fig. (4). d^' is the equivalent distance with distance penalty coefficient considering road damage. ξ_in is the distance penalty coefficient suitable for road damaged in the area with seismic intensity degree in, and ξ_in should satisfy the constraint that the ξ_in in higher seismic intensity degree is larger than that in lower seismic intensity degree, represented as:

(12)

Fig. (4) Illustration of road sections located in different seismic intensity degree areas.

The populations of the 22 counties are obtained by spatial analysis in Geographic Information System (GIS) using Gridded Population of the World, Version 3 (GPWv3) provided by NASA Socioeconomic Data and Applications Center (SEDAC) [20"Center for International Earth Science Information Network", Columbia University. Gridded Population of the World, Version 3 (GPWv3): Population Count. Grid.. (Accessed May 1st, 2015).
[http://dx.doi.org/10.7927/H4639MPP ] ]. The population distribution pattern of the study area is shown in Fig. (6) and the precise numbers are given in Table 5.

Fig. (6) The population distribution pattern in 22 counties.

Table 5
The populations of 22 counties.

The road distances between any two counties of the 22 counties are calculated by network analysis in GIS using the nation-wide road network layer (Fig. 7). The analysis result is presented as a symmetric matrix composed of elements d_pq(p,q=1,2,∙∙∙,22) whose value is given in Table 6. Furthermore, the lengths of road sections located in different seismic intensity degrees (dpqin) could be obtained through overlay analysis in GIS with the help of estimated seismic intensity map generated by the model proposed by Howell and Schultz [21B.F. Howell, and T.R. Schultz, "Attenuation of modified Mercalli intensity with distance from the epicenter", Bull. Seismol. Soc. Am., vol. 65, no. 3, pp. 651-665, 1975.] and modified by Dasheng Chen and Hanxing Liu [22D. Chen, and H. Liu, "Elliptical attenuation relationship of earthquake intensity", North China Earthquake Sci., vol. 7, no. 3, pp. 31-42, 1989.]. The value amount of d_pqⁱⁿ(in=7,8,∙∙∙,11) is large and not given in detail in this paper.

Fig. (7) Overlay analysis between road network and estimated seismic intensity map.

Table 6
The value of d_pq between country p and q.

The distance penalty coefficient is highly related with RDIin and it also complies with normal distribution. The normal distribution parameters in each seismic intensity degree area are analyzed and decided by experts. In the study area covering 22 counties, the roads are mainly located in the areas with seismic intensity degrees from 7 degree to 11 degree, so the parameters of are given in Table 7.

Table 7
The parameters of normal distributions for ξ_in in different seismic intensity degree areas.

Based on the input parameters specified in the above dataset, we solve the facility location problem for distribution center on the MATLAB R2012b platform. Besides, the deterministic model is also built and the parameters are given as the expectations of the stochastic parameters.

The computational results in Table 8 show that the objective values of deterministic model and the proposed stochastic model are almost the same, however, the solution sets between two models are not different. In the deterministic model, the road damage situation is set to be a known and explicit parameter, so the locations of DCs are optimally determined as county 1, 4, 8, 11 and 18, more visually indicated in the map in Fig. (8a). Since the stochastic model simulate the uncertain road damage which is similar to the actual situation, three solution sets are obtained through the repeated sampling of road damage in different seismic intensity degrees. The locations of DCs are varying among the following three solutions: (1) county 1, 4, 8, 11, 18; (2) county 1, 5, 8, 11, 18; (3) county 2, 5, 8, 11, 18 (Figs. 8b, c, d). The solutions are chosen according to the actual road situation after earthquake disaster. However, the pattern of probabilities of the above three solutions could be found through the iterations of model optimization. The probabilities of solutions are P1=0.51 for solution (1), P2=0.47 for solution (2), P3=0.02 for solution (3), indicating that solution (1) and (2) are generally adopted in most simulated road damage situations but solution (3) is just suitable for the extreme cases which rarely occurs in the simulations.

Table 8
The computational results.

Fig. (8) The solutions of numerical analysis: (a) the solution for the deterministic model; (b), (c), (d) the solutions for the stochastic model.

The deterministic model provides only one solution of DCs locations on the assumption that the road damage situations in different seismic intensity degrees are strictly equal to the expectations of the stochastic distributions. We have to admit that the model indeed give a location solution of DCs, however, the solution may not be suitable to all the road damage situations. The stochastic model has the advantage over the deterministic one for the more comprehensive and complete solution covering more situations, thus makes it possible to provide alternative solutions for decisions makers.