Journal of KONBiN 3(23)2012 ISSN 1895-8281

DOI 10.2478/jok-2013-0038

GLOBAL SENSITIVITY ANALYSIS OF MULTI-STATE MARKOV RELIABILITY MODELS OF POWER EQUIPMENT APPROXIMATED BY POLYNOMIAL CHAOS EXPANSION

ANALIZA GLOBALNEJ WRAZLIWOSCI WIELOSTANOWYCH MODELI NIEZAWODNOSCI MARKOWA DLA URZQDZEN ENERGETYCZNYCH APROKSYMOWANYCH ZA POMOC4 ROZWINI^CIA W CHAOS WIELOMIANOWY

1 2 Claudio M Rocco1, Enrico Zio2

:Universidad Central de Venezuela, Caracas, Venezuela 2Dipartimento di Energia, Politecnico di Milano, Milano, Italy

Abstract: Reliability and availability of electric power system equipment (e.g., generator units, transformers) are often evaluated by defining and solving Markov models. Transition rates among the identified equipment states are estimated from experimental andfield data, or expert judgment, with inevitable uncertainty. For model understanding and to guide validation and confidence building, it is of interest to investigate the effects of the uncertainty in the input transition rates on the output reliability and availability. To this aim, Global Sensitivity Analysis (GSA) can be used for defining importance (sensitivity) indexes that allow a ranking of the transition rates with respect to their influence on the uncertainty in the output. In general, GSA requires a large number of model evaluations. In this paper, a metamodel is defined to estimate the performance index of interest (e.g. reliability or availability). The metamodel is built based on polynomial chaos expansion (PCE), a multidimensionalpolynomial model approximation whose coefficients are determined by evaluating the model in a reduced set of predetermined values ofthe input. The proposed approach is illustrated on a power generating unit.

Keywords: Reliability; Availability; Multi-state systems; Homogeneous Continuous Time Markov Chain; Polynomial Chaos Expansion; Parameter Uncertainty; Sensitivity Analysis

Streszczenie: Niezawodnosc i gotowosc urzqdzen elektroenergetycznych jest czqsto oceniana poprzez definiowanie i rozwiqzywanie modeli tancuchow Markowa. Wspotczynniki prawdopodobienstwa przejscia pomiqdzy zdefiniowanymi stanami urzqdzen sq oceniane na podstawie badan doswiadczalnych i danych otrzymanych dla realnych systemow lub sq przedmiotem oceny ekspertow. W celu zrozumienia istoty modelu, kierowania procesem jego walidacji oraz budowania zaufania nalezy siq zainteresowac zbadaniem wptywu niepewnosci w okresleniu wejsciowych wspotczynnikow przejscia w modelu Markowa na uzyskiwane wyjsciowe wartosci niezawodnosci i gotowosci. W tym celu zostat zdefiniowany metamodel pozwalajqcy na okreslenie wspotczynnikow wptywu na parametry eksploatacyjne (np. niezawodnosc czy gotowosc). Ten metamodel zostat zbudowany w oparciu o rozwiniqcie w chaos wielomianowy, wielowymiarowy modelu aproksymacji wielomianowej, gdzie wspotczynniki modelu sq okreslane poprzez ewaluacjq modelu dla zredukowanego zbioru predefiniowanych wartosci wejsciowych. Zaproponowany sposob jest zilustrowany na przyktadzie bloku elektroenergetycznego.

Stowa kluczowe: Niezawodnosc; gotowosc; system wielostanowy; tancuchy Markowa z czasem ciqgtym; rozszerzenie w chaos wielomianowy; niepewnosc parametrow; analiza wrazliwosci.

1. Introduction

Power system equipment like generating units, transformers and others, show stochastic behaviors (e.g. due to failure occurrences, operational conditions) Markov models have been used to compute system performance indicators such as reliability, availability and performability, for power system equipment [1-6]. In general, the stochastic process of transition among the equipment states is modeled through a state transition rate matrix (the model input parameters) and a set of differential equations is obtained for the probability of state occupancy in time [1,9].

To understand what happens to the probability of a specific state if transition rates vary, different approaches have been suggested, such as the evaluation of partial derivatives [7-8], Monte Carlo Simulation (MC) [10], Interval Arithmetic [11-13], Affine Arithmetic [14], among others. These approaches fail to account for the possible interaction effects among transition rates due to the non-linear terms in the probabilities.

To account for this, in this paper we develop a Global Sensitivity Analysis (GSA) [15-16]. GSA is capable of providing (sensitivity) indexes of the importance of the input parameters with respect to their influence on the variability of the model output, by considering the effect of varying a model parameter (or factor, in SA terminology) while all others are varying as well, thus allowing the exploration of multi-dimensional parameter spaces [17].

The number of model evaluations required for a GSA depends on the number of factors in the model and the sample size. The evaluation of the importance (sensitivity) indexes depends on the complexity of the model and could require very large computational times. To speed up calculations,, surrogate or metamodels, like regression models, artificial neural networks or support vector machines, have been suggested to approximate the original model. In this paper, a particular technique named Polynomial Chaos Expansion (PCE) is adopted to create a metamodel represented by a polynomial of the parameters of the original model (i.e., the transition rates in our case) [18-21]. If the parameter uncertainties are modeled as random variables with known probability density functions, then the coefficients of the resulting polynomial contain "the complete response of the model" [18]. In this way, parameter importance indexes can be assessed directly from the coefficients of the polynomial decomposition. The presentation of the work is structured in the paper as follows. Sections 2 to 4 briefly review the basic concepts of Markov modeling, Sensitivity Analysis and

Polynomial Chaos expansion. Section 5 provides the results of the experimentation on a multi-state generating unit. Section 6 presents conclusions and future work.

2. Markov Modeling

The evolution in time of a system which changes state stochastically as its components change modes of operation at random times (due to degradation, reduced production, failure, repair, etc.) can be mathematically described by a discrete-state, continuous-time Markov chain (CTMC) Z={z(t), t > 0}, with finite state-space E ={1,2,..., n}. Let A(t) be the n x n infinitesimal generator matrix whose generic element ajt) represents the transition rate from state i to state j (i,j e E), aii(t)= Si#J ajt), and Pi(t) be the probability that the system is in state i at time t. The n-dimensional probability vector P(t) is obtained by solving [1]:

P (t) = P(t) A(t) (1)

given the initial conditions for P(0). The set of differential equations (1) can be solved, for example, using Ordinary Differential Equations methods. The numerical values of Pi(t) depend on the values of the elements of the transition rate matrix A(t), and it can be of interest to identify those most influential in terms of the variations that they induce on the state probabilities when their variability is considered, i.e., through a sensitivity analysis. To evaluate the influence of the k different transition rates on Pi(t) means to consider the variability in their values and see what the corresponding variability in the output is.

3. Sensitivity Analysis

Sensitivity analysis (SA) is "the study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input factors" X=(x;, x2, ...., xk), kbeing the number of parameters of the model (in our case, transition rates) [15].

In a model of the form Y = f(x\, x2, ..., xk), where Y is a scalar output and the xi are the k orthogonal input factors treated as random variables (i.e., uncertainty is characterized in terms of a probability density function (pdf)), the variance V(Y) of the output Y can be written as [15] :

V(Y) = +YLTV*+ (2)

i i j>i i j>i l>j

where:

V = V (¿(FX;)) is the main effect (or first-order term) due to xi,

V = V(¿"(fIx^ , Xj))— V _ V is the two-way interaction between xi and xj, and so on.

Main and total order sensitivity can then be defined as [15, 23]:

. VMi! and S„. (3)

' VY VY

where X.;- = (xi, xk) and E_i(7|xi) is the expected value of Y

conditioned on xi, and is thus a function of parameter xi alone. The main index S; is the fraction of variance of the output V(Y) that can be attributed to xi alone, while STi corresponds to the fraction of V(Y) that can be attributed to x; including all its interactions with the other parameters [15]. That means that "the computation of total sensitivity indices is a shortcut to measure the overall effect, or "main effect", and all higher order interactions" [23]. The main index Si is the measure to be used to determine which subset of parameters accounts for most of the output uncertainty (Factor Prioritization setting [16]), while STi is used for identifying the subset of non-influential parameters, i.e., those that, when fixed to any value within their uncertainty range, do not significantly reduce the output variance (Factor Fixing setting [16]). The estimation of S; and STi are approximated by: 1) assuming independence among parameters; 2) using random or special sampling techniques from the joint distribution of the space of input parameters; 3) evaluating the set of samples derived in 2) using the model under study. Techniques such as Sobol or extended FAST can be used to determine S; and STi [15].

4. Polynomial Chaos Expansion

The basic idea of PCE is to approximate the output of a system (i.e., the state probabilities in our case) through an orthonormal polynomial basis in the uncertain parameters (e.g., transition rates) spaces. Wiener [24] was the first to propose the use of Polynomial Chaos (PC), when parameters follow a Gaussian distribution. In this case, the basis is defined through Hermite polynomials. This was extended to generalized polynomial chaos (gPC), to consider other types of distributions [25]. Recently, Oladyshkin [22] defined the arbitrary polynomial chaos (aPC), which generalizes PC and is able to accommodate several arbitrary distributions or distributions of which only few statistical moments are known. Let 4«] define the vector of n random input factors and Q(£)=f(£)

a second-order model function with finite-variance output. PC theory [25] shows that the spectral expansion the model output Q(^) can be approximated as:

o©= X a ® i (4)

where:

ai are the coefficients of the expansion £ are random variables with a pdf defined 0,0 are polynomials of the selected basis

M +1 = [18]

n is the number of random variables of the model and d is the selected polynomial degree.

The polynomial basis is selected based on the pdf of If £ are uniform random variables in [-1,1], then the basis is defined as multivariate Legendre polynomials. Other distributions and the associated polynomials can be used, in the Wiener-Askey scheme [25].

The coefficients ai are unknown. In this paper, they are estimated using a non-intrusive approach, which exploits the solution of the set of differential equations (1), based on a set of transition rate deviates (i.e., the random vector and a numerical procedure for solving up to (M+1) multidimensional integrals. These integrals are approximated by evaluating in a set of selected ni integration values (points) £j (/=1,..,ni). This means that the computational cost for determining the polynomial coefficients are proportional to the number of integration values ni. In this paper, the Petras-technique [26], a special Smolyak-based method is used, as implemented in the Scilab toolbox NISP [27]. As suggested by Crestaux et al. [21], the Smolyak approach is, in general, a convenient method when the number of factors is less than 15, while Monte Carlo methods are "best suited for higher dimensional problems".

After determining the ai coefficient values, the polynomial metamodel is defined. As shown in [20-21], sensitivity indexes, including Si and STi, for all model parameters are evaluated analytically from the polynomial coefficients. So the computational cost for deriving them corresponds basically to the cost required to evaluate numerically a set of multidimensional integrals.

The PCE of dynamical systems, recently proposed in [18], is based on coefficients that also vary as a function of time. This means that at each time instant a set of polynomial coefficients are derived and sensitivity indexes are determined for each parameter in the model.

5. Example of Application

Recently, Lisniaski et al. [2] presented a multi-state model for a coal power generating unit. The proposed multi-state Markov model is based on a 4-space-state model, as shown in Figure 1. Table 1 shows the information on 10 transition rates among states (base case).

Table 1: Transition rates data

a,, 1 2 3 4

1 0.0800 0.0133 0

2 0.294 0.3235 0.0294

3 0 0.0288 0.3558

4 0.0002 0.0001 0.0007

The time-dependent behavior of Pi(t) is evaluated. The system of four first order differential equations is solved using the Runge-Kutta scheme of 4th order with a step size h=0.01. After nearly 80 hours the system reaches the steady state, with steady-state probabilities: n1 = 0.0018; n2 = 0.0008; n3 = 0.0025 and n4 = 0.9949. The sensitivity analysis is illustrated by arbitrarily assuming that the uncertainty of the transition rates can be represented by uniform pdfs with bounds at ± 10 % of the base case values. Sensitivity indexes are derived using the PCE approach, based on a Legendre polynomial basis. The PCE, as well as the sensitivity indexes, are evaluated using the Scilab toolbox NISP [27] and a Smolyak integration technique is selected, based on Petras integration points [26].

As suggested in [18], the selection of the appropriate degree d could be done iteratively, that is, trying different values for d and evaluating the average relative error between the model output and the PCE approximation. For polynomial degree d=2 the PCE is considered sufficient since the maximum error is less than 1.0 %, as shown in Figure 2 and the number of models evaluations is equal to 201 [28]. Figures 3 shows the time-dependent behavior of Si and STi for each parameter considered for each state probability. The title in each figure is coded as sx PETRAS yy z-ww, where: x is the state considered (1 to 4)

yy is the type of sensitivity index presented (Main Effect or Total Order) z is the polynomial degree

ww is the number of model evaluations

PETRAS stands for the option selected in NISP [28] for the determination of the ni integration points.

The behavior of Si or STi for state 1 shows that transition rates ai3 and ai2 are the most important and their importance ranking is maintained during the entire timespan. However, the behavior of Si or STi for states 2 to 4 shows that the ranking of the most important transition rate varies.

For example, consider the behavior of Si for state 2. Up to t=5 hours approximately, transition rate a42 has the highest main index value. For t>5, transition rate a23 accounts for most of the P2(t) uncertainty. From t=20 until the

steady-state, the importance of the transition rates does not change. In this example, the importance of each parameter is dependent not only on the state considered but also on the time instant at which it is evaluated. This last fact is important in the case of multi-state generating models, if short-term reliability assessment is required, as suggested in [2].

Plots of STi and Si of the most important parameters are almost equal and suggest that there are no interactions (i.e., STi«Si). On the other side, the behavior of the less important parameters suggests the presence of high interactions (i.e., Si « 0), but of low intensity.

6. Conclusions

The output of a model depends on the input parameters. In many situations, inputs are not precisely known because subject to uncertainty. It is important to know how uncertainty in the parameters affects the results of the model. Such assessment is possible using global sensitivity analysis, an approach able to analyze the variance of the output and determine the individual effects of each parameter (factor) as well as the possible interactions among them. In practice, global sensitivity analysis requires many evaluations of the model which may pose a problem of computation time.

In order to reduce the computational burden, a surrogate model (metamodel) based on Polynomial Chaos Expansion is proposed in this paper. GSA based on the PCE metamodel is performed for assessing the importance of parameters in a power system component modeled by a state-space Markov approach. The example presented relates to a multi-state generating unit described by a four-state Markov model with 10 uncertain transition rates. The metamodel is based on a PCE of second order whose coefficients are determined using only 201 model evaluations (with a maximum average relative error less than 1.0 %).

The results presented show that the importance of the transition rates depends not only on the state being analyzed but also on the time considered for the evaluation. Future works aim to test some approaches for: 1) evaluating the effects of different probability distributions associated to factors as well as to the initial state vector; 2) considering the convenience of using an integrated importance index able to reflect simultaneously the effects of one factor in several outputs; 3) reducing the number of evaluations without affecting the accuracy of the expansion, for example, through an initial screening approach, as the one suggested in [20].

Fig. 1. The State-space diagram [2]

Petras-2-201

vP 0.8

O 0.7-J

CL> 0.6

-i—i 03 0.5-

0) 0.4

o <w 0.2-

- S1 - S2

i —i- • -

10 20 30 40 50 60 70 80 90 100 t

Fig. 2. Error between the model output and the PCE approximation, using d=2

s2 PetrasTotal Qrder-2-201

SI 321

— __- - --

Fig. 3. Time-dependent behavior for Sh STi for each state

Acknowledgments

The participation of Enrico Zio to this research is partially within the China NSFC

project with grant number 71231001.

7. References

[1] Billinton R., Allan R.: "Reliability Evaluation of Engineering Systems", Second Edition, Plenum Press 1992

[2] Lisnianski A., Elmakias D., Laredo D., BenHaim H.: A multi-state Markov model for a short-term reliability analysis of a power generating unit, Reliability Engineering and System Safety, 2012, 98, ,1-6

[3] Billinton R, Fotuhi-Firuzabad M, Sidhu TS.: Determination of the optimum routine test and self checking intervals in protective relaying using a reliability model. IEEE Trans Power Syst 2002; 17(3).

[4] Seyedi H, Fotuhi M, Sanaye-Pasand M.: An extended Markov model to determine the reliability of protective system. In: Power India conference, IEEE; April 1, 2006. p. 10-12.

[5] Aminifar F, Firuzabad MF, Billinton R.: Extended reliability model of a unified power flow controller. Gener Transm Distrib IET 2007;1(6):896-903.

[6] Sefidgaran M., Mirzaie M., Ebrahimzadeh A.: Reliability model of the power transformer with ONAF cooling, Electrical Power and Energy Systems 35 (2012)97-104

[7] Do Van P., Barros A., Berenguer C.: Reliability importance analysis of Markovian systems at steady state using perturbation analysis, Reliability Engineering and System Safety, 2008, 93, 1605-1615.

[8] Do Van P., Barros A., Berenguer C.: From differential to difference importance measures for Markov reliability models, European Journal of Operational Research, 2010, 513-521

[9] Zio E.: Computational Methods For Reliability And Risk Analysis, Series on Quality, Reliability & Engineering Statistics Vol 14, World Scientific Publishing Company, 2009

[10] Marseguerra M., Padovani E., Zio E., Tarantola S., Saltelli A.: Sensitivity Analysis of a non-linear reliability model. Proceeding of the European Conference on Safety and Reliability ESREL 1998, Norway, 1998, pp. 13951401.

[11] Aperjis, D., White D.C., Schweppe F.C., Mettler M., Merrill H.M.: Energy Strategy Planning for Electric Utilities Part II, Smarte Methodology Case Study. IEEE Transactions on Power Apparatus and Systems 1982, PAS-101(2): 347-355.

[12] Merrill, H.M. and Schweppe F.C.: Energy Strategy Planning for Electric Utilities Part II, Smarte Methodology. IEEE Transactions on Power Apparatus and Systems, 1982, PAS-101(2): 340-346.

[13] Rocco C.M.: Effects of transition-rate uncertainty on steady-state probabilities of Markov models using Interval Arithmetic, Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability April 2012 vol. 226 no. 2 234-245

[14] Rocco C.M.: "Affine Arithmetic for assessing the uncertainty propagation on steady-state probabilities of Markov models due to uncertainties in transition rates", submitted to Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, 2012

[15] Saltelli A., Chan K., Scott E.M.: "Sensitivity Analysis", John Wiley & Sons, Probability and Statistics series, 2000

[16] Saltelli, A., Tarantola, S., Campolongo, F., Ratto, M.: "Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models", John Wiley & Sons, Probability and Statistics series, 2004

[17] Campolongo F., Saltelli A., Cariboni J.: From screening to quantitative sensitivity analysis. A unified approach Computer Physics Communications 182 (2011) 978-988

[18] Haro E., Anstett-Collin F., Basset M., Sensitivity study of dynamic systems using polynomial Chaos, Reliability Engineering and System Safety, 2012, 104 , pp. 15-26

[19] Ghanem R., Spanos P. D., Polynomial chaos in stochastic finite elements, Journal of Applied Mechanics 57 (1990) 197-202.

[20] Sudret B., Global sensitivity analysis using polynomial chaos expansions, Reliability Engineering and System Safety 93 (7) (2008) 964-979.

[21] Crestaux T., Maitre O. L, Martinez J., Polynomial chaos expansion for sensitivity analysis, Reliability Engineering & System Safety 94 (2009) 11611172.

[22] Oladyshkin S., Class H., Helmig R., Nowak W., A concept for data-driven uncertainty quantification and its application to carbon dioxide storage in geological formations, Advances in Water Resources 34 (2011) 1508-1518. doi:10.1016/j.advwatres.2011.08.005.

[23] Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D. Saisana, M., Tarantola, S.: "Global Sensitivity Analysis: The Primer", John Wiley & Sons, 2008

[24] Wiener N., The homogeneous chaos, American Journal of Mathematics 60 (4) (1938) 897-936.

[25] Xiu D., Karniadakis G., The Wiener-Askey polynomials chaos for stochastic Differential equations, Journal of Scientific Computing 26 Volume 24 Issue 2, 2002, 619 - 644.

[26] Petras K., Smolyak cubature of given polynomial degree with few nodes for increasing dimension, Numerische Mathematik, 2002, Volume 93, Number 4, Pages 729-753

[27] Baudin M., Martinez J., Polynômes de chaos sous Scilab via librairie NISP, in: 42 emes Journees de Statistique, 2010.

[28] Heiss F., Winschel V., Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, Volume 144, Issue 1, May 2008, Pages 62-80

C. M. Rocco received the Electrical Engineering and MSc. Electrical Engineering (Power System) degrees from Universidad Central de Venezuela (1980, 1982) and Ph.D. degree from The Robert Gordon University, Aberdeen, Scotland, UK (2000). He is a Full Professor at Universidad Central de Venezuela, currently at Operation Research post-graduate courses. His main areas of research interest are Statistics, Reliability, Evolutionary Multi-objective Optimization and Machine Learning techniques. He has published more than 170 refereed manuscripts related to these areas in technical journals, book chapters, conference proceedings and industry reports.

Enrico Zio (BSc in nuclear engng., Politécnico di Milano, 1991; MSc in mechanical engng., UCLA, 1995; PhD, in nuclear engng., Politécnico di Milano, 1995; PhD, in nuclear engng., MIT, 1998) is Director of the Chair in Complex Systems and the Energetic Challenge of the European Foundation for New Energy of Electricite' de France (EDF) at Ecole Centróle Paris and Supelec, full professor, President and Rector's delegate of the Alumni Association and past-Director of the Graduate School at Politécnico di Milano, adjunct professor at University of Sta\>anger. He is the Chairman of the European Safety and Reliability Association ESR4, member of the scientific committee of the Accidental Risks Department of the French National Institute for Industrial Environment and Risks, member of the Korean Nuclear society and China Prognostics and Health Management society, and past-Chairman of the Italian Chapter of the IEEE Reliability Society. He is serving as Associate Editor of IEEE Transactions on Reliability and as editorial board member in various international scientific journals, among which Reliability Engineering and System Safety, Journal of Risk and Reliability, International Journal of Performability Engineering, Environment, Systems and Engineering, International Journal of Computational Intelligence Systems. He has functioned as Scientific Chairman of three International Conferences and as Associate General Chairman of two others. His research focuses on the characterization and modeling of the failure/repair/maintenance behavior of components, complex systems and critical infrastructures for the study of their reliability, availability, maintainability, prognostics, safety, vulnerability and security, mostly using a computational approach based on advanced Monte Carlo simulation methods, soft computing techniques and optimization heuristics. He is author or co-author of five international books and more than 200 papers on international journals.