System Resilience Evaluation and Optimization Considering Epistemic Uncertainty

Abstract

Epistemic uncertainties, caused by data asymmetry and deficiencies, exist in resilience evaluation. Especially in the system design process, it is difficult to obtain enough data for system resilience evaluation and improvement. Mathematics methods, such as evidence theory and Bayesian theory, have been used in the resilience evaluation for systems with epistemic uncertainty. However, these methods are based on subjective information and may lead to an interval expansion problem in the calculation. Therefore, the problem of how to quantify epistemic uncertainty in the resilience evaluation is not well solved. In this paper, we propose a new resilience measure based on uncertainty theory, a new branch of mathematics that is viewed as appropriate for modeling epistemic uncertainty. In our method, resilience is defined as an uncertainty measure that is the belief degree of a system’s behavior after disruptions that can achieve the predetermined goal. Then, a resilience evaluation method is provided based on the operation law in uncertainty theory. To design a resilient system, an uncertain programming model is given, and a genetic algorithm is applied to find an optimal design to develop a resilient system with the minimal cost. Finally, road networks are used as a case study. The results show that our method can effectively reduce cost and ensure network resilience.

1. Introduction

In the engineering field, resilience reflects the ability of a system to withstand, adapt, and recover from disruptions (including both external disturbances and internal failures) [1]. If the impact of a disturbance cannot be effectively controlled, the disturbance may cause severe economic or other losses. For example, Hines et al. [2] pointed out that in the past 100 years, North America suffered power losses of as much as 186,000 MW due to natural disasters. Hurricane Isabel devastated the transportation system of the Hampton Roads, VA, region in 2003 and overwhelmed the emergency response system [3]. Besides these direct damages, MacKenzie et al. [4] analyzed the indirect production losses caused by disabled production facilities in the earthquake and tsunami that struck Japan on 11 March 2011. Hence, to reduce the direct and derived damages caused by disruption, more and more researchers begin to focus on system resilience evaluation and improvement. Resilience has become a hot research topic.

Many researchers have proposed a resilience definition from different views. Bruneau et al. [5] defined seismic resilience as the ability of systems to mitigate hazard, contain effects, and carry out recovery activities. Arcuri et al. [6] analyzed resilience on healthcare systems. They pointed out that healthcare system resilience has been linked to sustainability, safety, quality, and adaptive capacity. Cook and Long [7] studied the adaptive capacity in organization resilience. They proposed that the ability to borrow and adjust an adaptive capacity are features of resilience system and resilience engineering. Woods [8] proposed four senses of resilience including rebound, robustness, extensibility, and adaptability. Jain et al. [9] defined resilience as three abilities: avoidance, survival, and recovery. Pawar et al. [10] summarized resilience definitions and proposed that the resilience of industrial systems have three common themes, i.e., absorption, adaptation, and restoration. The same view can also be seen in Abbasnejadfard et al. [11].

To evaluate system resilience, researchers have proposed many measures, which can be categorized as deterministic and probabilistic measures. Deterministic resilience measures and compares a system’s behavior before and after a certain disruption, and typical deterministic resilience measure methods include the resilience loss [5], performance integral ratio-based measures [1,12,13,14], recovered performance ratio-based measures [15], and performance loss rate-based measures [16,17,18], etc. Probabilistic resilience measures consider the randomness of both a disturbance and the system response to it, and typical probabilistic resilience measures include deterministic resilience distribution-based measures [19], and performance degradation and recovery time-based measures [5,20,21], etc. Using the resilience measures described above, system resilience was evaluated using real data and simulations. Roach et al. (2018) [22] developed ten other metrics to evaluate the resilience of a water resources system with real data. Roberto and Patelli (2018) [23] simulated the perfor mance of a power grid system after disruptions and evaluated its resilience. Other researchers designed system resilience strategies such as protection strategies [24], reconfiguration strategies [25], and recovery strategies [26]. To design a resilient system, resilience optimization models and algorithms were also developed. For example, Li et al. [27] proposed a resilience optimization model to find the best recovery strategy for road networks after disruptions. Miller et al. [28] proposed a network resilience optimization model to find the optimal recovery strategy for a rail-based intermodal container network considering the uncertainty of traffic flow. Zhang et al. [29] proposed a resilience optimization model for road networks to find the best structure for a resilient road network. Salas and Yepes [30] proposed an improved planning support system to afford a planning alternative over the Spanish road network. Their method proved that decentralization optimization can improve the adaptive capacity of road networks and then improve resilience.

However, these resilience measures and their evaluation and optimization methods have problems when epistemic uncertainty existed. Epistemic uncertainty describes the subjective uncertainty caused by insufficient knowledge and information. In the system design process, no operation data can be obtained, and designers can only conduct tests or simulations to observe how disruptions affect a system. As is generally known, it is both time consuming and costly to conduct tests or simulations, especially for large-scale systems, so only limited data can be obtained. Due to this data insufficiency, the law of resilience cannot be obtained, which may lead to a poor knowledge of system resilience. In this situation, it is necessary to rely on empirical data to make decisions. When experts are investigated for empirical data, empirical data are asymmetrical due to the difference of knowledge between experts. To conclude, researchers face the problems of data insufficiency [31], poor knowledge of a system [32,33], and data asymmetry [34]. In this stage, the resilience evaluation for a new system is a typical problem with epistemic uncertainty that differs from traditional aleatory uncertainty [33,35]. Most existing resilience measures only consider aleatory uncertainty, i.e., the dispersion of both a disturbance and a system response to it. A few researchers have already studied epistemic uncertainty, and applied Bayesian theory [36] and evidence theory [37] to quantify the epistemic uncertainty of performance degradation and the disruption intensity in the resilience process, however, evidence theory will lead to the interval expansion problem for large systems [38], and Bayesian theory will bring subjective information into probability theory. There is no evidence to prove that subjective probability follows the operation law of probability theory [38,39,40]. Due to the reasons listed above, the problem of determining how to quantify epistemic uncertainty in resilience evaluation has not been well solved. The uncertainty theory proposed by Liu [41], a new branch of axiomatic mathematics, provides another solution for the epistemic uncertainty quantification problem, and this theory has been successfully applied in engineering. Hu et al. [42] applied uncertainty theory into risk assessment. Kang et al. [43] combined the reliability with uncertainty theory. Yang et al. [34] used uncertainty theory to measure epistemic uncertainty in spare parts transportation optimization. Li et al. [44] proposed a reliability analysis method for a muti-state deteriorating system based on belief reliability theory. These studies showed that compared with other methods for epistemic uncertainty, such as fuzzy theory, interval theory, and evidence theory, uncertainty theory is more applicable in epistemic quantification.

This paper proposes an uncertainty theory-based resilience measure and provides both resilience evaluation and optimization methods. Our main contributions are as follows:

• A new uncertainty theory-based resilience measure is proposed to quantify the epistemic uncertainty in resilience evaluation. Compared with other methods, our new resilience measure has a solid mathematical foundation in epistemic quantification;
• A resilience evaluation framework is provided for the new resilience measure. By building the performance model and the disruption response model, and obtaining the distribution function of uncertain variables, the distribution function of the disruption response can be calculated, and the system resilience can be evaluated;
• To build a resilient system with a minimum budget, an uncertain programming model is given and a genetic algorithm is applied to solve the optimization problem. A road network case verifies the effectiveness of our new model and algorithm.

The remainder of this paper is organized as follows. Section 2 gives a brief introduction of both resilience and uncertainty theory. Section 3 proposes our uncertainty theory-based resilience measure. Section 4 provides both resilience evaluation and optimization methods. Section 5 describes the application of two road networks as case studies to verify the effectiveness of our methods. Finally, concluding remarks are provided in Section 6.

2. Basic Concepts and Theories

2.1. System Response after Disruptions

After a disruption occurs, the system performance may degrade and then recover due to the implementation of recovery activities. This performance change process (i.e., the system response) after a disruption is shown in Figure 1. Some performance change processes are gradual and continuous, as Figure 1 shows, and some are discrete with sudden changes. In this figure, the system performed steadily at the beginning and suffered a disruption at time $t_0$. Then, the system performance started to decrease. After taking some recovery actions, the system performance gradually returns back. In some cases, the system may finally reach a better performance with resource reconfiguration, and also may get a worse performance as the system cannot recover completely. Considering that the ability of system to withstand, adapt, and recover from disruptions can be reflected in the performance changing process, many resilience measures are proposed based on the system’s response to the disruption.

2.2. Uncertainty Theory

Uncertainty theory is a new branch of axiomatic mathematics. In uncertainty theory, the belief degrees of events are quantified by defining the uncertain measure. Some basic concepts and results for uncertainty theory are presented as follows.

Definition 1

(Uncertain Measure [41]). Letting $\mathit{\Gamma }$ be a nonempty set and $ℒ$ be a $\sigma$-algebra over $\mathit{\Gamma }$, each element $\mathit{\Lambda }$ in $ℒ$ is an event. A set function $ℳ$ is called an uncertain measure if it satisfies the following axioms.

Axiom 1

(Normality Axiom). $ℳ(\mathit{\Gamma })=1$ for the universal set $\mathit{\Gamma }$.

Axiom 2

(Duality Axiom). $ℳ(\mathit{\Lambda })+ℳ\left({\mathit{\Lambda }}^{ c}\right)=1$ for any event $\mathit{\Lambda }$, where ${\mathit{\Lambda }}^{ c}$ is the complementary set of $\mathit{\Lambda }$.

Axiom 3

(Subadditivity Axiom). For every countable sequence of events ${\mathit{\Lambda }}_1,{\mathit{\Lambda }}_2,\cdots$, there is

$$ℳ\left\{{\displaystyle \underset{\infty }{\overset{i=1}{\cup }}{\mathit{\Lambda }}_{ i}}\right\}\le {\displaystyle \sum _{i=1}^{\infty }ℳ\left\{{\mathit{\Lambda }}_{ i}\right\}}.$$

(1)

Axiom 4

(Product Axiom [45]). Letting $\left({\mathit{\Gamma }}_k,{ℒ}_k,{ℳ}_k\right)$ be the uncertainty spaces for $k=1,2,\cdots$, there is

$$ℳ\left\{{\displaystyle \prod _{k=1}^{\infty }{\mathit{\Lambda }}_k}\right\}=\underset{k=1}{\overset{\infty }{{\displaystyle \wedge }}}ℳ\left\{{\mathit{\Lambda }}_k\right\}.$$

(2)

Definition 2

(Uncertain Variable [41]). An uncertain variable is a function $\xi$ from an uncertainty space $\left(\mathit{\Gamma },ℒ,ℳ\right)$ for the set of real numbers such that $\left\{\xi \in B\right\}$ is an event for any Borel set $B$ of real numbers.

Definition 3

(Uncertainty Distribution [41]). The uncertainty distribution $\mathsf{\Phi }$ of an uncertain variable $\xi$ is defined by

$$\mathsf{\Phi }\left(x\right)=ℳ\left\{\xi \le x\right\},$$

(3)

for any real number $x$.

Definition 4

(Regular Uncertainty Distribution [46]). An uncertainty distribution $\mathsf{\Phi }\left(x\right)$ is said to be regular if it is a continuous and strictly increasing function with respect to $x$ for which $0<\mathsf{\Phi }\left(x\right)<1$, and

$$\underset{x\to -\infty }{\lim }\mathsf{\Phi }\left(x\right)=0,\underset{x\to +\infty }{\lim }\mathsf{\Phi }\left(x\right)=1.$$

(4)

Definition 5

(Inverse Uncertainty Distribution [46]). Letting $\xi$ be an uncertain variable with a regular uncertainty distribution $\mathsf{\Phi }\left(x\right)$, then the inverse function ${\mathsf{\Phi }}^{-1}\left(x\right)$ is called the inverse uncertainty distribution of $\xi$.

Definition 6

(Independence [45]). The uncertain variables ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ are said to be independent if

$$ℳ\left\{{\displaystyle \underset{i=1}{\overset{n}{\cap }}\left({\xi }_i\in B_i\right)}\right\}=\underset{i=1}{\overset{n}{{\displaystyle \wedge }}}ℳ\left\{{\xi }_i\in B_i\right\},$$

(5)

for any Borel sets $B_1,B_2,\cdots ,B_n$ of real numbers.

Definition 7

(Strictly Monotone Function of Uncertain Variables [46]). A real-valued function $f\left(x_1,x_2,\cdots ,x_n\right)$ is said to be strictly monotone if it is strictly increasing with respect to $x_1,x_2,\cdots ,x_m$ and strictly decreasing with respect to $x_{m+1},x_{m+2},\cdots ,x_n$; that is

$$f\left(x_1,\cdots ,x_m,x_{m+1},\cdots x_n\right)\le f\left(y_1,\cdots ,y_m,y_{m+1},\cdots y_n\right),$$

(6)

whenever $x_i\le y_i$ for $i=1,2,\cdots ,m$ and $x_i\ge y_i$ for $i=m+1,m+2,\cdots ,n$ , and

$$f\left(x_1,\cdots ,x_m,x_{m+1},\cdots x_n\right)<f\left(y_1,\cdots ,y_m,y_{m+1},\cdots y_n\right),$$

(7)

whenever $x_i<y_i$ for $i=1,2,\cdots ,m$ and $x_i>y_i$ for $i=m+1,m+2,\cdots ,n$.

Theorem 1

(Operational Law: Inverse Distribution [46]). Letting ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ be independent uncertain variables with regular uncertainty distributions ${\mathsf{\Phi }}_1,{\mathsf{\Phi }}_2,\cdots ,{\mathsf{\Phi }}_n$ respectively, if $f\left({\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right)$ is a continuous and strictly increasing with respect to ${\xi }_1,{\xi }_2,\cdots ,{\xi }_m$ and strictly decreasing with respect to ${\xi }_{m+1},{\xi }_{m+2},\cdots ,{\xi }_n$, then

$$\xi =f\left({\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right),$$

(8)

has an inverse uncertainty distribution

$${\Psi }^{-1}\left(\alpha \right)=f\left({\mathsf{\Phi }}_1^{-1}\left(\alpha \right),\cdots ,{\mathsf{\Phi }}_m^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{m+1}^{-1}\left(1-\alpha \right),\cdots ,{\mathsf{\Phi }}_n^{-1}\left(1-\alpha \right)\right).$$

(9)

3. Uncertainty Theory-Based Resilience Measure

We first define the uncertainty theory-based resilience as follows.

Definition 8

(Resilience). Letting a system state variable $\xi$ be an uncertain variable that refers to the system’s state after disruption and letting $\Xi$ be the resilience domain of the system state, then the resilience is defined as the belief degree that the system’s state after disruption is within the resilience domain, i.e.,

$$R_E=ℳ\left\{\xi \in \Xi \right\}$$

(10)

where $\xi$ refers to the system’s state after disruption, and it can be described by the system’s behavior after the disruption.

This resilience measure is a general one, and it describes the chance that the system’s resilience is satisfactory for user requirements. Because most existing research uses a system’s performance change process to reflect the system’s behavior after disruption, we also define the state variable $\xi$ based on this. Let the system’s response to the disruption $r_e$ be a state variable and denote its threshold as $r_{eth}$. Then, the event $\xi \in \Xi$ can be recorded as $r_e>r_{eth}$, and the resilience can be measured as:

$$R_E=ℳ\left\{r_e>r_{eth}\right\}$$

(11)

In our paper, the average performance after disruption is selected to reflect the system’s response to the disruption event. It is called the disruption response and can be calculated as [1]

$$r_e=\frac{{\displaystyle {\int }_{t_0}^{t_0+T_a}Q\left(t\right)dt}}{T_a}$$

(12)

where $t_0$ is the time at which the disruption occurs, $T_a$ is the maximum allowable recovery time determined by users, and $Q\left(t\right)\left(0\le Q\left(t\right)\le 1\right)$ is the normalized performance of the system after the disruption. Using the min-max method [47], $Q\left(t\right)$ can be calculated as

$$Q\left(t\right)=\{\begin{array}{ll}\frac{P_{\max }-P\left(t\right)}{P_{\max }-P_0},\hfill & \mathrm{if} P \mathrm{is} \mathrm{STB}\hfill \\ \frac{P\left(t\right)-P_{\min }}{P_0-P_{\min }},\hfill & \mathrm{if} P \mathrm{is} \mathrm{LTB}\hfill \\ \min \left\{\frac{P\left(t\right)-P_{\min }}{P_0-P_{\min }},\frac{P_{\max }-P\left(t\right)}{P_{\max }-P_0}\right\},\hfill & \mathrm{if} P \mathrm{is} \mathrm{NTB}\hfill \end{array}$$

(13)

where $P_0$ is the performance of the system without disruption, and $P_{\min }$ and $P_{\max }$ are the minimum acceptable value and the maximum acceptable value of the system performance $P$, respectively. Once the performance exceeds $P_{\min }$ and $P_{\max }$, the system cannot be used completely. Here, LTB means the larger-the-better, and this is true for the system functions when $P>P_{\min }$. STB means the smaller-the-better, and this is true for the system functions when $P<P_{\max }$. NTB means the more nominal-the-better, and this is true for the system functions when $P_{\min }<P<P_{\max }$. The physical meaning of our normalized method is shown in Figure 2a. Using the performance normalization method in Equation (13), the performance in the perfect zone of Figure 2a are normalized as $Q\left(t\right)=100\%$, those in the fault zone are normalized as $Q\left(t\right)=0\%$, and those in the degraded zone are normalized as the percentage of the performance that the system can provide.

Figure 2b shows a normalized performance change process, and $r_e$ can be calculated using the ratio of the area of the dark blue zone to the areas of both the dark blue and the light blue zones. $r_e$ is a comprehensive index including robustness, adaption, and recovery capability. Its physical meaning is the average normalized performance in the maximum allowable recovery time after a disruption. The event $\left\{r_e>r_{eth}\right\}$ means that the system’s average performance after the disruption is within the resilience domain. In this way, our new resilience measure quantifies the belief degree of the system’s behavior meeting the requirements after disruptions.

According to uncertainty theory, $r_e$ is an uncertain variable. Letting $\mathsf{\Phi }$ be the distribution of $r_e$, our new resilience measure can be calculated as

$$R_E=ℳ\left\{r_e>r_{eth}\right\}=1-ℳ\left\{r_e\le r_{eth}\right\}=1-\mathsf{\Phi }\left(r_{eth}\right).$$

(14)

After the actual or experimental data are obtained, researchers can use uncertainty statistics to obtain the distribution function $\mathsf{\Phi }$ and calculate the resilience within the threshold $r_{eth}$.

4. Resilience Evaluation and Optimization

In this section, we provide a resilience evaluation framework first, and then we propose a resilience optimization method to assist with the system’s resilience design.

4.1. Resilience Evaluation

Based on the resilience definition given above, a resilience evaluation framework is proposed as shown in Figure 3. The resilience evaluation steps are as follows.

Step 1: The performance model is built to describe the system’s performance change over time. The general performance model is as follows:

$$P\left(t\right)=f\left(t,{\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right),$$

(15)

where $P\left(t\right)$ is the system performance at time $t$, and the function $f(\cdot )$ is the performance model. It is worth noting that in this paper we focus on the performance changing process of a system after disruption, and disruption modeling and generation will be our future work. Therefore, ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ are parameters that describe the performance characteristics such as the performance degradation and recovery time. These parameters differ in different resilience processes. In this paper, parameters ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ are uncertain variables. With performance changing process data after several disruption events, the distribution functions of ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ can be obtained by an uncertain statistic method. Here, we assume that ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ are independent uncertain variables. This limitation will be overcome in future work.

Step 2: The disruption response model is built to calculate $r_e$ that is defined in Equation (12). The general form of the disruption response model can be denoted as:

$$r_e=g\left({\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right),$$

(16)

where $g(\cdot )$ is the disruption response model. It can be calculated by substituting Equation (15) into Equations (12) and (13). According to the definition of $r_e$ in Equation (12), this model can be used to calculate the average value of a system’s normalized performance during a performance changing process. Here the disruption response model is a deterministic model. As described in Step 1, parameters ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ are uncertain variables. Therefore, the disruption response $r_e$ is also an uncertain variable. Its distribution function can be obtained based on the distribution functions of ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ and the operation law of uncertain theory.

Step 3: The uncertainty quantification is determined. As described in Section 3, resilience is defined as an uncertain measure, and it can be explained as the quantification of the uncertainty of a system’s disruption response. After building the disruption response model, the system resilience can be obtained based on the distribution function of the parameters in the model. Therefore, the uncertainty quantification can be divided into two sub-steps:

• The distribution functions of uncertain variables ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ are obtained. The distribution function of the parameters can be obtained by using the uncertain statistics method [46,48] or data form experts [49];
• The distribution function of $r_e$ is calculated. Using the distribution function of the uncertain variables ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$, the resilience can be evaluated using Equation (14). According to Theorem 1, the distribution function of $r_e$ can be calculated as follows.

It is assumed that (1) the disruption response function $g\left({\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right)$ is continuous and strictly increasing with respect to ${\xi }_1,{\xi }_2,\cdots ,{\xi }_m$ and strictly decreasing with respect to ${\xi }_{m+1},{\xi }_{m+2},\cdots ,{\xi }_n$, and (2) the independent uncertain variables ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ have the regular uncertainty distributions ${\mathsf{\Phi }}_1,{\mathsf{\Phi }}_2,\cdots ,{\mathsf{\Phi }}_n$. Then the disruption response $r_e$ has the inverse uncertainty distribution

$${\Psi }_{r_e}^{-1}\left(\alpha \right)=g\left({\mathsf{\Phi }}_1^{-1}\left(\alpha \right),\cdots ,{\mathsf{\Phi }}_m^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{m+1}^{-1}\left(1-\alpha \right),\cdots ,{\mathsf{\Phi }}_n^{-1}\left(1-\alpha \right)\right),$$

(17)

where $\alpha$ is the belief degree, and ${\Psi }_{r_e}^{-1}$ is the inverse distribution of the disruption response. The monotonicity of $g$ can be analyzed according to the physical meaning of parameters ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ and $r_e$.

Based on the inverse distribution of the disruption response in Equation (17), the 99-method [46] can be used to obtain the distribution function of the disruption response, and then the system resilience can be calculated.

4.2. Resilience Optimization

A resilience optimization method is proposed to find the optimal system design. In this research, resilience is defined as an uncertain measure, and the resilience-based optimization model is an uncertain programming model. In 2009, Liu [50] proposed the theory of uncertain programming, and provided the basic models for the uncertain programming and optimization methods. Based on conclusions from that research, the steps of our optimization method are as follows.

4.2.1. Optimization Model

Many resilience-based models have been provided in previous studies. For example, Zhang et al. [29] proposed a resilience-based network optimization model with the aims of finding the minimum cost design with the resilience constraint. Based on the same idea, our optimization model is as follows:

$$\begin{array}{l}\underset{\mathit{x}}{\min } cost\left(\mathit{x}\right)\\ \mathrm{s}.\mathrm{t}. ℳ\left\{r_e\left(\mathit{x},{\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right)>r_{eth}\right\}\ge \alpha \end{array}$$

(18)

where the function $cost\left(\mathit{x}\right)$ is the system design cost, $ℳ\{\cdot \}$ is the system resilience, $\alpha$ is the resilience requirement, and $\mathit{x}$ is the vector of the optimization variable, which is determined by the system design. Equation (18) only accounts for the cost and resilience of the system, and researchers can add other constraints as needed for practical problems.

4.2.2. Optimization Model Transformation

Liu [50] proved that the uncertain programming model could be transformed into a crisp mathematical programming model according to the operation laws in uncertainty theory. According to Theorem 1, the constraint function $ℳ\left\{r_e\left(\mathit{x},{\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right)>r_{eth}\right\}$ in Equation (18) is equivalent to crisp mathematical programming

$$\begin{array}{l}\underset{\mathit{x}}{\min } cost\left(\mathit{x}\right)\\ \mathrm{s}.\mathrm{t}. r_e\left({\mathsf{\Phi }}_1^{-1}\left(1-\alpha \right),\cdots ,{\mathsf{\Phi }}_m^{-1}\left(1-\alpha \right),{\mathsf{\Phi }}_{m+1}^{-1}\left(\alpha \right),\cdots ,{\mathsf{\Phi }}_n^{-1}\left(\alpha \right)\right)\ge r_{eth}\end{array}$$

(19)

where ${\mathsf{\Phi }}_1^{-1},{\mathsf{\Phi }}_2^{-1},\cdots ,{\mathsf{\Phi }}_n^{-1}$ are inverse uncertain distributions of ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$.

4.2.3. Optimization Method

To solve the crisp mathematical programming shown in Equation (19), common optimization algorithms such as the gradient descent method and the Hessian matrix method can be used. For some complex optimization problems, intelligent algorithms such as genetic algorithms and particle swarm optimization algorithms are better choices.

5. Case Study

We use two road network cases to demonstrate the application of our resilience evaluation and optimization method.

5.1. Case Introduction

Case 1: A road network with a simple structure is shown in Figure 4. This case is taken from Hiller et al. [51]. Nodes S and T are the source and destination nodes, respectively. In the figure, the numbers in the brackets represent the link number and the capacity. The maximum flow between Nodes S and T is chosen as the key performance parameter of the system. This model represents a determinate situation, in which the network topology is deterministic.

Case 2: Figure 5 illustrates another road network with a slightly more complex structure. This case is taken from Dai and Poh [52]. In this problem, all of the links represent potential roads. Designers need to find an optimal road network design that minimizes the costs within the constraint of the system resilience. Similar to Case 1, the maximum flow between Nodes S and T is selected as the key performance parameter. The numbers in the brackets represent the number, capacity, and cost of the links.

To simplify the problem, we make the following assumptions:

• Only one link is affected during each disruption;
• The capacity of the road links degrades suddenly after disruptions, and the capacity recovery processes are linear processes. This is a widely used assumption in resilience research [53]. In most disruptions, especially natural disasters and traffic accidents, the road link capacity degradation time is very short, and the capacity degradation time can be regarded as zero;
• The uncertain variable $C_d$ follows a linear uncertain distribution as follows

  $${\mathsf{\Phi }}_{C_d}\left(x\right)=\{\begin{array}{ll}0,\hfill & x\le 0\hfill \\ \frac{x}{C_i}\hfill & 0<x\le C_i\hfill \\ 1\hfill & x>C_i\hfill \end{array}$$

  (20)

  where ${\mathsf{\Phi }}_{C_d}$ is the distribution function of the maximum performance degradation $C_d$;
• The uncertain variable $T_r$ follows a lognormal uncertain distribution as follows

  $${\mathsf{\Phi }}_{T_r}\left(x\right)={\left(1+\exp \left(\frac{\pi \left(e-\ln x\right)}{\sqrt{3}\sigma }\right)\right)}^{-1},x\ge 0$$

  (21)

  where ${\mathsf{\Phi }}_{T_r}$ is the distribution function of the recovery time $T_r$, and $e$ and $\sigma$ are the expected value and the standard deviation of the normal uncertainty distribution, respectively.

5.2. Case 1: Road Network Resilience Evaluation

5.2.1. Resilience Evaluation

Using the resilience evaluation method proposed in Section 4.1, the resilience of the road network in Case 1 is evaluated as follows.

(1) The system performance model is built. In this case, the maximum flow between Nodes S and T is determined by the capacity of the road link. Therefore, the system performance model should be built based on the road link capacity disruption model and the correlation between the road link capacity and the maximum flow.

After a certain disruption, the capacity of road link i at time t can be expressed as

$$c_i\left(t\right)=h\left(C_i,C_d,T_r,t\right),$$

(22)

where $C_i$ is the capacity of road link $i$ in the normal state, $C_d$ is the maximum capacity degradation of link $i$, and $T_r$ is the recovery time of link $i$.

For linear capacity recovery processes, $g(\cdot )$ can be written as follows:

$$c_i\left(t\right)=\{\begin{array}{ll}\frac{t-T_r}{T_r}{C}_d+C_i,\hfill & 0\le t\le T_r\hfill \\ C_i,\hfill & t>T_r\hfill \end{array}$$

(23)

After the disruption, the maximum flow can be calculated according to the capacity of the road links as

$$F_{\max }\left(t\right)=f\left(C_1,\cdots ,C_{i-1},c_i\left(t\right),C_{i+1},\cdots ,C_n\right),$$

(24)

where $c_i\left(t\right)\le C_i$.

Obviously, the function in Equation (22) is a piecewise linear function with two sections, as shown in Figure 6. In this case, the disruption of link $i$ affects the maximum flow as follows:

$$F_{\max }\left(t\right)=\{\begin{array}{ll}{c}_i\left(t\right)+14-C_i^P,\hfill & 0\le c_i\left(t\right)<C_i^P\hfill \\ 14,\hfill & C_i^P\le c_i\left(t\right)\le C_i\hfill \end{array}$$

(25)

where $C_i^P$ is the inflection point of the piecewise linear function for link $i$, which is determined by the network structure. Figure 6 shows how the capacity of each link affects the maximum flow, and the values for the parameters $C_i^P$ and $C_i$ are shown in

Table 1. Parameter values in Figure 6.

$\mathbf{Link} \mathit{i}$	${\mathit{C}}_{\mathit{i}}^{\mathit{P}}$	${\mathit{C}}_{\mathit{i}}$
1	3	5
2	6	7
3	3	4
4	0	1
5	3	3
6	0	2
7	4	4
8	3	5
9	3	4
10	8	9
11	1	1
12	6	6

.

(2) The disruption response model is built. In this case, the maximum flow is an LTB parameter. According to the performance normalization method given in Equation (13), the normalized maximum flow can be calculated as:

$$Q\left(t\right)=\{\begin{array}{ll}\frac{c_i\left(t\right)-C_i^P}{14}+1,\hfill & 0\le c_i\left(t\right)<C_i^P\hfill \\ 1,\hfill & C_i^P\le c_i\left(t\right)\le C_i\hfill \end{array}$$

(26)

where $P_0=14$, $P_{\min }=0$.

Combining Equations (23) and (26), the normalized maximum flow of the network is as follows:

$$Q\left(t\right)=\{\begin{array}{ll}1,\hfill & C_d\le C_i-C_i^P\hfill \\ \frac{C_{d}t-\left(C_d-C_i-14+C_i^P\right)T_r}{14T_r},\hfill & 0\le t<\frac{C_i^P-C_i+C_d}{C_d}{T}_r \mathrm{and} C_d>C_i-C_i^P\hfill \\ 1,\hfill & \frac{C_i^P-C_i+C_d}{C_d}T\le t_r \mathrm{and} C_d>C_i-C_i^P\hfill \end{array}$$

(27)

Taking Link 1 as an example, the normalized maximum flow of the network when Link 1 is disrupted is as follows:

$$Q_1\left(t\right)=\{\begin{array}{ll}1,\hfill & C_d\le 2\hfill \\ \frac{C_d}{14T_r}t-\frac{C_d}{14}+\frac{8}{7},\hfill & 0\le t<\frac{C_d-2}{C_d}{T}_r \mathrm{and} C_d>2\hfill \\ 1,\hfill & \frac{C_d-2}{C_d}T\le t_r \mathrm{and} C_d>2\hfill \end{array}$$

(28)

According to Equation (12), the disruption response of the network can be calculated as follows:

$$r_e^1=\{\begin{array}{ll}1,\hfill & C_d\le 2\hfill \\ \frac{C_d{T}_a}{28T_r}-\frac{C_d}{14}+\frac{8}{7},\hfill & 0\le T_a<\frac{C_d-2}{C_d}{T}_r \mathrm{and} C_d>2\hfill \\ \frac{\left(C_d-2\right)T_r}{7C_d{T}_a}\left(\frac{1}{2}-\frac{C_d}{4}\right)+1,\hfill & T_a\ge \frac{C_d-2}{C_d}{T}_r \mathrm{and} C_d>2\hfill \end{array}$$

(29)

(3) The uncertainty is quantified. In this case, $C_d$ and $T_r$ are uncertain variables. Taking Link 1 as an example, let $C_i=5$, $e=2$, and $\sigma =1$, so the inverse distributions of the variables are as follows:

$${\mathsf{\Phi }}_{C_d}^{-1}\left(\alpha \right)=5\alpha ,$$

(30)

$${\mathsf{\Phi }}_{T_r}^{-1}\left(\alpha \right)=\exp \left(2-\frac{\sqrt{3}}{\pi }\ln \left(\frac{1}{\alpha }-1\right)\right),$$

(31)

where $\alpha$ is a belief degree between 0 and 1.

After the distribution functions of the uncertain variables are obtained, the inverse distribution function of the network disruption response $r_e$ can be computed using Equation (17). According to the physical meanings of $C_d$, $T_r$, and $r_e$, it is obvious that $r_e$ is continuous and strictly decreasing with respect to $C_d$ and $T_r$. Letting $T_a=10$, the inverse distribution of $r_e$ when Link 1 is disrupted can be expressed as follows:

$$\begin{array}{cc}\hfill {\Psi }_{r_e^1}^{-1}\left(\alpha \right)& =r_e^1\left({\mathsf{\Phi }}_{C_d}^{-1}\left(1-\alpha \right),{\mathsf{\Phi }}_{T_r}^{-1}\left(1-\alpha \right)\right)\hfill \\ & =\{\begin{array}{ll}1,\hfill & \alpha \ge 0.6\hfill \\ \frac{5T_a\left(1-\alpha \right)}{28\exp \left(2-\frac{\sqrt{3}}{\pi }\ln \left(\frac{1}{\left(1-\alpha \right)}-1\right)\right)}-\frac{5\left(1-\alpha \right)}{14}+\frac{8}{7},\hfill & \alpha <0.27\hfill \\ \frac{\left(3-5\alpha \right)\exp \left(2-\frac{\sqrt{3}}{\pi }\ln \left(\frac{1}{\left(1-\alpha \right)}-1\right)\right)}{35\left(1-\alpha \right)T_a}\left(\frac{1}{2}-\frac{5\left(1-\alpha \right)}{4}\right)+1,\hfill & 0.27\le \alpha <0.6\hfill \end{array}\end{array}$$

(32)

As Equation (32) shows, the inverse distribution function of $r_e^1$ is complex. Using the linear interpolation method, the distribution function of $r_e$ when Link 1 is disrupted can be obtained, as shown in Figure 7. According to Equation (14), letting $r_{eth}=0.9$, the resilience when Link 1 is disrupted can be calculated as

$$R_E^1=ℳ\left\{r_e^1>r_{eth}\right\}=1-ℳ\left\{r_e^1\le 0.9\right\}=1-\mathsf{\Phi }\left(0.9\right)=0.88.$$

(33)

From Figure 7, one can see that:

• The belief degree of the event $\left\{r_e^1\le 0.7857\right\}$ is 0. It is obvious that the value of the network resilience is at a minimum when Link 1 is totally interrupted and cannot recover. In this situation, the maximum flow after the disruption has a constant value of 11. $r_e^1$ can also be calculated using the ratio of 11 to 14, and is equal to 0.7857;
• The belief degree of the event $\left\{r_e^1<1\right\}$ is 0.6. According to the network structure, the maximum flow of the road network starts to decrease only when the capacity degradation of Link 1 exceeds two. According to Equation (20), the belief degree of the event $\left\{C_d>2\right\}$ is 0.6.

The two results are intuitive when Link 1 is disrupted, and the effectiveness of our method is verified.

5.2.2. Results and Discussion

Using the same method, the network resilience when each link is disturbed is evaluated. Letting $r_{eth}=0.9$, the system resiliencies when different links are disturbed are shown in

Table 2. Resilience of system.

Link	Resilience	Link	Resilience
1	0.879	7	0.66
2	0.644	8	0.879
3	0.833	9	0.833
4	1	10	0.565
5	0.746	11	1
6	1	12	0.557

.

As shown in Table 2, the network resilience is equal to one when Links 4, 6, and 11 are disturbed. Links 4 and 6 are redundant road links, and the system operates normally even if Links 4 and 6 fail completely. If Link 11 completely fails and cannot recover, $r_e^1$ is 0.9286 (ratio of maximum flow when Link 11 completely fails to the initial max-flow), i.e., $r_e^1$ is always larger than 0.9, and it is always in the resilience domain.

Figure 8 further illustrates how the system resilience $R_E$ varies with $r_{eth}$ when different links are disturbed.

To highlight the novelty of our method, we compared our resilience evaluation result with a probability-based measure. The probability-based measure is defined as

$$R_{E^P}^1=P\left\{r_e^1>r_{eth}\right\},$$

where $P$ is the probability measure, $R_{E^P}^1$ is the probability-based resilience measure when Link 1 is disrupted. $R_{E^P}^1$ can be calculated using the Monte Carlo method.

The resilience evaluation result is shown in Figure 9.

As is shown in Figure 9, the uncertainty theory-based resilience evaluation result is always less then probability-based resilience. This result shows the advantage of our new resilience measure. When epistemic uncertainty exists, our measure can provide a more conservative result, which is beneficial to decision makers. Similar results can also be found in [44].

5.3. Case 2: Road Network Resilience Optimization

5.3.1. Network Optimization

According to the resilience optimization method in Section 4.2, the road network in Case 2 is optimized as follows.

(1) The optimization model is built. In this case, the network cost is regarded as the optimization goal, and the cost function can be expressed as

$$cost={\displaystyle \sum _{i=1}^{30}{x}_i{l}_i},$$

(34)

where $cost$ is the cost of the network, $l_i$ is the cost of link $i$, $x_i$ is a binary decision variable, and

$$x_i=\{\begin{array}{ll}1,\hfill & \mathrm{if} \mathrm{link} i \mathrm{is} \mathrm{included} \mathrm{in} \mathrm{network}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}$$

(35)

Our network should provide enough transportation and network resilience. The requirement of the transportation amount can be expressed by an inequality constraint as follows:

$$F\left(\mathit{x}\right)\ge p,$$

(36)

where $F$ is the maximum flow of the network with a link composition vector $\mathit{x}$, which is a vector of binary decision variables, and $p$ is the minimum requirement of the maximum flow.

In this case, the system resilience is regarded as the minimum value of the system resilience when a single link is disrupted. This means that designers care about the worst situation. Ahmadian [54] applied the same resilience measure. The system resilience requirement can be expressed as follows:

$$\min \left\{R_E^1,R_E^2,\cdots ,R_E^{30}\right\}\ge \alpha ,$$

(37)

where $R_E^1,R_E^2,\cdots ,R_E^{30}$ are the system resilience when Links $1,2,\cdots ,30$ are disrupted. Letting $r_{eth}=0.9$, $R_E^i$ can be evaluated by using the resilience evaluation methods described in Section 4.1.

It is worth noting that it is probable that not all of the links are included in the optimal solution. However, Equation (37) is still applicable. The reason for this is that if link $i$ is not included in the optimal network, the network performance is not be influenced by the disruption to link $i$. $R_E^i$ is equal to one and does not affect the result of Equation (37).

The optimization model of this case is as follows:

$$\begin{array}{l}\underset{\mathit{x}}{\min }{\displaystyle \sum _{i=1}^{30}{x}_i{l}_i}\\ \mathrm{s}.\mathrm{t}. \min \left\{R_E^1,R_E^2,\cdots ,R_E^{30}\right\}\ge \alpha \\ F\left(\mathit{x}\right)\ge p\end{array}$$

(38)

(2) The optimization model is transformed. The resilience constraint in the optimization model of Equation (38) represents the fact that the minimum resilience of the system with different disruptions should be greater than the requirement value. This is equivalent to the event of the resilience of the system with each disruption being greater than the requirement value. Additionally, the resilience constraint is an uncertainty constraint. According to the optimization model transformation method discussed in Section 4.2, the resilience constraint can be transformed as follows

$$\begin{array}{c}{r}_e^1\left(\mathit{x},{\mathsf{\Phi }}_{C_d^1}^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{T_r}^{-1}\left(\alpha \right)\right)\ge r_{eth}\\ r_e^2\left(\mathit{x},{\mathsf{\Phi }}_{C_d^2}^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{T_r}^{-1}\left(\alpha \right)\right)\ge r_{eth}\\ ⋮\\ r_e^{30}\left(\mathit{x},{\mathsf{\Phi }}_{C_d^{30}}^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{T_r}^{-1}\left(\alpha \right)\right)\ge r_{eth}\end{array}$$

(39)

where $r_e^1,r_e^2,\cdots ,r_e^{30}$ are the disruption responses of the network when Links $1,2,\cdots ,30$ are disrupted, $C_d^1,C_d^2,\cdots ,C_d^{30}$ are the maximum performance degradations of Links $1,2,\cdots ,30$, and ${\mathsf{\Phi }}_{C_d^1}^{-1},{\mathsf{\Phi }}_{C_d^2}^{-1},\cdots ,{\mathsf{\Phi }}_{C_d^{30}}^{-1},{\mathsf{\Phi }}_{T_r}^{-1}$ are their inverse uncertainty distributions.

In this way, the resilience optimization model given in Equation (38) is equivalent to the crisp mathematical programming, as follows:

$$\begin{array}{cc}\hfill \underset{\mathit{x}}{\min }& {\displaystyle \sum _{i=1}^{30}{x}_i{l}_i}\hfill \\ \hfill \mathrm{s}.\mathrm{t}. & r_e^1\left(\mathit{x},{\mathsf{\Phi }}_{C_d^1}^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{T_r}^{-1}\left(\alpha \right)\right)\ge r_{eth}\hfill \\ & r_e^2\left(\mathit{x},{\mathsf{\Phi }}_{C_d^2}^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{T_r}^{-1}\left(\alpha \right)\right)\ge r_{eth}\hfill \\ & \hfill ⋮\hfill \\ & r_e^{30}\left(\mathit{x},{\mathsf{\Phi }}_{C_d^{30}}^{-1}\left(\alpha \right),{\mathsf{\Phi }}_{T_r}^{-1}\left(\alpha \right)\right)\ge r_{eth}\hfill \\ & F\left(\mathit{x}\right)\ge p\hfill \end{array}$$

(40)

(3) The optimization problem is solved. The network structure optimization problem is a 0–1 integer optimization problem, and the genetic algorithm, which is widely used in network resilience research [28,55,56], is used to solve the optimization problem given in Equation(40). In the genetic algorithm, every individual in the population is judged to be a feasible solution or not. However, in practice, $r_e$ varies with the network structure vector $\mathit{x}$ and the disruptions to different links. It is inefficient to derive the expressions for $r_e$ with different disruptions for every individual during optimization. To solve the problem, we propose a numerical algorithm to calculate $r_e$ in Equation (40). The Algorithm 1 is as follows.

Algorithm 1:$r_e^i$calculation

Input: The sample number N The link capacity $C_i$ The resilience constraint $\alpha$ The maximum allowable recovery time $T_a$ Output: The disruption response $r_e^i$ Step 1 let $T_r^i=\exp \left(2-\frac{\sqrt{3}}{\pi }\ln \left(\frac{1}{\alpha }-1\right)\right)$ $C_d^i=\alpha C_i$ $\Delta t=\frac{T_a}{N}$ Calculate the maximum flow $F_0$ when all links are in normal state. Step 2 For $j=0$ to $N$     $t_j=j\Delta t$     $c_i\left(t_j\right)=\{\begin{array}{ll}\frac{t_j-T_r}{T_r}{C}_d+C_i,\hfill & 0\le t_j\le T_r\hfill \\ C_i,\hfill & t_j>T_r\hfill \end{array}$     $C_i=c_i\left(t_j\right)$, and calculate the maximum flow $F\left(t_j\right)$ at time $t_j$.     $Q\left(t_j\right)=\frac{F\left(t_j\right)}{F_0}$ Endfor Step 3 $r_e^i=\frac{1}{T_a}{\displaystyle \sum _{j=1}^N\left[\frac{Q\left(t_{j-1}\right)+Q\left(t_j\right)}{2}\Delta t\right]}$

Using the genetic algorithm, the optimization problem can be solved. The optimization results are discussed in Section 5.3.2.

5.3.2. Results and Discussion

To solve the optimization problem, the related parameters of the optimization model and the genetic algorithm are set as shown in

Table 3. Parameters of optimization model and algorithm.

Parameter		Value
Optimization model	$\mathrm{Disruption} \mathrm{response} \mathrm{threshold} r_{eth}$	0.9
	$\mathrm{Resilience} \mathrm{constraint} \alpha$	0.5
	Max-flow constraint $p$	20
Genetic algorithm	Population NIIND	100
	Generation MAXGEN	100
	Generation gap GGAP	0.4
	Cross rate XOVR	1
	Mutation rate MUTR	0.1

.

Figure 10 presents the optimization process using the genetic algorithm. With the increase in the generations, both the average and the best fitness of the population increase and are close to a stable value. The evolution process of the optimal network structure in each generation is shown in Figure 11. The resiliencies, maximum flows, and costs of the networks in Figure 11 are shown in Figure 12. The results in Figure 11 and Figure 12 verify that our method can effectively minimize the network cost for both the maximum flow and the resilience constraints.

The individual with the best fitness in the last generation is the optimization result. The network structure of this individual is shown in Figure 11e. All of the information about the optimal network is shown in Figure 13.

The resilience evaluation results for the optimal network are provided in

Table 4. Resilience of optimal network for different disruptions.

Link	Resilience	Link	Resilience
1	0.5184	18	0.9486
3	0.5184	20	0.9744
6	0.9281	23	0.6188
8	0.9827	28	0.6864
13	0.5741	29	0.5184

. Figure 14 shows the resilience for different $r_{eth}$. The results prove that the resilience of the optimal network meets the requirements. The lowest resilience can be obtained when Link 1, Link 3, or Link 29 is disturbed.

6. Conclusions

In this study, the uncertainty theory is introduced to quantify the epistemic uncertainties in resilience research. A new uncertainty theory-based resilience measure is proposed based on a system’s response to disruption. Our new resilience measure has a solid mathematical foundation and a clear physical significance. It can be used to evaluate resilience and design a resilient system with insufficient data, such as the system design stage. To evaluate the system resilience, a resilience evaluation framework is proposed, and the framework contains three steps: (1) The system performance model is built; (2) the system’s disruption response model is built; and (3) the uncertainty of the variables in the system performance model and the uncertainty of the disruption response are quantified, and finally the system resilience is calculated. Furthermore, a resilience optimization method is proposed to design a resilient system with minimal costs. Based on the uncertain programming theory, a resilience-based optimization model and a model transmit method are proposed, and intelligent algorithms (e.g., genetic algorithms) can be used to solve this type of optimization problem. Two road networks are used in the case study to verify the effectiveness of our resilience evaluation and optimization method.

In the future, we will apply the uncertain statistics methods to the resilience evaluation and optimization problem. Without assumptions about the distributions of uncertain variables, these methods can be used in real data cases. This method will also be used in multi-scale systems to build a bottom-up resilience evaluation framework.