An Interval Mathematic-Based Methodology for Reliable Resilience Analysis of Power Systems in the Presence of Data Uncertainties

: Prevention and mitigation of low probability, high impact events is becoming a priority for power system operators, as natural disasters are hitting critical infrastructures with increased frequency all over the world. Protecting power networks against these events means improving their resilience in planning, operation and restoration phases. This paper introduces a framework based on time-varying interval Markov Chains to assess system’s resilience to catastrophic events. After recognizing the difﬁculties in accurately deﬁning transition probabilities, due to the presence of data uncertainty, this paper proposes a novel approach based on interval mathematics, which allows representing the elements of the transition matrices by intervals, and computing reliable enclosures of the transient state probabilities. The proposed framework is validated on a case study, which is based on the resilience analysis of a power system in the presence of multiple contemporary faults. The results show how the proposed framework can successfully enclose all the possible outcomes obtained through Monte Carlo simulation. The main advantages are the low computational burden and high scalability achieved.


Introduction
Power systems were traditionally planned and operated according to security and reliability principles, which involved the analysis of credible events, such as N-1 or N-2 contingencies.
However, due to the increasing frequency of natural disasters [1,2] and the greater vulnerability of power systems to cyber-attacks [3], planning studies and operational paradigms are changing. Hence, researchers and power system operators are increasingly interested to the concept of resilience [4].
Resilience is defined by the IEEE as [5]: the ability to withstand and reduce the magnitude and/or duration of disruptive events, which include the capability to anticipate, absorb, adapt to, and/or rapidly recover from such an event.
On the basis of this definition, resilience studies can be classified in three main branches [6]: • The ability to anticipate High Impact Low Probability (HILP) events is addressed in resilience-based planning, which aims at identifying long-term strategies to strengthen power networks. • The capacity to absorb and adapt is analysed in resilience-based response, which aims at improving power system operation in terms of reserves, network topology control, vulnerability analyses. • Finally, the ability to recover from a catastrophic event, is explored in resilience-based restoration analysis, which aims at minimizing the impacts of power outage after a catastrophic event.
Although the scientific literature proposes a wide range of quantitative and qualitative methods and metrics to assess resilience, it is still lacking a global consensus over a unified framework [7]. Qualitative evaluation is based on questionnaires, investigations, ratings from Transmission System Operators, and overall considerations regarding correlated systems, such as fuel supply chain or telecommunication systems. The results are analyzed through qualitative planning methods, such as Analytic Hierarchy Process (AHP), in decision making tools [8].
On the other hand, quantitative methods introduce metrics to compare the level of resilience of different systems. They can be simulation-based, analytical or based on statistical analysis of historical data.
The main metrics aim at evaluating the difference between the expected system performance and the real system performance under a HILP event. Some examples are restoration costs, restoration time or loss of load frequency and loss of load expectation [9].
In particular, the authors of Reference [10], show how the definition of resilience implies a dependency on time-a broad classification can be done between short-term and long-term resilience. Short-term resilience refers to the characteristics that the power system should have before, during and after a HILP contingency, while long-term resilience refers to the adaptability of the system to unknown future threats. The latter can be achieved by performing multiple risk and reliability studies, and planning ahead for unexpected events.
Resilience is a concept that is strictly related to mitigation, which is the process of reducing the unavoidable damage suffered by a system due to a HILP event. Mitigation strategies vary depending on the considered time-scale. They can be evaluated during the planning and operation phase, and they are mainly classified in physical interventions on the network to improve its overall robustness, or implementation of more advanced monitoring and control capabilities [11]. Examples of actions aimed at reducing the impact of unexpected events in the long run are [12]: In all these contexts, the analysis of the potential impacts of renewable energy sources (RES) is an interesting open problem to address. In particular, although the use of different, decentralized RES can be beneficial for the power system in mitigating the impacts of HILP events, their intermittent power profiles, the limited black start capability, and the misoperation of traditional coordinated protections, are some of the most critical issues hindering the RES contributions to resiliency enhancement, especially in long-term planning scenarios. In this domain, the research for reliable computing paradigms aimed at managing the uncertainties induced by both RES and other relevant sources, for example, weather dependent fault/reparation rates, represents a relevant issue to address.
A model often used to deal with resilience-based planning problems is based on Markov processes [13,14], which are employed in the task of modelling the system states evolution consequent to HILP events, once the deterministic transition probabilities are assessed. Anyway, the reliable estimation of these parameters is a challenging issue to address, especially in predictive resilience studies, which are mainly related to low frequency events. In this context, the available historical information are often not enough to reliably describe the state transition rates by deterministic variables. To face this issue, sampling-based approaches, such as Monte Carlo simulations, could be a useful solution strategy, but the associated computational burden, especially for large-scale systems, makes them unfeasible for practical purposes. Furthermore, stochastic models often need assumptions about the probability distributions of the uncertain variables, which are typically assumed Gaussian, although this assumption is often not supported by experimental evidence [15].
To overcome these limitations, stochastic models could be enhanced by advanced modelling approaches, such as convex models [16]. A convex model is a set-theoretic approach, which aims at modelling the data uncertainty by defining subsets over the set of allowed events. No probability measure is defined over the set, such as in the case of probabilistic approaches. An instance of convex models for handling uncertainty is interval mathematics [17], which considers hypercubes as convex sets to describe uncertainty.
Interval Mathematics (IM) was first introduced to handle the uncertainty raising from numeric computations [18], however it spread over different disciplines as an alternative tool to deal with stochastic behaviours or inaccurate models. It was deployed in different power systems and power electronics applications, such as load flow [19] or power converter's characteristics [20]. IM was also applied to model uncertain Markov Chains, typically referred as Interval Markov Chains, where the state transition probabilities are uncertain but bounded in fixed real-value intervals [21]. This methodology could play an important role in resilience-based long term planning and mitigation, where the long time-horizon, the scarcity of data characterizing HILP events, and the large number of analyzed scenarios, require the conceptualization of reliable computing techniques, which allow overcoming the intrinsic limitations of homogeneous, deterministic Markov Chains.
Armed with such a vision, in this paper Interval Markov Chains will be adopted for modelling uncertain time-varying models applied to resilience-based long term planning and mitigation. The main advantage deriving by the application of this computing paradigm is its low computational burden, which grants large scalability levels, and its robustness, since it allows to reliably compute an outer estimation of all possible outcomes of a given set of random scenarios, considering all the possible combinations of the uncertain variables.
The paper is organized as follows: Section 2 describes the mathematical methods employed for the analysis, such as Markov Chains, Interval Mathematics and finally Interval Markov Chains, Section 3 presents the case study aimed at validating the proposed methodology, Section Section 4 suggests future directions for research on this topic and Section Section 5 draws the main conclusions from the results.

Markov Chains
A Markov Chain is a stochastic dynamic process X(t) : x ≥ 0 which satisfies the Markov property: Hence X(t), which is a variable representing the state s occupied over the time t in a random occurrence, evolves based only on its present state. Markov Chains can model both discrete and continuous processes and the number of states can be infinite or finite.

Discrete-Time Markov Chains
In this paper, the analysis will be focused on discrete time Markov Chains. A Discrete Time Markov Chain describes the evolution of a system over a set S of N states, S = {s 1 , ..., s N }. At each time t, the system can be only in a single state. The probability that the system is in the state s i at time t is given by the state probability π i (t). Since each state is characterized by a state probability and these probabilities vary over time, a vector Π(t) = {π 1 (t), ..., π N (t)} of state probabilities is defined for each step t of the evolution over time [22].
State probabilities are strictly related to transition probabilities p ij (t), which describe the probability that the system evolves from state s i to state s j at time t.
Transition probabilities can be constant or vary over time. In the first case the Markov Chain is called homogeneous. In this paper, we will deal with time-varying transition probabilities [23].
The transition matrix P(t) is a square matrix of dimension n, where n is the number of possible states. It contains the set of transition probabilities from a state to another. The transition matrix is stochastic, since each of its elements belongs to the interval [0, 1], and each row of the matrix P(t) satisfies the condition: The trajectory of the state probabilities from an initial step to a certain time T, can be obtained by solving a finite difference equation for each time step t: Markov Chains find different applications in power system engineering, such as the analysis of voltage profiles [24], characterization of the harmonics [25] or analysis of the electrical vehicles [26] to name some.

Effects of Data Uncertainties
Transition matrices, which rule the evolution of Markov Chains, are commonly described by deterministic real-valued coefficients. Such a description does not consider the large amount of uncertainty sources which affect transition rates. For example, the potential sources of uncertainty within power systems include natural temporal variations, monitoring and forecasting errors, wrong estimation of physical parameters [15]. Uncertainty can be either aleatory or epistemic: aleatory uncertainty is irreducible, since it is inherited from the random behaviour of the considered process. Epistemic uncertainty can be reduced by improving the accuracy of the models or the estimation of the parameters.
An example of epistemic uncertainty is the estimation of transition rates that results from a limited amount of available data. For example some influencing variables might be known with low accuracy or they might not even be measurable [27]. Another source of epistemic uncertainty is due to the classification of continuous variables into discrete classes: the distinction between the classes might not always be trivial, especially for boundary values.
Uncertainty can be represented with different methods. The most common is the probabilistic model, however it might not always be correctly quantifiable. Furthermore, exogenous variables, which influence transition rates, might not be easily represented in probabilistic form.
A wiser approach, in case there is no accurate probabilistic representation of a process, is to represent the transition matrix as an interval matrix: that is, a matrix whose elements are interval variables [28]. The following section presents the mathematical framework under which interval Markov Chains can be defined.

Interval Mathematics
Interval mathematics has been widely used in the scientific literature to address uncertain computations in different engineering branches [29][30][31]. With respect to power system engineering, some examples of its application have been described by References [32][33][34].
Given the field of real numbers R, an interval is a subset of R of the form [35]: The numbers a andā represent the lower and upper values of the interval. The set of all closed intervals in R can be called I(R). A graphical representation is shown in Figure 1. Real numbers can be considered a particular instance of I(R), in which a =ā. Binary operations can be defined for I(R). Given two intervals: A and B, basic operations are: The interval can also be characterized by a central value x 0 , a width w(X) and a radius ∆x: In the analysis of uncertain problems, the central value is usually taken as the value considered in case of a deterministic model. The radius depends on the degree of uncertainty of the parameter considered: greater uncertainty involves a larger radius.
An issue to deal with in interval analysis is the dependency problem [36]. When an interval variable appears more than once in computations, it is treated as a different and independent variable at each step of the computation [37]. Namely, in case of the product X · X, given x as an unknown number known to lie within the interval X, the product x · x does not imply that the first factor is the same as the second: they vary independently and this causes the widening of the resulting interval. For example, given the interval X = [−1, 1], the correct value for X 2 would be X 2 = [0, 1]. However, due to the dependency problem, by applying Equation (5), X 2 = [−1, 1], which contains the correct interval solution, but is wider.
A simple means to reduce the dependency problem is to rewrite, when possible, the expression to be computed in a way that equal variables appear only once. If overconservatism must be avoided, it is possible to apply affine arithmetic [38], at the cost of increased computational burden. Affine arithmetic maintains the dependencies between variables even after computations which would normally lead to wider intervals.
Interval mathematics can be extended to more than one dimension. An interval vector is a vector whose components are interval elements. For example, in case of a two dimensional interval vector X: The geometrical representation of that vector is a rectangle in the plane, as shown in Figure 2 that is, it is the set of all the points (x 1 , x 2 ) such that: For an interval vector, the following quantities can be defined: • The width of the interval vector w(X), that is the largest of the widths of its component intervals: • The norm, which is a generalization of the absolute value: • The midpoint, which is the point corresponding to the central values: Similarly, a matrix whose elements are intervals is called interval matrix. An interval matrix can be considered as a set of point matrices.

Interval Markov Chains
Transition probabilities and initial state probabilities of a Markov Chain might not be known with the desired degree of accuracy. The traditional approach to handle the uncertainty is to take the best estimates of the parameters, however the reliability of this approach is questionable. When considering high risk situations, such as in resilience studies including HILP events, such unreliability should be minimized, in order to avoid huge losses.
A solution to overcome the unreliability problem is the incorporation of imprecision into the model. This can be achieved in two ways [39]: • The assumption of time homogeneity, that is, transition probabilities considered constant over time, can be preserved, but they are considered imprecisely known [40,41]. • Classical models can be further improved by omitting the assumption of time homogeneity, as in the case of Markov set-chains theory [42].
In general, the described approaches involve the substitution of single and precisely known transition probabilities with sets of feasible probabilities. Modelling capabilities depend on the form used to represent such sets, which is connected to the constraints used to determine them.
The most common form of sets are intervals of probabilities: the transition probability p ij from state i to j is given in the form of interval [p ij ,p ij ], which is supposed to contain the unknown true transition probability. Although true transition probabilities are unknown, the values within intervals must be taken such that they sum to 1. Furthermore, a fundamental assumption to work with intervals is that all values of the interval must be reachable, that is, the property of coherence is satisfied.
Interval Markov chains map elementary events, which are the target states, to intervals of probabilities [40]. The theory of interval probabilities was developed by Weichselberger [43], and the main elements of this theory are illustrated below, in order to contextualize the theory of interval valued Markov Chains.
Consider a non-empty set Ω and a σ-algebra A of its subsets, which represent random events in Ω. A set function p : A → R is called probability if it satisfies Kolmogorov's axioms.
The introduction of interval probabilities requires the definitions of other two set functions on A, such as 0 ≤ L ≤ U ≤ 1, L(∅) = U(∅) = 0 and L(Ω) = U(Ω) = 1. The interval probability P(·) is an interval valued set function, containing the values p(·), which follows two axioms: In this context, Interval Valued Markov Chains can be defined as Markov Chains whose transition matrix's elements are interval probabilities.
A Markov Chain is a particular case of interval Markov Chain in which all interval probabilities are closed point intervals [44], that is, L(A) = U(A) ∀A ∈ A . The trajectories of the interval state probabilities can be computed by solving the equation: π I (t + 1) = P I (t) · π I (t) (16) iteratively for each t. By analyzing (16), it is worth noting as the main complexity characterizing this computing paradigm is basically related to the iterative product between uncertain interval variables. This non-linear interval operator add further uncertainties to the overall solution process, here referred as endogenous uncertainties, which could induce unrealistic outer estimations of the solution bounds, especially in the presence of highly correlated uncertainties.

Case Study
Interval valued and time varying discrete Markov Chains can be applied to resilience-based planning for a network in which N critical components are identified. The fault and repair rates of each component, λ i (t) and µ i (t), are taken as time varying and uncertain.
Considering time varying and uncertain transition probabilities means taking interval vectors of µ I (t) and λ I (t) for each component, whose elements are the interval parameters for each line, at each time step. In the following analysis, the superscript I to indicate an interval quantity will be dropped, since every value considered will be an interval.
Each analyzed component can be either in the working (1) or faulty (0) state, Figure 3, hence the total number of states of the system for each time step t is 2 N .
For example, in case of two lines and 3 time steps, the graph showing the possible evolution of the system is shown in Figure 4.
In such case, the space S of the possible states is: Assuming that the fault and repair rates of each line are independent probabilities, the corresponding interval transition matrix P, in case of N lines, is made up by the following interval elements p ij (t), which are the transition interval probabilities from state i to state j, when passing from time t = T to time t = T + 1: In this paper, a case study which considers four critical components and their evolution over 8 time steps is analyzed. The central values λ i0 (t) and µ i0 (t), for each component i, are computed for each time step, based on the fragility curve and on the repair process involved.
The radius of the intervals for the transition probabilities at each step is taken as 10% of their deterministic value.
The study has been conducted using the INTLAB toolbox of Matlab [45]. This toolbox allows the definition of interval valued elements and interval matrices and provides a set of operations on such elements.
The interval calculations lead to the results shown in Figure 5 for State 2, taken as an example. The results for all the other states are similar, with the interval mathematics robustly enclosing all possible scenarios. Actually, if the constraint ∑ n i=1 π i = 1 is not imposed, the obtained results are sets wider than the feasible set of solutions. Reducing the level of conservatism, by introducing the constraint, would come at the cost of larger computational time.
To validate the results and to check the robustness of interval computations, given the uncertainty previously defined, 10,000 sample trajectories are obtained from Monte Carlo simulations. In this case all the uncertainty sources are independent. The trajectories cover all the possible values of the state probabilities over time. Again, Figure 5 shows how interval mathematics can perfectly bound all the possible scenarios for state probabilities. Other Monte Carlo simulations are performed in case of dependent uncertainty sources. In this case, as shown in Figure 6, interval mathematics can still enclose all the possible trajectories, however it results in an outer bound: more robust, but also more conservative. Given this result, Monte Carlo simulations may seem better to estimate the real possible trajectories. However, such simulations require a prohibitive computational burden in case of large systems. On the contrary, interval mathematics, at the cost of greater conservatism, can reduce the computational times by different orders of magnitude. To give an idea, in the considered case with 4 components and 8 time steps, using an Intel Core i7 CPU, interval mathematics took 0.05 s to provide the results, while Monte Carlo simulation 33.2 s.
The limited number of steps is due to the increasing uncertainty bounds over time in non-homogeneous interval Markov Chains. The resulting intervals become wider as time steps increase, so a threshold on the radius of intervals should be defined. In this case, given the threshold ∆π max = 0.25, the maximum number of time steps chosen is 8.
The maximum number of time steps should be chosen by inspecting each state over time. Indeed, different states might have a different propagation of the error over time: as shown in Figure 7 for state 14, the final interval is much narrower than state 2, hence it should not be chosen to predict the maximum number of time step required to maintain the desired level of accuracy.  Table 1 shows the lower and upper bounds of state probabilities at time step t = 8. The width of the uncertainty intervals is not uniform over the states, since it depends by the uncertainty assigned to the each element of the interval transition matrix. Among the states, state s 1 is the one that crosses the defined threshold at time step t = 8.  In the scope of resilience-based analysis, using interval mathematics can help handling cases in which fault and repair rates are expected to change greatly and rapidly due to an external HILP event. In such a case, this tool would be necessary to take the necessary preventive, control and mitigating measures.
Furthermore, in case of resilience analysis, the conservatism is not a big concern such as in case of reliability studies, because the HILP events can cause enormous damage. Thus a faster, less accurate but more robust and conservative solution, can be accepted.
In this case, the results show that the system is more likely to be in state s 1 , that corresponds to each component working properly: s1 = {1, 1, 1, 1}. However there are some states, such as s 2 , s 3 , s 5 , whose probabilities' upper bound has a non negligible value. Such states refer to the fault of a single component. The analysis suggests that the state s 4 , corresponding to the single fault of the third element s 4 = {1, 1, 0, 1} is much less probable than the other possible single faults. Operatively, this means that the system operator can focus on components other than the third, which appears to be robust. Finally, the state s 9 , which is: {1, 0, 1, 0} also has a non-negligible upper bound probability, although the fault in this case is double. This result means that reinforcing strategies on components two and four are even more urgent.

Future Development
The analysis of the obtained results, and, in particular, the solution bounds reported in Figure 6, demonstrates as the performance of the proposed algorithm may degenerate in the presence of correlated uncertain sources. This phenomenon is mainly related to the inability of the interval-based operators in managing correlated uncertain variables. To face this issue, in this paper different mathematical operators aimed at computing reliable enclosures of the products between uncertain interval variables have been tested, and a novel product operator has been conceptualized. The main idea is to express each uncertain interval variable by an affine formx, which is a first order polynomial defined as follows: where x 0 , [x 1 . . . x n ] and [ 1 . . . n ] ∈ [0, 1] are the central value, the partial deviations and the noise symbols of the affine formx, respectively. The noise symbols are symbolic variables, which represent independent uncertainty sources affecting the uncertain variable x, and whose value is only known to be in the interval [−1 ,1] [38]. The connection between affine forms and interval variables is straightforward, once the radius of the affine form R(x) is defined: The multiplication between affine forms can be defined as: This non-affine operation is traditionally approximated by bounding the second order terms: Starting from this result a new noise symbol describing the approximation error (i.e., endogenous uncertainty) can be defined as follows: Unfortunately, the adoption of this operator seems inadequate in reliably solving the problem formalized in (16), since it lacks in describing the propagation of the approximation errors, and the number of noise symbols describing the endogenous uncertainty increases at each iteration. Hence, a new method for computing multiplication between affine forms is here proposed, which is based on the following robust formalization: Compared to the traditional approach, this operator does not take into account the approximation error by introducing a new noise symbol, rather than spreading the corresponding endogenous uncertainty on the partial deviation of the primitive noise symbols, which are no longer numbers but affine forms. These variables, which are here referred as second order affine forms, can be properly managed in order to reliably solve the uncertain problem (16), considering all the possible combinations of the exogenous uncertainty, and the worst case instance of the approximation errors. The robustness of the computed solutions, which guarantee the rigorous inclusion of the real solution sets, is currently under investigation by the Authors.

Conclusions
Resilience-based planning must be integrated to reliability studies of power systems, in order to assure their security and robustness. Markov Chains can be an efficient tool to analyze power system's reliability and resilience over time, however their accurate description requires a huge amount of data, that is not always present, especially in case of HILP events.
When studying the impact of HILP events, since there are not enough data available, transition matrices should not be defined as point values, but as intervals of uncertainty. This paper employed interval Markov Chains to find robust bounds of the state probabilities over time. The validity of the method was assessed through Monte Carlo simulations. The results show that interval Markov Chains always result in outer bounds, enclosing all possible scenarios of the evolution of state probabilities. The conservatism of such bounds is related to the dependency between the uncertainty sources affecting the transition probabilities.
As the case study showed, incorporating uncertainties into resiliency planning tools can improve the confidence associated to the chosen mitigating actions. The results showed that an uncertainty of 10% on the parameters affecting state transitions can lead to largely different state probabilities, which affect the calculation of risk and resiliency indices and the resulting mitigating actions.
Furthermore, the computational burden of this technique, compared to Monte Carlo, allows the simulation of large scale systems.
The future research direction, which is being explored by the authors, is the application of the described technique, enhanced by affine-based operators, to real world data-sets provided by transmission system operators.
Author Contributions: A.P. and A.V. conceived and designed the methodology; A.P. performed the experiments; A.P. and A.V. analyzed the data; A.P., A.V., L.L.L. wrote the paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.