Zhiguo Zeng

Shijia Du
email: dushijia111@gmail.com

Yi Ding

Resilience Analysis of Multi-state Systems with Time-dependent Behaviors

Keywords: Resilience, extreme events, multi-state system, non-homogeneous semi-Markov process, accumulated reward

Resilience analysis of multi-state systems with time-dependent behaviors

distribution system. Li et al. [START_REF] Li | Minimax-regret robust co-optimization for enhancing the resilience of integrated power distribution and natural gas systems[END_REF] developed a minimax-regret robust co-optimization model for resilience assessment and design of integrated power distribution and natural gas systems. A Bayesian optimization method based on ranking and selection was proposed in [START_REF] Giahi | Design optimization for resilience for risk-averse firms[END_REF] for resilience assessment and design for riskaverse decision-makers.

In topology-based methods, resilience is modeled and analyzed based on topological models of the systems (usually in terms of network models). This type of model is often used in vulnerability analyses, which is related to the capability of the system to resist the disruptive events and remain resilient. For example, Bose et al. [START_REF] Bose | Vulnerability assessment of a power transmission network employing complex network theory in a resilience framework[END_REF] conducted a resilience and vulnerability analysis for an electrical network using the topology-based methods.

Topological-based evaluation, physics-based performance simulation, and multi-criteria decision analysis were combined into a unified framework to assess pre-event community resilience against extreme rainfall [START_REF] Yang | A federated pre-event community resilience approach for assessing physical and social sub-systems: An extreme rainfall case in hong kong[END_REF]. In Liu et al. [START_REF] Liu | Resilience analysis framework for interconnected critical infrastructures[END_REF], the topological model has combined with system dynamic models for resilience evaluation of interconnected gas and electricity networks.

Simulation-based methods use simulation to capture the uncertain behaviors involved in the resilience quantification. For example, To investigate the resilience of power systems against extreme weather events, Panteli et al. [START_REF] Panteli | Modeling and evaluating the resilience of critical electrical power infrastructure to extreme weather events[END_REF] developed a time-series Monte Carlo simulation method. Cardini et al. [START_REF] Cadini | A modeling and simulation framework for the reliability/availability assessment of a power transmission grid subject to cascading failures under extreme weather conditions[END_REF] developed a simulation model for resilience quantification, which allows analyzing simultaneously the normal system behavior and extreme weather events. Rocchetta et al. [START_REF] Rocchetta | A power-flow emulator approach for resilience assessment of repairable power grids subject to weather-induced failures and data deficiency[END_REF] proposed a simulation-based framework to evaluate the resilience of power grids subject to extreme weather-induced failures.

Most of the existing methods for resilience modeling and analysis, as reviewed above, assume that the performance of the system is continuous. A lot of practical engineering systems, however, are multi-state in nature or needs to be modeled by multi-state models to control the modeling/computational complexity [START_REF] Bao | A multi-state model for reliability assessment of integrated gas and power systems utilizing universal generating function techniques[END_REF]. A typical application of multi-state models is to use them for modeling the demands and capacities of energy generation systems [START_REF] Jiang | Optimal loading strategy for multi-state systems: Cumulative performance perspective[END_REF]. Tsai et al. [START_REF] Tsai | A multiple-state discrete-time markov chain model for estimating suspended sediment concentrations in open channel flow[END_REF] developed a multiple-state discrete-time Markov chain model to describe the suspended sediment concentration distributions in water with different depths. There are only a few works regarding quantifying resilience using multi-state models. Thekdi and Chatterjee [START_REF] Thekdi | Toward adaptive decision support for assessing infrastructure system resilience using hidden performance measures[END_REF] developed a multi-state hidden Markov model to quantify the resilience of infrastructure systems and discussed how to make decisions based on the model. Nuss et al. [START_REF] Nuss | Toward resilience as a tradable parameter during conceptual trade studies[END_REF] used a continuous-time discrete-state Markov model to quantify the resilience against disruptions and support system design selections considering the multi-state engineering system behaviors of uncrewed surface vessels. In previous research of the authors [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF], a Markov reward process-based model is developed for quantifying resistance, absorption, recovery and overall resilience. This work only applies to time-static systems, as the model is Markovian. The major difference between [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF] and the present work is that, the present work is based on non-homogeneous semi-Markov process model, and, therefore, can be applied to problems with time-dependents bahaviors.

All the works mentioned above assume Markov (or hidden Markov) models. A significant drawback of these Markov model-based methods is that they have to assume that the transition times among states follow exponential distributions. This assumption indicates that, if the state jump depends on the performance of some safety systems, these safety systems must have time-static behaviors, i.e., their failure probabilities cannot change over time. Such models, however, cannot adequately describe the time-dependent behaviors of a large variety of real-world engineering systems. Ferrario et al. [START_REF] Ferrario | Goal tree success tree-dynamic master logic diagram and monte carlo simulation for the safety and resilience assessment of a multistate system of systems[END_REF] used goal-tree-success-tree as a tool for quantifying August 17, 2020 DRAFT seismic resilience of infrastructure systems considering their multi-state behaviors and developed a Monte Carlo simulation framework to consider the non-exponentially distributed recovery times. However, their work did not consider the time-dependent behaviors of the safety barriers in the accident protection phase. In this paper, we develop a non-homogeneous Semi-Markov Reward Process (NHSMRP) model to quantify the resilience of multi-state systems, considering the time-dependent behaviors of both the pre-disruption protection processes and the post-disruption recovery processes.

Compared to the existing works, the major contributions of this paper can be summarized as follows:

1) A multi-state resilience model is developed based on NHSMRP.

2) Time-dependent behaviors in both pre-and post-disruption phases are considered.

3) An efficient algorithm is developed for resilience analysis.

The rest of the paper is organized as follows. Section II presents the developed NHSMRP-based resilience model. Section III develops an efficient simulation method for resilience analysis. Section IV and V applies the developed methods on real-world case studies of oil storage tank farms and reconfigurable computing systems, respectively. Finally, section VI concludes this paper with a discussion of future works.

II. A NON-HOMOGENEOUS SEMI-MARKOV REWARD PROCESS-BASED RESILIENCE MODEL

A. The model

Let S be a discrete (finite or countable) space (called state space hereafter) in which a discrete random variable

X n takes values. Let θ n be a continuous random variable taking values in [0, ∞) and let τ n = n i=0 θ n , n = 1, 2, • • • Physically, θ n is often used to model the inter-arrival time between the nth and (n -1)th state transition in X n , and τ n represents the accumulated time up to the nth state transition. A bivariate process {X n , τ n }, n = 0, 1, 2, • • • is called a non-homogeneous Renewal Process (NHRP) when the following assumptions hold [START_REF] Grabski | Semi-Markov processes: applications in system reliability and maintenance[END_REF]:

P r(X n+1 = j, τ n+1 -τ n ≤ t | X n = i, τ n = τ, X n-1 , τ n-1 , • • • , X 0 , τ 0 ) = P r(X n+1 = j, τ n+1 -τ n ≤ t | X n = i, τ n = τ ).
(1) and P r(X 0 = i, τ 0 = 0) = P r(X 0 = i).

(

Let

N (t) = sup n {τ n ≤ t}, n = 0, 1, • • • Then, a stochastic process {X(t), t ≥ 0} with piecewise constant
and right continuous sample paths given by X(t) = X N (t) is called a non-homogeneous Semi-Markov Process (NHSMP) associated with the NHRP {X n , τ n } [START_REF] Grabski | Semi-Markov processes: applications in system reliability and maintenance[END_REF]. It is called NHSMP as it is a semi-Markov process with reward structure. A semi-Markov process is called "semi", as it has Markov property only at instants when state transitions occur. In an NHSMP model, future system behaviors depend not only on the current state but also on the age of the system, which makes it an ideal tool to capture time-dependent system behaviors.

Let us assume that the performance of the system can be characterized by an NHSMP X(t) with n + 1 performance levels: X(t) ∈ S = {0, 1, 2, • • • , n}, as shown in Figure 1. Without loss of generality, we assume that the system performance decrease as the value of X(t) increases (i.e., X(t) = 0 indicates perfect performance while X(t) = n indicates the worst performance). Disruptive events might impair system performances and cause the system to jump from a higher performance level to a lower performance level. After the disruptions, August 17, 2020 DRAFT recovery measures can be taken to restore the system performances, resulting in backward jumps from lower performance levels to higher ones. The NHSMP can be determined by initial distribution and renewal kernel [START_REF] Grabski | Semi-Markov processes: applications in system reliability and maintenance[END_REF]. Initial distribution is a probability vector representing the distribution of states at t = 0. Let π 0 = [P r(X(0) = i) : i ∈ S] represent the initial probability vector of the NHSMP. In this paper, we assume that at t = 0, the system starts operation from the perfect performance state. Hence, we have

π 0 = [1, 0, • • • , 0]. Let Q(t, τ ) = [Q i,j (t, τ ) : i, j ∈ S]
represent the renewal kernel, where

Q i,j (t, τ ) = P r(X n+1 = j, θ n+1 ≤ t | X n = i, τ n = τ ). (3)
In Eq. ( 3), τ n represents the time when the nth state transition occurs; θ n+1 is the inter arrival time between the (n + 1)th and the nth state transition: θ n+1 = τ n+1 -τ n , where τ n+1 is the time when the (n + 1)th state transition occurs. The physical meaning of the renewal kernel Q i,j (t, τ ) is the joint probability that the next state is j and the time to next transition is no greater than t, given that the current state is i and that the age up to the nth transition is τ .

When rewards can be accumulated when system jumps between states and/or stays in a state, a NHSMP becomes a NHSMRP. In this paper, we use a NHSMRP to characterize the resilience of the system of interest, as shown in Figure 1, where the rewards represent the direct and indirect losses caused by disruptive events.

Direct losses are generated directly by disruptive events and do not depend on the length of the disruption.

The degraded system performances cause indirect losses before the system is fully recovered (e.g., revenue losses) [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]. Hence, indirect loss depends on the length of the recovery process. In this paper, reward rates d i,j

represents the direct loss suffered by the system when it degrades from state i to state j. It is easy to verify

that      d i,j > 0, if i < j, d i,j = 0, if i ≥ j, (4)
as i < j indicates performance degradation, while i ≥ j indicates recovery. Similarly, the indirect losses per unit of time of sojourn in the performance degraded state i, i = 1, 2, • • • , n are modeled by reward rate l i .

Please note that l 0 = 0 as state 0 represent perfect performance and there is no indirect loss for this state.

B. Resilience metrics

Achieving resilience requires the system to possess various capabilities. In the previous works, different papers might emphasize different aspects of these capabilities. For example, in the classical definition by Holling [START_REF] Holling | Resilience and stability of ecological systems[END_REF],

resilience is defined as "the persistence of systems and of their ability to absorb change and disturbance and still maintain the same relationships between populations or state variables". This definition mainly focuses on the system capability to resist and absorb the potential damages caused by disruptive events. Other definitions of resilience concentrate more on the recovery aspect. For example, Iervolino and Giorgio [START_REF] Iervolino | Stochastic modeling of recovery from seismic shocks[END_REF] defined seismic resilience as the characteristic of a system which "measures its capability to recover from a shock rapidly." Du et al. summarized the different capabilities required for resilience into three categories [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF]:

• resistant capability, which refers to the capability of the system to resist the impact of the disruptive event and remain normal operations [START_REF] Holling | Resilience and stability of ecological systems[END_REF];

• absorption capability, which refers to the capability of the system to absorb the influence of the disruptive event (possibly by degrading its performance) and remains resilient, so that the system can return to normal operation states when the disruptive event disappears [START_REF] Hollnagel | Resilience engineering: Concepts and precepts[END_REF];

• recovery capability, which refers to the capability of the system to quickly restore normal operation after the disruptive event disappears [START_REF] Pfefferbaum | Building resilience to mass trauma events[END_REF].

Four numerical metrics have been developed in [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF] to measure resilience based on the three resilience-related capabilities. However, as they were defined based on Markov reward process model, they only apply to timestatic problems. In this paper, we extended these numerical metrics to consider time-dependent system behaviors based on the NHSMRP model.

Suppose we are interested in the resilience behavior in (0, t). The first resilience metric is called resistant probability and is defined as [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF]:

p Rs (t) = P r (X(τ ) = X 0 , ∀τ ∈ (0, t)) , (5)
where X 0 is the system state with perfect performance. The physical meaning of p Rs (t) is the probability that the system successfully resisted the attack of disruptive events, so that its performance remains unaffected in (0, t). As shown in Figure 2, the system might encounter several disruptive events in (0, t). To resist these attacks, safety barriers need to be designed and the resistant probability depends on the performances of these safety systems. The resistant probability takes values in [0, 1] and a larger value of p Rs indicates better resistant capability.

The second resilience metric is related to the absorption capacity. It is defined based on the resilience limit of the system beyond which the system collapses since it can no longer absorb the damages. Let us first divide the state space S into two groups B 0 and B 1 . Group B 0 contains all the resilient states in which the damages of the disruptive event can be absorbed and the system is recoverable after the disruptive events disappears;

while B 1 contains all the non-resilient states where the system loses its resilience so that it can no longer be recovered after the disruptions. An example of the non-resilient states is the Fukushima nuclear accident. As the safety barriers fail to stop the initiating damages promptly, the damages propagate and become so severe that the power plant cannot be repaired and has to be abandoned. It is easy to see that in the NHSMRP model, the states in B 1 are absorbing states. The second resilience metric, called resilient probability (denoted by p Re (t)),

August 17, 2020 DRAFT Fig. 2: An illustration of the resilience metrics is defined as [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF]:

p Re (t) = P r (X(t) ∈ B 0 | X(τ ) > 0, ∃τ ∈ (0, t)) . (6)
The physical meaning of p Re (t) is the conditional probability that the system can absorb the damages from the disruptive events (performance degradation might be caused) and remain in resilient states in (0, t), given that initial damages were caused by the disruptive events, as shown in Figure 2. It takes values in [0, 1] and a larger value indicates better absorption capability.

The third resilience metric is related to recovery capability and defined based on recovery time. Let us introduce a random variable T Rc (t) to represent the accumulated time spent in recovering system performances for a system that is operating in (0, t). Let T th,Rc represent the maximal tolerable recovery time. Then, the recovery capability can be measured by the recovery probability p Rc (t) [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF]:

p Rc (t) = P r (T Rc (t) ≤ T th,Rc | X(τ ) > 0, ∃τ ∈ (0, t)) . (7)
It can be seen from Eq. ( 7) that the physical meaning of the recovery probability is that, given that damages have been caused by the disruptive events, the conditional probability that the system recovery time is within a acceptable threshold. The recovery probability takes values in [0, 1] and a larger value indicates better recovery capability.

Apart from measuring the three aspects of resilience separately, an integrated resilience metric can also be defined based on the total loss caused by disruptive event. Let L D (t), L ID (t) and L(t) denote the direct, indirect and total loss in (0, t), respectively. For the NHSMRP model defined in Figure 1, it is easy to show that the total loss is

L(t) = L D (t) + L ID (t) = n i=0 n j=0 d i,j • N i,j (t) + n i=0 l i • T i (t), (8)
where N i,j (t) denotes the number of state transition from state i to j that occur in (0, t); d i,j is defined in Eq. ( 4) and represents the direct loss associated with each state transition; l i represents the indirect loss per unit of sojourn time in state i (see Figure 1); while T i (t) represents the accumulated sojourn time in state i for a system operating in (0, t). The L(t) in Eq. ( 8) naturally integrates the resistant, absorption and recovery capabilities: the resistant and absorption capabilities affect the direct loss, while the recovery capability affects the indirect loss. A numerical metric of overall resilience is, then, defined based on L(t) [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF]:

Re(t) = P r(L(t) < L tol ), (9)
where Re(t) is called overall resilience and L tol is a predefined threshold value representing the maximal tolerable losses. An underlying assumption behind this definition is that the resilience objective of the system is that the total loss caused by the disruptive event should not exceed L tol . It is easy to see that Re(t) takes values in [0, 1] and that a larger value of Re(t) indicates better overall resilience. It has been shown in [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF] that the overall resilience defined in Eq. ( 9) is equivalent to the Cumulative Distribution Function (CDF) of the resilience triangle defined in [START_REF] Bruneau | A framework to quantitatively assess and enhance the seismic resilience of communities[END_REF] (see also Figure 2 for illustration). Out of the four resilience metrics, resistant, resilient, and recovery probability are individual measures of resistant, absorption and recovery capabilities, respectively. They are independent from one another. Overall resilience is an integrated resilience metric, which considers the joint contribution of resistant, absorption and recovery capabilities. Therefore, overall resilience depends on the three aforementioned resilience metrics.

III. RESILIENCE ANALYSIS BASED ON NON-HOMOGENEOUS SEMI-MARKOV REWARD PROCESSES

A. Procedures of applying the model Figure 3 summarizes the procedures for applying the developed model for resilience modeling and analysis.

The analysis starts with identifying the disruptive events that might threaten the system of interest. Typical disruptive events considered in the resilience analysis of energy systems include extreme weather, failure of components, natural disasters, etc. In this paper, we limit our analysis to the case where the system is subject to the threat of only one disruptive event. The developed methods can also be extended naturally to systems subject to multiple disruptive events. Homogeneous Poisson processes are widely used for modeling the occurrence of disruptive events like earthquakes, floods, hurricanes, etc [START_REF] Katz | Stochastic modeling of hurricane damage[END_REF]. As these existing works, we also assume that the occurrence of the disruptive event follows a homogeneous Poisson process with a rate λ D . The value of λ D can be estimated from historical data or based on expert judgment [START_REF] Baker | An introduction to probabilistic seismic hazard analysis[END_REF].

Various safety barriers protect the system from the disruptive event. Depending on the performance of these safety barriers, different consequences can result following the occurrence of a disruptive event. For example, when an earthquake occurs, safety barriers like NPP structure strengthening, emergency trip system, emergency cooling systems, etc., are activated to protect a NPP from severe consequence like core meltdown. Depending on whether these safety barriers successfully perform their designed function, different consequences can result, e.g., (in ascending order of severity), unaffected operation, operation interruption, core meltdown. In the second step of the analysis, a consequence analysis determines the possible consequences that might be caused by the disruptive event and the losses associated with each consequence. The losses include both direct and indirect losses: the former are determined directly by estimating the financial losses caused to the system by the disruptive event, while the latter is determined by estimating the cost per unit of time of sojourn in the performance degraded states. Each consequence is mapped into a state in the NHSMRP model (see Figure 1). The values of d i,j and l i , i, j = 0, 1, • • • , n are determined based on the estimated losses caused by the corresponding consequences.

Then, a probabilistic analysis calculates the occurrence probability of each consequence. When the performance of the safety barriers does not change with time, the analysis can be easily done by combining the Poisson process model with probabilistic combinational models like event tree [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF]. In this paper, we consider safety barriers with time-dependent performances, for which the existing models cannot be directly applied. In Sect. III-B, we develop an analytical approach to calculate the Probability Density Function (PDF) and CDF of the occurrence time for each consequence.

The next step is estimating the recovery time. The estimation can be done by collecting field recovery time data or recovery exercise data and estimate the CDF of the time needed to recover to a given performance level.

These CDFs, together with the distribution of the occurrence time for each consequence, are used to derive the renewal kernel matrix of the NHSMRP.

B. Determining renewal kernels based on time-dependent probabilistic analysis

Event tree is a widely used model in risk analysis. In an event tree, possible event combinations are enumerated through a tree structure to calculate the occurrence probability of the consequences, as shown in Figure 4. In Figure 4, the initiating event is the occurrence of the disruptive event, and each intermediate event corresponds

to the state of one safety barrier. In the traditional event tree analysis, the occurrence probabilities are constants.

In this paper, we extend the traditional methods to consider time-dependent behaviors arising from the arrival process of the disruptive events and the degradation/repair behavior of the safety barriers.

Fig. 4: An illustration of an event tree

Let us consider a system with age τ ; that is, at t = 0, the system has been operating for τ units of time.

Let T C,i represent the time when the ith consequence occurs for the first time after τ . As T C,i is a random variable, we are interested in its PDF and CDF. First, let us consider an arbitrary time x, 0 < x < t. From the event tree model, it is easy to see that for dx → 0, we have:

P r(x < T C,i ≤ x + dx) = P r(x < T D ≤ x + ∆x) • n C i j=1 n E,i,j k=1 P r(X SB,k (τ + x) = q) ( 10
)
where T D represents the occurrence time of the disruptive event; n Ci is the number of event sequences leading to the ith consequence; n E,i,j is the number of events in the jth event sequence leading to the ith consequence;

X SB,k (τ + x) represents the state of the kth safety barrier at the current time and the value of q depends on the state of the corresponding safety barrier in the event sequence (from the event tree), where q = 0 indicates that the safety barrier cannot perform its intended function while q = 1 indicates otherwise. For non-repairable safety barriers, P r(X SB,k (τ + x) = 1) refers to its reliability at τ + x; while for repairable safety barriers, P r(X SB,k (τ + x) = 1) is its instantaneous availability at τ + x. Reliability/availability analysis can be used to calculate these values [35].

As we assumed that the arrival of the disruptive event follows a homogeneous Poisson process with rate λ D , T D is exponentially distributed [START_REF] Grabski | Semi-Markov processes: applications in system reliability and maintenance[END_REF]:

P r(x < T D ≤ x + dx) = λ D e -λ D x dx. (11)
Eqs. ( 10) and ( 11) define the PDF of T C,i , from which its CDF can be calculated:

P r(T C,i < t) = t 0 λ D • e -λ D x • n C i j=1 n E,i,j k=1 P r(X SB,k (τ + x) = q)dx = n C i j=1 t 0 λ D • e -λ D x • n E,i,j k=1 P r(X SB,k (τ + x) = q)dx. (12)
August 17, 2020 DRAFT Let η i,j (τ ) denote the time for a system aged τ to make a transition from state i to state j, irrespective of the transition times to the other states. Let f ηi,j (t, τ ) and F ηi,j (t, τ ) represent the PDF and CDF of η i,j (τ ), respectively. It is easy to see that the state transitions i → j (i < j) are results of disruptive events. Hence, the distribution of η i,j (τ ), i < j can be determined based on Eqs. [START_REF] Zhang | Improving the resilience of postdisaster water distribution systems using dynamic optimization framework[END_REF] and [START_REF] Li | Minimax-regret robust co-optimization for enhancing the resilience of integrated power distribution and natural gas systems[END_REF]. The state transitions i →, i > j are caused by system recovery. The distributions of η i,j (τ ), i > j can, therefore, be determined based on the associated time-to-repair distributions. Then, from Eq. ( 3), the elements of the renewal kernel of the NHSMP can be derived:

Q i,j (t, τ ) = t 0 f ηi,j (x, τ ) k∈S,k =j P r (η i,k (τ ) > x | η i,j = x) dx = t 0 f ηi,j (x, τ ) k∈S,k =j 1 -F η i,k | ηi,j (x, τ ) dx. (13)
In Eq. ( 13),

F η i,k | ηi,j (x, τ ) is the conditional CDF of η i,k (x, τ ) given the value of η i,j . When η i,k is independent from η i,j , F η i,k | ηi,j (x, τ ) = F η i,k (x, τ ).

C. Efficient Monte Carlo simulation for resilience analysis

Resilience analysis of the developed NHSMRP model requires analyzing its time-dependent behavior. Although there are a few existing approaches for transient analysis of an NHSMRP model, a common challenge is the high computational burden, especially when the state space of the NHSMRP gets large. In this section, we develop an efficient Monte Carlo simulation algorithm to improve the computational performance of resilience analysis. The performance of this algorithm will be discussed through both theoretical analyses and numerical experiments in the next section.

The algorithm is based on the embedded chain and holding time distributions of the NHSMP. Let τ i , i = 1, 2, • • • be the time instant when the ith system state transition (X i ) occurs. For an NHSMP described in Figure 1, it is well-known that X i follows a non-homogeneous discrete time discrete state Markov chain with a initial distribution π 0 and a transition probability matrix P (τ ) = [p i,j (τ ) : i, j ∈ S], where p i,j (τ ) is given by [START_REF] Grabski | Semi-Markov processes: applications in system reliability and maintenance[END_REF]:

p i,j (τ ) = lim t→∞ Q i,j (t, τ ). (14)
The discrete time Markov chain

X i , i = 0, 1, • • • is called the embedded chain of the NHSMP. Let F i,j (t, τ )
denote the CDF of the holding time T i,j , which is defined by:

F i,j (t, τ ) = P r(T i,j ≤ t) = P r(θ n+1 ≤ t | X n = i, X n+1 = j, τ n = τ ), i, j ∈ S, t ≥ 0. (15)
In Eq. ( 15), θ n+1 denote the interarrival time between the n and (n + 1)th state transition, and τ n is the accumulated time up to the nth state transition. It is easy to see that [START_REF] Grabski | Semi-Markov processes: applications in system reliability and maintenance[END_REF]

F i,j (t, τ ) = P r(θ n+1 ≤ t, X n = i, X n+1 = j, τ n = τ ) P r(X n = i, X n+1 = j, τ n = τ ) = P r(θ n+1 ≤ t, X n = i | X n+1 = j, τ n = τ ) P r(X n = i | X n+1 = j, τ n = τ ) = Q i,j (t, τ ) p i,j (τ ) . ( 16
)
The PDF of T i,j can also be derived:

f i,j (t, τ ) = dF i,j (t, τ ) dt = 1 p i,j (τ ) dQ i,j (t, τ ) dt = 1 p i,j (τ ) f ηi,j (x, τ ) k∈S,k =j 1 -F η i,k | ηi,j (x, τ ) . (17)
Having derived p i,j (τ ) and F i,j (t, τ ), The NHSMRP can be easily simulated for resilience analysis. In practice, however, calculating p i,j (τ ) through Eq. ( 14) often requires numerical integration, as the integration in Eq. ( 13) is often too complicated to have analytical solutions. The numerical integration could bring very high computational burdens, as it has to be evaluated for each generated sample. In the developed algorithm, we attempt to solve this problem by introducing a linear interpolation model (Algorithm 1 -Algorithm 3). The developed algorithms comprise of two phases. In the training phase, a linear interpolation model is trained based on n tr training samples and used to approximate the transition probability matrix P (τ ) of the embedded chain

X i , i = 0, 1, • • •
The trained linear interpolation model is, then, used in the simulation phase to generate state jumps. In this way, numerical integration of Eq. ( 14) can be avoided in the simulation phase. When n tr is large, the linear interpolation model can approximate P (τ ) fairly accurately. At the same time, the computational costs of running the linear interpolation model are much less than directly doing a numerical integration. Since P (τ )

needs to be evaluated for each generated state jump, the linear interpolation model can greatly improve the computational efficiency of the algorithm. In general, a larger value of n tr would have better approximation accuracy, but paying the price of increasing computational costs. In practice, a trade-off needs to be made to balance the accuracy and computational costs. Another major difference between the developed algorithms and the traditional Monte Carlo simulation is that, in the developed algorithms, samples are generated using the idea of vectorization. Instead of generating state jumps interactively using loops, vectorized functions are designed that take vector inputs and generate vectors of the state jumps directly. Vectorization can greatly reduce the overhead cost of the program, and, therefore, dramatically improve the computational efficiency of the algorithm. Algorithm 1 -Algorithm 3 also apply to Markov reward models. In this case, only p i,j (τ ) and F i,j (t, τ ) need to be replaced by their counterparts in the Markov reward models, while the algorithms can remain unchanged.

In Algorithm 1, T represent the evaluation horizon; x cur and x next represent the current and the next state, respectively; t cur and t next present the time instant for the current and the next jump, respectively; θ represents the holding times. The logical operations on vectors, e.g., θ > 0, returns a logical index vector whose element is 1 if the the logical operation on the corresponding element of θ is true (and 0 otherwise). In programming languages like Matlab, the logical index vector can be directly used to access the corresponding elements in the original vector. For example, θ( θ > 0) returns the positive elements in θ.

Algorithm 2 is vectorized generations of the inverse transform sampling method for generating random numbers [START_REF] Ross | Introduction to probability models[END_REF]. The inputs x cur is the current state; τ is a column vector that contains the ages of the samples; pxcur,j , j = 1, 2

t next ( x cur == x cur ) ← τ + θ( x cur == x cur ); t next ( t next > T ) ← T ; x next ( t next > T ) ← x cur ( t next > T ) ; L( x cur == x cur ) ←
P = [p i,j ] ← pxcur,j ( τ (i)), i = 1, 2, • • • , length( τ ), j = 0, 1, • • • , n ;
C ← Calculate the column-wised cumulative sum of the matrix P ;

u ← Generate a column vector of random numbers with the same size as τ from U [0, 1];

x next ← For each row in C, find the index of the first element larger than the element in the same row of u;

Algorithm 3 is a vectorized generation of the acceptance rejection sampling method [START_REF] Ross | Introduction to probability models[END_REF]. The inputs f i,j (x, τ )

is the PDF of the holding time, as defined in Eq. ( 17); g i,j (x, τ ) is a proposal density function which is easier to simulate; c(τ ) is a function that satisfies f i,j (x, τ ) ≤ c(τ ) • g i,j (x, τ ), ∀x and τ. The operators . * and ./ represent multiplication and division operators for vectors, i.e., apply the corresponding operations on each element of the vectors.

August 17, 2020 DRAFT Algorithm 3: Generate holding times using vectorized acceptance-rejection method. input : f i,j (x, τ ), g i,j (x, τ ), c i,j (τ ), τ output: θ while τ is not empty do u ← Generate a column vector of random numbers with the same size as τ from U [0, 1];

θ ← Generate a column vector of random numbers with the same size as τ from g i,j (x, τ );

Keep
O d = O sim + O oc = O(n s nt n st )O p + O(n s nt )O ht + O(n t )O oh , (18)
where n s is the sample size of the simulation, nt is the average number of state jumps in (0, t), n st is the number of states, O p is the computational complexity of evaluating one element of Eq. ( 14), O ht is the computational complexity of generating one sample from the hitting time distribution, and O oh represents the overhead cost of per loop.

In the literature, there are two types of methods for analyzing the behaviors of Markov/semi-Markov reward models: Monte Carlo simulations [START_REF] Masala | Survival probabilities for hiv infected patients through semi-markov processes[END_REF] and numerical integrations [START_REF] Janssen | Numerical solution of non-homogeneous semi-markov processes in transient case[END_REF][START_REF] Tijms | A fast algorithm for the transient reward distribution in continuous-time markov chains[END_REF]. Table I lists some of the most widely used methods in the literature and discusses their computational complexities. It can be seen from the Table that, compared to the numerical integration methods, a significant benefit of the developed method is that its computational complexity only grows linearly with the state numbers, while the numerical integration methods have at least quadratic growth rates. Hence, the computational performance of the developed method would be significantly better than the numerical integration methods for problems with large state space. Compared to the standard Monte Carlo simulation, the computational performance of the developed method outperforms in the two aspects: first, thanks to the linear interpolation model, the O p of the developed methods is much lower; second, the overhead cost is also much lower in the developed method.

2) Numerical experiment -I: Two numerical studies are conducted to compares the computational performances. The first is based on a homogeneous Markov reward process model (a special case of the NHSMP) described in Sect. 4.2 of [START_REF] Cloth | Five performability algorithms: A comparison[END_REF]. The data used in this experiment are from Table 4.2 of [START_REF] Cloth | Five performability algorithms: A comparison[END_REF]. Let Y t denote the accumulated reward at t. As in [START_REF] Cloth | Five performability algorithms: A comparison[END_REF], we calculate P (Y 10 ≤ 10) under different sizes of state spaces (n st ). Three methods are tested: the developed method, Tijms' method [START_REF] Tijms | A fast algorithm for the transient reward distribution in continuous-time markov chains[END_REF] and standard Monte Carlo simulation [START_REF] Masala | Survival probabilities for hiv infected patients through semi-markov processes[END_REF]. The parameter values, i.e., the step size ∆ in the Tijms' method and the sample size n s in the other two methods, are set to be ∆ = 0.01, n s = 10 6 . The reason for setting parameters like this is that it can ensure the three methods achieve the same degree of errors.

The computational experiment has been carried out on a computer with a CPU of 2.59 GHz (12 cores) and 64 GB RAM. Table II presents the results of the numerical experiment. In Table II, M is a parameter that defines the problem and controls the size of the state space (see [START_REF] Cloth | Five performability algorithms: A comparison[END_REF] for details), and n st is the number of states of the derived model. Theoretical values for M = 4, 5, 6, 7 are provided in [START_REF] Cloth | Five performability algorithms: A comparison[END_REF] through high-precision numerical integrations and are used in this study as reference values. As in [START_REF] Cloth | Five performability algorithms: A comparison[END_REF], we keep increasing the state space until the running time of the algorithm exceeds 300 (seconds) with ∆ = 0.01, n s = 10 6 . Figure 5 shows the comparisons of the running times in log-scales. It can be seen from the comparisons that when the state space is small, the Tijms' method performs better than the developed methods. While as the state space increases, the developed method outperforms. This is because, as discussed in the theoretical analysis, the computation time of the numerical integration methods increases quadratically with the state space, while the developed method has a linear increasing rate. The developed method performs always better than the standard Monte Carlo simulations. This is due to the benefits from the reduced overhead costs through vectorization, and also from the reduced cost of simulating next jumps through the use of linear interpolation. 3) Numerical experiment -II: A second numerical experiment is conducted to compare the performances of developed method with Janssen's method [START_REF] Janssen | Numerical solution of non-homogeneous semi-markov processes in transient case[END_REF] and standard Monte Carlo simulation [START_REF] Masala | Survival probabilities for hiv infected patients through semi-markov processes[END_REF] on an NHSMRP model. We design a simple case study following the same protocol as [START_REF] Janssen | Numerical solution of non-homogeneous semi-markov processes in transient case[END_REF]. A three-state system is defined to describe the event tree model in Figure 4 and used in this case study. The arrival of the initiating event is assumed to be a Poisson process with rate λ and p II is assumed to be a time-dependent function whose values are given by the inverse CDF of a Weibull distribution with parameters η f and β f . The recovery time August 17, 2020 DRAFT III. The three methods are used to calculate the state probabilities at t = 1000 : P r(X(1000) = i), i = 0, 1, 2.

The computational experiment has been carried out on a computer with a CPU of 2.59 GHz (12 cores) and 64 GB RAM. The results and computation times are compared in Table IV. For the two Monte Carlo simulation methods, the two-sided confidence intervals with a confidence level α = 0.05 are also given. It can be seen from the comparison that even for a relatively small-scale NHSMP, the developed efficient Monte Carlo simulation performs much better than the numerical integration method. As shown in the previous section, when the sample size grows larger, the computational benefits of the developed method over the numerical integration methods would become more significant. This is because, since the system is nonhomogeneous, considering a given time t requires considering all the history up to t, which makes numerical integration method computationally expensive. Compared to the standard Monte Carlo simulation, the developed method can achieve significant performance improvement. This confirms our theoretical analysis that the developed method improves the standard Monte Carlo simulation by reducing the computational costs in evaluating the transition probability matrix of the embedded chain through linear interpolation, and also by reducing the overhead cost through vectorization.

IV. APPLICATION ON CRUDE OIL STORAGE TANKS

A. System descriptions

A crude oil storage tank farm, initially reported in [START_REF] Wu | Quantitative risk assessment of fire accidents of large-scale oil tanks triggered by lightning[END_REF] all the failure probabilities of safety barriers were assumed to be constants; while in the present work, timedependent failure probabilities are considered, thanks to the non-homogeneous semi-Markov reward model.

The tank farm comprises of ten External Floating Roof (EFR) storage tanks. Each tank has a storage capacity of 1 × 10 5 (m 3 ), and the tank farm has a total storage capacity of 1 × 10 6 (m 3 ) [START_REF] Wu | Quantitative risk assessment of fire accidents of large-scale oil tanks triggered by lightning[END_REF]. Lightning is the most frequently encountered threat to oil storage tanks [START_REF] Chang | A study of storage tank accidents[END_REF]. Therefore, we consider the resilience of the oil storage tank against the threat of lightning.

When lightning hits a storage tank, it might catch fire due to two possible mechanisms [START_REF] Chang | A study of storage tank accidents[END_REF]. The first is the ignition of the leaked oil due to the puncturing of the tank wall by the lightning. The second is the ignition of flammable vapors that exists inside the tank by the lightning. Safety barriers are often designed to protect the oil tank and contain the propagation of the fire if it occurs. In this paper, we consider four safety barriers, i.e.,

Lightning Protection Mast (LPM), Automatic Rim seal fire Extinguishing System (ARES), Fixed Foam System (FFS) and fire brigade. Detailed descriptions can be found in Table 1 of [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]. If a fire is caused by lightning but not contained promptly by the fire extinguishing system or the fire brigade, it could spread to neighboring tanks (domino effect) and cause severe damages to the tank farm and also social well beings [START_REF] Necci | Accident scenarios triggered by lightning strike on atmospheric storage tanks[END_REF].

B. Resilience modelling 1) Occurrence of the disruptive event:

As in [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF], we assume that the lightning that hits one storage tank arrives following a homogeneous Poisson process with a rate λ L . In this paper, we use the estimated value of λ L from [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]: λ L = 9.842 × 10 -4 (y -1 ). Let n tank = 10 be the number of storage tanks and let us consider an initiating event as at least one tank hit by the lightning. The initiating event can, then, be modeled as a homogeneous Poisson process with a rate λ D = n tank • λ L [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF].

2) Consequences and losses: Depending on the performance of the safety barriers, different consequences might be caused. An event tree model is developed in Figure 7 of [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF] to describe the possible consequences.

Depending on the severity of the consequences, they can be classified into three categories, as shown in Table V. The event occurrence probabilities in the event tree model are summarized in Table VI. Unlike most of the existing works that assume time-invariant event occurrence probabilities (e.g., see [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]), in this paper, we consider the time-dependent behaviors of the LPM, ARES, FFS to better reflect the reality, as in practice, a large number of factors degrade these safety barriers physically, making their reliability decrease over time [START_REF] Yu | Stochastic chemo-physical-mechanical degradation analysis on hydrated cement under acidic environments[END_REF]. For example, the lightning protection mast is subject to corrosion failure that degrades its reliability significantly over time [START_REF] Ghavamian | Concerns of corrosive effects with respect to lightning protection systems[END_REF]. To simplify the analysis, let us assume that the time to failure of the three safety barriers follow Weibull distributions. Then, the time-dependent failure probabilities characterize their time-dependent failure August 17, 2020 DRAFT behavior:

p L (t + τ ) = P r(X LP M (t + τ ) = 1) = e -t+τ η L β L , (19)
where η L and β L represent the characteristic life and shape parameter of the LPM, respectively. The same expression also applies to the ARES and FFS. The values of these parameters are summarized in Table VI.

In this case study, the performance of the tank farm can be represented by the number of damaged tanks, as if a tank is damaged, it has to be shut down for repair and cannot be used to generate profits [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]. Let As shown in Table V, when local damages occur, X(t) jumps from i to i + 1, 0 ≤ i ≤ 5. When global damages occur, more than half of the tanks would be affected, i.e., the system directly enters state 6. We also assume that due to limited recovery resource, only one storage tank can be repaired at a time and minimal repair policy is taken [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]. The possible state transitions of the NHSMRP is summarized in Figure 6. For the direct and indirect losses associated with the state transitions, we assume the values based on those in [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF] (in arbitrary units): Fire spread to other tanks and caused large-scale damages that affect more than half of the tanks.

X(t) = 0, 1
d i,j =            6, 0 ≤ i ≤ 5, j = i + 1, 60, 0 ≤ i ≤ 4, j = 6, 0, otherwise. ( 20
)
l i = 0.6 • i (arbitrary unit per day), 0 ≤ i ≤ 5. (21)
3) Construct the renewal kernel: The renewal kernel matrix can be determined based on Eq. ( 13). To show how to apply this equation, we take Q i,i+1 (t, τ ), 1 ≤ i ≤ 5 as an example. As can be seen from Figure 6, from state i, there are four possible destinations: state i -1, i, i + 1 and 6. Also, from the event tree model, it can be seen that the transitions to states i, i + 1 and 6 are mutually exclusive, given the occurrence of the initiating event. Hence, for Q i,i+1 (t, τ ), Eq. ( 13) becomes The transition from state i to state i + 1 is caused by local damages. From the event tree model (Figure 7 of [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]) and Table V, f T,i,i+1 (x, τ ) can be calculated by

Q i,i+1 (t, τ ) = t 0 f T,i,i+1 (x, τ ) (1 -F T,i,i-1 (x, τ )) dx. (22)
f T,i,i+1 (x, τ ) = λ D • e -λ D x (1 -p L (x + τ ))• (p DD p B + (1 -p DD )(p A (x + τ ) + (1 -p A (x + τ ))(p F (x + τ ) + (1 -p F (x + τ ))(1 -p B )))) , (23)
where p DD represents the probability that a direct damage is caused by the lightning strike; p B represents the probability that the fire brigade responds correctly and promptly; p L (x + τ ), p A (x + τ ) and p F (x + τ ) represent respectively the probability that the LPM, ARES and FFS work properly at time (x + τ ).

The repair of the damaged tank causes the transition from state i to i-1. As in [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF], we assume that only one tank can be repaired at a time, and the recovery time between adjacent states follows a LogNormal distribution.

Therefore, we have

f T,i,i-1 (x, τ ) = 1 √ 2πσ r x e - (ln(x)-µr ) 2 2σ 2 r , t > 0, i = 1, 2, • • • , 5. (24)
where µ r and σ r are the mean and standard deviation of the associated normal distribution, whose values are given in Table VI.

Similarly, we can derive

Q i,i-1 (t, τ ), Q i,i (t, τ ) and Q i,6 (t, τ ): Q i,i (t, τ ) = t 0 f T,i,i (x, τ ) (1 -F T,i,i-1 (x, τ )) dx, (25)
Q i,i-1 (t, τ ) = t 0 f T,i,i-1 (x, τ ) (1 -F T,i,i (x, τ ) -F T,i,i+1 (x, τ ) -F T,i,6 (x, τ )) dx = t 0 f T,i,i-1 (x, τ ) • e λ D •x dx, (26)
Q i,6 (t, τ ) = t 0 f T,i,6 (x, τ ) (1 -F T,i,i-1 (x, τ )) dx, (27)
where f T,i,i (t + τ ) and f T,i,6 (t + τ ) can be derived based on the event tree model:

f T,i,i (t + τ ) = λ D e -λ D t • p L (t + τ ). (28)
f T,i,6 (t + τ ) = λ D e -λ D t • (1 -p L (t + τ ))• (p DD (1 -p B (t + τ )) + (1 -p DD )(1 -p A (t + τ ))(1 -p F (t + τ ))(1 -p B )) . (29)
For state 0, there are three possible destinations: state 0, 1 and state 6. We can derive the corresponding renewal kernels in a similar way:

Q 0,0 (t, τ ) = t 0 f T,0,0 (x, τ )dx = t 0 f T,i,i (x, τ )dx, (30)
Q 0,1 (t, τ ) = t 0 f T,0,1 (x, τ )dx = t 0 f T,i,i+1 (x, τ )dx, (31)
Q 0,6 (t, τ ) = t 0 f T,0,6 (x, τ )dx = t 0 f T,i,6 (x, τ )dx. (32)
Eqs. ( 22) -( 32) define the renewal kernel of the developed NHSMRP. The values of the model parameters are summarized in Table VI. It should be noted that in our data source [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF], the failure probabilities of LPM, ARES and FFS are assumed to be constant values. As we consider the time-dependent failure behavior of three safety barriers in this paper, the parameters of their lifetime distributions need to be set based on the original values in [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]. To do this, we assume that the failure probabilities of these safety barriers at 2 (years) equal the original values in [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF], and that the failure probabilities at T = 50 (years) increase to 0.99. The distribution parameters are obtained by solving these two equations, given that the failure probabilities have been defined

in Eq. [START_REF] Rocchetta | A power-flow emulator approach for resilience assessment of repairable power grids subject to weather-induced failures and data deficiency[END_REF]. p F B 0.693 [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF] T th,Rc 45 (d) [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF] L tol 70 (arbitrary unit) Assumed based on [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF] C. Numerical experiment

We conduct an experiment to test the computational efficiency of the developed resilience analysis method (Algorithm 1). Ten groups of experiments are designed, and the sample size of each group increases evenly from 10 4 to 10 6 . In each group, both the developed method and a standard Monte Carlo simulation are used to evaluate the resilience of the tank farm in (0, 50) (years). The experiment is repeated five times, as elapsed time of the program depends on the number of lightning that hit the tank farm in the evaluation horizon, which is uncertain and varies in each run of the simulation. Average CPU times over the five runs are used to compare the computational efficiency of the two methods in Figure 7. As the computational time of both the two methods depend on the number of events that need to be generated in the evaluation horizon. In Figure 7, we present two cases: in Figure 7 (a), the arrival rate is set to a low value λ D = 9.42 × 10 -3 (year -1 ), which indicates an average number of 2.12 initiating events would occur per simulation run; while in Figure 7 (b), we select a high value of λ D = 9.42 × 10 -2 (year -1 ), which makes the average number of initiating events in the evaluation horizon to be 21.2. The computer that runs the experiment has a CPU of 2.59 GHz (12 cores) and 64 GB RAM.

Figure 7 shows that in both cases, the computational times of both the two methods increase linearly with sample size, but the developed method has a significantly smaller slope. The smaller slope means that the computational time of the developed method increases much slower compared to standard Monte Carlo methods.

The numerical values of the slopes are given in Figure 7. Based on these values, a detailed comparison can be made: the developed method is 10.89 times faster compared to the standard Monte Carlo methods in case (a), while in case (b), the advantage increases to 29.74. The main reason for the benefit of the developed method is that implements vectorization in simulation and uses a linear interpolation model to reduce the computational burden in evaluating the renewal kernels. The comparison between case (a) and (b) also reveals that as the number of needed simulations increases, the efficiency of the developed simulation method also improves.

Figure 8 investigates the accuracy of the developed methods by comparing it to a benchmark of 10 7 simulations using the standard Monte Carlo simulation method. Average relative errors over the five repetitions are used as a measure of simulation accuracy. As a comparison, we also present the result of the standard Monte Carlo simulation under the same conditions. We can see that the accuracy of the developed method is comparable to the standard Monte Carlo simulation. We should note that the relative errors for the resilient and recovery probabilities are higher compared to the other two metrics. The reason is that these two metrics are conditional probabilities conditioned on the occurrence of initial damages. Hence, the sample sizes for these two metrics (years), respectively. The sample size of the simulation is 10 6 . The results of the analysis are given in Figure 9. We test the correctness of the developed method by comparing its results to those from a standard Monte Carlo simulation (without vectorization and interpolation). In the standard Monte Carlo simulation, one sample path is generated per simulation run by simulating sequentially occurrences of disruptive events, consequences of disruptive events (based on the event tree model in Figure 7 of [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]) and recovery times (Eq. ( 24)). The comparisons show that developed methods are able to estimate the resilience metrics accurately.

Figure 9 also reveals that as the evaluation horizon (t) gets longer, all the four resilience metrics decrease.

Two reasons mainly cause this. First, as shown in Figure 10, the reliabilities of the three key safety barriers, i.e., LPM, ARES and FFS, degrade over time. The degradation of these key safety barriers impairs the system's ability to stay resilient and, as a result, worsens the evaluated resilience metrics. Besides, more disruptive events (lightning) are expected in a longer interval with a larger t, as predicted by the Poisson process model.

Since lightning triggers the initiating event in the resilience analysis, more lightning arrivals could decrease the resilience, even though failure probabilities of safety barriers do not degrade over time. This finding suggests the necessity to monitor the degradation of resilience to ensure that it can operate with high degree of resilience.

Also, intervention measures like preventive maintenance can be planned based on the measured resilience. When the resilience metrics drop below tolerable values, the systems need to stop for maintenance to regain resilience.

At a first glance, the overall resilience is Figure 9 (d) seems quite satisfactory (above 0.99) in (0, 50) (years), (ii) High arrival rate: λD = 9.42 × 10 -2 (year -1 ).

Fig. 8: Accuracy of the developed method.

indicating that we have high confidence (above 99%) that the total costs of operating the tank farm for 50 years would not exceed the tolerance. However, a closer examination on Figure 9 (b) and (c) reveals that the system is actually quite vulnerable. Once the lightning arrives and the protection mask fails, it is highly likely that the tank farm will enter non-resilient states and/or fail to satisfy recovery requirements. This is caused by the fact that the safety barriers degrade over time, which increases significantly the risks when the protection barriers fail. The overall resilience (and similar ones like the resilience triangle) still yield satisfactory results mainly because its value also depends on the occurrence probability of the initiating events (which is usually very low).

The low occurrence probability of the initiating events provides us with "fake confidence" in this case. Resilient and recovery probability, on the other hand, can complement overall resilience as they are defined conditional probabilities that draw our attention directly to the performance of safety barriers and recovery processes. It should be noted that it seems the uncertainty bounds in the beginning range of Figure 9 (c) is larger than the other three metrics. This is because we have relatively small number of samples in this range. Figure 11 compares the developed time-dependent resilience model to the time-static ones [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]. In the timestatic model, failure probabilities of LPM, ARES and FFS are assumed to be constants and takes the values in [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]. In Figure 11, we also present the result of a worst-case time-static resilience analysis, where the failure probabilities of LPM, ARES and FFS are all assumed to be one. The comparison shows that depending on the parameter values chosen, time-dependent resilience analysis might overestimate (in the first case) or underestimate the resilience (in the second case). The difference is mainly determined by whether the assumed constant-value failure probabilities of safety barriers overestimate or underestimate the true failure probabilities.

As the developed method allows capturing the realistic time-dependent behaviors of safety barriers, it is more accurate than the existing time-static resilience analysis method in [START_REF] Du | Resilience analysis of multistate systems based on markov reward processes[END_REF]. 2) Improving the resilience: There are different design measures that can improve resilience. Depending on which phase in the accident evolution process these measures take place, we consider three types of design measures as summarized in Table VII. In the analysis, we improve the associated parameters by 25%

to investigate their effectiveness for improving different resilience metrics. Note that η LP M , η ARES , η F F S are characteristic life of the corresponding safety barriers. Improving them improves the reliabilities of the corresponding safety barriers. The results of the analysis are presented in Figure 12. It can be seen that only improving protection measures (i.e., increase the reliability of the lightning protection mask) can improve the resistant probability. This can be explained from the event tree model in Figure 7 of [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]: only the protection mask can protect the tank farm from the damages caused by lightning. The resilient probability can only be improved through improving safety barriers (i.e., the ARES and FFS). The most efficient way of improving recovery probability is to improve the recovery measures, as the recovery measures directly affect the recovery time distribution. However, improving the protection and safety barriers can also help improve the recovery probability. The result is caused by the fact taht the two improvement measures can reduce the probability that the system enters non-resilient states. From the definition of recovery probability in Eq. ( 7), we can see that entering non-resilient states could decrease the resilient probability, as non-resilient states are not recoverable. All three measures can improve overall resilience. Among them, improving protection measures is the most one. This is because improving protection measures reduces the occurrence probability of the initiating events in the event tree model (see Figure 7 of [START_REF] Zeng | An integrated modeling framework for quantitative business continuity assessment, Process August 17[END_REF]). Without initiating events, no subsequent damages and losses would occur. However, being more efficient does not mean that improving protection measures alone is enough. From Figure 12, we can see that protection measures cannot improve the resilient and recovery probability. Hence, different improvement measures need to be combined to achieve better resilience.

Fig. 12: Effectiveness of different improvement measures.

V. APPLICATION ON A RE-CONFIGURABLE COMPUTING SYSTEM

A. System descriptions

In this section, we apply the developed model on a second case study, a re-configurable computing system from [START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF]. The system has 16 processors, 16 memories an interconnection network that allows each processor to access any memory. For the system to be operational, a minimal number of 4 processors, 4 memories and one interconnection network are required. Otherwise, the system enters a non-operational state. Above these threshold values, more available processors and memories indicate better performance. Let c, m and s represent the number of working processors, memories and interconnection networks, respectively. Then, the performance of the computing system could be quantified by [START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF]

r =        max(c, m) 1 -1 - 1 max(c, m) min(c,m) , if c, m ≥ 4, s = 1, 0, otherwise. (33)
August 17, 2020 DRAFT All the three types of components could fail. According to [START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF], the failure times can be assumed to follow exponential distributions. Once failure occurs, the computing system can be reconfigured: there are automatic mechanisms for failure detection and isolation, which locates the failed component and remove it from the normal ones. Therefore, the system can still operate, but with a degraded performance level. The reconfiguration might also fail with a probability of (1 -p c ), where p c is called coverage probability. Once a reconfiguration fails, the computing system goes into a non-operational state. The reconfigured system can be recovered by repairing the isolated failed component. In this paper, we take the same assumption as [START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF] that only one failed component can be repaired at a time and the repair times are also exponentially distributed. The failure rates and repair rates are given in Sect. II of [START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF].

B. Resilience modeling

The system state can be defined by a tuple (c, m, s). The state space comprises all the combinations of (c, m, s), with 4 ≤ c, m ≤ 16, s = 0, 1, and an additional state F, which represents the non-operational state.

Unlike the original case study in [START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF], we assume in this paper that the computing system provides essential service so that its breakdown is unacceptable. Once the system enters the non-operational state, it cannot be recovered and remains in the non-operational state. Hence, the non-operational state is an absorbing state. The behavior of the system can be characterized by deriving the renewal kernels. Let's take the jump from (16, 16, 1)

to (16, 15, 1) as an example. From state (16, 16, 1), we can reach three states, i.e., (16, 15, 1), (15, 16, 1) and F. From Eq. ( 13), the corresponding element in the renewal kernel can be derived:

Q (16,16 ,
= λ m λ m + λ c + λ s 1 -e -(λm+λc+λs)t ,
where λ m , λ c and λ s are the failure rates of the memory, processor and network switch, respectively. It should be noted that in this case study, the state transitions are homogeneous. Therefore, the τ in Eq. ( 13) can be dropped. From Eq. ( 14) and Eq. ( 15), we can further derive the transition probability of the embedded chain, and the associated holding time distributions: = 1 -e -(λm+λc+λs)t .

P (16,
It can be seen that the holding time distribution is an exponential distribution. Hence, the system degenerates to a homogeneous Markov reward process. The system has 170 states and the other elements in the renewal kernel can be determined in a similar way.

August 17, 2020 DRAFT Both direct and indirect losses can be observed in the computing system. In this case study, we assume the values of the direct losses to be where n f is the total number of failed memories and processors. The indirect cost per unit time l i is assumed to be proportional to the performance loss:

d i,j =           
l i = 10
• max i r i -r i (arbitrary unit per hour), for all the working states, where r i is the performance rate for state i and can be calculated by Eq. [START_REF] Katz | Stochastic modeling of hurricane damage[END_REF].

C. Results and discussions

A time-dependent resilience analysis is conducted based on Algorithm 1 to calculate the four resilience metrics. The threshold for maximal allowable recovery time is taken to be T th,Rc = 10 (hours) and the maximal allowable loss is taken to be L tol = 50 (arbitrary unit). The evaluation is done every 36 (hours), up to 720 (hours). The other parameters take the same values as those in [START_REF] Smith | Performability analysis: measures, an algorithm, and a case study[END_REF]. The sample size of the simulation is 10 6 .

The results are given in Figure 13. From Figure 13, it can be seen that all the four resilience metrics decrease as the evaluation horizon (t) increases. This is mainly because, as t increases, the reliabilities of the components drop. Thus, more component failures are expected, causing severe performance loss and lowering the resilience.

Such an analysis can help decision makers make operational decisions considering resilience. For example, it is not recommended to operate the computing system continuously for the whole evaluation horizon (720 hours), as the analyses reveal that we would have very high risks to suffer resilience problems. For example, from Figure 13 (d), if the system is operating for 720 hours, the risk that the loss is beyond tolerance is 66.64%.

Fig. 13: Time-dependent resilience analysis for the computing system.

VI. CONCLUSION

In this paper, we developed a non-homogeneous Semi-Markov reward process-based resilience model for multi-state systems. The time-dependent behavior of system resilience can be captured in the developed model thanks to its non-homogeneous nature. Four resilience metrics, i.e., resistant probability, resilient probability, recovery probability, and overall resilience, were used to quantify system resilience. A simulation algorithm was developed for resilience analysis based on vectorization and linear interpolation. Computational experiments were designed to compare the computational performance of the developed algorithm to existing methods.

The results showed that the developed outperforms the existing methods, thanks to the reduced overhead costs through vectorization and the reduced computational costs of generating samples through linear interpolation.

The developed methods were applied on two real-world case studies. The main findings from th case studies include: first, compared to existing time-state analysis methods, the developed method can more accurately describe the actual resilience behavior, due to its capability to capture time-dependent features. Second, existing methods based on global resilience metrics like resilience triangle might be overconfident in some cases, especially when the occurrence probability of disruptive events is small. Applying the developed four-metric framework can alleviate this overconfidence as resilient probability and recovery probability focus on system behavior after the disruptive event occurs. Third, if we want to improve the resilience of the system, the most efficient way is to improve the reliability of protection measures. However, doing this alone is not enough, as it does not affect resilient and recovery probability. Thus, decision-makers need to coordinate the different improvement measures for designing a system toward comprehensive resilience.

In the future, some improvements might be considered to the developed methods. First, the developed resilience analysis methods are still based on simulation. Analytical methods could be investigated to reduce the computational burden of the analysis further. Second, the developed model is multi-state in the sense that the system performance takes discrete values. In practice, some systems might have discrete and continuous performance levels at the same time. The developed model can be extended to these systems. Third, maintenance, especially preventive maintenance and predictive maintenance, can be considered in the developed model, so that optimal maintenance strategies can be planned to maximize the system resilience. Also, in this paper, we assumed that the parameter values of the NHSMRP model is precisely known to us. In practice, however, we often have epistemic uncertainty regarding these parameter values. How to handle epistemic uncertainty in the developed model, then, becomes an interesting question that deserves further investigations. Another limitation of the current model is that, cascading failures, which are commonly observed in an actual network system [START_REF] Dobson | Complex systems analysis of series of blackouts: Cascading failure, critical points, and self-organization[END_REF],

are not considered. However, this can be done by extending the NHSMRP model to include state transitions corresponding to common cause failures (simultaneous failure of multiple components due to a common cause).

The added state transition can be naturally used to describe cascading failures. Lastly, in a lot of real-world systems, the reliability of safety barriers might depend on the operation levels of the system. This fact can also be considered in an extension of the developed model so that joint optimization of operation and resilience can be considered.

Fig. 1 :

 1 Fig. 1: An illustration of the NHSMRP-based resilience model.

Following

  the previous steps, a NHSMRP is constructed for resilience analysis. Compactly, let us denote the developed model by a tuple < π 0 , Q(t, τ ), D, l >, where D = [d i,j : i, j ∈ S] is a matrix whose elements represent the corresponding direct losses and l = [l 0 , l 1 , • • • , l n ] is a vector that contains the unit indirect loss for each state. Based on the developed model, the four resilience metrics defined in Sect. II-B are calculated for resilience analysis. In Sect. III-C, we present an efficient Monte Carlo simulation algorithm to calculate the resilience metrics.

Fig. 3 :

 3 Fig. 3: Procedures of applying the developed model for resilience analysis

  the elements in θ which satisfy u. * c i,j ( τ ) ≤ f ij ( θ, τ )./g ij ( θ, τ ) and discard the other elements; Delete the elements in τ whose holding time has been sampled in θ; end D. Performance analysis and numerical experiments 1) Theoretical analysis: The computational complexity of the developed method, denoted by O d , depends on the complexity of the simulation algorithm O sim , and the overhead cost O oc :

Fig. 5 :

 5 Fig. 5: Running times of the first numerical experiment.

Fig. 6 :

 6 Fig. 6: State transition diagram of the NHSMRP.

  are smaller than the other two metrics, which results in high errors. Compared to the standard Monte Carlo simulation, the only source for systematic error is the approximation error of the linear interpolation model, which is controlled by the training sample size. A larger training sample size would reduce the evaluation error, but requires longer execution time. D. Resilience analyses 1) Time-dependent resilience analysis: The parameters in Table VI are used to investigate the resilience behavior of the tank farm in [0, t] using Algorithm 1. In the training phase of Algorithm 1, 10 5 training (a) Low arrival rate: λD = 9.42 × 10 -3 (year -1 ). (b) High arrival rate: λD = 9.42 × 10 -2 (year -1 ).

Fig. 7 :

 7 Fig. 7: Comparison of computational efficiency.

  (i) Low arrival rate: λD = 9.42 × 10 -3 (year -1 ).

Fig. 9 :

 9 Fig. 9: Time-dependent resilience analysis (Sample size: 10 6 ).

Fig. 10 :

 10 Fig. 10: Reliability of the three major safety barriers.

Fig. 11 :

 11 Fig. 11: Comparison to time-static resilience analysis.

0 λ 0 λ

 00 1),(16,15,1) (t) = P r(X n+1 = (16, 15, 1), θ n+1 ≤ t | X n = (16, 16, 1)) = t m e -λmx P r(T c > x)P r(T s > x)dx = t m e -λmx e -λcx e -λsx dx

  Resilience analysis based on the NHSMRP model. Training phase: τ tr ← Equally insert n tr points into [0, T ];for i, j ∈ S do p tr ← Evaluate Eq. (14) with τ = τ tr , using numerical integration methods; pi,j (τ ) ← Train a linear interpolation model based on training data ( τ tr , p tr ); Set x cur , x next , θ, t cur , t next to column vectors with N S elements of zeros;while At least one element in t cur < T do for x cur = Unique states in x cur do τ ← t cur ( x cur == x cur );

	end
	Simulation phase:

, • • • , n state are the linear interpolation models to approximate the transition probabilities p xcur,j , j = 1, 2, • • • , n state , respectively. The output x next is a vector of the same size as τ , which contains the generated next states for samples with current state x cur and age τ . August 17, 2020 DRAFT Algorithm 1: x next ( x cur == x cur ) ← Starting from state x cur , generate the next state with parameter τ (Algorithm 2); θ( x cur == x cur ) ← Starting from state x cur and ending with states x next ( x cur == x cur ), generate holding times with parameter τ based on Algorithm 3;

TABLE I :

 I Similar methods in the literature.

	Method	Category	Applicable models	Computational complexity
	Tijms' method [39]	Numerical integration	Markov reward model	O	n 2 st ty ∆ 2	[40]
	Janssen's method [38]	Numerical integration	Semi-Markov reward model	O n 2 st	t 2 ∆ 2 + t ∆
	Standard Monte Carlo [37]	Monte Carlo simulation	Both Markov and semi-	O(ns ntnst)Op +O(ns nt)Oht +
			Markov reward model	O(ns)O oh

TABLE II :

 II P r(Y 10 ≤ 10) calculated by different methods.

	M	nst	Developed method Tijms' method [39] Standard Monte Carlo [37] Theoretical value from [40]
	4	5	0.445951	0.444693	0.444687	0.434068
	5	13	0.226824	0.226727	0.226679	0.214623
	6	25	0.171730	0.170841	0.171006	0.160671
	7	41	0.142809	0.142358	0.142764	0.137375
	9	81	0.121781	0.121266	0.121818	-
	30	1513	0.085572	-	0.085364	-
	60	6613	0.080194	-	-	-

-: Running time exceeds 300 (seconds).

TABLE III :

 III Parameter values of the second numerical experiment (in arbitrary units).

	Parameters	λ	p I	η f	β f	ηr	βr
	Values	10 -3	1.6 × 10 -3	100	2	1	2

TABLE IV :

 IV Results of the second numerical experiment.

	Methods	Developed method	Standard Monte Carlo	Janssen's method
	P r(X(1000) = 0)	0.996768 ± 0.000883 0.996701 ± 0.000915	0.995583
	P r(X(1000) = 1)	0.002349 ± 0.000111 0.002384 ± 0.000112	0.002022
	P r(X(1000) = 2)	0.000058 ± 0.000095 0.000059 ± 0.000096	0.002394
	Computation time (s)	6.80	328.78	172.12

  , • • • 5 represent the number of damaged tanks when less than half of the tanks are damaged (X(t) ≤ 5). When more than half of the tanks are damaged, it is considered an unrecoverable state because the costs become too high to make a recovery economically reasonable. Let X(t) = 6 represent this unrecoverable state. Therefore, the state space of the NHSMRP model can be defined byS = [0, 1, • • • , 6],where the first five states are resilience state while the last one is non-resilient state (B 0 = [0, 1, 2, 3, 4, 5] and B 1 =[START_REF] Caputo | Computing resilience of process plants under na-tech events: Methodology and application to sesmic loading scenarios[END_REF]).

TABLE V :

 V Classification of the consequences. Fire was caused by lightning but contained promptly. Only one tank is affected.

	Category	Event sequences Descriptions
	Normal functioning (C 0 )	S 0	No damage is caused by the lightning.
	Local damage (C 1 )	S 1 , S 3 , S 4 , S 5
	Global damage (C 2 )	S 2 , S 6

TABLE VI :

 VI Parameter values of the developed model

	Parameter	Value	Source
	λ D	9.842 × 10 -3 (year -1 )	[41]
	η LP M	24.889 (year)	Assumed based on [28]
	β LP M	2.189	Assumed based on [28]
	η ARES	19.616 (year)	Assumed based on [28]
	β ARES	1.632	Assumed based on [28]
	η F F S	23.419 (year)	Assumed based on [28]
	β F F S	2.014	Assumed based on [28]
	µr	3.388 (d)	[28]
	σr	0.166 (d)	[28]
	p DD	1.6 × 10 -2	[41]

TABLE VII :

 VII Improvement measures.

	Category	Target	Affected parameters
	Improve protection measure Improve the capability to resist the damages of disruptive events.	η LP M
	Improve safety barrier	Reduce the likelihood of system collapse.	η ARES , η F F S
	Improve recovery measure	Reduce the recovery time.	µr

  10 • n f (arbitrary unit), if remaining in operational states, 60 (arbitrary unit), if entering the non-operational state,

	0, otherwise,

August 17, 2020 DRAFT

ACKNOWLEDGMENT

The authors would like express their deepest gratitude to the associate editor, domain editor, and two anonymous reviewers, for their insightful comments and suggestions that help to significantly improve this paper. This work was finished during Dr. Shijia Du's visit to Université Paris-Saclay. She would like to express her deepest gratitude to Professor Enrico Zio for hosting her visit. Dr. Du's research is supported by Natural Science Foundation of China (NSFC) under grant No. 71601010. The work in this paper is also partially supported by High-end Foreign Experts Recruitment Program of Shanghai university, China (July -August, 2019). Lastly, Shijia Du and Zhiguo Zeng would like to thank their baby, Zhuochen (Duke) Zeng, who arrives during the revision of this manuscript, for making the revisions process much longer and more difficult.