<figure><img src="image_1.png" /></figure>

HAL
open science

Two-Stage Stochastic Resilience Optimization of
Converter Stations Under Uncertain
Mainshock-Aftershock Sequences

Kai Wang, Zhihang Xue, Di Cao, Yu Liu, Yi-Ping Fang

<figure><img src="image_2.png" /></figure>

► To cite this version:

Kai Wang, Zhihang Xue, Di Cao, Yu Liu, Yi-Ping Fang. Two-Stage Stochastic Resilience Optimization of Converter Stations Under Uncertain Mainshock-Aftershock Sequences. IEEE Transactions on Power Systems, 2025, pp.1-13. 10.1109/TPWRS.2025.3561890. hal-05050530

HAL Id: hal-05050530

https://hal.science/hal-05050530v1

Submitted on 29 Apr 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Two-Stage Stochastic Resilience Optimization of
Converter Stations Under Uncertain
Mainshock-Aftershock Sequences

Kai Wang, Zhihang Xue, Di Cao, Yu Liu, Senior Member, IEEE, and Yi-Ping Fang

**Abstract**—Reducing vulnerability and enhancing resilience of power systems, particularly converter stations, against earthquakes are crucial for maintaining their safe operation. Nevertheless, earthquakes are often accompanied by aftershocks and encompass plenty of uncertainties, posing significant challenges for the development of a joint pre-earthquake preparation and post-earthquake restoration strategy. In this article, a novel two-stage stochastic programming model is put forth to bolster the resilience of converter stations under uncertainties associated with mainshock-aftershock sequences. The model is bifurcated into two pillars: the first stage focuses on designing equipment hardening strategy (EHS) and spare parts strategy (SPS) for converter stations before earthquakes. The second stage is dedicated to the optimization of recovery scheduling (RS) in the wake of mainshock-aftershock sequences. To mitigate the endogenous uncertainties of the model, a tolerated random number generation method is employed to generate deterministic scenarios in the second stage. The resulting optimization model is, then, efficiently solved using a combination of the sample average approximation and progressive hedging algorithm. To evaluate the performance of the proposed resilience optimization method, comparative studies are carried out on a real ±800 kV converter station located in the city of Yibin, Sichuan Province, P.R. China. Results from this application demonstrate that the EHS pre-dominantly focuses on hardening the equipment of the same circuit to guarantee the low transmission capacity under diverse earthquake scenarios, while the SPS is proven to be more effective than the EHS in scenarios with high values of peak ground acceleration (PGA) of mainshock and occurrence intensity of aftershocks.

**Index Terms**—Converter station, endogenous uncertainty, mainshock-aftershock sequences, resilience, stochastic programming.

NOTATIONS

<table><tr><td>E</td><td>Set of earthquake scenarios</td></tr><tr><td>E<sup>ξ</sup></td><td>Set of shocks in earthquake scenario ξ</td></tr><tr><td>T</td><td>Set of recovery time point</td></tr></table>

This work was supported by the National Natural Science Foundation of China under Grant 72331002, and in part by Sichuan Science and Technology Program under Grant 2023YFG0089. The work of Yi-Ping Fang was supported by Agence Nationale de la Recherche under the project RoScaResilience ANR-22-CE39-0003. (Corresponding author: Yu Liu.)

Kai Wang, Zhihang Xue, and Di Cao are with the School of Mechanical and Electrical Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China (e-mail: wangkai3517@163.com; xuezhihang2@163.com; e-mail: caodi@uestc.edu.cn).

Y. Liu is with the School of Mechanical and Electrical Engineering, and the Center for System Reliability and Safety, University of Electronic Science and Technology of China, Chengdu 611731, China (e-mail: yuliu@uestc.edu.cn).

Yi-Ping Fang is with the Laboratoire Génie Industriel, CentraleSupélec, Université Paris-Saclay, 91190 Gif-sur-Yvette, France (e-mail: yiping.fang@centralesupelec.fr).

<table><tr><td>T<sup>ξ</sup></td><td>Set of shock occurrence time in earthquake scenario ξ</td></tr><tr><td>N</td><td>Set of equipment in converter station function logic diagram</td></tr><tr><td>L</td><td>Set of lines in converter station function logic diagram</td></tr><tr><td>S</td><td>Source node in converter station function logic diagram</td></tr><tr><td>D</td><td>Sink node in converter station function logic diagram</td></tr><tr><td>p<sub>i</sub><sup>T<sup>ξ</sup></sup></td><td>Failure probability of equipment i under k th shock in earthquake scenario ξ</td></tr><tr><td>δ<sub>i</sub><sup>T<sup>ξ</sup></sup></td><td>Reduced failure probability of equipment i under k th shock in earthquake scenario ξ</td></tr><tr><td>c<sub>ij</sub></td><td>Capacity of line (i, j) in converter station function logic diagram</td></tr><tr><td>R<sub>i</sub><sup>T<sup>ξ</sup></sup></td><td>Generated random number of equipment i under k th shock in earthquake scenario ξ</td></tr><tr><td>T<sub>k</sub><sup>ξ</sup></td><td>Occurrence time of k th shock in earthquake scenario ξ</td></tr><tr><td>w<sub>i</sub></td><td>Recovery time of equipment i</td></tr><tr><td>s<sub>i</sub><sup>ξ</sup>(t)</td><td>State of equipment i at time t in earthquake scenario ξ</td></tr><tr><td>g<sub>i</sub><sup>ξ</sup>(t)</td><td>Recovery schedule of equipment i at time t in earthquake scenario ξ</td></tr><tr><td>τ<sub>i</sub><sup>T<sup>ξ</sup></sup></td><td>0 if the equipment i is disrupted under the k th shock in earthquake scenario ξ, and 1, otherwise</td></tr><tr><td>n<sub>i</sub><sup>T<sup>ξ</sup></sup></td><td>1 if the hardening strategy of equipment i is effective under the k th shock in earthquake scenario ξ, and 0, otherwise</td></tr></table>

Decision variables

<table><tr><td>x<sub>ij</sub><sup>ξ</sup>(t)</td><td>Flow of line (i, j) at time t in earthquake scenario ξ</td></tr><tr><td>m<sub>i</sub><sup>ξ</sup>(t)</td><td>1 if the equipment i is chosen to be recovered at time t in earthquake scenario ξ, and 0, otherwise</td></tr><tr><td>n<sub>i</sub></td><td>1 if the equipment i is chosen to be hardened, and 0, otherwise</td></tr><tr><td>r<sub>i</sub></td><td>Number of spare parts of equipment i</td></tr></table>

I. INTRODUCTION

THE safe operation of power systems is paramount for the stability and advancement of society [1]. Recent years have
witnessed many large-scale blackouts due to earthquakes, lead-
ing to not only tremendous economic losses but also delayed
earthquake rescue processes [2]. These blackouts are primarily
caused by the broken or collapsed critical power equipment,
such as transformers and transmission lines. For instance, a se-
ries of eleven earthquakes, with a maximum magnitude of 6,
struck Malcolm, Aba Prefecture, Sichuan Province, China, on
June 10, 2022, resulting in the shutdown of 4 substations, 11
lines, and the blackout of 2,098 customers. The power restora-
tion has been lasting for three days after the earthquakes [3]. To
mitigate such huge losses, it necessitates a holistic framework
that encompasses both pre-earthquake preparation measures
and post-earthquake restoration measures to maintain and en-
sure the continuity of power supply.

Ultra-high voltage (UHV) power transmission, which refers to power lines operating at voltages greater than 800 kilovolts (kV), is of paramount importance in power systems for its advantages of high capacity and transmitting electricity over long distances with low power losses [4]. Thanks to the UHV power transmission technology, the energy supply-demand disparity issue between western and eastern China can be effectively solved [5]. Hence, reliable operation of UHV power transmission serve as a critical role for the safe and economic operation of power systems [6]. The converter stations, as core hubs, have a wide application in the UHV power transmission. The core function of converter stations is to convert low-voltage alternating current (AC) power into high-voltage direct current (DC) power for transmitting electricity over long distances. However, due to the towering structure of the converter station equipment and the complex physical connections between equipment, converter stations are exposed to high risks of natural disasters [7]. Once converter stations are damaged by earthquakes, it may not only cause large outages to power systems but also lead to huge economic losses. The ability of converter stations to resist earthquakes and quickly recover after earthquakes are, therefore, essential for the safe operation of power systems.

High-impact low-probability (HILP) events, such as earth-
quakes, hurricanes, and storms, pose significant challenges to
the safe operation of power systems due to their unpredictability
and catastrophic consequences [8]-[10]. For example, earth-
quakes can cause extensive damage to transmission lines, sub-
stations, and critical facilities, severely disrupting power trans-
mission and grid stability [11]. Hurricanes may devastate power
infrastructure in coastal regions, including cable towers and un-
derground conduits, resulting in widespread blackouts [12].
Sandstorms can contaminate equipment surfaces, increase the
failure rates of transmission lines, and pose long-term threats to
the performance of the power network [13]. These events often
disrupt the functionality of power systems, leading to prolonged
restoration processes and substantial economic losses. In recent
years, resilience has been extensively studied to against HILP
events [14],[15]. In general, resilience is defined as the system's
ability to withstand, respond, and adapt to various disruptive
events [16]. In this regard, resilience is composed of two critical
dimensions: vulnerability and recoverability. The former refers
to the inability of a system to sustain the consequence of dis-
ruptions, whereas the latter represents the ability of a system to

recover from disruptions promptly. Building on this foundation,
resilience optimization has been intensively implemented to
power systems under disasters [17]. Regarding vulnerability,
tremendous efforts have been made to address the pre-disrup-
tion preparation problems of power systems via protection, ca-
pacity expansion, backup component construction, voltage con-
trol, resource allocation, and facility location [18]-[20]. These
researches aimed to take preventive measures to reduce the im-
pact of disasters directly based on the disaster damage model
before the occurrence of disasters. Tremendous studies have
also focused on the recoverability to enhance resilience of
power systems [21]. These methods can be mainly categorized
by different strategies, i.e., network reconstruction, emergency
power supply, and emergency resource scheduling [22]-[24].
These research attempts aimed to take restoration measures to
recover the power supply as soon as possible after disasters.

For converter stations, resilience is typically defined as their
comprehensive ability to maintain operational functionality, re-
cover quickly, and restore stable power transmission under nat-
ural disasters. Enhancing the resilience of converter stations can
provide crucial support for the long-term stability of the power
networks [25]. However, existing research on converter station
resilience primarily focuses on resilience assessment method.
For instance, Liu et al. [26] proposed a multi-model probabilis-
tic framework to evaluate the resilience of converter stations,
but effective models and methods for resilience enhancement
remain lacking. A new probability-based method was proposed
to assess the seismic resilience of complex systems, such as
substations and converter stations [27]. Even though Wang et
al. [28] proposed a novel converter station structure to improve
AC and DC fault ride-through performance of the multiterminal
high-voltage direct current (HVDC) grid, effective methods for
enhancing converter station resilience are still lacking.

Additionally, natural disasters, as a typical type of HILP
event, are often fraught with uncertainties, such as the magni-
tude and timing of occurrence [29]. Most research efforts on the
resilience enhancement of power systems treat the pre-disrup-
tion preparation and post-disruption restoration independently.
For instance, system hardening strategies were made without
taking account of their impacts on the restoration process [30]-
[32]. The pre-disruption preparation strategies of power sys-
tems are, however, usually designed to prevent system failures
or mitigate disaster-induced losses, and therefore, influence the
effectiveness of post-disruption restoration strategies. Hence, it
necessitates the integration of the pre-disruption preparation
and post-disruption restoration for the effective allocation of
limited resources to enhance the resilience of power systems.
Nevertheless, only a few studies have explored resilience opti-
mization of power systems from the perspective of pre-disrup-
tion preparation and post-disruption restoration considering the
uncertainty of earthquakes [33]-[35]. Stochastic programming
has been implemented to enhance resilience of various types of
infrastructure networks [36], such as power systems, transpor-
tation networks, supply chains, and energy systems. Lagos et al.
[37] developed a resilience enhancement framework based on a
hierarchical approach, where the network investment and resil-
ience assessment were considered in the first and second layers,
<table><caption>TABLE I
COMPARISON OF FEATURES CONSIDERED IN RESILIENCE OPTIMIZATION
STUDIES</caption><thead><tr><th>Study</th><th>Optimization model</th><th>Both preparation and restoration</th><th>Disaster uncertainties</th><th>After-shocks</th><th>Endogenous uncertainty</th></tr></thead><tbody><tr><td>[19]-[24]</td><td>√</td><td></td><td>√</td><td></td><td></td></tr><tr><td>[17][26][27]</td><td></td><td>√</td><td>√</td><td></td><td></td></tr><tr><td>[31]-[38]</td><td>√</td><td>√</td><td>√</td><td></td><td></td></tr><tr><td>[44]</td><td></td><td></td><td>√</td><td>√</td><td></td></tr><tr><td>[41]</td><td>√</td><td></td><td>√</td><td></td><td>√</td></tr><tr><td>This study</td><td>√</td><td>√</td><td>√</td><td>√</td><td>√</td></tr></tbody></table>

respectively. Wang et al. [38] introduced a two-stage mixed-integer planning-based resilience optimization method to develop line hardening, distributed generator, and switches placement strategies. However, earthquakes, as a class of special natural disasters and HILP events, are often accompanied by after-shocks after the mainshock [29]. Aftershocks can not only cause secondary damage to the power facilities but also hinder the emergency restoration of power systems [34]. To the best of our knowledge, only a few resilience assessment studies accounted for the impacts of aftershocks [44], and all earthquake resilience enhancement methods have, however, overlooked the impacts of aftershocks.

Furthermore, existing studies on earthquake resilience enhancement methods primarily concentrates on exogenous uncertainty, namely, the uncertainty that is independent of decisions [39]. Put another way, the pre-disruption preparation strategies would not influence the probabilities of disaster damages. In reality, the decisions can also affect the uncertainty of damage degree, leading to endogenous or decision-dependent uncertainty [40]. For instance, Peeta et al. [41] formulated a two-stage stochastic model for the pre-disruption investment problems to enhance infrastructure resilience under endogenous uncertainty. In their work, the failure probabilities of links could be decreased by implementing protective actions. Likewise, the above studies only focused on the initial shock damage caused by disasters and failed to consider the impacts of secondary disasters or multiple disasters, say aftershocks. It was observed that the impacts of sequential disasters on the resilience of infrastructures or networks are much more complicated than a single disaster [29].

Table I summarizes the differences between existing literature and this study in terms of features considered in power system resilience optimization research. A checkmark (√) indicates that the feature is considered in the respective study, while the absence of a checkmark indicates otherwise. In summary, the presence of aftershocks and the plethora of uncertainties pose significant challenges to the development of a joint pre-earthquake preparation and post-earthquake restoration strategy. The endogenous uncertainties associated with the pre-earthquake preparation are also difficult to be addressed. To bridge this gap, this article put forth a three-level framework for the resilience optimization of converter stations, considering both pre-earthquake preparation and post-earthquake restoration under the uncertainties of mainshock-aftershock sequences. Based on the timing of earthquakes, the decision problems are divided into pre-earthquake preparation decision problem, and recovery

decision problems under mainshock and aftershocks, corresponding to the three levels of the framework. Differing from the recovery decision problem under mainshock, the recovery decision problem under aftershocks considers the impacts of secondary shocks which can hinder the recovery process. The equipment hardening strategy (EHS) is utilized to reduce the vulnerability of converter stations. The spare parts strategy (SPS) and recovery scheduling (RS) are utilized together to minimize the recovery time of converter stations. The primary contributions of this article are trifold: 1) Considering the effects of the mainshock-aftershock sequences in resilience enhancement from a joint pre-earthquake preparation and post-earthquake restoration perspective; 2) Developing a novel two-stage stochastic programming model for enhancing the resilience of converter stations, and employing a tolerated random number generation method to generate the deterministic second-stage scenarios that approximate the elimination of endogenous uncertainties in the model; 3) Proposing a combined approach with sample average approximation (SAA) algorithm and progressive hedging (PH) algorithm to efficiently solve the resulting optimization model.

The remainder of this article is rolled out as follows. The three-level modeling framework is elaborated in Section II. Section III details the proposed two-stage stochastic programming model, followed by the solution methodology in Section IV. Section V applies the model to a ±800 kV converter station. A summary of concluding remarks and future work are provided in Section VI.

## II. MODELING FRAMEWORK

Converter stations are critical components in UHV transmission and smart grids. Various electrical equipment, such as valve halls, converter transformers and DC filters, are placed in the converter stations, and they are highly susceptible to damage under earthquakes. In this article, an optimization framework is put forth to maximize the resilience of converter stations. The resilience of converter station is defined as the ability to resist earthquakes and quickly recover from earthquakes. Two pre-earthquake preparation strategies, i.e., EHS and SPS are considered in the framework. The EHS is utilized to reduce the converter station vulnerability, while the SPS is used to replace damaged equipment after the earthquake. The framework also considers the uncertainties associated with the mainshock-aftershocks, viz., magnitude, occurrence time, and occurrence intensity. The optimization process is structured into a three-level framework, as illustrated in Fig. 1. The first level is the pre-earthquake preparation model, the second level is the recovery model under the mainshock, and the third level is the recovery model under a sequence of aftershocks. The framework is detailed as follows.

### A. First Level

An effective and reasonable pre-earthquake preparation strategy is crucial for converter stations to withstand various earthquakes. In the first level, the pre-earthquake preparation model is given and the magnitudes and occurrence times of earthquakes are characterized via the earthquake prediction model.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Structure of the three-level framework.</figcaption></figure>

In this article, the earthquake prediction model is carried out
from two aspects of the mainshock and aftershocks. For
mainshock magnitude prediction, the well-known Gutenberg-
Richter law (G-R law) is employed. This law elucidates the re-
lationship between earthquake magnitude and occurrence inten-
sity in a specific area and time frame, and it is given by [42]:

$$
\lg N_M = a - bM . \tag{1}
$$

where $N_M$ is the number of earthquakes with magnitude greater than or equal to $M$. $a$ represents the earthquake activity level. A larger value of $a$ indicates a higher potential magnitude of the earthquake. $b$ is a parameter for the number of earthquakes, and it indicates that the number of earthquakes increases by a factor of $10^b$ if the magnitude decreases by one.

To analyze aftershocks, it is assumed that the magnitude and
occurrence time of aftershocks are independent. The aftershock
magnitude and occurrence time are described by triangular and
exponential distributions respectively [44]. The parameters of
triangular distribution $m_{min}$, $m_{mid}$, and $m_{max}$ are the minimal,
middle, and maximal magnitude of aftershocks, respectively. In
this article, it is assumed that $m_{min} = m - 1$, $m_{mid} = (2m - 1)/2$, and $m_{max} = m$, that is, the minimum magnitude of aftershock is
considered to be the mainshock magnitude minus 1 and the
maximum magnitude of aftershocks is considered to be the
same as the mainshock magnitude [44]. The occurrence inten-
sity of aftershock is assumed to be $\lambda$, then the probability den-
sity function $\lambda e^{-\lambda\tau}$ of exponential distribution is utilized to ob-
tain the aftershock occurrence time $\tau$ [45].

The maximum likelihood estimation (MLE) method is uti-
lized to estimate the parameters *a*, *b* and *λ* [43]. Combining
the mainshock magnitude simulation model with the G-R law,
the aftershock magnitude simulation model with triangular dis-
tributions, and the aftershock occurrence time simulation model
with exponential distributions, the mainshock-aftershock se-
quence can be simulated.

Subsequently, the vulnerability of converter stations in the
event of earthquakes needs to be assessed. The converter station
has numerous equipment, and all of them can be damaged by
the earthquakes. The vulnerability of converter stations is de-

fined as the failure probability of equipment under various mag-
nitudes of earthquakes, and can be determined using some fit-
ting models. In engineering practice, the lognormal cumulative
distribution is often utilized to assess the vulnerability of elec-
trical equipment subject to earthquakes [14]. The failure prob-
ability of equipment *i* under earthquake magnitude *m* is,
therefore, given by:

$$
p_i(m) = \Phi \left( \frac{\ln(PGA / \mu_i)}{\beta_i} \right) = \Phi \left( \frac{\ln(2^{1.29m-9.76} / \mu_i)}{\beta_i} \right). \quad (2)
$$

where *PGA* is the peak ground acceleration (PGA) of the earthquake and represents the degree of the impact of the earthquake on electrical equipment. The relationship between *PGA* and the earthquake magnitude *m* is expressed as
*PGA* = 2<sup>1.29*m*−9.76</sup>. *μ*<sub>i</sub> and *β*<sub>i</sub> are the mean value and standard deviation of the lognormal cumulative distribution of equipment *i*, respectively. Both coefficients can be estimated by MLE through history data [14]. Given the magnitude of any earthquake, the corresponding PGA and the failure probability of equipment can be calculated.

The EHS and SPS are adopted for the converter station to
mitigate the impact of earthquakes. The EHS can reduce the
failure probability of equipment under earthquake, and it is of-
tentimes modeled by increasing mean value of the vulnerability
of equipment subject to earthquakes. Therefore, the mean value
of the vulnerability for equipment *i* with hardening strategy is
*μ*<sub>*i*</sub>*k*<sub>*i*</sub> and *k*<sub>*i*</sub> is the hardening factor of equipment *i*. The SPS is
to prepare spare parts for different equipment in the converter
station, and the spare parts will be used to replace the damaged
equipment in the post-earthquake recovery process. In this arti-
cle, it is assumed that if the spare parts are not prepared for a
specific type of equipment, the damaged equipment cannot be
recovered after the earthquake. The spare parts can only be used
sequentially according to the RS and only one failure equipment
can be recovered at a time during the recovery process. The
EHS and SPS have distinct costs for different equipment, and
both of two strategies are subject to a budget.

B. Second Level

A highly resilient converter station is characterized not only
by its capability to resist the damage caused by earthquakes but
also by its ability to enact effective recovery actions quickly af-
ter earthquakes. In this article, it is assumed that the mainshock
is the most harmful disruptive event and the magnitudes of af-
tershocks are lower than that of the main earthquake. In this
level, the recovery model pertaining to the mainshock is de-
scribed as follows.

A typical ±800 kV converter station can be divided into two
circuits, pole I and pole II, and each pole circuit contains a high
and low valve hall, serving the core function of converting AC
to DC, and each valve hall is connected to multiple transformers.
According to the different voltage levels of converter trans-
formers (e.g., ±800 kV, ±600 kV, ±400 kV and ±200 kV) within
each valve hall, the system functioning state of the converter
station is defined as the functional state of the four bipolar
<figure><img src="image_5.png" /><figcaption>Fig. 2. The function logic diagram model of the converter station.</figcaption></figure>

branch circuits. Each branch circuit contains a voltage level of converter transformers and valve halls, and the functional state of each branch circuit is determined by the functional state of different equipment. Specifically, each branch contributes 25% to the total transmission capacity of the converter station. In this article, the maximum power output from the converter station after earthquakes is selected as the index to assess the station's performance. Therefore, the function of the converter station $F(t)$ is divided into the following five transmission capacities, namely 100%, 75%, 50%, 25%, 0%.

Based on the composition of the converter station, the function logic diagram model of the converter station is established to analyze the performance of the converter station, which is shown in Fig. 2. The converter station is composed of one power distribution equipment (ACDE), ten station transformers (TF), two filter fields (FLTR), three relay rooms (RM), four control buildings (CB), four valve halls (VH), four bypass switches (BS), four smoothing reactors (SR), two damping filters (DF), two DC filters (DCF), two earthing switches (ES), and one DC output equipment (DCOE). The different color squares indicate different equipment, of which states could be either 0 or 1, representing normal working and disrupted by earthquakes respectively. The overall system function of the converter station, in the event of an earthquake, can be effectively determined and analyzed using this logic diagram model. It is worth noting that while the functionality of incoming and outgoing feeders can influence the system function of the converter station, this study does not include incoming and outgoing feeders as part of the research scope. This is because the focus of this study is on the converter station as a single-node unit within the power system, rather than on the entire power network. Regardless of whether the external incoming or outgoing feeders are damaged, the converter station should restore the damaged equipment within the station as soon as possible to ensure normal power transmission capability. Typically, the

post-earthquake recovery task of external feeders is not handled by the personnel of the converter station and does not require the allocation of the converter station's recovery resources. In other words, the research object is the converter station under known states of incoming and outgoing feeders

Based on the definition of $F(t)$, the resilience of converter station can be assessed based on the system's ability to withstand, respond, and adapt to earthquake damage. Therefore, the resilience of converter station after the earthquake $\xi$, denoted as $R(\xi)$, is described as the area surrounded by $F(t)$ and $t$ in time period $T$, and it is formulated as:

$$R(\xi) = \frac{\int_{T_1}^{T} F(t) dt}{\Delta t} = \frac{\sum_{t=T_1}^{T} \Delta t F(t)}{T} \quad (3)$$

where $\Delta t$ represents the duration of the time step. $T_1^\xi$ denotes the occurrence time of the mainshock.

By combining the vulnerability curve of the converter station and the mainshock magnitude model, the equipment damage condition of the converter station after the mainshock can be simulated. The recovery model under the mainshock to maximize the resilience of the converter station is established to develop the RS. In this article, the RS is defined as the sequential recovery order for the damaged equipment.

### C. Third Level

In the third level, the recovery model under aftershocks is constructed, and it contains the following two cases: If the equipment has not been recovered at the time of an aftershock, the ongoing recovery process is halted, and the disrupted equipment should be recovered again after the aftershock. If the equipment has already been recovered when an aftershock occurs, the condition of equipment should be re-assessed based on the aftershock magnitude and vulnerability curve. Additionally, the EHS is viewed as failed if the equipment is damaged. After the occurrence of any aftershock, the recovery optimization model will be repeated.

The recovery process of the converter station after the mainshock is dynamically adjusted in response to multiple aftershocks. This article assumes that the time limit of post-earthquake recovery process is 60 days. Specifically, the post-earthquake stage in this article concentrates solely on the recovery process of the converter station within 60 days, relying solely on existing on-site resources, such as personnel and spare parts [46]. The time limit is set because the equipment in converter station are typically large, heavy, and delicate, making it challenging to be manufactured and transported in a short period. After reaching the time limit, the recovery process ceases, and further repairs would require external resources, which are beyond the scope of this study. Moreover, the 60-day limit can be customized based on the specific situation of different converter stations. Note that the equipment damage condition of the converter station is not only affected by the earthquakes, but also influenced by the EHS which should be implemented prior to earthquakes. This phenomenon underscores the challenge of
formulating effective EHS strategies before earthquakes, and it
is referred to the endogenous uncertainties in this article.

In summary, the uncertainties of mainshock-aftershock se-
quences and the endogenous uncertainties within the optimiza-
tion model pose great challenges in developing a joint pre-
earthquake preparation (EHS and SPS) and post-earthquake
restoration (RS) strategy. To this end, a two-stage stochastic
programming model and a tolerated scenario generation method
are proposed to deal with the endogenous uncertainties in the
following section.

### III. PROPOSED STOCHASTIC PROGRAMMING MODEL

To address the above three-level problem, a two-stage sto-
chastic programming model is put forth. The pre-earthquake
preparation model constitutes the first-stage problem, while the
recovery models for the mainshock and aftershocks are inte-
grated into the second-stage problem. Hereafter, the
mainshock-aftershock sequences are referred to as earthquakes
and the k th shock represents the k th earthquake. The pro-
posed model can be formulated as follows:

**First-stage:**

$$
z = \max \mathbb{E}[R(\xi)]. \quad (4)
$$

$$
\sum_{i \in \mathbf{N}} n_i C_i^h + \sum_{i \in \mathbf{N}} r_i C_i^s \leq C. \quad (5)
$$

$$
n_i \in \{0,1\}, r_i \text{ is integer, } \forall i \in \mathbf{N}. \quad (6)
$$

The first-stage problem focuses on maximizing the expecta-
tion of converter station resilience under various scenarios,
where $\mathbb{E}[R(\xi)] = \sum_{\xi \in \mathbb{E}} R(\xi) / |\mathbb{E}|$, and $|\mathbb{E}|$ represents the num-
ber of earthquake scenarios. Constraint (5) specifies that the to-
tal budget allocated to equipment hardening and spare parts is
limited. Constraint (6) demonstrates that the decision variables
consist of binary or integer variables.

**Second-stage:**

$$
\max R(\xi) = \int_{T_1^{\xi}}^{T_2^{\xi}+T} F^{\xi}(t)/T dt = \sum_{t=T_1^{\xi}}^{T_2^{\xi}+T} \Delta t F^{\xi}(t)/T . \quad (7)
$$

$$
\sum_{\forall j \in \mathbf{N}:(i,j) \in \mathbf{L}} x_{ij}^{\xi}(t) = F^{\xi}(t), \forall i \in \mathcal{D}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E}. \quad (8)
$$

$$
\sum_{\forall j \in \mathbf{N}:(v,j) \in \mathbf{L}} x_{jv}^{\xi}(t) = \sum_{\forall j \in \mathbf{N}:(v,j) \in \mathbf{L}} x_{vj}^{\xi}(t), \quad (9)
$$

$$
\forall v \in \mathbf{N} / \{O, D\}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E}
$$

$$
0 \le x_{ij}^{\xi}(t) \le c_{ij}, \forall (i, j) \in \mathbf{L}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E}. \quad (10)
$$

$$
\sum_{i \in \mathbf{N}} m_i^{\xi}(t) \le 1, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E}. \quad (11)
$$

$$
\sum_{t \in \mathbf{T}} m_i^{\xi}(t) / w_i \le r_i, \forall i \in \mathbf{N}, \forall \xi \in \mathbf{E}. \quad (12)
$$

$$
\begin{align}
R_i^{T_i^\xi} - (p_i^{T_i^\xi} - \delta_i^{T_i^\xi} n_i^{T_i^\xi}) &\le \tau_i^{T_i^\xi} \le 1 + R_i^{T_i^\xi} - (p_i^{T_i^\xi} - \delta_i^{T_i^\xi} n_i^{T_i^\xi}), \quad (13) \\
\forall i \in \mathbf{N}, \forall k \in \mathbf{E}^\xi, \forall \xi \in \mathbf{E}^\xi & \nonumber
\end{align}
$$

$$
n_i^{T_i^\xi} = \max \{0, n_i^{T_i^\xi} + \tau_i^{T_i^\xi} - 1\}, \forall i \in \mathbf{N}, \forall k \in \mathbf{E}^\xi, \forall \xi \in \mathbf{E}. \quad (14)
$$

$$
n_i^{T_i^\xi} = n_i, \forall i \in \mathbf{N}, \forall \xi \in \mathbf{E}. \quad (15)
$$

$$
g_i^{\xi}(t+1) = \begin{cases}
g_i^{\xi}(t) + m_i^{\xi}(t)/w_i, & \forall i \in \mathbf{N}, \forall t \in \mathbf{T}/T_k^{\xi}, \forall \xi \in \mathbf{E} \\
\min\{g_i^{\xi}(t) + m_i^{\xi}(t)/w_i, \tau_i^{T_k^{\xi}}\}, & \forall t \in T_k^{\xi}, \forall i \in \mathbf{N}, \forall \xi \in \mathbf{E}
\end{cases} \quad (16)
$$

$$
\begin{gather}
s_i^{\xi}(t) = \max \left\{ s_i^{\xi}(t-1), \mathbf{1}_{g_i^{\xi}(t-1) \ge 1} \right\}, \tag{17} \\
\forall i \in \mathbf{N}, \forall t \in \mathbf{T} / T_k^{\xi}, \forall \xi \in \mathbf{E}. \nonumber
\end{gather}
$$

$$
\begin{equation}
\begin{split}
s_i^{\xi}(T_k^{\xi} + 1) = \min \left\{ s_i^{\xi}(T_k^{\xi}) + g_i^{\xi}(T_k^{\xi}), \tau_i^{T_k^{\xi}} \right\}, & \\
\forall i \in \mathbf{N}, \forall k \in \mathbf{E}^{\xi}, \forall \xi \in \mathbf{E}
\end{split}
\tag{18}
\end{equation}
$$

$$
g_i^{\xi}(T_1^{\xi}) = s_i^{\xi}(T_1^{\xi}) = \tau_i^{T_1^{\xi}}, \forall i \in \mathbf{N}, \forall \xi \in \mathbf{E}. \quad (19)
$$

$$
\begin{equation}
\begin{aligned}
& 1 \ge s_i^{\xi}(t+1) + m_i^{\xi}(t+1) \ge m_i^{\xi}(t), \\
& \quad \forall i \in \mathbf{N}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E}.
\end{aligned}
\tag{20}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
x_{ij}^{\xi}(t) \le M \min \{s_i(t), s_j(t)\}, & \\
\forall i \in \mathbf{N}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E} &
\end{split}
\tag{21}
\end{equation}
$$

$$
\left\{
\begin{array}{@{}l@{}}
0 \le g_i^{\xi}(t) \le 1, \forall i \in \mathbf{N}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E} \\
s_i^{\xi}(t), m_i^{\xi}(t) \in \{0, 1\}, \forall i \in \mathbf{N}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E}. \\
x_{ij}^{\xi}(t) \ge 0, \forall (i, j) \in \mathbf{L}, \forall t \in \mathbf{T}, \forall \xi \in \mathbf{E}
\end{array}
\right.
\tag{22}
$$

The objective of the second-stage optimization is to maxim-
ize the resilience of converter stations under each earthquake
scenario $\xi$, where constraint (8) defines the formula for as-
sessing converter station performance at a given time $t$. Con-
straints (9) and (10) state that the flow through line $(i, j)$ is lim-
ited, and the inflow towards the equipment $i$ is equivalent to
its outflow. According to constraints (11) and (12), only one
equipment can be recovered at a time, and the recovery times
for equipment $i$ must be not greater than the available number
of spare parts. Constraint (13) defines the state of equipment $i$
following the $k$ th shock in earthquake scenario $\xi$, employing
a tolerated random number generation method to address en-
dogenous uncertainties. Specifically, $(p_i^{T_i^\xi} - \delta_i^{T_i^\xi} n_i^{T_i^\xi})$ repre-
sents the failure probability of equipment $i$ under the $k$ th
shock in earthquake scenario $\xi$, where $n_i^{T_i^\xi}$ is determined by
constraints (14) and (15). If the generated random number $R_i^{T_i^\xi}$
is less than the failure probability $(p_i^{T_i^\xi} - \delta_i^{T_i^\xi} n_i^{T_i^\xi})$, the equip-
ment $i$ is disrupted and $\tau_i^{T_i^\xi} = 0$, otherwise, $\tau_i^{T_i^\xi} = 1$. $g_i^{\xi}(t)$
serves as an instrumental variable which is defined by con-
straint (16) and represents the recovery schedule of equipment
$i$. Once an aftershock occurs during the recovery process, the
recovery schedule will be reset. Constraint (17) stipulates that
$s_i^{\xi}(t)$ equates to 1 only if $g_i^{\xi}(t)=1$. The recursive definition of
each equipment state at time $t$ is governed by constraints (17)
and (18), where $\mathbf{1}_{g_i^{\xi}(t-1) \ge 1}$ is equal to 1 only if $g_i^{\xi}(t-1)$ is true,
and 0, otherwise. To sum up, the state of the equipment is determined by its disruption status, its prior state, and the recovery schedule. The initial conditions of $g_i^\xi(t)$ and $s_i^\xi(t)$ are detailed in constraint (19), and $T_{1,i}^\xi$ is equal to 0 in the optimization model. Put another way, the mainshock is assumed to occur at the outset of the second-stage problem. Constraint (20) indicates that recovery of other equipment cannot commence until one equipment has been fully recovered. Constraint (21) specifies that the flow of line (i, j) is restricted by the states of equipment i and j. The decision variables for the second-stage problem are constrained to be nonnegative, as stipulated in constraint (22). To convert the optimization model into an integer programming model, two linearization techniques for the max/min and signum functions are introduced. Constraints incorporating the max/min function, e.g., constraints (14), (16), (17), (18), and (21), can be linearized through the introduction of new binary variables. For instance, through the introduction of new binary variables $u_{1,i}^{T_{1,i}^\xi}$ and $v_{1,i}^{T_{1,i}^\xi}$, constraint (14) can be transformed into constraint (23), defined as following:

$$ egin{cases} n_i^{T_{k+1}^\xi} \ge 0, n_i^{T_{k+1}^\xi} \ge n_i^{T_k^\xi} + \tau_i^{T_k^\xi} - 1, n_i^{T_{k+1}^\xi} \le 0 + M(1 - u_{1,i}^{T_{k+1}^\xi}) \\ n_i^{T_{k+1}^\xi} \le n_i^{T_k^\xi} + \tau_i^{T_k^\xi} - 1 + M(1 - v_{1,i}^{T_{k+1}^\xi}), u_{1,i}^{T_{k+1}^\xi} + v_{1,i}^{T_{k+1}^\xi} \ge 1 \quad \forall i \in \mathbf{N}, \forall k \in \mathbf{E}^\xi/\{1\}, \forall \xi \in \mathbf{E} \end{cases} \quad (23) $$

The signum function $1_{g_i^\xi(t-1) \ge 1}$ is equated to the new instrumental variables $sig_i^\xi(t-1)$ via the addition of the new constraint (24), defined as:

$$ \begin{cases} g_i^\xi(t-1) \ge sig_i^\xi(t-1) \ge g_i^\xi(t-1) - 1 + 1/M \\ \forall i \in \mathbf{N}, \forall t \in \mathbf{T} / T_k^\xi, \forall \xi \in \mathbf{E} \end{cases} \quad (24) $$

Linearization techniques are utilized to transform the proposed two-stage stochastic programming into an integer linear programming model.

## IV. SOLUTION METHODOLOGY

To address the proposed two-stage stochastic programming model, it is necessary to generate a finite number of deterministic second-stage scenarios. In the article, a deterministic post-earthquake scenario is characterized as a damage series of the converter station, comprising both a failure equipment sequence and an occurrence time sequence of earthquake. The uncertainties associated with earthquake scenarios, such as magnitude, occurrence time, and occurrence intensity, pose challenges for the determination of the second-stage scenario. SAA [47] is employed to generate the second stage scenarios. Once the second-stage scenarios are determined, the proposed two-stage stochastic programming model becomes essentially a linear integer programming model. However, the presence of a large number of second-stage scenarios results in numerous variables and constraints in the proposed model, thereby hindering the efficiency of traditional solvers, such as CPLEX and LINGO. To this end, the PH algorithm [48] is utilized to effectively solve the model with a finite number of deterministic second-stage scenarios.

### A. Tolerated Scenario Generation

The second-stage scenarios essentially represent the disruption scenarios of converter station equipment, determined by earthquake scenarios and the EHS. Earthquake scenarios, comprising mainshock-aftershock sequences, can be generated using the earthquake prediction model introduced in Section II. Subsequently, the failure probabilities of equipment, both with and without the hardening strategy, can be determined under these sequences. Consequently, the second-stage scenarios in the stochastic programming model are converted into the random sets of all equipment under mainshock-aftershock sequences, and the endogenous uncertainties in the model can be approximately eliminated. The details are given in Eq. (13). Thereby, a finite number of deterministic second-stage scenarios is generated using the SAA algorithm.

The fundamental concept of the SAA algorithm involves generating a small number of second-stage scenarios and calculating the optimal solution based on these scenarios. First, several scenarios $N$ are generated, and the objective function of the original optimization model under $N$ scenarios can be approximated as:

$$ \min \left\{ \nu_N = 1 - \frac{1}{N} \sum_{n=1}^{N} R(\zeta_n) \right\}. \quad (25) $$

The optimal solution $X_N$ and objective function $\nu_N$ will converge to the values of the original problem with the increase of $N$. In this article, $X_N = (n_N, r_N)$ where $n_N$ and $r_N$ represent the set of the EHS and SPS in the first stage, respectively. The lower bound of the SAA is objective function calculated using $X_N$ across a large number of two-stage scenarios, and the upper bound is the maximum objective function obtained by repeatedly solving the optimization model $M$ times when the number of scenarios is $N$. $M$ represents the number of times the optimization model is repeatedly calculated, typically taking a value between 5 and 20 in engineering applications. $N$ represents the number of generated second-stage scenarios, with its value usually determined based on the problem's optimization performance. The gap of the solution is the difference between the upper bound and the lower bound. In practice, the number of scenarios $N$ will directly affect the gap of the solution and the computational efficiency. In this article, the sensitivity analysis of $N$ is presented in Section V. SAA achieves a set of optimal solutions by repeatedly sampling second-stage scenarios and subsequently selects the highest quality solution as the final solution.

### B. Progressive Hedging Algorithm

As the number of scenarios increases, both the scale and computational demand of the two-stage stochastic programming models grow. Consequently, the PH algorithm is employed to enhance computational efficiency. The method addresses a scenario-wise decomposed version of the problem, effectively managing the computational challenges posed by large-scale issues. The PH algorithm is a well-known scenario-
based decomposition method in stochastic programming framework proposed by Rockafellar and Wets [48]. The basic idea of PH algorithm is similar to the exterior point method, which solves the problem by adding a penalty term to the objective function and iterating continuously [48].

In this article, let $(\mathbf{n}, \mathbf{r})$ define the first-stage decision variables which contains the EHS and SPS, $p(\xi)$ is the occurrence probability of scenario $\xi$ and equal to $1/|\mathbf{E}|$, $z_{\xi}^k$ is the objective function for the $k$ th iteration under scenario $\xi$ and shown in Eq. (26), $\rho$ is the algorithm control parameters, the PH algorithm is elaborated below:

**Step 1:** Initialize the problem and set the number of iterations $k=0$.

**Step 2:** For each scenario $\xi \in \mathbf{E}$, solve the following optimization problem:

$$z_{\xi}^{k} = \min[1 - R(\xi)]. \quad (26)$$

**Step 3:** The optimal solution of the first-stage problem is obtained as $(\mathbf{n}_{\xi}^k, \mathbf{r}_{\xi}^k)$ in the $k$ th iteration of scenario $\xi$, and the mean value of the first-stage decision variables for each scenario is calculated as $(\bar{\mathbf{n}}^k, \bar{\mathbf{r}}^k) = \sum_{\xi \in \mathbf{E}} p(\xi)(\mathbf{n}_{\xi}^k, \mathbf{r}_{\xi}^k)$.

**Step 4:** The penalty term of the objective function is calculated as: $\mathbf{w}_{\xi}^k = \rho[(\mathbf{n}_{\xi}^k, \mathbf{r}_{\xi}^k) - (\bar{\mathbf{n}}^k, \bar{\mathbf{r}}^k)]$, $\forall \xi \in \mathbf{E}$.

**Step 5:** For each scenario $\xi \in \mathbf{E}$, solving the following augmented generalized Lagrangian optimization problem and set $k=k+1$.

$$z_{\xi}^{k} = \min \left\{ \begin{array}{l} 1 - R(\xi) + \left[ \mathbf{w}_{\xi}^{(k-1)} \right]^{\mathrm{T}} (\mathbf{n}^{k}, \mathbf{r}^{k}) + \\ \frac{\rho}{2} \left\| (\mathbf{n}^{k}, \mathbf{r}^{k}) - \left[ \bar{\mathbf{n}}^{(k-1)}, \bar{\mathbf{r}}^{(k-1)} \right] \right\|^{2} \end{array} \right\}. \quad (27)$$

where $\| (\mathbf{n}^k, \mathbf{r}^k) - [\bar{\mathbf{n}}^{(k-1)}, \bar{\mathbf{r}}^{(k-1)}] \|^2$ represents the Euclidean distance between $(\mathbf{n}^k, \mathbf{r}^k)$ and $[\bar{\mathbf{n}}^{(k-1)}, \bar{\mathbf{r}}^{(k-1)}]$.

**Step 6:** Updating the mean values of the first-stage decision variables $(\bar{\mathbf{n}}^k, \bar{\mathbf{r}}^k) = \sum_{\xi \in \mathbf{E}} p(\xi)(\mathbf{n}_{\xi}^k, \mathbf{r}_{\xi}^k)$.

**Step 7:** Updating the penalty term under each scenario as $\mathbf{w}_{\xi}^k = \mathbf{w}_{\xi}^{(k-1)} + \rho[(\mathbf{n}_{\xi}^k, \mathbf{r}_{\xi}^k) - (\bar{\mathbf{n}}^k, \bar{\mathbf{r}}^k)]$, $\forall \xi \in \mathbf{E}$.

**Step 8:** Calculating the bias weighted values for each scenario as $g^k = \sum_{\xi \in \mathbf{E}} p(\xi) \| (\mathbf{n}_{\xi}^k, \mathbf{r}_{\xi}^k) - (\bar{\mathbf{n}}^k, \bar{\mathbf{r}}^k) \|^2$.

**Step 9:** If $g^k < \varepsilon$, the optimal solution is reached; otherwise, return to Step 5.

The PH algorithm is capable of reducing the size of the problem when solving stochastic programming problems. As the iterative process depends on the penalty term, the ratio of the penalty term to the original objective function's magnitude dictates the convergence speed. It is necessary to reasonably control the size of the penalty term $\{[\mathbf{w}_{\xi}^{(k-1)}]^{\mathrm{T}}(\mathbf{n}^k, \mathbf{r}^k) + \frac{\rho}{2}\| (\mathbf{n}^k, \mathbf{r}^k) - [\bar{\mathbf{n}}^{(k-1)}, \bar{\mathbf{r}}^{(k-1)}] \|^2\}$, i.e., the size of the algorithm parameters $\rho$. The size of the penalty term $\rho$ is usually set as an empirical constant. In this article, $\rho$ is set to vary with the number of iterations $k$, i.e.,

$$\rho_{k+1} = \gamma \rho_k. \quad (28)$$

where $\gamma$ is an iterative multiplier, and one has $1 \le \gamma \le 2$. The

<figure><img src="image_6.png" /><figcaption>Fig. 3. The basic information of the ±800 kV converter station located in Yibin, Sichuan Province, China.</figcaption></figure>

<table><caption>TABLE II
THE PARAMETERS OF DIFFERENT EQUIPMENT</caption><thead><tr><th>Equipment</th><th>μ<sub>i</sub></th><th>β<sub>i</sub></th><th>C<sup>h</sup><sub>i</sub> (10<sup>6</sup> dollars)</th><th>C<sup>s</sup><sub>i</sub> (10<sup>6</sup> dollars)</th><th>w<sub>i</sub> (days)</th><th>k<sub>i</sub></th></tr></thead><tbody><tr><td>ACDE</td><td>0.7</td><td>0.3</td><td>1</td><td>3</td><td>8</td><td>0.65</td></tr><tr><td>TF</td><td>0.45</td><td>0.5</td><td>0.5</td><td>1</td><td>4</td><td>0.65</td></tr><tr><td>RM</td><td>0.8</td><td>0.3</td><td>1</td><td>—</td><td>8</td><td>0.65</td></tr><tr><td>CB</td><td>0.8</td><td>0.3</td><td>1</td><td>—</td><td>8</td><td>0.65</td></tr><tr><td>FLTR</td><td>0.52</td><td>0.5</td><td>1</td><td>1.5</td><td>8</td><td>0.65</td></tr><tr><td>VH</td><td>0.5</td><td>0.5</td><td>2</td><td>3</td><td>8</td><td>0.65</td></tr><tr><td>BS</td><td>0.6</td><td>0.3</td><td>—</td><td>0.5</td><td>4</td><td>—</td></tr><tr><td>SR</td><td>0.5</td><td>0.5</td><td>0.6</td><td>1</td><td>4</td><td>0.65</td></tr><tr><td>DF</td><td>0.55</td><td>0.5</td><td>0.6</td><td>1</td><td>4</td><td>0.65</td></tr><tr><td>DCF</td><td>0.5</td><td>0.5</td><td>0.6</td><td>1</td><td>4</td><td>0.65</td></tr><tr><td>ES</td><td>0.7</td><td>0.3</td><td>—</td><td>0.5</td><td>8</td><td>—</td></tr><tr><td>DCOE</td><td>0.7</td><td>0.3</td><td>1</td><td>3</td><td>8</td><td>0.65</td></tr></tbody></table>

initial value of $\rho$ is determined as $\rho_0 = 0.005 \sim 0.3$.

## V. CASE STUDY

### A. Experimental Setting

The proposed method was implemented to a ±800 kV converter station located in the city of Yibin, Sichuan Province, P.R. China and the converter station is close to the major earthquake zone. The geographic information of the converter station is shown in Fig. 3, which is close to the major earthquake zone. Figure 3 also displays the function logic diagram model of the converter station and a segment of the historical earthquake data pertaining to the converter station. The detailed functional logic diagram of the converter station is depicted in Fig. 2. It is noted that the converter station has four inputs and four outputs. Based on data from the China Earthquake Networks Center (CENC), the historical earthquake records from the past ten years for the geographic area corresponding to the latitude and longitude of the converter station are collected [49]. The data underwent initial cleaning and processing, including extracting earthquake occurrence times and magnitudes. The least squares and MLE methods are utilized to estimate the earthquake activity level $a$, the number of earthquakes $b$, and aftershock occurrence intensity $\lambda$. They are set to 19, 3, and 0.025 day-1, respectively. The mean values $\mu_i$ and standard deviations $\beta_i$ of equipment vulnerability curve, the cost of EHS $C_i^h$, the cost of SPS $C_i^s$, the equipment hardening factor $k_i$ and the equipment recovery time $w_i$ for different equipment are tabulated in Table II. The computation process is implemented in the Python API of Gurobi 10.0.1 and performed on a desktop PC with 3.00 GHz CPU and 24 GB RAM.
<table><caption>TABLE III
THE SENSITIVE ANALYSIS OF THE NUMBER OF SCENARIOS</caption><thead><tr><th>Number of scenarios <i>N</i></th><th>5</th><th>10</th><th>20</th><th>30</th><th>50</th></tr></thead><tbody><tr><th>Number of calculations <i>M</i></th><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><th>Lower bound of the solution</th><td>0.374</td><td>0.428</td><td>0.457</td><td>0.458</td><td>0.481</td></tr><tr><th>Upper bound of the solution</th><td>0.653</td><td>0.602</td><td>0.599</td><td>0.582</td><td>0.600</td></tr><tr><th>Gap of the solution</th><td>0.279</td><td>0.174</td><td>0.142</td><td>0.124</td><td>0.119</td></tr><tr><th>Computation time (hours)</th><td>0.494</td><td>1.463</td><td>3.518</td><td>5.162</td><td>6.144</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 4. The results of SAA.</figcaption></figure>

## B. Pre-earthquake Preparation Strategy

The time limit $T$ for post-earthquake recovery process is set to be 60 days and the total budget $C$ is 15 million dollars. Thirty mainshock-aftershock sequences are generated by the earthquake simulation model. The tolerated scenario generation method is employed to create 30 sets of random numbers for all equipment, serving as scenarios in the two-stage stochastic programming model. Subsequently, the PH algorithm is utilized to efficiently solve the optimization model. To enhance the solution quality, the SAA method is implemented to calculate the model with ten times and the results are depicted in Fig. 4. The upper bound of each calculation denotes the objective function derived from the optimization model, whereas the lower bound reflects the resilience expectation, determined by applying the optimal solution in ten thousand simulation scenarios. The gap of the solution is the difference between the upper bound and the lower bound. The solution exhibiting the greatest lower bound is selected as the optimal solution. The gap of the optimal solution is 0.124 and the variance of the gap is 0.0015. Additionally, the sensitivity analysis regarding the number of scenarios, i.e., $N$, is presented in Table III. It should be noted that the increment of the number of scenarios results in a smaller gap in the solution, indicating a higher quality of the solution. Elevating the number of scenarios from 30 to 50 results in a 4% reduction in the solution gap, yet it escalates the computational time by 19%. Consequently, 30 is chosen as the number of scenarios to keep a balance between computational burden and quality of solution.

Figure 5 illustrates the optimal EHS and SPS, with the resilience expectation of the converter station is 0.582. Red and white colored blocks represent “hardening” and “no hardening” of the equipment, respectively, and the numbers of spare parts are displayed on the right side. It is noted that equipment such as TF, VH, FLTR, SR, DF, and DCF exhibit great vulnerability, and the EHS tends to harden the equipment of the same circuit

to ensure the low transmission capacity of the converter station against various earthquakes. It is evident that the SPS prioritizes to stock spare parts for equipment that has not been hardened, such as BS.

## C. Experimental Comparison

Three additional cases under the same experimental settings: no EHS (Case II), no SPS (Case III), and neither EHS nor SPS (Case IV), are provided to examine the effectiveness of the proposed pre-earthquake preparation strategies (Case I). Figure 6 illustrates the resilience expectation of the converter station under different pre-earthquake preparation strategies. When no pre-earthquake preparation strategies are taken, the converter station can hardly cope with the impact of earthquakes, leading to a poor resilience expectation. However, the resilience expectation of the converter station can be significantly improved if an additional EHS or SPS is taken. It indicates that the implementation of the EHS and SPS contributes to achieving a higher resilience expectation (improved by 0.409, 0.369, and 0.127 compared with the results from Cases I, II, and III, respectively). It also reveals that the SPS is more effective than the EHS, because the lack of the SPS leads to a compromised recoverability of the converter station. A modest increase in the SPS budget allocation can yield significant improvements in the resilience of the converter station.

<figure><img src="image_8.png" /><figcaption>Fig. 5. The results of the optimal EHS and SPS.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Resilience expectations in different cases.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Resilience and budget proportions with increasing total budget.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Resilience and budget proportions of different experiments.</figcaption></figure>

Moreover, existing methods often adopt resilience enhance-
ment strategies based on Monte Carlo-based resilience assess-
ment method [50],[51]. The general process of these methods
involves generating vulnerability curves for converter station
equipment and earthquake simulation scenarios, simulating
post-earthquake damage and recovery process of the converter
station, and producing resilience assessment results. Through
sensitivity analysis or similar techniques, these methods iden-
tify earthquake vulnerable equipment in the converter station
and allocate budgets to the equipment to improve resilience.
When Monte Carlo-based resilience assessment approach is
used, the resilience assessment of the converter station is deter-
mined to be 0.104. The identified top four earthquake vulnera-
ble equipment are FLTR, VH, TF, and SR. By evenly distrib-
uting a budget of 15 million dollars across the four equipment,
the resulting resilience increased to 0.283. However, the result
is significantly lower than 0.554 of the proposed method in this
article. The comparative results demonstrate that the proposed
method can yield a better solution than that of the Monte Carlo-
based resilience assessment method in terms of improving con-
verter station resilience.

Figure 7 depicts the variation of the resilience expectation and budget proportion when the total budget increases in Ex- periment I. It is observed that the resilience expectation of the converter station increases as the total budget increases, while an excessive budget no longer yields further improvements. It offers insights to decision-makers to determine an appropriate budget. Additionally, the budget allocation for the SPS tends to increase alongside the total budget. This is because that there is

<table><caption>TABLE IV
THE PARAMETERS OF THE EARTHQUAKES IN DIFFERENT EXPERIMENTS</caption>
   <thead>
    <tr>
     <td>
      Experiment ID
     </td>
     <td>
      λ (day
      <sup>
       -1
      </sup>
      )
     </td>
     <td>
      Mean value of mainshock PGA
     </td>
    </tr>
   </thead>
   <tbody>
    <tr>
     <td>
      Experiment I
     </td>
     <td>
      0.025
     </td>
     <td>
      0.365
     </td>
    </tr>
    <tr>
     <td>
      Experiment II
     </td>
     <td>
      0.01
     </td>
     <td>
      0.315
     </td>
    </tr>
    <tr>
     <td>
      Experiment III
     </td>
     <td>
      0.025
     </td>
     <td>
      0.365
     </td>
    </tr>
    <tr>
     <td>
      Experiment IV
     </td>
     <td>
      0.01
     </td>
     <td>
      0.315
     </td>
    </tr>
   </tbody>
  </table>

<table><caption>TABLE V
THE COMPARISON OF MODEL I AND MODEL II</caption>
   <thead>
    <tr>
     <td>
      Index
     </td>
     <td>
      Model I
     </td>
     <td>
      Model II
     </td>
    </tr>
   </thead>
   <tbody>
    <tr>
     <td>
      Resilience expectation
     </td>
     <td>
      0.683
     </td>
     <td>
      0.745
     </td>
    </tr>
    <tr>
     <td>
      Resilience expectation after the implementation of SPS and EHS
     </td>
     <td>
      0.646
     </td>
     <td>
      0.452
     </td>
    </tr>
   </tbody>
  </table>

still a probability that the equipment will be damaged by earth-
quakes, so that the effectiveness of the EHS is limited. In other
words, if the budget is adequate, a substantial SPS investment
with a suitable EHS investment is suitable for supporting the
converter station resilience enhancement. Otherwise, moving
the focus to the EHS may be an efficient alternative.

D. *Analysis of Aftershock Impact*

By setting different parameters of *a*, *b*, and *λ*, four experiments are undertaken to further evaluate the effectiveness of the optimization model and analyze the impact of aftershocks. The mean values of the mainshock PGA are generated by the earthquake simulation model and the parameters in different experiments are shown in Table IV. Within this study, a mean value of the mainshock PGA of 0.365g is classified as high, while 0.315g is low. Similarly, an aftershock occurrence intensity of 0.025/ day is classified as high, and 0.01/day as low. Figure 8 illustrates the resilience expectations and budget proportions of different experiments. It indicates that: (1) the mean values of the mainshock PGA and aftershock occurrence intensity have huge impacts on the resilience expectation of the converter station. Consequently, the impact of aftershocks cannot be overlooked; (2) the results from Experiments I and II reveal that in scenarios characterized by increased aftershock occurrence intensity and high value of the mainshock PGA, directing additional budget towards the EHS proves more beneficial than allocating them to the SPS. Conversely, Experiments III and IV suggest that in scenarios with increased aftershock occurrence intensity and low value of the mainshock PGA, prioritizing budget allocation towards the SPS is more advantageous than that towards the EHS.

Furthermore, the impact of aftershocks is examined by comparing the proposed model (Model I) with the model without aftershocks (Model II), as detailed in Appendix A. The results, as given in Tables V and VI, show that, with a total budget limit C of 20 million dollars, solving Model II yields a resilience expectation value of 0.745. This value suggests an overestimation comparing with the result of Model I. Furthermore, upon the implementation of the SPS and EHS derived from Models I and II, the resilience expectation values under a large number N' of the second-stage scenarios of Models I and II are 0.646 and 0.452, respectively. Apparently, the SPS and EHS obtained by Model II are less effective than those by Model I, underscoring the importance of taking account of aftershocks. Similar to
the findings in Experiment II, the SPS is allocated a greater budget proportion than the EHS. Additionally, as evident from the data presented in Table VI, for equipment with higher vulnerability and repair costs, such as VH, the model without aftershocks leans towards the SPS, while the model with aftershocks gives more weights on the EHS. This analysis offers decision-makers a novel perspective on the budget distribution.

## VI. CONCLUSION AND FUTURE WORKS

This article aimed at maximizing the converter station resilience under the mainshock-aftershock sequence. A three-level modeling framework was developed to deal with the uncertainty of mainshock-aftershock sequences. Three resilience enhancement strategies, i.e., EHS, SPS, and RS, were considered in the three-level model. A tolerated random number generation method and a novel two-stage stochastic programming model were proposed to deal with uncertainties associated with earthquakes, including magnitude, occurrence time, occurrence intensity, and endogenous uncertainties. Subsequently, the optimization model is efficiently solved using SAA and PH. The proposed framework was examined on a ±800 kV converter station, with results indicating that 1) the implementation of the optimal pre-earthquake preparation strategy could greatly improve the resilience expectation of the converter station from 0.145 to 0.554, and the effect of aftershocks could not be ignored; 2) The EHS tended to harden the equipment of the same circuit to ensure the low transmission capacity of the converter station against various earthquakes, while the SPS tended to prepare spare parts for the equipment which has not been hardened; 3) if the budget is adequate, prioritizing the SPS can yield substantial improvements in the resilience of the converter station, otherwise, shifting focus to the EHS may be a valid alternative; 4) The SPS was more effective than the EHS in improving the converter station resilience with high values of mainshock PGA and occurrence intensity of aftershock.

It is noteworthy that the proposed method focuses on converter stations, which have simpler composition and operation compared to large-scale power networks. This work has not yet considered factors such as generation capacity, load demand, and the vulnerability of incoming and outgoing feeders in large-scale power networks. It can be further extended to more realistic studies of large-scale power networks. However, for large-scale power networks, the resilience assessment and optimization method might be time-consuming, necessitating some assumptions and approximation methods. Therefore, future research will focus on the earthquake resilience assessment and optimization for large-scale power networks. Furthermore, this work assumed that damaged equipment is repaired sequentially and that accurate information regarding the damage of each equipment is collected post-earthquake. In reality, it is very challenging to quickly and accurately identify the damage severity of each piece of equipment in the converter station. The optimization method for earthquake resilience of the converter station by taking account of the inspection time and errors will be conducted.

# APPENDIX A

## Two-stage stochastic programming model
### without aftershocks

In this model, the earthquake scenario $\xi$ is characterized by the mainshock occurring at $T_1^\xi$. The two-stage stochastic programming model without aftershocks is shown as following:

#### First-stage:

$$z = \max \mathbb{E}[R(\xi)]. \qquad (A.1)$$

$$\sum_{i \in \mathbf{N}} n_i C_i^h + \sum_{i \in \mathbf{N}} r_i C_i^s \leq C. \qquad (A.2)$$

$$n_i \in \{0,1\}, r_i \text{ is integer}, \forall i \in \mathbf{N}. \qquad (A.3)$$

#### Second-stage:

$$\max R(\xi) = \sum_{t=T_1^\xi}^{T_1^\xi+T} \Delta t F^\xi(t)/T. \qquad (A.4)$$

Eqs. (8)-(12), (15), and (19)-(22).

$$\begin{cases} \tau_i^{T_1^\xi} \geq R_i^{T_1^\xi} - (p_i^{T_1^\xi} - \delta_i^{T_1^\xi} n_i^{T_1^\xi}) \\ \tau_i^{T_1^\xi} \leq 1 + R_i^{T_1^\xi} - (p_i^{T_1^\xi} - \delta_i^{T_1^\xi} n_i^{T_1^\xi}) \end{cases}, \forall i \in \mathbf{N}, \forall \xi \in \mathbf{E}. \qquad (A.5)$$

$$g_i^\xi(t+1) = g_i^\xi(t) + \frac{m_i^\xi(t)}{w_i}, \qquad (A.6) \\ \forall i \in \mathbf{N}, \forall t \in \mathbf{T} / T_1^\xi, \forall \xi \in \mathbf{E}$$

$$s_i^\xi(t) = \max \left\{ s_i^\xi(t-1), \mathbf{1}_{g_i^\xi(t-1) \ge 1} \right\}, \qquad (A.7) \\ \forall i \in \mathbf{N}, \forall t \in \mathbf{T} / T_1^\xi, \forall \xi \in \mathbf{E}$$

## REFERENCES

[1] Y. Guo, J. Li, Z. Hu, et al, "Impacts of information dissemination framework on public compliance behaviour: Evidence from sudden power outages," in Proc. 3rd Int. Conf. Public Manag. Digital Economy Internet Technol, Jinan, China, Sep. 6–8, 2024.

[2] Y. Tan, F. Qiu, A. K. Das, D. S. Kirschen, P. Arabshahi, J. Wang, "Scheduling post-disaster repairs in electricity distribution networks," IEEE Trans. Power Syst., vol. 34, no. 4, pp. 2611-2621, 2019.

[3] X. Yin, X. Zhang, R. Cai, et al, "Fast 1-D velocity optimization inversion to 3D velocity imaging: A case study of Sichuan Maerkang earthquake swarm in 2022," Sustainability., vol. 14, no. 23, p. 15909, 2022.

[4] G. L. Ma, Q. Xie, A. Whittaker, "Seismic performance assessment of an ultra-high-voltage power transformer," Earthq. Spectra., vol. 35, no. 1, pp. 423-445, 2019.

[5] W. Yu, Y. Xue, J. Luo, et al, "An UHV grid security and stability defense system: Considering the risk of power system communication," IEEE Trans. Smart Grid., vol. 7, no. 1, pp. 491-500, 2015.

[6] G. L. Ma, Q. Xie, A. Whittaker, "Seismic performance assessment of an ultra-high-voltage power transformer," Earthq. Spectra., vol. 35, no. 1, pp. 423-445, 2019.

[7] S. Wang, K. H. Ahmed, G. P. Adam, et al, "A novel converter station structure for improving multiterminal HVDC system resiliency against AC and DC faults," IEEE Trans. Ind. Electron., vol. 67, no. 6, pp. 4270-4280, 2019.

[8] S. Mousavizadeh, T. G. Bolandi, M. R. Haghifam, et al, "Resiliency analysis of electric distribution networks: A new approach based on modularity concept," Int. J. Electr. Power Energy Syst., vol. 117, p. 105669, 2020.
[9] S. A. Modaberi, S. Tohidi, S. G. Zadeh, et al, "A review of power system resilience assessment and enhancement approaches by focusing on wind farms and wind turbines," IET Renew. Power Gener., vol. 17, no. 9, pp. 2391-2410, 2023.

[10] Q. Sun, Z. Wu, W. Gu, et al, "Seismic-resilient planning for integrated energy system: A risk-economic coordination perspective," IEEE Trans. Power Syst., 2024.

[11] Y. Zhang, L. Xu, C. Deng, et al, "Resilience improvement strategy of distribution network based on network reconfiguration in earthquake disaster scenario," in Proc. 2023 Panda Forum Power Energy, Chengdu, China, Apr. 2023, pp. 2193-2197.

[12] Z. Yang, A. Martí, Y. Chen, et al, "Optimal resource allocation to enhance power grid resilience against hurricanes," IEEE Trans. Power Syst., vol. 38, no. 3, pp. 2621-2629, 2022.

[13] A. V. Díaz, I. H. Moya, "Urban wind energy with resilience approach for sustainable cities in tropical regions: A review," Renew. Sustain. Energy Rev., vol. 199, p. 114525, 2024.

[14] A. M. Stanković, K. L. Tomsovic, F. De Caro, et al, "Methods for analysis and quantification of power system resilience," IEEE Trans. Power Syst., vol. 38, no. 5, pp. 4774-4787, 2022.

[15] G. Huang, J. Wang, C. Chen, et al, "Integration of preventive and emergency responses for power grid resilience enhancement," IEEE Trans. Power Syst., vol. 32, no. 6, pp. 4451-4463, 2017.

[16] S. Hosseini, K. Barker, J. E. Ramirez-Marquez, "A review of definitions and measures of system resilience," Rel. Eng. Syst. Saf., vol. 145, pp. 47-61, 2016.

[17] M. Panteli, C. Pickering, S. Wilkinson, "Power system resilience to extreme weather: Fragility modeling, probabilistic impact assessment, and adaptation measures," IEEE Trans. Power Syst., vol. 32, no. 5, pp. 3747-3757, 2016.

[18] D. Cao, J. Zhao, J. Hu, et al, "Physics-informed graphical representation-enabled deep reinforcement learning for robust distribution system voltage control," IEEE Trans. Smart Grid, vol. 15, no. 1, pp. 233-246, 2023.

[19] J. Hou, F. Teng, W. Yin, Y. Song, Y. Hou, "Preventive-Corrective Cyber-Defense: Attack-Induced Region Minimization and Cybersecurity Margin Maximization," IEEE Trans. Power Syst., vol. 39, no. 3, pp. 5324-5336, 2023.

[20] D. Alvarado, R. Moreno, A. Street, et al, "Co-optimizing substation hardening and transmission expansion against earthquakes: A decision-dependent probability approach," IEEE Trans. Power Syst., vol. 38, no. 3, pp. 2058-2070, 2022.

[21] Y. Huang, X. Fan, X. Xu, et al, "Fault recovery of integrated energy system considering electric-gas coupling characteristic and robust chance-constrained power balance," Electr. Power Syst. Res., vol. 227, p. 109979, 2024.

[22] Z. Xu, J. E. Ramirez-Marquez, Y. Liu, et al, "A new resilience-based component importance measure for multi-state networks," Rel. Eng. Syst. Saf., vol. 193, p. 106591, 2020.

[23] Y. Wang, Y. Xu, J. He, et al, "Coordinating Multiple Sources for Service Restoration to Enhance Resilience of Distribution Systems," IEEE Trans. Smart Grid., vol. 10, no. 5, pp. 5781-5793, 2019.

[24] G. Zhang, F. Zhang, X. Zhang, et al, "Sequential disaster recovery model for distribution systems with co-optimization of maintenance and restoration crew dispatch," IEEE Trans. Smart Grid., vol. 11, no. 6, pp. 4700-4713, 2020.

[25] C. Luo, L. Liang, X. Bai, "Review of HVDC systems for enhancing resilience after extreme natural events," J. Eng., vol. 2021, no. 11, pp. 653-664, 2021.

[26] X. Liu, Q. Xie, "A multi-model probabilistic framework to evaluate seismic resilience of UHV converter stations," Eng. Struct., vol. 300, p. 117153, 2024.

[27] J. Li, T. Wang, Q. Shang, "Probability-based seismic resilience assessment method for substation systems," Struct. Infrastruct. Eng., vol. 18, no. 1, pp. 71-83, 2021.

[28] S. Wang, K. H. Ahmed, G. P. Adam, et al, "A novel converter station structure for improving multiterminal HVDC system resiliency against AC and DC faults," IEEE Trans. Ind. Electron., vol. 67, no. 6, pp. 4270-4280, 2019.

[29] W. Wen, C. Zhai, D. Ji, et al, "Framework for the vulnerability assessment of structure under mainshock-aftershock sequences," Soil Dyn. Earthq. Eng., vol. 101, pp. 41-52, 2017.

[30] K. Wang, Z. P. Xu, Y. Liu, et al, "Resilience enhancement for multi-state interdependent infrastructure networks: From a preparedness perspective," IEEE Trans. Rel., vol. 72, no. 1, pp. 190-203, 2022.

[31] C. Zou, K. Wang, T. Xiahou, et al, "Two-stage distributionally robust optimization for infrastructure resilience enhancement: A case study of 220 kV power substations under earthquake disasters," IEEE Trans. Rel., 2024.

[32] C. Zou, C. Deng, Z. Xue, et al, "Distributionally robust optimal strategy for resilience enhancement of power distribution networks under uncertain earthquake failures," in Proc. IEEE 2nd Int. Conf. Power Sci. Technol, Chengdu, China, May 2024, pp. 1200-1205.

[33] Y. P. Fang, E. Zio, "An adaptive robust framework for the optimization of the resilience of interdependent infrastructures under natural hazards," Eur. J. Oper. Res., vol. 276, no. 3, pp. 1119-1136, 2019.

[34] C. Zhao, N. Yu, T. Peng, et al, "Vulnerability assessment of AP1000 NPP under mainshock-aftershock sequences," Eng. Struct., vol. 208, p. 110348, 2020.

[35] M. Ouyang, Y. Fang, "A mathematical framework to optimize critical infrastructure resilience against intentional attacks," Comput.-Aided Civ. Infrastruct. Eng., vol. 32, no. 11, pp. 909-929, 2017.

[36] Q. P. Zheng, J. Wang, A. L. Liu, "Stochastic optimization for unit commitment-A review," IEEE Trans. Power Syst., vol. 30, no. 4, pp. 1913-1924, 2014.

[37] T. Lagos, R. Moreno, A. N. Espinosa, et al, "Identifying optimal portfolios of resilient network investments against natural hazards, with applications to earthquakes," IEEE Trans. Power Syst., vol. 35, no. 2, pp. 1411-1421, 2019.

[38] Z. Wang, B. Chen, J. Wang, et al, "Robust optimization based optimal DG placement in microgrids," IEEE Trans. Smart Grid., vol. 5, no. 5, pp. 2173-2182, 2014.

[39] X. V. Doan, "Distributionally robust optimization under endogenous uncertainty with an application in retrofitting planning," Eur. J. Oper. Res., vol. 300, no. 1, pp. 73-84, 2022.

[40] Y. Li, S. Lei, W. Sun, C. Hu, Y. Hou, "A distributionally robust resilience enhancement strategy for distribution networks considering decision-dependent con-tingencies," IEEE Trans. Smart Grid., vol. 15, no. 2, pp. 1450-1465, 2024.

[41] S. Peeta, F. S. Salman, D. Gunnec, et al, "Pre-disaster investment decisions for strengthening a highway network," Comput. Oper. Res., vol. 37, pp. 1708-1719, 2010.

[42] A. H. Darooneh, A. Mehri, "A nonextensive modification of the Gutenberg-Richter law: q-stretched exponential form," Physica A., vol. 389, no. 3, pp. 509-514, 2010.

[43] M. Taroni, "Back to the future: Old methods for new estimation and test of the Gutenberg-Richter b-value for catalogues with variable completeness," Geophys. J. Int., vol. 224, no. 1, pp. 337-339, 2021.

[44] Y. Dong, D. M. Frangopol, "Risk and resilience assessment of bridges under mainshock and aftershocks incorporating uncertainties," Eng. Struct., vol. 83, pp. 198-208, 2015.

[45] M. Taroni, A. Akinci, "Good practices in PSHA: Declustering b-value estimation, foreshocks and aftershocks inclusion; A case study in Italy," Geophys. J. Int., vol. 224, no. 2, pp. 1174-1187, 2021.

[46] X. Liu, Q. Xie, "Resilience-based post-earthquake recovery strategies for substation systems," Int. J. Disaster Risk Reduct., vol. 96, p. 104000, 2023.

[47] A. J. Kleywegt, A. Shapiro, T. Homem-de-Mello, "The sample average approximation method for stochastic discrete optimization," SIAM J. Optim., vol. 12, no. 2, pp. 479-502, 2002.

[48] J. P. Watson, D. L. Woodruff, "Progressive hedging innovations for a class of stochastic mixed-integer resource allocation problems," Comput. Manag. Sci., vol. 8, no. 4, pp. 355, 2011.

[49] C. M. Wei, G. J. Meng, W. W. Wu, et al, "Time-independent grid-based forecast model for $M \ge 6.0$ earthquakes in Southeastern Tibetan Plateau using GNSS strain rates and seismicity," Bull. Seismol. Soc. Am., vol. 114, no. 1, pp. 258-277, 2024.

[50] A. Dadashzade, H. Bagherzadeh, M. Mottaghizadeh, et al, "Electricity distribution networks resilience in area exposed to salt dust: Fragility curve modeling of insulators, Monte Carlo-based resilience assessment, and enhancement measures," Heliyon., vol. 10, no. 16, 2024.

[51] M. Panteli, C. Pickering, S. Wilkinson, et al, "Power system resilience to extreme weather: Fragility modeling, probabilistic impact assessment, and adaptation measures," IEEE Trans. Power Syst., vol. 32, no. 5, pp. 3747-3757, 2016.
<figure><img src="image_12.png" /></figure>

**Kai Wang** received the B.S. degree in industrial engineering and M.S. degree in mechanical engineering from the University of Electronic Science and Technology of China, Chengdu, China, in 2020 and 2023, respectively. He is currently a Ph.D. student in industrial engineering at CentraleSupélec, Université Paris-Saclay, France.

His research interest includes the assessment and optimization of resilience in infrastructure networks.

<figure><img src="image_13.png" /></figure>

**Zhihang Xue** received the M.S. degree in measurement and control technology and automation from the University of Electronic Science and Technology of China, Chengdu, China, in 2012. He is currently with the State Grid Sichuan Electric Power Research Institute and is also a Ph.D. student at the University of Electronic Science and Technology of China.

His research interest includes the seismic resilience and disaster prevention in power systems.

<figure><img src="image_14.png" /></figure>

**Di Cao** received the Ph.D. degree in control science and engineering from the University of Electronic Science and Technology of China in 2021. He is currently an associate professor at this university.

His research interest includes optimization of distribution network and applications of machine learning in power systems.

<figure><img src="image_15.png" /></figure>

**Yu Liu** (Senior Member, IEEE) received the Ph.D. degree in mechatronics engineering from the University of Electronic Science and Technology of China, Chengdu, China, in 2010.

He is currently a Full Professor with the School of Mechanical and Electrical Engineering, University of Electronic Science and Technology of China. He was a Visiting Predoctoral Fellow with the Department of Mechanical Engineering, Northwestern University, Evanston, IL, USA, from 2008 to 2010, and a Postdoctoral Research Fellow with the Department of Mechanical Engineering, University of Alberta, Edmonton, AB, Canada, from 2012 to 2013. He has authored/co-authored more than 100 peer-reviewed papers in international journals and conferences. His research interests include system reliability modeling and analysis, maintenance decisions, prognostics and health management, and design under uncertainty.

Prof. Liu is an editorial board of Reliability Engineering & System Safety, Quality and Reliability Engineering International, and Engineering Optimization. He is also an Associate Editor for *IISE Transactions* and IEEE TRANSACTIONS ON RELIABILITY. He is an ISEAM Fellow.

<figure><img src="image_16.png" /></figure>

**Yi-Ping Fang** received his Ph.D. degree in Industrial Engineering from École Centrale Paris (ECP), France, and was a Research Fellow at ETH Zurich, Switzerland from 2015 to 2017. He is currently a Chair professor at the Laboratoire Génie Industriel, Centrale-Supélec, Université Paris-Saclay, France.

His primary research focuses on the study of advanced computational methods for risk, reliability and resilience in complex industrial systems and infrastructure networks. His specific interests include data-

driven stochastic and robust optimization, risk and decision analysis, and reinforcement learning.