Resilience-based seismic retrofit of urban infrastructure systems

INTRODUCTION

The battle against natural hazards and disasters has been a persistent theme throughout human history. Over time, mankind has made remarkable progress in eliminating or alleviating the threat posed by natural hazards and their resulting disasters. Earthquakes, as one of the common natural disasters, present significant risks to urban infrastructure, human life, and property. For instance, the 2008 Wenchuan earthquake collapsed 20 power plants and 170 substations, destroyed 80% of the buildings in Beichuan County, and caused approximately 150 billion dollars in economic losses and 69, 180 known deaths^[1]. Urban infrastructure systems (UISs), including electric power, water supply, and gas systems, are quite vulnerable to earthquakes as their physical facilities, such as substations, pipes, and roads, are sensitive to ground shaking^[2]. The collapse of these physical facilities not only threatens the lives of people but also affects post-disaster humanitarian relief due to the essential roles of UISs in the functioning of modern society. A successful solution to mitigate the impact of disaster events is to build more resilient UISs^[3,4]. Compared to other concepts such as reliability, risk, and safety, resilience emphasizes the comprehensive ability of UISs to resist and absorb negative impacts, recover rapidly from disasters, and adapt to better cope with future events^[5,6].

Numerous researchers committed to the resilience enhancement problem of UISs against earthquakes and put forward a variety of constructive strategies across the four phases of disaster relief, including mitigation, preparedness, response, and recovery phases. The mitigation phase entails identifying risks and hazards to either substantially reduce or eliminate the impact of an incident usually through structural measures. The preparedness phase intends to reduce the system failure probability, and the relevant actions include deploying backup systems^[7], extending system topology^[8], and retrofitting components^[9-11]. The response phase attempts to resist the diffusion of failures and takes emergency actions to ensure the normal functioning of critical facilities. For example, operators can adjust the topology of electric power networks by the pre-installed switches on transmission lines to isolate the faulty section^[12]. The recovery phase aims to make an effective plan to repair damaged facilities and restore the services of UISs^[13,14]. Researchers established many mathematical models to describe the recovery process of UISs, in which various factors are taken into account, such as available resources, routes of repair crews, preferences of stakeholders, decision environments, and interdependencies across UISs^[15-17]. This article mainly focuses on the strategy of component retrofitting, which can decrease the failure probabilities of components facing disturbances and has been widely adopted in the literature and practice.

Generally, a UIS consists of thousands of components with different types, indicating that retrofitting each component is impracticable and costly. Hence, only partial components that are essential to system resilience are selected to be retrofitted. Researchers in the literature have proposed numerous methods to explore the critical components of UISs, and a common approach is based on the component importance index, which describes the importance of a component to the whole system^[18-20]. The component importance index could be measured in accordance with component types, topological characteristics (e.g., degree and betweenness), and geographical location^[21-23]. In the context of electric power systems, plants are the most critical facilities, followed by transmission substations and lines and distribution substations and lines. Also, a component with a large degree value has a high priority to be retrofitted as it connects many components in the system, and its failure might cause a large impact. Moreover, some studies measure the component importance index from the perspective of reliability and vulnerability. Espiritu et al. introduced several existing reliability criticality measures in the literature, including Birnbaum importance, criticality importance, reliability reduction worth, and reliability achievement worth, and then developed a novel measure for electric power systems^[24]. Salman et al. utilized risk achievement worth to evaluate the importance of components in electric power distribution systems subjected to hurricanes^[25]. Here, risk achievement worth describes the "worth" of a component in achieving system reliability. Li et al. developed a probability-based method to evaluate the seismic reliability of substations and identify the critical components in an electric power system^[26]. Rocco et al. described a vulnerability analysis method to identify critical components for protection from the perspective of improving network performance^[27].

The aforementioned component importance-based methods are straightforward and easily implemented, but they ignore the synergistic effect between components, resulting in a low improvement performance to system resilience. Selecting a limited set of components for reinforcement is a typical combinatorial optimization problem, and researchers have established corresponding optimization models to seek the set of components that can bring the largest benefit to system resilience. Yuan et al. focused on the resilience enhancement problem of electric power distribution systems and established a two-stage robust optimization model with the consideration of the uncertain occurrence of disasters due to natural hazards^[28]. Two resilience enhancement strategies, retrofitting components and deploying distributed generation resources, are incorporated into this model. Yan et al. studied the seismic retrofit problem of a railway system and formulated a stochastic model to seek the optimal railway stations and tracks to be strengthened under a limited budget^[29]. Ma et al. developed a tri-level optimization model to enhance the resilience of power distribution networks against extreme weather events, whose objective is to minimize the hardening investment and load shedding cost^[30]. Lu et al. developed a mean-risk two-stage stochastic programming model to investigate the transportation network protection problem against extreme events^[31]. The first stage intends to seek the retrofitting strategy of highway bridges with the minimum retrofitting cost, whereas the second stage minimizes the travel cost given retrofitting decisions and hazard scenarios. Liu et al. established a two-stage stochastic model to study the seismic retrofit problem of UISs and developed a novel heuristic method to solve this problem^[32]. Numerical experiments were implemented on three electric power systems under seismic scenarios to illustrate the validity and superiority of this method. Several studies integrated the post-disaster restoration decision problem into the pre-disaster retrofit problem. Miller-Hooks et al. studied the resilience enhancement of freight transportation networks and developed a two-stage stochastic model which simultaneously considers pre-disaster preparedness and post-disaster recovery actions^[33]. Gomez and Baker also developed a two-stage stochastic model to address the coupled pre-disaster and post-disaster decision problem in a transportation network under seismic hazards^[34].

The retrofit problem of UISs includes massive uncertainties, such as seismic hazard occurrence probabilities, components failure probabilities, and restoration time, which make those developed optimization models difficult to be exactly solved.Researchers utilized various methods to reduce the computation complexity, including robust programming, Monte Carlo simulation, and importance sampling. Miller-Hooks et al. proposed a solution methodology with the incorporation of the integer L-shaped method and Monte Carlo simulation^[33]. Romero et al. developed a knapsack-based heuristic method to optimize the selection of seismic retrofit strategies^[35]. Also, several studies put forward methods to generate damage scenarios. Adachi and Ellingwood used seismic hazard maps to determine component failure probabilities and generate component damage scenarios^[7]. Gomez and Baker utilized a probabilistic risk assessment of transportation networks to generate hazard-consistent scenarios^[34].

This article studies the seismic retrofit problem of UISs and formulates a two-stage stochastic model. The first stage attempts to seek an optimum seismic retrofit strategy under a limited budget that takes its future benefit into account, whereas the future benefit is quantified by the expected system resilience to all possible seismic scenarios in the second stage. The restoration decision model for each seismic scenario is established in the second stage, and the system resilience describes the cumulative system functionality during the whole restoration process (i.e., from the initial time of a disaster event to the completion time of all repair actions). Then, a resilience-based heuristic method is introduced to solve this mathematical model. This heuristic method first generates a limited set of component damage scenarios, and then the original stochastic model is reformulated into an approximated model, which is solved by the integer L-shaped method. Also, the sample average approximation method is adopted to enhance the solution accuracy. Finally, this heuristic method is applied to the resilience-based seismic retrofit problem of the Qujing power transmission system to demonstrate its validity.

The main contributions of this article include (1) establishing a resilience-based seismic retrofit optimization model for UISs with the incorporation of post-disaster repair actions; (2) proposing an efficient heuristic method to solve the stochastic model; and (3) validating the heuristic method on an electric power transmission system. The proposed retrofit optimization model and efficient heuristic algorithm can be integrated into a decision support system and help the government officials with the seismic investment decision-making. The remainder of this article is organized as follows. Section 2 presents the resilience metric and models the seismic hazard scenarios. Section 3 formulates the resilience-based seismic retrofit problem. Section 4 introduces a heuristic solution method. Taking the Qujng power transmission system as an example, Section 5 presents the numerical results. Section 6 concludes and discusses the findings and provides directions for future research.

MODEL DESCRIPTION

This article studies the seismic retrofit problem of UISs to maximize system resilience under a limited budget. A general UIS can be modeled as an undirected network $$ G $$($$ N $$, $$ E $$), where $$ N $$ and $$ E $$ denote the set of nodes (e.g., power/water plants and substations) and edges (i.e., transmission lines and pipes), respectively. Nodes are roughly classified into source nodes (e.g., power/water plants) and demand nodes (i.e., substations). The set of source nodes and demand nodes are denoted by $$ N^S $$ and $$ N^D $$. Let $$ s_n $$ and $$ \bar{s}_n $$ be the real and maximal output of a source node $$ n \in N^{S} $$ and $$ d_n $$ and $$ \bar{d}_n $$ be the real and required demand of a demand node $$ n \in N^{D} $$.

As the investment budget is limited, only some components can be retrofitted. Compared with un-retrofitted components, retrofitted components have lower failure probabilities of being damaged when facing a seismic event. For simplicity, this article only considers node retrofit, and the candidate of nodes to be retrofitted is denoted by $$ N^R $$. Let binary variable $$ w_n $$ denote retrofit decision, with 1 indicating that node $$ n \in N^{R} $$ is retrofitted, and 0 otherwise. The cost of retrofitting node $$ n $$ is represented by $$ c_n $$, and the total investment budget is expressed by the parameter $$ B^R $$.

Note that UISs provide essential services (i.e., electricity, water, and transportation) to customers (i.e., residents, factories, and other UISs); the evaluation of UIS resilience needs to consider the expectations of customers, i.e., whether system service can be restored within expected critical times after a disruptive event. So, it is more reasonable and practical to quantify UIS resilience based on those time points that are critical to customers. As shown in Figure 1, $$ P_R $$($$ t $$) denotes the real restoration curve of a UIS, and time series $$ \left\{t_{c_{1}}, t_{c_{2}}, t_{c_{3}}, t_{c_{4}} \right\} $$ are four critical time points of concern. To capture this characteristic, this article quantifies UIS resilience based on a series of time points that are critical to customers. Here, critical time points could be determined by the preference of customers or expert opinions. Hence, the resilience $$ R $$ of a UIS under a disruptive event is measured as follows^[36]:

Figure 1. A typical restoration curve and four critical time points of concern.

where $$ \left\{t_{c1}, t_{c2}, \ldots, t_{cm}\right\} $$ represent the critical time points of concern, $$ P_R(t{_{ci}i}) $$ denotes the system functionality level at the time point $$ t_({ci}) $$, and $$ w_i $$ describes the weight of system functionality at time point $$ t_({ci}) $$.

To calculate the system functionality level at each time point, two factors should be known: the state of each component and the operation mechanisms of the UIS of concern. The state of each component depends on the initial component damage scenario $$ \xi \in \Xi $$ and the restoration decision. Here, the component damage scenario $$ \xi $$ is a vector consisting of the state of each component. The definitions of damage states are based on the $$ Hazus's technical manual $$, a total of five damage states are defined for UIS components, and they are none, slight, moderate, extensive, and complete^[37]. This article assumes that a component will lose its functionality if it falls into the damage limit state exceeding "extensive" (i.e., extensive and complete). Let binary variables $$ \xi_{n}^{no} $$, $$ \xi_{n}^{ex} $$, and $$ \xi_{n}^{co} $$ denote the normal, extensive, and complete damage states, respectively. For a given seismic scenario, the failure probability of each component can be estimated through its fragility curves for the "extensive" and "complete" limit states. As for the restoration decision, this article assumes that the available repair resources are characterized by the number of the repair crews, and the maximum amount of available repair resources is represented by $$ RR $$. This article assumes that the working efficiency of repair crews is identical, and each repair crew repairs the assigned damaged components independently, with each damaged component needing only one repair crew. The repair time of damaged component $$ n $$ under the component damage scenario $$ \xi $$ is expressed by $$ \tau_{n}(\xi) $$.

Each type of UISs, transporting commodities (i.e., electricity, water, and gas) from the supply side (i.e., power and water plants) to the demand side (i.e., factories and residential districts) through lines (i.e., transmission lines and pipes), has a particular operating mechanism. The network flow model has been frequently used to simulate the operating mechanisms of UISs. However, for the electric power transmission system to be studied in the case study, the direct current power flow (DCPF) model is a better alternative and has been frequently used in the field of electrical engineering^[38].

MATHEMATICAL FORMULATION

This section proposes a two-stage stochastic optimization model for the resilience-based seismic retrofit to maximize the seismic resilience of UISs under a limited retrofit budget. A graphical representation of the optimization model is shown in Figure 2, the first stage is for the system planner to make retrofit decisions, and the second stage is for the system operator to accelerate the restoration process and optimize system operation under the given retrofit strategy. The first stage attempts to seek an optimum seismic retrofit strategy under a limited budget that takes its future benefit into account, whereas the future benefit is quantified by the expected resilience to all possible seismic scenarios in the second stage. The second stage searches for the best repair decision given component damage scenarios and retrofit strategy.

Figure 2. A graphical representation of the two-stage stochastic optimization model.

Denote a binary variable $$ x_n $$($$ \xi $$, $$ {t_{c_i}} $$) by the damage state of node $$ n $$ at time point $$ {t_{c_i}} $$ under damage scenario $$ \xi $$, with its value 1 indicating normal operation, and 0 otherwise. Denote a binary variable $$ {r_{kn}} $$($$ \xi $$, $$ {t_{c_i}} $$), which represents whether node $$ n $$ is repaired by repair team $$ k $$ at the time point $$ {t_{c_i}} $$, with 1 for repaired and 0 otherwise. The proposed mathematic model is introduced as follows:

The objective function (2) is to maximize the expected resilience under all generated component damage scenarios, where $$ P(\xi|\mathit{\boldsymbol{w}}, q) $$ is a retrofit decision-dependent probability, capturing the fact that the occurrence probability of component damage scenario $$ \xi $$ is affected by retrofit strategy $$ \mathit{\boldsymbol{w}} $$ and seismic scenario $$ q $$. The expansion equation of $$ P(\xi|\mathit{\boldsymbol{w}}, q) $$ is shown as follows:

where $$ p_{n, q}^{n o, b}, p_{n, q}^{e x, b}, p_{n, q}^{c o, b}, p_{n, q}^{n o, a}, p_{n, q}^{e x, a}, p_{n, q}^{c o, a} $$ represent the probability of node $$ n $$ falling into three different damage states before and after being retrofitted under a seismic scenario $$ q $$. The product term on the right side of Equation (17) is the occurrence probability that nodes have the functionality states defined by $$ \xi $$, where $$ \left[\left(1-w_{n}\right) p_{n, q}^{n o, b}+w_{n} p_{n, q}^{n o, a}\right] $$ gives the probability that node $$ n $$ works normally ($$ \tilde{\xi}_{n}^{no}=1 $$) under the retrofit decision $$ {w_{n}^{N}} $$, $$ \left[\left(1-w_{n}\right) p_{n, q}^{e x, b}+w_{n} p_{n, q}^{e x, a}\right] $$ and $$ \left[\left(1-w_{n}\right) p_{n, q}^{c o, b}+w_{n} p_{n, q}^{c o, a}\right] $$ are the probabilities that node $$ n $$ falls into "extensive" ($$ \tilde{\xi}_{n}^{e x}=1 $$) and "complete" ($$ \tilde{\xi}_{n}^{co}=1 $$) damage states, respectively, under $$ {w_{n}^{N}} $$.

Constraint (3) limits the retrofit budget, and Constraint (4) enforces binary retrofit decision variables. Constraints (5)-(10) describe the DCPF model. Constraint (5) sets the phase angle of the reference node as zero. Constraint (6) states that the flow of each edge is determined by its susceptance and the phase angles and operation states of its origin and destination nodes. Constraint (7) ensures flow conservation, and Constraint (8) states the flow capacity of each edge. Constraints (9)-(10) state the maximum output of each source node and the target demand of each demand node at different critical time points. The constraints for the recovery decision variables are described by Constraints (11)-(16). Constraint (11) ensures that if node $$ n $$ is operational in the network at the beginning of time period $$ i $$, it must have been repaired by some repair group at the beginning of that period. Constraint (12) states that non-damageable nodes are always operational. At the beginning of the restoration process, no damaged component is repaired, as stated by constraint (13). Constraint (14) ensures that the total elapsing time for those components which have been repaired from period 1 to $$ i $$ does not exceed time point $$ t_{c_i} $$. Constraint (15) states that a damaged node is only repaired one time at most. Constraint (16) enforces the recovery decision variables as binary.

SOLUTION ALGORITHM

This section introduces an efficient heuristic method that takes advantage of several existing methods. The proposed method mainly includes the following three steps: (1) generates limited component damage scenarios to reformulate the original problem as an approximated model; (2) adopts a retrofit efficacy-based method to reduce the solution space and applies the integer L-shaped method to solve the approximated model; (3) employs the sample average approximation method to enhance the solution quality.

For the first step, this article adopts the following procedures to generate the limited component damage scenarios: (1) randomly generates a large number $$ K $$ of pre-retrofit and post-retrofit component damage scenarios $$ \left\{r_{n, k, q}^{N, b}, r_{n, k, q}^{N, a}\right\} $$ (1 if node $$ n $$ fails, and 0 otherwise, $$ k $$=1,2,…,$$ K $$) for each seismic scenario $$ q $$ using Monte Carlo simulations ("post-retrofit" refers to the case that all candidate components are retrofitted); (2) selects $$ H $$ from $$ K $$ scenarios to minimize the total gap error for both component failure probability and system resilience. The gap error, in terms of component failure probability, is the sum of the overestimating and underestimating errors for component failure probabilities before and after being retrofitted. The gap error, in terms of system resilience, is the sum of the absolute difference between the estimated system resilience calculated using those limited component damage scenarios and the "true" system resilience calculated using large-scale component damage scenarios in the case with no component being retrofitted, plus that difference for the case with all the components being retrofitted. The two types of gap errors are respectively normalized by the "true" component failure probability and "true" system resilience to make them comparable, and the weight coefficient $$ \alpha $$ is initially set as 0.5.

Define a binary decision variable $$ y_{k, q} $$ which is 1 if pre-generated scenario $$ k\in{1,2,,K} $$ under seismic scenario $$ q $$ is selected; and define the occurrence probability of this scenario by a continuous decision variable $$ \rho_{k, q} $$. Denote the gap errors resulting from overestimating and underestimating the pre-retrofit and post-retrofit failure probabilities of node $$ n $$ under seismic scenario $$ q $$ by $$ {err}_{n, q}^{b+}, {err}_{n, q}^{b-}, {err}_{n, q}^{a+}, {err}_{n, q}^{a-} $$; the gap errors between the estimated and the "true" system resilience for the pre-retrofit and post-retrofit systems by $$ {err}_{r, q}^{b+}, {err}_{r, q}^{b-}, {err}_{r, q}^{a+}, {err}_{r, q}^{a-} $$; the system resilience under pre-retrofit component damage scenario $$ {r}_{n, k, q}^{N, b} $$ and post-retrofit scenario $$ {r}_{n, k, q}^{N, a} $$ by $$ {f}_{k, q}^{b} $$, $$ {f}_{k, q}^{a} $$; the "true" system resilience for the pre-retrofit and post-retrofit systems by $$ {F}_{r, q}^{b} $$, $$ {F}_{r, q}^{a} $$, respectively. The mathematical model for identifying a limited set of component damage scenarios together with their occurrence probabilities under seismic scenario $$ q $$ is formulated as follows:

The objective function (18) minimizes the sum of gap errors for component failure probability and system resilience. Constraints (19)-(20) define the gap errors with respect to pre-retrofit and post-retrofit failure probabilities of components. Constraints (21)-(22) define the gap errors with respect to pre-retrofit and post-retrofit system resilience. Constraint (23) ensures the number of selected scenarios under seismic scenario $$ q $$ not larger than a pre-set size $$ h_q $$. The value of $$ h_q $$ is proportional to the product of the occurrence probability of seismic scenario $$ q $$ and the expected functionality loss of the pre-retrofit system under seismic scenario $$ q $$. Note that a rigorous restriction of $$ H $$ (the total number of required component damage scenarios under all seismic scenarios) might cause $$ h_q $$ decimal, so to the algorithm rounds $$ h_q $$. This may make the number of final generated component damage scenarios slightly more or less than $$ H $$. Constraint (24) forces the occurrence probability of a component damage scenario to zero if it is not selected. Constraint (25) makes the sum of the occurrence probabilities of all scenarios as one. Constraints (26)-(28) ensure that each error term and occurrence probability is nonnegative, and each scenario selection variable is binary.

Based on limited component damage scenarios, the original seismic retrofit optimization problem can be reformulated as an approximated model, which is a standard integer program with linear constraints, with the objective function (2) updated as follows:

For the second step, this article adopts a retrofit efficacy-based method to reduce the solution space of the seismic retrofit optimization problem and then solves the updated problem using the integer L-shaped method. In the retrofit efficacy-based method, the original objective function is replaced by some leading terms of its Taylor series expansion. Peeta et al. approximated the objective function of a pre-disaster investment decision problem by using the first-order term of its Taylor series expansion^[9] and then reformulated the problem as a knapsack problem, with its solution being a local optimum of the original problem. Denote the objective function in Equation (29) by $$ F(\mathit{\boldsymbol{w}}) $$, which is defined only at the vertices of the unit hypercube, $$ U=\left\{\boldsymbol{w} \mid 0 \leq w_{n} \leq 1, n \in N^{R}\right\} $$. Relaxing the integrality restrictions on the components of $$ \mathit{\boldsymbol{w}} $$ allows $$ F(\mathit{\boldsymbol{w}}) $$ to be continuously differentiable in the domain $$ U $$, and hence enables the consideration of its Taylor series expansion in the neighborhood of some $$ \mathit{\boldsymbol{w}}_0 \in U $$:

where $$ g_{n}\left(\boldsymbol{w}_{0}\right)=\left.\frac{\partial F(\boldsymbol{w})}{\partial w_{n}}\right|_{\boldsymbol{w}=\boldsymbol{w}_{\bf{0}}} $$ is the first-order derivative with respect to the investment decision for node $$ n $$ at $$ \boldsymbol{w}_{0} $$, $$ g_{n_{1} n_{2}}\left(\boldsymbol{w}_{0}\right)=\left.\frac{\partial^{2} F(\boldsymbol{w})}{\partial w_{n_{1}} \partial w_{n_{2}}}\right|_{\boldsymbol{w}=\boldsymbol{w}_{0}} $$ is the second-order derivative with respect to the investment decision for node $$ n_1 $$ and node $$ n_2 $$ at $$ \boldsymbol{w}_{0} $$, and so forth. Applying the retrofit efficacy-based method with the first order approximation, the original seismic retrofit problem is reformulated as follows:

Where $$ g_n (0)=F(u_n)-F(0) $$, $$ u_n $$ is the unit vector of dimension $$ |N| $$ with 1 at node $$ n $$ and 0 at the remaining nodes. $$ F(\mathit{\boldsymbol{w}}) $$ could be calculated through the Monte Carlo simulation method or based on the importance sampling. Once $$ g_{n}(0) $$ is calculated, the reformulated program is a 0-1 knapsack problem, which can be solved efficiently by the dynamic programming algorithm^[39] or by the branch-and-bound algorithm.

The solution space reduction is realized through the following procedures: (1) enlarges the retrofit budget to be $$ cp \ast B^R $$, where $$ cp \ge 1 $$ is a control parameter; (2) applies a retrofit efficacy-based method under the retrofit budget $$ cp \ast B^R $$ to get an initial solution $$ {w_n^{initial}} $$, which defines the reduced solution space. The idea is based on the assumption that most components in the optimum solution of the original problem are included in the solution obtained from the retrofit efficacy-based method under an enlarged budget $$ cp \ast B^r $$. If $$ cp \ge \sum_{n \in N^{r}} c_n/B^R $$, all components are potential candidates to be retrofitted, so the solution space is not reduced; if $$ cp $$ is slightly larger than 1, the solution space will be largely reduced. After reducing the solution space, the following constraints should be added to the approximated model:

Constraint (32) ensures the solution space is limited to the space solved by the retrofit efficacy-based method with $$ cp > 1 $$. Other constraints for the approximated model include Constraints (3)-(16). Together with the objective function (29), the approximated model can be easily solved using the integer L-shaped method.

The third step applies the sample average approximation method to enhance the solution quality. Kleywegt et al. stated that if the computational complexity of solving the sample average approximation problem increases faster than linearly with the sample size, it is more efficient to choose a smaller sample size and to replicate generating and solving several sample average approximation problems^[40]. Hence, the proposed method replicates generating and solving several sample average approximation problems with middle sample sizes (a small number of component damage scenarios) through the above two steps to return several candidate retrofit strategies. Among those strategies, the optimum retrofit strategy is determined by comparing their performance gain using large-scale Monte Carlo simulations.

RESULTS

System test data and seismic hazard simulation

This article adopts the electric power transmission system in the central area of Qujing, Yunnan province, China, for a case study, which contains 8 gate stations, 35 substations (twenty-six 110 kV substations and nine 35 kV substations), and 56 transmission lines (forty-five 110 kV transmission lines and eleven 35 kV transmission lines), as shown in Figure 3. The cost of retrofitting a component is estimated by multiplying the cost of anchoring a transformer by the number of transformers in the component^[41].

Figure 3. The electric power transmission system in the central area of Qujing, Yunnan province, China.

Yunnan Province (21-29$$ ^{\circ} $$ N, 97-106$$ ^{\circ} $$ E), situated in the southeastern region of the Qinghai-Tibetan Plateau, experiences high levels of crustal activity as a result of being extruded by the Indian and Eurasian plates. Qujing is positioned on the edge of the Qiaojia-Dongchuan seismic zone in eastern Yunnan Province, where earthquakes occur frequently. In accordance with the Seismic Motion Parameter Zonation Map of China GB18306-2015, there are three distinct ground motion levels: basic ground motion (50-year exceedance probability of 10%), frequent ground motion (50-year exceedance probability of 63%), and rare ground motion (50-year exceedance probability of 2%). This article adopts the rare ground motion to demonstrate the proposed approach and the spatial distribution of peak ground acceleration (PGA) is shown in Figure 4.

Figure 4. Peak ground acceleration (PGA) distribution for rare ground motion in the central area of Qujing.

The comparison of the "true" fragilities and estimated fragilities of components in the electric power transmission system located in the central area of Qujing under the rarely occurred earthquake is shown in Figure 5. The upper sub-figure is fragilities for all components before retrofitting, and the bottom sub-figure is fragilities for all components retrofitting. The circle represents the "true" fragility of each component, which is calculated by using large-scale Monte Carlo simulations with sufficient component damage scenarios, and the dot represents the estimated fragility of each component with the generated limited component damage scenarios by the proposed heuristic method (PHM). Results show that the positions of the circle and dot basically coincide and the average error between the "true" and the estimated fragilities is 0.0004, which is quite small.

Figure 5. Comparison of the "true" (circle) and the estimated (dot) fragilities of components in the electric power transmission system located in the central area of Qujing under the rarely occurred earthquake.

Impact of restoration resources on retrofit strategy

This article further analyzes the impact of restoration resources on retrofit strategy. Table 3 shows the retrofit strategies under different retrofit budgets and different numbers of restoration resources. Results show that under the same retrofit budget, the retrofit strategies are varied with different restoration resources, which means that the number of post-earthquake restoration resources would impact the pre-earthquake retrofit strategies. For example, when the retrofit budget is set to ＄2.4M, one retrofitted component is node 34 with $$ RR $$ = 1 and 3, and this component is replaced by node 38. Also, nodes 3, 5, 6, 26, 32, and 36 are always retrofitted under different retrofit budgets and $$ RR $$, which indicates that these nodes are the most critical components of this system.

Table 3

Retrofit strategies under varied retrofit budgets and different amounts of restoration resources

Retrofit budget	RR = 1	RR = 2	RR = 3
1, 800	3, 5, 6, 26, 32, 36	3, 5, 6, 26, 32, 36	3, 5, 6, 26, 32, 36
2, 400	3, 5, 6, 26, 28, 32, 34, 36	3, 5, 6, 26, 28, 32, 36, 38	3, 5, 6, 26, 28, 32, 34, 36
3, 000	3, 5, 6, 26, 28, 32, 33, 34, 36, 38	3, 5, 6, 7, 26, 28, 32, 34, 36, 38	3, 5, 6, 7, 26, 28, 32, 34, 36, 38
3, 600	3, 5, 6, 7, 26, 28, 30, 32, 33, 34, 36, 38	3, 5, 6, 7, 26, 28, 30, 32, 33, 34, 36, 38	3, 5, 6, 7, 26, 28, 30, 32, 33, 34, 36, 38
4, 200	3, 5, 6, 7, 19, 26, 28, 30, 32, 33, 34, 35, 36, 38	3, 5, 6, 7, 19, 26, 28, 30, 32, 33, 34, 35, 36, 38	3, 5, 6, 7, 19, 21, 26, 28, 30, 32, 33, 34, 36, 38

Figure 6 further shows the resilience curves under varied retrofit budgets and different amounts of restoration resources. Results show that when the amount of restoration resources is equal to 1, and the retrofit budget increases from 1.8 to 4.2 million dollars, the resilience of the electric power transmission system in the central area of Qujing is linearly improved from 0.7123 to 0.7448. Moreover, with the number of restoration resources increasing from 1 to 2 and then to 3, the resilience value enlarges with a similar extent under different retrofit budgets, and the resilience curves are approximately parallel.

6. Resilience curve under varied retrofit budgets and different amounts of restoration resources.

In addition, to illustrate the regional differentiation of resilience, Figure 7 shows the spatial distribution of resilience for each street district in the central area of Qujing under varied retrofit budgets and different amounts of restoration resources. Results show that the retrofit budget and the amount of restoration resources both influence the resilience of different districts. On the one hand, when the amount of restoration resources $$ RR=1 $$, the area of the green zone gets larger and larger when the retrofit budget increases from $$ B^R=0 $$ to $$ B^R=\text{＄}4.2M $$. On the other hand, when the retrofit budget $$ B^R=\text{＄}4.2M $$, the area of green zone continues increasing when the amount of restoration resources increases from $$ RR=1 $$ to $$ RR=3 $$. In sum, with the increase of the retrofit budget and the amount of restoration resources, the resilience values of different street districts show an overall rising trend.

Figure 7. Spatial distribution of resilience for each street district in the central area of Qujing under varied retrofit budgets and different amounts of restoration resources.

CONCLUSION AND FUTURE WORK

This article proposes a resilience-based seismic retrofit optimization model for UISs under a limited retrofit budget and an efficient heuristic algorithm for its solution and also analyzes the impact of post-earthquake restoration resources on pre-earthquake retrofit strategies. Results show that the PHM performs better than the existing SHM. In addition, the amount of post-earthquake restoration resources not only influences the calculated resilience level but also affects the pre-earthquake retrofit strategies when the retrofit budgets are identical. Also, the retrofit budget and the amount of restoration resources influence the spatial distribution of the resilience at the street district levels served by the UIS. The proposed model and the solution algorithm can be used by local and central government agencies to aid investment decisions to upgrade UISs for disaster response.

However, this study still has certain limitations, and there are several areas that can be explored in future research directions. First, apply the proposed method to the seismic retrofit optimization of interdependent infrastructure systems and collect hazard scenario data with a higher resolution. Second, formulate and solve the seismic retrofit optimization problem from a life-cycle perspective with the consideration of different types of hazards. Third, as the disruptions outside the system of concern significantly affect the system performance, the system boundary issue needs to be integrated into the problem formulation and taken into account for its solution.

DECLARATIONS