Resilience-Based Restoration Sequence Optimization for Metro Networks: A Case Study in China

Abstract

Metro station restoration sequence optimization is crucial during post-disaster recovery. Taking both budget limitations and repair time uncertainty into account, this paper proposes a resilience-based optimization model for choosing an optimal restoration sequence scheme, maximizing the global average efficiency, under the condition that the network accessibility meets given resilience requirements. Evolutionary algorithm NSGA-II is applied to solve the model. A Case study in Nanjing and Zhengzhou gives insights into restoration sequence strategies for decision-makers. Results show that a ring network is more robust than a radial network under the same scale attack. Under limited budget, the optimal restoration sequence is closely related to the damaged stations’ location and repair time. Specifically, if damaged stations’ distribution is relatively centralized and transfer stations need more repair time, giving repair priority to transfer stations is not always the best strategy. If damaged stations’ distribution is relatively scattered and all stations’ repair time is the same, the station with a bigger node degree should be repaired earlier. However, this conclusion may be invalid if transfer stations repair time is far longer than others. Sensitivity analysis show that the total budget is more sensitive than one day’s budget in the entire restoration phase. However, in the emergency phase, increasing one day’s budget is more significant for shortening recovery time. The proposed model can contribute to effective and flexible decision-making for metro network restorations.

This work is licensed under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

1. Introduction

As one of the fastest and high-capacity means of transportation, metros have become an effective way of solving congestion problems. However, metro networks can be vulnerable to certain natural disasters, terrorist attacks or equipment failures, and any damages can cause severe socio-economic losses [1–3]. In order to mitigate losses caused by these types of disruptions, damages need to be repaired immediately. Working with a limited budget, decision-makers need to determine a set of criteria for a restoration sequence to ensure the system’s rapid and effective recovery.

Some research has been done in the pursuit of the optimal restoration sequence scheme for transportation networks, which can be conceptualised as a network design problem (NDP). Throughout the last decade, researchers have combined resilience evaluation and the NDP to explore the methodology of the resilience-based optimal restoration plan. Infrastructure resilience is an ability to absorb, adapt to, and rapidly recover from potentially disruptive or destructive events [4, 5]. The advantage of resilience analyses is to reduce the disruption’s risks and improve the system’s recovery ability. The resilience index has been widely used to quantify the system’s recovery ability [6–8].

Resilience-based restoration optimization was widely used in road network or road-bridge network, but few studies applied to the metro network. Previous research considered various network metrics, including (1) connectivity [2, 9, 10]; (2) travel delay cost [11]; and (3) network capacity [12, 13]. However, as these studies focus on a single metric, ignoring the trade-offs or connections between different metrics, they cannot meet the needs for multi-objective decision processes, especially for long-term restoration. Focusing on network topology characteristics, this paper proposes an optimal metro restoration sequence model to guide multi-objective decision-making. In existing research, node degree, network accessibility, and global average efficiency are the main network topology characteristics for metro networks [14–16]. Those metrics have all been considered in our model.

Some researchers considered multi-objective optimization, such as Matisziw’s [17], Liu, Zhai, and Dong [18], Somy, Shafaei, and Ramezanian [19]. They all adopted a weighting method to deal with a multi-objective problem. However, different weighting factors will lead to contradictory results. The proposed model focuses on metro topology characteristics, and also can transfer a multiobjective problem to a single-objective combinatorial problem with additional constraints that respect the other goals. This method can reduce the subjectivity caused by improper weighting factor setting. Furthermore, the proposed model also can present a flexible framework. The objectives and constraints can be set flexibly by decision-makers to meet multi-stage objectives.

The majority of existing restoration optimization models based on resilience only considered the resilience index of the network performance, but ignored the recovery trajectories of restoration schemes. For example, in Figure 1, the resilience index of curve 1 and curve 2 are equal, but with significant differences in trajectories. Thus, we add an extra objective that makes the network performance recover to the pre-given value at a specific time. This method can avoid the above question validly.

[figure omitted; refer to PDF]

Genetic algorithm [10, 18, 20], tabu search algorithm [13], and greedy algorithm [21] were widely applied to deal with discrete NDP problem. Compared with the greedy algorithm and tabu search algorithm, the Genetic Algorithm has better global search ability and convergence control. NSGA-II introduces the elitist strategy to expand the sampling space, prevent the loss of the best individual, and improve the operation speed and robustness of the algorithm [22]. Thus, NSGA-II is adopted to deal with optimization questions.

Most metro stations and lines are underground. The restoration process is relatively complex and may continue for a long time. Multi-stage decision-making is necessary for large-scale metro network restoration. Focusing on network topology characteristics, this paper provides a valuable and flexible tool for optimizing the metro restoration sequences, while taking limited budget restraints and repair time uncertainty into account, giving a guideline for multi-objective decision-making. Unlike most bi-objective research, we transfer a multi-objective problem to a single-objective combinatorial problem with additional constraints that respect the other goals. The research can enhance the effective allocation of disaster relief resources and improve post-disaster recovery, which improves the metro’s resilience following a disruption.

The remainder of the paper is organised as follows. Firstly, the framework and methodology of the resilience-based, multiobjective optimal restoration sequence model are introduced. Then, three empirical case studies using a real-world metro network are presented. Finally, the conclusion and recommendations are summarised.

2. Methodology

2.1. Problem Statement

As metro networks can be particularly vulnerable during extreme events, potentially causing damage to stations, actions should be adopted to ensure the network is able to recover as quickly and effectively as possible. In reality, postdisaster recovery is generally a long-term project. As decision-makers may have different objectives at different stages. Based on the need for practice and the key points of this research, some assumptions are made in order to simplify the complexity of the problem.

(1) As alternative paths often exist in a damaged network, passengers can be evacuated from affected zones to safer areas without the help of other means of transportation, proving that network accessibility is more important than global average efficiency in the short term. Networks with different restoration schemes may have the same network accessibility. However, the global average efficiency may vary. Hence, the proposed model aims to seek an optimal restoration sequence scheme to maximize the global average efficiency, under the condition that network accessibility meets predetermined resilience requirements (such as the recovery time, the resilience index, and recovery rate).

(2) Consideration of uncertainties associated with restoration projects is critical for effective decision-making. Uncertainties can include the network damage status, the residual network capacity [12], and the effects of aftershocks [23], among others. This paper assumes that postdisaster network damaged conditions are known, and ignores the effects of any potential aftershocks. Besides, the restoration duration for a damaged metro station obeys a certain probability distribution.

(3) Fund allocation is crucial for network recovery and any restoration scheme will be particularly sensitive to budget changes [24]. This paper assumes that the total budget is limited, and the utilisation of funds can be different throughout different stages, meaning that some stations can be repaired at the same time.

2.2. Network Topology Performance Metrics

Based on graph theory, a physical metro network can be abstracted as an undirected graph G (N, L) composed of N nodes and L links. In this graph, stations are abstracted as nodes, and links between two adjacent stations are abstracted as links. The most common network modelling methods include the Space L, Space P, and the Space R method [15]. Because the Space L method can directly reflect a metro network’s physical structure, using it to calculate network topology metrics can better reflect passengers’ travel behaviour. Thus, the Space L method is adopted to build a network topology model.

2.2.1. Network Accessibility

A graph G can be completely represented by an adjacency matrix $AN\times N$ with element a_ij = 1 if nodes i and j connect directly, otherwise, a_ij = 0. By default, when i = j, the a_ij = 0.

Accessibility matrix B (N × N) can describe the network accessibility. If node i to node j have a path (i ≠ j), b_ij = 1, otherwise, b_ij = 0. By default, when i = j, the b_ij = 1. In the operational state, the value of b_ij (i ≠ j) is equal to 1. When some nodes are damaged and there are no alternative paths from node i to node j, the b_ij = 0.

To depict the node importance in the model, the node degree is considered. The degree of node i is defined as the number of the nodes that directly connect with node i [25]. The degree matrix D is a 1 × N matrix. Because we constructed an undirected graph, so the D does not consider the incoming and outgoing links. The degree of node i equals to the sum of the value in matrix A’s ith row. Matrix C = D × B is used to consider the node importance and accessibility comprehensively. This method also can help to decrease the matrix dimension and built a mathematical model. Matrix C is a 1 × N matrix, the ith element of the matrix represents the product of the degree of node i and the total number of node j which connect with node i by at least one path. To generate a network accessibility resilience index, the dimension of the matrix C should be reduced further. Z (t) = Sum (C_i) (i ∈ N) is adopted to depict the network accessibility at time t.$$\begin{array}{}\text{(1)}& Z\text{t}={\displaystyle \sum _{i=1}^N{C}_i}\text{t}={\displaystyle \sum _{i=1}^{N}D\times B_t},\end{array}$$where Z (t) is the network accessibility at time t. C_i (t) is the value of ith element of matrix C at time t. D is the degree matrix $N\times N$. B_t is the accessibility matrix $1\times N$ at time t. N is the total number of nodes.

2.2.2. Global Average Efficiency

The global average efficiency E quantifies the network efficiency in the global scope, it can be calculated by equation (2) [26]. The shortest path between two nodes can be found by the Dijkstra algorithm [27].$$\begin{array}{}\text{(2)}& E=\frac{1}{NN-1}{\displaystyle \sum _{i\ne j}^N\frac{1}{d_{ij}}},\end{array}$$where E is the global average efficiency. d_ij is the shortest path length between node i and node j.

2.3. Resilience Index

As shown in Figure 1, F (t₀) represents the system’s performance value in the original stable state. F (t) represents the system’s performance curve, whose value continues to change with time. t_e represents the time of disruption occurring. t_f represents the time for the system to recover to its original state. The resilience of the system at time t_r can be quantified as resilience index R_F, shown as equation (3). It is equal to the ratio of the area covered by the performance curve F (t) and the area covered by the curve F (t₀) during the period of t_e to t_r. In some research, the shadow part is named the “resilience triangle” [6].$$\begin{array}{}\text{(3)}& R_F{t}_r=\frac{{\displaystyle {\int }_{t_e}^{t_r}Ft\text{d}t}}{F{t}_0{t}_r-t_e}.\end{array}$$

2.4. Resilience-Based Multiobjective Optimization for Restoration Sequence

Taking budget limitations and repair time uncertainty into account, the optimal restoration sequence scheme aims to maximize global average efficiency, when the network accessibility satisfies pregiven resilience requirements (including the resilience index and recovery rate). A multiobjective problem can be transferred to a single-objective combinatorial problem with additional constraints that respect the other goals of the problem. Equations (14) and (15) represent the other goals of the problem. The optimal sequence scheme should make the network accessibility’s resilience index and recovery rate reach the threshold at least and maximize as much as possible under the budget limitation. Then, from the partial solutions, the solution with the maximum global average efficiency is selected as the optimal scheme.$$\begin{array}{}\text{(4)}& N_d=s_1,s_2,s_3,\dots ,s_m,\text{(5)}& t_i\sim N\overline{t_l},{\delta }_i.\end{array}$$

Objective:$$\begin{array}{}\text{(6)}& \text{Max}E{t}_r=\frac{1}{NN-1}{\displaystyle \sum _{i=1,i\ne j}^N\frac{1}{d_{ij}},}\end{array}$$

Decision variables:$$\begin{array}{}\text{(7)}& \begin{array}{cc}{\beta }_{ik}^t\in \mathrm{0,1},& \forall i\in N_d\end{array},\text{(8)}& t_i^s,\text{(9)}& t_i^e=t_i^s+t_i,\text{(10)}& N_a^t=s_a,s_b,s_c,\dots ,s_p.\end{array}$$

Constraints:$$\begin{array}{}\text{(11)}& \theta T_Z\le U,\text{(12)}& K\le \text{INT}\frac{u}{\theta },\text{(13)}& {\displaystyle \sum _{i=1,i\in N_a^{t_r}}^{N_a^{t_r}}{t}_i}\le T_Z,\text{(14)}& R_Z=\frac{{\displaystyle {\int }_{t_0}^{t_r}Zt\text{d}t}}{Z{t}_0{t}_r-t_0}\ge f_1,\text{(15)}& Q_Z=\frac{Z{t}_r}{Z{t}_0}\ge f_2,\text{(16)}& t_i^s\ge 0,\text{(17)}& \max t_i^e\le t_r.\end{array}$$

Equation (4) is the set of damaged stations after disaster. Equation (5) reveals the repair time of node i obeys normal distribution whose mean is $\overline{t_i}$ and variance is δ_i. It reflects the uncertainty of repair time. Equation (6) is the objective of the model, which maximizes the global average efficiency. Equations (7)–(10) is the decision variables. Equation (7) represents whether the station i is repaired at time t by team k. Equation (8) is the repair start time of node i. Equation (9) is the repair finished time of node i. It equals to the node i’s repair start time plus the i’s repair time. Equation (10) is the set of repaired stations at time t, and it can get by equation (7). Equation (11) demonstrates the total budget limitation. Equation (12) demonstrates the limitation of the number of maintenance teams in one day, which can reflect the abundance of resources. Equation (13) indicates the limitation of total repair time for repaired nodes. Equation (14) is the resilience requirement of the network accessibility, which requires the resilience index to reach f₁ at least at time t_r. Equation (15) is the other requirement of the resilience, which requires the network accessibility to get back into Z (t₀) $\ast$ f₂ at least at time t_r. The range of f₁ and f₂ is $0,1$. The former reflects the recovery degree of Z with continuous time. The latter reflects the recovery degree of Z at discrete time. In practice, decision-makers can preset the different f₁ and f₂ to reflect the different resilience requirements. Equation (16) indicates the repair start time of node i should be bigger than 0. Equation (17) indicates the maximized repair finished time of node i should not be bigger than the pre-given deadline. The notations are shown in Table 1.

Table 1

Notations.

$N_d$	The set of damaged nodes at time t0	t0 (day)	The time of the disruption occurring
tr (day)	The pregiven deadline	E	Global average efficiency
N	The total nodes in the network	dij	The length of the shortest path from node i to node j
${\beta }_{ik}^t$	Binary variable – if node i is repaired by team k at time t, then ${\beta }_{ik}^t=1$; otherwise, ${\beta }_{ik}^t=0$	$N_a^t$	The set of repaired nodes at time t
$t_i^s$	Repair start time of node i	$t_i^e$	Repair finished time of node i
$Z$	Network accessibility	Bt	The accessibility matrix (N × N) at time t
Rz	The resilience index of network accessibility at time tr	Qz	The ratio of network accessibility at time tr and t0
f1, f2	Constants that represent the resilient requirement	θ ($)	A constant that represents the cost of one day’s work for a maintenance team
$U$ ($)	The total budget	Tz (day)	The total restoration time of one team under the total budget limitation
$u$ ($)	The maximized budget which can be paid in one day	K	The number of the maintenance teams in one day
ti (day)	The repair time of damaged node i	D	The degree matrix (1 × N) of the network
$\overline{t_i}\hspace{0.17em}$ (day)	The average repair time of node i	δi (day)	The variance of repair time of node i

2.5. Solution Process and Algorithm

The solution process of the proposed model for restoration sequence is shown as Figure 2. The solution process is divided into two-stage. The first stage aims to make the network accessibility’s resilience index R_Z and recovery rate Q_Z reach to the threshold and maximize as much as possible. The NSGA-II algorithm is adopted to deal with this problem. Under the budget limitation, the Pareto front solutions can be available. The second stage aims to find out the optimal restoration sequence scheme to maximize the network efficiency E from the Pareto front solutions.

[figure omitted; refer to PDF]

The solution process is as follows:

(1) Parameters input. All parameters can be found in scenario descriptions. Because the damaged station’s repair time obeys normal distribution, the repair time of each node t_i is randomly allocated according to the distribution function.

(2) According to equations (11) and (12) to calculate the value of T_Z, K_max. k = 1, 2,…, K_max

(3) Stage One: using NAGA-II to deal with multi-objective optimization.

  Objective 1: maximize the network accessibility’s resilience index R_Z

  Objective 2: maximize the network accessibility’s recovery rate Q_Z

The NSGA-II has four important parameters: the number of individuals in the population (Population), the maximum number of generations (Maximum generations), cross-over rate, mutation rate, and the number of individuals in each tournament. The cross-over rate is used to decide the probability of two individuals need to crossover. The mutation rate is used to decide the probability that any individual needs mutation. The termination criterion is the number of iterations is equal to the maximum generations.

(1) Initialization—Define the basic parameters to control the NSGA-II: number of individuals in the population (n_p), the maximum number of generations ($n_{g\text{max}}$), cross-over rate (c), mutation rate (m_r), number of individuals in each tournament (n_tour), and generation counter (Gen).

(2) Generation of the initial population—The illustration of chromosome coding is shown in Figure 3. x₁ (i) represents the order priority value of generation 1, 0 < x (i) < 1. The smaller the value, the higher the priority. y₁ (i) represents the resource allocation of generation 1, k = 1, 2, …, K_max. Generate 3n_p individuals and among them select n_p with distinct characteristics. With one individual, calculate the $\hspace{0.17em}{t}_i^s$, $t_i^e,$R_Z, Q_Z, and select the individual that meets the constraints (13)–(17). Otherwise, an individual is generated repeatedly until it meets all the constraints. Repair scheduling by chromosome decoding is shown in Figure 4.

  Because the R_Z is integral form, it should be replaced by the sum of discrete forms. The repair finished time of node i can be got by decoding. Sort the repair finished time of node i in ascending order, the time axis can be divided into p segments. The calculation of R_Z can be achieved by equation.$$\begin{array}{}\text{(18)}& R_Z=\frac{{\displaystyle {\int }_{t_0}^{t_r}Zt\text{d}t}}{Z{t}_0{t}_r-t_0}=\frac{{\displaystyle {\sum }_0^{p-1}Z{t}_{p+1}+Z{t}_p}{t}_{p+1}-t_p}{2Z{t}_0{t}_r-t_0}.\end{array}$$

(3) Classification according to the dominance criterion (nondominated sorting)—classify the individuals in fronts based on the value of R_z, Q_z. The nondominated sorting criterion states that for two solutions, X_a and X_b, X_a dominates X_b if X_a is better than or equal to X_b in all objectives and X_a is better than X_b in at least one of the objectives. If there are no other decision variables that can dominate a decision variable, the decision variable is called a nondominated solution. All nondominated solutions in chromosome set X is defined as the first nondominated layer.

(4) Iterative process—while Gen < $n_{g\text{max}}$, repeat the following loop:

(a) Selection—the individuals selected for cross-over and mutation are chosen by tournaments which then select each individual in the following way: among the individuals of the population, n_tour are randomly chosen and the best classified is selected according to the dominance criteria.

(b) Crossover—using operators for the crossover, resulting in c $\ast$ n_p descendants.

(c) Mutation—using operators for the mutation, resulting in m $\ast$ n_p individuals.

(d) Generation of the temporary population—with the n_p individuals obtained through crossover and mutation, a new population is generated and added to the previous temporary population of 2n_p individuals.

(e) Classification of the temporary population—the 2n_p individuals of the temporary population are classified again according to the criteria of dominance and diversity.

  For the solution $X_q\in X$ in the i^st dominating layer, the crowding distance $C{X}_q$ can be calculated as equation.$$\begin{array}{}\text{(19)}& C{X}_q=\frac{Q_Z{X}_q+1-Q_Z{X}_q-1}{Max{Q}_{Z}X-Min{Q}_{Z}X}+\frac{R_Z{X}_q+1-R_Z{X}_q-1}{Max{R}_{Z}X-Min{R}_{Z}X},\end{array}$$

  where $\text{Max}{Q}_{z}X$ and $\text{Min}{Q}_{z}X$ are the maximum and minimum values of $Q_{z}X$ in the i^st dominant layer, respectively. $\text{Max}{R}_{z}X$ and $\text{Min}{R}_{z}X$ are the maximum and minimum values of $R_{z}X$ in the i^st dominant layer, respectively. $Q_z{X}_q+1$ and $Q_z{X}_q-1$ are the previous and next values of $Q_z{X}_q$ after the ascending order of $Q_{z}X$ in the i^st dominant layer. $R_z{X}_q+1$ and $R_z{X}_q-1$ are the previous and next values of $R_z{X}_q$ in the i^st dominant layer. The illustration can be shown in Figure 5. The farther each individual is located form the remaining solutions, the better its classification.

(f) Selection of a new population—A new population is obtained from the temporary population.

(g) Updating of the generation counter—Gen = Gen + 1.

(h) Output the Pareto front solutions.

(4) Stage two: output the optimal sequence scheme that maximizes the network efficiency E (t_r).

[figure omitted; refer to PDF]

Calculate the network efficiency E (t_r) of each Pareto front solution and find out the one that makes the E (t_r) maximized.

The output of model is $X=i,k,t_i^s,\hspace{0.17em}{t}_i^e,\dots ,$. i represents the node ID. k represents the node is repaired by team k. $t_i^s$ is the repair start time of node i. $t_i^e\hspace{0.17em}$ is the repair finished time of node i.

3. Case Study

The performance indicators of restoration sequence scheme include “The percent of repaired nodes,” “The total budget,” “The total recovery time,” “The resilience index of the Z,” “The recovery ratio of Z,” and “The recovery ratio of E.” Indicators 2 to 5 also appeared in the scenario setup. However, scenario setup only provides those indicator’s upper or lower limits. The optimization results can provide the specific values of those indicators. Those values can be used to compare the different restoration sequence schemes.

3.1. Scenario I: The Location of the Damaged Stations Is Aggregated in Nanjing Metro

3.1.1. Scenario Descriptions

Until May 2020, Nanjing’s metro operated ten lines, with a total length of 378 km, and 174 stations, including 13 transfer stations. The Space L method is used to construct network topology. Some stations of Nanjing metro in the CBD are closed because of a natural disaster, which can be shown in Figure 6. There are twenty stations and four lines affected. Table 2 shows that each damaged station’s repair time obeys normal distribution. It reflects the uncertainty of repair time. Because the transfer stations often serve multiple lines, we assume they have a longer repair time than the remaining stations.

[figures omitted; refer to PDF]

Table 2

The distribution of repair time for damaged stations in scenario I.

ID	Mean, $\overline{t_i}$	Variance, δi	ID	Mean, $\overline{t_i}$	Variance, δi
7, 10, 24, 97, 117, 118	3	1	8, 13, 25, 99, 100, 116	4	1
12, 27, 28, 101	5	1	11, 98	6	1
9, 26	7	1

The objective of the base scenarios is to maximize the E under the condition that the Z’s resilience index and recovery rate reach 0.5 and 90% at least, respectively. The parameters of the model are shown in Table 3.

Table 3

The value of parameters of the model in scenario II and III.

The meaning of the parameter	Parameter	Value
The pregiven deadline time (day)	tr	19
The limitation of the total budget ($)	$U$	55000
The lower limit of Z’s resilience index	f1	0.50
The lower limit of Z’s recovery ratio	f2	90%
The cost of each day’s work for a maintenance team ($)	θ	1000
The maximized budget for one day ($)	$u$	3200

3.1.2. The Optimal Restoration Sequence Scheme in Scenario I

The NSGA-II and GA are used to search the optimal solutions for scenario I. The numerical application is performed on an ordinary desktop with Intel(R) Core i7-6700K CPU @ 8.00 GHz CPU and 16 GB memory. All algorithms are solved by Python.

According to Section 2.5, the solution process is divided into two-stage. The first stage aims to make the network accessibility’s resilience index R_Z and recovery rate Q_Z reach to the threshold and maximize as much as possible. The NSGA-II algorithm is adopted to find the Pareto front solutions. The second stage aims to find out the optimal restoration sequence scheme to maximize the network efficiency E from the Pareto front solutions.

Before using the GA algorithm, the model’s objective should be changed as equation (20). The sum of l₁, l₂, and l₃ is one. In this case, l₁ = 0.5, l₂ = l₃ = 0.25.$$\begin{array}{}\text{(20)}& \text{SAT}=\text{Max}{l}_{1}E+l_2{R}_Z+l_3{Q}_Z.\end{array}$$

NSGA-II and GA algorithm are both based on the damaged stations’ repair time equal to the mean value. The most appropriate parameter settings in NSGA-II and GA are shown in Table 4.

Table 4

The most appropriate parameter settings of NSGA-II and GA.

Parameter	NSGA-II	GA
Population (np)	100	100
Maximum generations ($n_{g\text{max}}$)	100	100
Cross-over rate (c)	0.9	0.9
Mutation rate (mr)	0.7	0.6
Number of individuals in each tournament (ntour)	10	—

Figures 7 and 8 illustrate the convergence of NSGA-II and GA over generations. The average computing time of NSGA-II is 6.9 min, and the GA is 3.3 min. Although the GA can achieve convergence about the 30th generation, NSGA-II also has a faster convergence speed, indicating fast convergence and good stability. However, GA for multi-objective optimization is limited by weight. Different weights lead to different results. NSGA-II can fill this gap and generate the Pareto frontier. NSGA-II can obtain some non-dominated solutions. The optimal sequence solution maximizes the network efficiency E from the Pareto frontier. The solutions can be shown in Figure 9.

[figure omitted; refer to PDF]

Table 5 and Figure 10 illustrate the optimal restoration sequence in scenario I. Because every node i’s repair time follows a specific normal distribution, the optimal restoration sequence scheme is not the only one. When each station’s repair time is shorter, the best results of the optimal scheme can be achieved. When each station’s repair time is longer, the worst results of the optimal scheme are realised. Through comparison of the best and the worst, the total budget only increases by 7.84%, the number of repaired stations increases by 36.36%, the resilience index of the Z increases by 74.21%, the recovery ratio of the Z and the E improve by 13.26% and 8.38% respectively, but the total recovery time only extends one day. In this case, the recovery time is mainly limited by the total budget. If the total budget is adequate and can support all nodes repaired, the best scheme’s total recovery time must be smaller than the worst. After all, the former repaired 75% of the nodes in 18 days, but the latter only repaired 55% of the nodes in 19 days. It is worth mentioning that the worst scheme does not satisfy the resilience constraints of objective of Table 3. It aims to reveal the sensitivity of the repair time’s uncertainty.

Table 5

The optimal restoration sequence in scenario I.

Indicators	The best of the optimal scheme	The worst of the optimal scheme	Average value	Ignore uncertainty
The sequences of the repaired nodes	A: 24, 117, 7, 117, 10, 116	A: 7, 24, 8, 13	—	A: 24, 98, 9, 26
	B: 97, 12, 9, 11, 12,	B: 97, 98, 99, 25		B: 97, 25, 99, 11
	C: 118, 25, 99, 26, 27	C: 118, 117, 9		C: 118, 7, 117, 10
The percent of repaired nodes	75%	55%	65%	60%
The total budget ($)	51000	55000	53000	52000
The total recovery time (day)	18	19	18.6	18
The resilience index of the Z	0.608	0.349	0.486	0.503
The recovery ratio of Z	98.46%	86.93%	93.74%	92.80%
The recovery ratio of E	87.91%	81.11%	84.60%	82.25%

[figure omitted; refer to PDF]

3.1.3. The Comparison of Proposed Optimal Scheme with Other Schemes

This paper also analyses the advantage of optimal schemes compared to random alternatives, which can be obtained by the way of free combination. In Figure 11, the blue line is the optimal restoration scheme based on the model. The red dots are Z’s value of random schemes on different days. The green line is the average Z of random schemes, which can be called an average random scheme. From the 1^st to the 11^th day, Z’s recovery is very slow. However, after the 11^th day, the advantage of the optimal scheme becomes more obvious, especially between days 12 and 13 where the Z is increases by approximately 113%. Although the average random scheme’s Z is higher than the optimal scheme on some days, its resilience index dropped by 42.8%.

[figure omitted; refer to PDF]

This paper also compares the objective values (the blue line in Figure 11) and KPIs with a solution that repairs the transfer stations first (the yellow line in Figure 11). The yellow line gradually increases and has a significant upward trend after the 13th day. The blue one gradually increases for the first 12 days, then goes up suddenly and always stays at a higher value. The resilience index of the yellow is smaller than the blue. The reasons are as follows: the repair time of transfer stations is longer than nontransfer stations. They need more time to be repaired, so the network’s recovery is slow in the early time. By analysing the transfer stations’ location, even they are repaired, they cannot connect with other no-damaged stations. Thus, giving them priority has little influence on the whole network recovery. For the optimal scheme, the number of repaired nodes are increasing as time goes on. Alternative paths will exist in the damaged network. Other nodes may replace the function of transfer station. Thus, even without giving priority to transfer stations, Z can continue to grow.

3.1.4. Sensitivity Analysis

(1) The Influence of Changing the Restoration Objective in Scenario I. The optimization of this part is based on the results of scenario I without considering uncertainty. After the restoration of scenario I, there are 8 stations are left behind and need to be repaired in Nanjing metro (ID: 8, 12, 13, 27, 28, 100, 101, 116). In this stage, the new objective is to make the Z recover to 100% and maximize the E. Besides, the pregiven deadline time and the budget limitation change to 35 days and 90000$, respectively. Other parameters are same as Table 3. The meaning of part is to illustrate the flexibility of the model. Decision-makers can make different objectives according to the recovery situation and remaining resources. Figure 12 demonstrates the optimal restoration schemes with different maintenance teams. All the nodes are repaired in those conditions, however with 3 teams compared to with 2 teams. The Z sees little improvement after the 22^nd day, the total recovery time shortens by 3 days, the resilience index improves by 3.51%, and the Z can recover to 100% three days earlier.

[figure omitted; refer to PDF]

(2) The Influence of the Maximized Budget for One Day in Scenario I. This part aims to analyze the influence of the maximized budget for one day (u) on the optimal restoration scheme. Assume there are three different $u$ ($2200, $3200, $4000) in scenario I, and the optimal sequence schemes of them be shown in Figure 13. Because the total budget is the same, except for the total recovery time and resilience index of the Z, other indicators do not change. The maximized budget for one day directly affects the number of synchronous maintenance teams, and when the teams increase from 2 to 3, the total recovery time is reduced by eight days. However, the resilience index has little change, with an increase of 6.34%. When the number of synchronous teams increases to 4, the total recovery time only shortens by three days and the resilience index only increases by 1.19%. Therefore, we can imply that increasing the synchronous number of maintenance teams blindly does not achieve the desired effects and may lead to the wasting of resources.

[figure omitted; refer to PDF]

(3) The Influence of the Total Budget in Scenario I. This part aims to analyze the influence of the total budget ($U$) on optimal restoration scheme. Assume there are three different total budget limitations ($40000, $55000, $70000) in scenario I, and the optimal sequence schemes can be shown in Table 6. All the indicators grow with total budget increase. When the number of repaired nodes equals 80% of the total number, the Z can recover to 100%. Examining budgets of $55,000 and $70,000, the E’s recovery ratio and Z’s resilience index increase by 8.57% and 21.07%, respectively, and the total recovery time is reduced by six days. It is worth mentioning that the scheme with $40,000 does not satisfy the resilience constraints, where the purpose is to reveal the sensitivity of the total budget. As the total budget affects the number of repaired nodes directly, when 50% of the total number nodes are repaired, the resilience index only is 0.362. However, when 60% of the total number nodes are repaired, the resilience index and recovery ratio of Z can meet the resilience constraints of objective, showing that the total budget is vital for recovery.

Table 6

The sensitivity analysis of the limitation of the total budget in scenario I.

Total budget limitation $U$	40000 ($)	55000 ($)	70000 ($)
The sequences of the repaired nodes	A: 24, 98, 9, 10	A: 24, 98, 9, 26	A: 24, 98, 9, 26, 27, 100
	B: 97, 25, 99	B: 97, 25, 99, 11	B: 97, 25, 99, 11, 13
	C: 118, 7, 117	C: 118, 7, 117, 10	C: 118, 7, 117, 10, 116
The total recovery time (day)	14	18	24
The percent of repaired nodes	50%	60%	80%
The total budget ($)	39000	52000	69000
The resilience index of the Z	0.362	0.503	0.609
The recovery ratio of Z	86.9%	92.8%	100%
The recovery ratio of E	76.91%	82.25%	89.3%

3.2. Scenario II: The Location of the Damaged Stations Is Aggregated in Zhengzhou Metro

3.2.1. Scenario Descriptions

We also set up a similar scenario to compare the results with different network structure. Until March 2021, Zhengzhou metro has six lines and 114 stations, shown in Figure 14. Zhengzhou Metro is a typical radioactive network with a ring line. Line 5 is a ring line. Twenty stations and four lines also are affected because of a natural disaster. There are eight stations damaged in ring line 5. The repair time distribution of damaged stations is similar to the Nanjing metro, shown in Table 7. It reflects the uncertainty of repair time.

[figure omitted; refer to PDF]

(2) The Influence of Different Network Structure Based on Zhengzhou Metro. This part aims to analyze the influence of the damaged stations’ structure. Comparing the damaged stations’ locations of scenario I and II, we find that scenario I’s transfer stations are closer to the unaffected stations. This scenario changes the damaged stations’ location based on the Zhengzhou metro, shown in Figure 18. It is more similar to scenario I than scenario II. Except for the network structure, other model setups are the same as scenario II.

[figure omitted; refer to PDF]

We contrast the optimal restoration sequence with transfer station priority and no transfer station priority in scenario II-I. The results can be shown in Figures 19 and 20. The red line’s resilience index is larger than the green line, increasing by 3.6%. This result is entirely different from scenario II but similar to the scenario I. In this scenario, giving transfer station priority is not the optimal strategy.

[figure omitted; refer to PDF]

Because the optimization model also considered the recovery process (resilience index), so the repair sequence is closely related to each node’s repair time. After many simulations, we can get the result of how much the variability of restoration time affects the variability of network accessibility. Taking the repair time uncertainty into account, Figure 24 shows the 150 optimal restoration sequence schemes’ network accessibility probability distribution when four teams work synchronously. It can be obtained by executing the program 150 times after randomly giving each node’s repair time according to the distribution function. The network accessibility’s resilience index range is (0.734–0.754). The mean value is 0.745, and the variance is 0.011. Under the budget limitation, the Z’s resilience index and node’s repair time are negatively correlated. The shorter the repair time of nodes, the more nodes are repaired, and the greater the resilience index of network accessibility. The resilience index’s mean value is approximately equal to the resilience index that does not consider repair uncertainty (0.744).

[figure omitted; refer to PDF]

Figure 25 illustrates the probability distribution of total recovery time. The shortest recovery time is 8 days, and the longest recovery time is 12 days. The average of recovery time is 10.2 days. Because each station’s repair time obeys normal distribution, the total recovery time obeys normal distribution too.

[figure omitted; refer to PDF]

3.3.3. Sensitivity Analysis

(1) The Influence of the Maximized Budget for One Day in Scenario III. If prolong the pre-given deadline time to 20 days, Figure 22 also illustrates the different network performance curves of schemes with different synchronous maintenance teams. This comparison is based on the damaged station repair time equal to their mean value, which means it does not consider repair time uncertainty. By contrasting the indicators of Table 10, it can be found that with the increase in the number of maintenance teams, the resilience index becomes bigger. By comparing the 4 teams and 2 teams, the total recovery time shortens by 50% and the resilience index improves by 36.51%. But the total recovery time of 3 teams and 4-2 teams is the same because the multiple teams repairing one station at the same time is not allowed. Interestingly, the stations’ repair sequence does not change. The station with a bigger node degree is repaired earlier. Firstly, the node 44 with degree 5 is repaired, then the set of nodes whose degree are 4 are repaired, the nodes with degree 3 and 2 are subsequently repaired, which means transfer stations with bigger node degree should be given priority in restoration to make network recover quickly and efficiently.

Table 10

The sensitive analysis of synchronous maintenance teams in scenario III.

Indicators	2 teams	3 teams	4-2 teams	4 teams
The sequences of the repaired nodes	A: 44, 14, 2, 23	A: 44, 98, 23	A: 44, 2, 23	A: 44, 2
		B: 26, 2, 108	B: 26, 5, 108	B: 26, 5
	B: 26, 98, 5, 108	C: 14, 5	C: 14	C: 14, 23
			D: 98	D: 98, 108
The percent of repaired nodes	100%	100%	100%	100%
The total budget ($)	40000	40000	40000	40000
The total recovery time (day)	20	15	15	10
The resilience index of Z	0.545	0.697	0.718	0.744
The recovery ratio of Z	100%	100%	100%	100%
The recovery ratio of E	100%	100%	100%	100%

(2) The Influence of Damaged Stations’ Repair Time in Scenario III. This part aims to check the sensitivity of damaged stations’ repair time distribution. Assume damage transfer stations’ repair time equal 6 days in scenario III, and other damaged stations’ repair time equal 3 days, the results shown in Figure 26. Those schemes’ resilience index is very similar, but the optimal sequence is not based on transfer station priority. The optimal restoration plan shown in Figure 27. In this case, giving transfer station priority is not the optimal strategy. It can prove that damaged stations’ repair time significantly influences the optimal sequence scheme. Repairing the transfer station firstly is unsuitable for all situations. The reason is that the transfer stations need more time to be repaired. If all transfer stations are repaired first, the early stage’s recovery will be prolonged and affect the whole stages’ recovery.

[figure omitted; refer to PDF][figure omitted; refer to PDF]

4. Conclusions and Recommendations

Given the ever-increasing concerns for metro network safety, based on resilience evaluation, a multi-objective restoration sequence optimization method is proposed. Taking budget limitations and repair time uncertainty into account, combined with the network topology characteristics, the optimal restoration sequence scheme aims to maximize E, when the Z satisfies pregiven resilience requirements. The resilience requirements include the recovery ratio and the resilience index that can reflect the network performance change in process. The proposed method is applied to the Nanjing and Zhengzhou Metro network. The benefits of the proposed model are as follows:

(1) The proposed method focuses on network topology characteristics and has a good handle on the node degree, network accessibility, and global average efficiency. The model presents a flexible framework. It can integrate multiple goals. The objectives can be set flexibly to the decision-makers preference. In the early emergency stage, decision-makers often want the network performance to recover to an acceptable level as early as possible, and they can change objectives based on the current recovery results after this objective is achieved.

(2) Unlike other traditional restoration sequence optimization models, the proposed method combines the resilience evaluation and optimization problem, ensuring that network performance is optimal in both process and result. Assuming the damaged stations’ repair time obeys different normal distributions and comparing the optimal sequence under different repair times, we can find that damaged stations’ repair time will impact the optimal restoration sequence. In some cases, the result is entirely different from the result based on importance ranking. More detailed conclusions are as follows:

(1) By comparing the Zhengzhou metro and Nanjing metro, the performance of a ring network is more excellent than a radial network under the same network attack. After a similar disaster, Zhengzhou’s Z decreased by 36.13%, but the Nanjing metro’s Z decreased by 74.97%. It proves that a ring network is more robust than a radial network.

(2) Many studies prove that transfer stations should be given priority in the restoration process. This conclusion may be invalid in the following scenario.

  When damaged stations’ location is relatively centralized, giving transfer station priority is not always the optimal strategy, especially when most transfer stations are far away from unaffected stations and they need more time to be repaired. It can be proved by Sections 3.1.2 and 3.2.3(2). Because even if they are repaired, they cannot play the role of network connection. The network accessibility recovery is prolonged in the early stage. It will affect the network accessibility resilience index because the metro resilience focuses on final results and the recovery process.

  In scenario I and II-I, the repair time of transfer stations is longer than non-transfer stations. They need more time to be repaired, so the network’s recovery is slow in the early time. By analysing the transfer stations’ location, even they are repaired, they cannot connect with other no-damaged stations. Thus, giving them priority has little influence on the whole network recovery. As time goes on, alternative paths will exist in the damaged network. Other nodes may replace the function of transfer station. Thus, even without giving priority to transfer stations, Z can continue to grow with accepted level. It can further prove that giving transfer station priority is not always the optimal strategy.

  In the case of the location of the damaged stations is scattered, for example in the scenario III, if the damaged stations’ repair time are the same, station with a bigger node degree should be given priority to be repaired. Because those damaged transfer stations are very close to unaffected stations, so if they are repaired, they can immediately play the role of connecting other nodes. However, if the damaged transfer stations repair time are far longer than other damaged stations, the advantage of giving the transfer stations priority becomes smaller. It can be proved by Section 3.3.3(2).

  Thus, the transfer station’s location and repair time together affect the optimal sequence scheme. In practice, we should shorten the repair time of the transfer station as much as possible. It can make the network recovery faster. But if the transfer stations’ repair time cannot change, the stations connected to the unaffected station need to be repaired firstly, even those stations have smaller node degree, then the stations with bigger node degree and connect with unaffected or repaired stations should be given priority, which can make the recovery more effectively.

(3) According to the sensitivity analysis, the total budget is more sensitive than the input of the one-day budget in the entire restoration phase. However, in the emergency phase, increasing the number of synchronous maintenance teams can effectively shorten recovery time. For example, in Section 3.3.2, when the number of maintenance teams increases from 2 to 4, the recovery time was cut by 42.31%. After the Z recovered to an improved level, increasing the number of synchronous maintenance teams does not result in significant change. For example, in Section 3.3.1, when the number of maintenance teams rises from 2 to 3, the resilience index only improves by 3.51%, and the total recovery time shorten by 8.57%. In summary, according to the available resources, a reasonable resource input plan should be set up to avoid delays or the wasting of resources.

The scheduling restoration is a complex process that may involve many factors, the paper only focuses on network topology characteristics, and the global average efficiency ignores the difference of the edge’s length. Except for the topology characteristics, there are other critical operational metrics, such as passengers demand, network flow, travel cost, etc. Those metrics can be integrated into the model in future work. For example, passengers demand and flow can be used to reflect the network performance. Due to the variability of passenger flow with time, the optimal restoration scheme may change with time and different form this paper proposed. Although the proposed model can provide a near-optimal scheme in the case study, not an exact solver, the restoration strategy can be used in similar scenarios. Finally, the applicability of restoration strategy needs to be further analyzed. For example, the damaged station’s distribution is relative scatter. Using the proposed method in different scenarios to generate corresponding restoration strategies and extend them to other similar scenarios, which is very meaningful to enhance the network resilience facing disruptions.

Acknowledgments

The authors would like to acknowledge the financial support for this study provided by the National Key Research and Development Program of China (Grant no. 2019YFB1600200) and the Postgraduate Research and Practice Innovation Program of Jiangsu Province, China (Grant no. KYCX20_0129).

References

[1] M. D’Lima, F. Medda, "A new measure of resilience: an application to the London underground," Transportation Research Part A: Policy and Practice, vol. 81, pp. 35-46, DOI: 10.1016/j.tra.2015.05.017, 2015.

[2] T.-Y. Liao, T.-Y. Hu, Y.-N. Ko, "A resilience optimization model for transportation networks under disasters," Natural Hazards, vol. 93 no. 1, pp. 469-489, DOI: 10.1007/s11069-018-3310-3, 2018.

[3] X. Wang, Y. Koç, S. Derrible, S. N. Ahmad, W. J. A. Pino, R. E. Kooij, "Multi-criteria robustness analysis of metro networks," Physica A: Statistical Mechanics and its Applications, vol. 474, pp. 19-31, DOI: 10.1016/j.physa.2017.01.072, 2017.

[4] S. Hosseini, K. Barker, J. E. Ramirez-Marquez, "A review of definitions and measures of system resilience," Reliability Engineering & System Safety, vol. 145, pp. 47-61, DOI: 10.1016/j.ress.2015.08.006, 2016.

[5] L.-G. Mattsson, E. Jenelius, "Vulnerability and resilience of transport systems-a discussion of recent research," Transportation Research Part A: Policy and Practice, vol. 81, pp. 16-34, DOI: 10.1016/j.tra.2015.06.002, 2015.

[6] A. W. Righi, T. A. Saurin, P. Wachs, "A systematic literature review of resilience engineering: research areas and a research agenda proposal," Reliability Engineering & System Safety, vol. 141, pp. 142-152, DOI: 10.1016/j.ress.2015.03.007, 2015.

[7] R. Das, "Approach for measuring transportation network resiliency: a case study on Dhaka, Bangladesh," Case Studies on Transport Policy, vol. 8 no. 2, pp. 586-592, DOI: 10.1016/j.cstp.2020.04.001, 2020.

[8] N. Y. Aydin, H. S. Duzgun, F. Wenzel, H. R. Heinimann, "Integration of stress testing with graph theory to assess the resilience of urban road networks under seismic hazards," Natural Hazards, vol. 1 no. 91, pp. 37-68, DOI: 10.1007/s11069-017-3112-z, 2018.

[9] W. Zhang, N. Wang, "Resilience-based risk mitigation for road networks," Structural Safety, vol. 62, pp. 57-65, DOI: 10.1016/j.strusafe.2016.06.003, 2016.

[10] W. Zhang, N. Wang, C. Nicholson, "Resilience-based post-disaster recovery strategies for road-bridge networks," Structure and Infrastructure Engineering, vol. 13 no. 11, pp. 1404-1413, DOI: 10.1080/15732479.2016.1271813, 2017.

[11] R. Faturechi, E. Miller-Hooks, "Travel time resilience of roadway networks under disaster," Transportation Research Part B: Methodological, vol. 70, pp. 47-64, DOI: 10.1016/j.trb.2014.08.007, 2014.

[12] L. Chang, F. Peng, Y. Ouyang, A. S. Elnashai, B. F. Spencer, "Bridge seismic retrofit program planning to maximize postearthquake transportation network capacity," Journal of Infrastructure Systems, vol. 2 no. 18, pp. 75-88, DOI: 10.1061/(asce)is.1943-555x.0000082, 2012.

[13] Q. Ye, S. V. Ukkusuri, "Resilience as an objective in the optimal reconstruction sequence for transportation networks," Journal of Transportation Safety & Security, vol. 7 no. 1, pp. 91-105, DOI: 10.1080/19439962.2014.907384, 2015.

[14] L. Zhang, J. Lu, B.-b. Fu, S.-b. Li, "A review and prospect for the complexity and resilience of urban public transit network based on complex network theory," Complexity, vol. 2018,DOI: 10.1155/2018/2156309, 2018.

[15] D.-m. Zhang, F. Du, H. Huang, F. Zhang, B. M. Ayyub, M. Beer, "Resiliency assessment of urban rail transit networks: shanghai metro as an example," Safety Science, vol. 106, pp. 230-243, DOI: 10.1016/j.ssci.2018.03.023, 2018.

[16] M. Zanin, X. Q. Sun, S. Wandelt, "Studying the topology of transportation systems through complex networks: handle with care," Journal of Advanced Transportation, vol. 8 no. 9,DOI: 10.1155/2018/3156137, 2018.

[17] T. C. Matisziw, A. T. Murray, T. H. Grubesic, "Strategic network restoration," Networks and Spatial Economics, vol. 10 no. 3, pp. 345-361, DOI: 10.1007/s11067-009-9123-x, 2010.

[18] K. Liu, C. Zhai, Y. Dong, "Optimal restoration schedules of transportation network considering resilience," Structure and Infrastructure Engineering, vol. 17 no. 8, pp. 1141-1154, DOI: 10.1080/15732479.2020.1801764, 2021.

[19] S. Somy, R. Shafaei, R. Ramezanian, "Resilience-based mathematical model to restore disrupted road-bridge transportation networks," Structure and Infrastructure Engineering,DOI: 10.1080/15732479.2021.1906711, 2021.

[20] A. Karamlou, P. Bocchini, "Optimal bridge restoration sequence for resilient transportation networks," Structures Congress, vol. 2014, pp. 1437-1447, DOI: 10.1061/9780784413357.127, 2014.

[21] G. Lu, Y. Xiong, C. Ding, Y. Wang, "An optimal schedule for urban road network repair based on the greedy algorithm," PLos One, vol. 11 no. 10,DOI: 10.1371/journal.pone.0164780, 2016.

[22] L. A. Pereira, S. Haffner, G. Nicol, T. F. Dias, "Multiobjective optimization of five-phase induction machines based on NSGA-II," IEEE Transactions on Industrial Electronics, vol. 64 no. 12, pp. 9844-9853, DOI: 10.1109/tie.2017.2701768, 2017.

[23] A. Decò, P. Bocchini, D. M. Frangopol, "A probabilistic approach for the prediction of seismic resilience of bridges," Earthquake Engineering & Structural Dynamics, vol. 42 no. 10, pp. 1469-1487, DOI: 10.1002/eqe.2282, 2013.

[24] Z. Li, C. Jin, P. Hu, C. Wang, "Resilience-based recovery strategy optimization in emergency recovery phase for transportation networks," System Engineering Theory and Practice, vol. 39 no. 11, pp. 2828-2841, 2019.

[25] Y. Yang, Y. Liu, M. Zhou, F. Li, C. Sun, "Robustness assessment of urban rail transit based on complex network theory: a case study of the Beijing subway," Safety Science, vol. 79, pp. 149-162, DOI: 10.1016/j.ssci.2015.06.006, 2015.

[26] V. Latora, M. Marchiori, "Efficient behavior of small-world networks," Physical Review Letters, vol. 87 no. 19,DOI: 10.1103/PhysRevLett.87.198701, 2001.

[27] T. P. Mo, H. H. Zhao, W. Mo, "Design and implementation of shortest travel path searching based on improved Dijkstra algorithm," Applied Mechanics and Materials, vol. 157-158, pp. 390-394, DOI: 10.4028/www.scientific.net/amm.157-158.390, 2012.