Modeling reliability of power systems substations by using stochastic automata networks

Reliability Engineering and System Safety 157 (2017) 13–22

Contents lists available at ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

Modeling reliability of power systems substations by using stochastic automata networks Mindaugas Šnipas a,n, Virginijus Radziukynas b, Eimutis Valakevičius a a b

Kaunas University of Technology, Department of Mathematical Modeling, Lithuania Lithuanian Energy Institute, Laboratory of Systems Control and Automation, Lithuania

art ic l e i nf o

a b s t r a c t

Article history: Received 30 September 2015 Received in revised form 22 April 2016 Accepted 12 August 2016 Available online 18 August 2016

In this paper, stochastic automata networks (SANs) formalism to model reliability of power systems substations is applied. The proposed strategy allows reducing the size of state space of Markov chain model and simplifying system specification. Two case studies of standard configurations of substations are considered in detail. SAN models with different assumptions were created. SAN approach is compared with exact reliability calculation by using a minimal path set method. Modeling results showed that total independence of automata can be assumed for relatively small power systems substations with reliable equipment. In this case, the implementation of Markov chain model by a using SAN method is a relatively easy task. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Reliability modeling Markov chain Stochastic automata network Power system Substation

1. Introduction Markov chain modeling has been applied in power system reliability modeling for a long time [1,2]. It has a well-developed mathematical apparatus and can describe complex behavior of a system [3]. Two alternatives to Markov models are common in power system reliability modeling: methods based on the independence of system components (e.g., fault tree analysis, reliability block diagrams, etc.) and Monte Carlo simulation. Methods which assumes total independence of system components are simpler, but even these techniques can be cumbersome, and classical methods, such as fault tree analysis [4] or generation of minimal set paths [5], require computer assistance and use of special algorithms. These methods can also be combined in modeling of Markov chains [6,7]. A Monte Carlo simulation method is often applied in reliability modeling and has been used to simulate various types of objects, such as substations [8], power plants [9] or standard and composite generation and transmission systems [10,11]. However, simulation of rare events (which is common in reliability modeling) requires special attention [12,13], and Markov chain models tend to be more efficient if a system is not too large and highly reliable [14]. Moreover, since by nature Monte Carlo method is based on n

Corresponding author. E-mail addresses: [email protected] (M. Šnipas), [email protected] (V. Radziukynas), [email protected] (E. Valakevičius). http://dx.doi.org/10.1016/j.ress.2016.08.006 0951-8320/& 2016 Elsevier Ltd. All rights reserved.

random experiments, it is beneficial to be able to compare its results to other methods, for which Markov chain modeling can be used [15,16]. Markov chain modeling allows evaluating complex system behavior [17,18]; however, the main drawback is rapid growth of system states and, consequentially, rising complexity. Therefore, most examples in Markov chain reliability modeling deal with relatively small state space [17,19]. For a more complex system, special software can be applied [20,21], but there is still a need for systematic approach, which would allow for simplified description of complex systems. Some proposed strategies are based on system decomposition into smaller independent subsystems [22,23], but the problems might arise if subsystems are interdependent. One of the methods suitable for complex Markov chain model creation is stochastic automata networks (SANs) formalism [24]. SAN formalism applies Kronecker algebra operations, which enables to store infinitesimal generator matrix in compact format; therefore, it is especially suitable for solving the problem of dimensionality. It also allows for systematic description of interaction between smaller subsystems which provides an exact solution. The SAN method was applied to different areas of research. For example, SAN formalism was used to evaluate availability of large-scale computer networks [25], in system theory, SANs were applied in creating the influence model [26]; in cell biology, SANs can be used to model ion channels [27], etc. Kronecker algebra approach is not new in system reliability, but most examples deal with systems of independent components

14

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

[28]. A more sophisticated use of Kronecker algebra with functional transitions in respect to system reliability is considered in [29], though it uses different terminologies than that of SAN formalism. In [30], theoretical k out of n system was specified by using SAN formalism. However, practical application of SANs in solving real-world reliability problems is still uncommon. In [31], SAN formalism is used to specify reliability of cogeneration power plant substation. In this paper, it has been elaborated on these ideas to propose a methodology to model system reliability of power systems substations. System division into different automata and the use of arrowhead matrices are addressed in detail. We also consider the use of SANs together with classic reliability modeling techniques like minimal path sets, failure modes and effect analysis (FMEA) and structure functions. Reliability models of two standard configurations of substation under different model assumptions are created using the proposed methodology. The calculation of system measures and evaluation of independence of individual automata are also considered. The main principles of Markov chain numerical modeling and SANs are introduced in Section 2, while the theoretical background on SAN formalism is presented in Section 3. In Section 4, the application of SANs of real world power systems is addressed. In this chapter, the formation of individual automaton is considered due to preventative circuit breaker actions of power system. The proposed technique leads to automata, whose infinitesimal generators are arrowhead matrices. We demonstrated that arrowhead matrices and SAN formalism allow specifying various scenarios, which are common in reliability modeling. This methodology is applied in two case studies. In Section 5 we present the model of sectionalized bus and in Section 6 – a ring bus substation configuration. Both reliability models are considered in detail under different assumptions, which can be easily implemented by the use of functional transition rates. Functional transition rates are used to specify preventative circuit breaker operation, shared load and repair capacity. In Section 7, it is shown how SAN reliability modeling can be combined with minimal path set methods and structure functions for estimation of system availability. In Section 8, the modeling results are presented. In Section 9, we discuss the implication of modeling results and the conditions, under which the independence of automata can be assumed. In this case, the implementation of SAN modeling becomes a much easier task.

2. Development of Markov chain model In this paper we assume stationary analysis of homogenous irreducible continuous time Markov chain reliability models, with finite number of system states. Development of Markov chain model can be divided into three main stages: 1) Defining the set of states of the system and possible transitions amongst them. 2) Computation of steady-state probabilities. 3) Computation of necessary probabilistic characteristics of the system, using steady-state probabilities. The first step is model specification. For reliability model this could mean defining failure and repair rates, describing possible system reconfiguration after failure, etc. Memoryless property is satisfied under assumption that transition rates among system states are those of exponential distribution. Basically, system specification results to an infinitesimal generator Q of continuous time Markov chain.

The row vector of steady-state probabilities π is the solution of the system of linear equations

π⋅Q = 0.

(1)

If Markov chain is irreducible, Q is singular – its rank is equal to if n is the number of system states. Thus, an additional condition is used

( n − 1),

(2)

π⋅e = 1;

where e denotes a column vector (of size n) consisting of 1. Necessary system measures can be calculated from steady-state probabilities. E.g., suppose that S is the set of all possible states s( i) of a system reliability model, while SF ⊂ S is a subset of states in which the system is failed. Then probability that the system is failed Pr( F ) can be estimated as follows

Pr( F ) =

∑ s ( i )∈ SF

( )

π s ( i) ; (3)

( ) is a steady-state probability that the system is in the

where π s

( i)

state s( i). One of the main problems in Markov chain modeling is a rapid growth of the number of system states. E.g., Markov chain model of a system consisting of k different items (each of them can be at fault or operating) has a state-space of size 2k. Thus a twenty-item system would have more than one million states. Moreover, the number of states doubles if the number of items is increased by one. This phenomenon – the number of states of the Markov chain grows exponentially when the number of system components grows linearly – is called a state space explosion. The state space explosion makes the practical application of Markov chains a difficult task. If the number of system states is large, the model creation and solution of (2) might be challenging. Therefore, the use of efficient Markov chain modeling techniques and numerical methods is crucial.

3. Stochastic automata networks One of the methods which allow creating large Markov chain models is stochastic automata networks (SANs) formalism. SANs allows storing infinitesimal generator Q in a compact form by using tensor (Kronecker) algebra operations. Using SAN formalism the system is described as a few different automata which can interact among themselves. Each automaton is represented by a Markov chain, i.e., a set of states and possible transitions among them. If two automata interact, transition in one automaton may depend on the state in another. The state of the system (global state) is compositional state of all automata. For more information about SANs we refer to [24]. Infinitesimal generator matrix of the whole system (global generator matrix) can be expressed by infinitesimal generators of individual automata, using tensor algebra operations. We recall basic definitions of tensor algebra. Tensor (Kronecker) product A ⊗ B of two matrices A ∈ Rm × n and B ∈ Rp × q is given by

⎛ a11B ⋯ a1nB ⎞ ⎟ ⎜ A ⊗ B = ⎜ ⋮ ⋱ ⋮ ⎟ ∈ R mp × nq. ⎟ ⎜ ⎝ a m1B ⋯ a mnB ⎠

(4)

Tensor (Kronecker) sum A ⊕ B of two squared matrices A ∈ Rm × m and B ∈ Rn × n is given by

A ⊕ B = A ⊗ In + Im ⊗ B ∈ R mn × mn.

(5)

More about Kronecker product with regard to SANs can be

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

found in [32]. If the network consists of k independent stochastic automata A( i) with infinitesimal generators Q( i) , i = 1... k , when an infinitesimal generator of the whole system (so called global infinitesimal generator) Q can be expressed as a tensor sum k

Q = ⊕ Q ( i). i=1

(6)

Expression (6) is called SAN descriptor of the system. If all automata are completely independent, steady-state probability vector π of the whole system is given by k

π = ⊗ π ( i). i=1

(7)

Here π( i) is steady-state probability vector of an individual automaton. The more difficult task is dealing with interaction among automata, if the system cannot be divided into completely independent modules. Plateau and Atif [24] expressed two different ways to describe the interaction among automata: 1) Functional transition rates – a transition rate in a single automaton may depend on the state of the other automata. Transition rates which are independent on the state of other automata are called constant transition rates. 2) Synchronizing events – transition in one automaton can cause a transition in other automata. Transition rates are called local, if they are not transition rates of synchronizing events. Synchronizing transitions may also be functional. The main problem – the steady-state solution cannot be expressed as a simple product form (7). This means that much larger linear system of equations must be solved, in order to find π. However, system descriptor can still be expressed using tensor products and extended tensor algebra concepts for functional transition rates [33]. Special software can also be used for specification of SAN models [34]. Steady-state probabilities can be found either from solving (1) after Q is built from SAN descriptor, or explicitly from SAN descriptor. I.e., building and storing of Q is not necessary, if special numerical methods are applied. These problems are considered in detail in [35].

4. Power system reliability modeling by the use of SANs The reliability models are created under the following three assumptions: 1) each item can be in one of two possible states – operating or failed; 2) a failed item is detected and repair is initiated immediately; 3) the repaired item is as good as new. 4.1. System division into distinct automata Application of SAN formalism requires choosing a proper system division into distinct automata. E.g., consider substation with a sectionalized bus (see Fig. 1): Even if one assumes absolutely reliable line segments, it still consists of 17 different items. Modeling each item as an individual automaton with 2 possible states would result in Markov chain with 217 states. Moreover, specification of all possible interactions among different automata would require a lot of functional transition rates and/or synchronizing events. Therefore derivation of SAN descriptor would be a very complicated task in this case. We propose the following strategy for division into different automata which is more convenient for power system reliability

15

modeling – all items, which are disconnected together in case any of them is at fault, must be modeled by the same individual automaton. In order to choose this division at first we perform FMEA analysis – i.e., we assume preventative circuit breaker action, which is described in Tables 1 and 2. The first column of Table 1 is very important, because it presents items which are disconnected together under repair. Therefore it is beneficial to model them by the same individual automaton. The effect of circuit breaker failure is assumed to be as presented in Table 2. The proposed model specification technique leads to a network consisting of 6 automata. We enumerate them according to Table 1, i.e. L1 and S1 are modeled by the first automaton A( 1) , L2 and S2 are modeled by A( 2) etc. Assigning a circuit breaker to the automaton is somewhat arbitrary – e.g., C1 could be assigned either to A( 1) (together with L1, S1) or to A( 3), which also includes S3, S5, B1 and S7. 4.2. The use of arrowhead matrices The main advantage of proposed strategy – it allows reducing the size of system state space and also the number of interactions among different automata. This follows from the simple assumption that an item cannot fail if it is disconnected under repair. Therefore multiple failures are impossible in a single automaton. In this case an automaton modeling ni different items has infinitesimal generator which is given by T1

S9

C4

S7 B1

S3

C1

S4

C2

S1

L1

S5 C3 S6

T2

S10

C5

S8

B2

S2

Fig. 1. Sectionalized bus configuration of the substation.

Table 1 Preventative disconnection due to the failure of other items. Failed

Disconnected

L1, S1 L2, S2 B1, S3, S5, S7 B2, S4, S6, S8 T1, S9 T2, S10

C1-S3 C2-S4 C1-S1, C3-S6, C4-S9 C2-S2, C3-S5, C5-S10 C4-S7 C5-S8

Table 2 Preventative disconnection due to the failure of circuit breakers. Failed Disconnected C1 C2 C3 C4 C5

C3-S6, C4-S9 C3-S5, C5-S10 C1-S1,C2-S2,C4-S9,C5-S10 C1-S1, C3-S6 C2-S2, C3-S5

L2

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

16

⎛ * ⎜ ⎜μ Q ( i) = ⎜ 1 ⎜ ⋮ ⎜μ ⎝ ni

λ1 … λ n i ⎞ ⎟ * 0 0⎟ ⎟. 0 * 0⎟ 0 0 * ⎟⎠

(8)

The first line in (8) denotes an operating state, while the rest denotes a failure of items number 1, 2, …, ni respectively. Eq. (8) can also be written in the block form as follows

⎛ d Λ1, ni ⎞ 1,1 ⎟; Q ( i) = ⎜⎜ ⎟ ⎝ Μ ni,1 Dni ⎠

(9)

where Λ1, ni – a row-vector consisting of failure rates; Μ ni,1 – a column-vector consisting of repair rates;

Dni = − diag (Mni,1);

d1,1 = − Λ1, nieni,1. The structure of non-zero entries in matrices (8)–(9) has a distinct shape, which is known as arrowhead matrix [36]. In this paper we want to demonstrate that the use of arrowhead matrices and SAN formalism allows describing variety of events common in system reliability. 1) Repairable item – I.e. an item which can be in two possible states: operating and failed. If λ and μ denotes a failure and repair rates respectively, its' Markov chain has an infinitesimal generator

Q=

⎛ −λ λ ⎞ ⎜ ⎟. ⎝ μ −μ ⎠

(10)

It is obvious that (10) is a particular case of (8). 2) Planned maintenance – Suppose that a repairable item can be repaired not only because of a failure, but also to ensure its long-run availability. If λ p and μp denotes planned “failure” and “repair” rates, an item has an infinitesimal generator

⎛− λ − λ p ⎜ Q=⎜ μp ⎜ ⎜ μ ⎝

(

)

λp −μ p 0

λ ⎞⎟ 0 ⎟. ⎟ ⎟ −μ ⎠

(11)

3) System with shared load and repair capacity – Consider a twoitem system with a shared load. I.e. failure and repair rates can vary due to the failure of another item (because of increased load, limited repair capabilities, etc.). Each item can be described by the infinitesimal generator with functional transition rates:

⎛ −f λ f2 λ1 ⎞ 1 ⎟⎟; Q ( 2) = Q ( 1) = ⎜⎜ 2 ⎝ g2μ1 −g2μ1⎠

⎛ −f λ f1λ2 ⎞ ⎟⎟. ⎜⎜ 1 2 ⎝ g1μ2 −g1μ2 ⎠

(12)

Suppose that failure rate increases 50% due to the shared load, and repair rate decreases 20% due to the shared repair capabilities. In that case functional transition rates in (12) can be defined as

⎧ ⎧ ⎪ 1, ⎪ 1, s ( j) = O ; s ( j) = O ; ; gi = ⎨ fi = ⎨ . ⎪ ⎪ j ⎩ 1.5, s ( ) = F; ⎩ 0.8, s ( j) = F;

(13)

where s( j) denotes a state of automaton A( j); i, j = 1, 2, i ≠ j . SAN descriptor of the two-item system is given by

Q = Q ( 1) ⊕ Q ( 2);

systems with shared load. The main difference – functional transition rates are defined as indicators of failure. E.g., imagine a two-item system and suppose that the failure of the second item prevents the first item of failure, while the failure of the first item does not have an effect on the second item. This can be modeled by changing (13) into

⎧ 2 ⎪ 0, s ( ) = F ; f1 = ⎨ f2 = g1 = g2 = 1. ⎪ ( ⎩ 1, s 2) = O;

(15)

5) Generalized arrowhead matrix – Putting together all aforementioned reliability scenarios one can define an arrowhead matrix with functional transition rates (in respect to generalized tensor algebra we call it generalized arrowhead matrix): i ⎞ … fn ⋅λ n(i ) ⎟ i ⎟ 0 0 ⎟ ⎟. 0 ⎟ * ⎟ 0 * ⎟⎠

⎛ i f1⋅λ1( ) ⎜ * ⎜ ⎜ g ⋅μ ( i) * Q ( i) = ⎜ 1 1 0 ⎜ ⋮ ⎜⎜ ( i) 0 ⎝ gni⋅μ ni

(16)

Matrix (16) can also be written in a block form (9), if Λ1, ni , Μ ni,1 and Dni include functional transition rates. SAN consisting of k automata, each of them having infinitesimal generator (16), can be written as k

Q = ⊕ Q ( i). i=1 g

(17)

In the next chapter we present reliability models of standard substation configurations described by SAN with (17) type of descriptor.

5. Reliability modeling of sectionalized bus configuration In this chapter we present the reliability model of substation with sectionalized bus, presented in Fig. 1. We use the proposed methodology for system division into different automata, when each automaton is described by an arrowhead matrix. Three different models were considered: one with independent automata and two models with preventative failures and shared load. In both later cases functional transition rates will be used. 5.1. SAN descriptor of sectionalized bus configuration The repair and failure rates of an incoming line, a disconnector, a circuit breaker, a bus bar and transformer are denoted as λl and μl , λs and μs , λc and μc , λb and μb , λt and μt respectively. We assume that the first automaton A( 1) describes items L1, S1 and circuit breaker C , while the second automaton A( 2) represents 1

items L2, S2 and circuit breaker C2. Since both of these parts consist of identical items, its infinitesimal generators are given by:

⎛* f λ f λ f λ ⎞ Li l Si s Ci c ⎟ ⎜ ⎜μ * 0 0 ⎟ ⎟, i = 1, 2. Q ( i) = ⎜ l ⎜ μs 0 * 0 ⎟ ⎜⎜ ⎟ * ⎟⎠ 0 ⎝ μc 0

(18)

(14)

In (18) fL , fS and fC denotes the functions which model pre-

where ⊕g denotes generalized Kronecker sum, which is used when dealing with functional transition rates [33]. 4) System where failure of one item prevents other from failure – These systems can be described in a similar manner as the

ventative disconnection under the repair. These and other functional transition rates are considered in detail in Sections 5.2–5.4. For the rest of infinitesimal generators Q( i) we use block form (9). In that case Q( i) is defined by vectors Λ( i) and Μ( i) . E.g., the

g

i

i

i

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

third automaton A( 3), which models behavior of the items S3, S5, S7, B1 and C3, is given by

Λ( 3) =

(

fB λb fS λ s fS λ s fS λ s fC λ c 1

3

5

7

3

)

( = (s = (s = (s

1

fS fS

5

T

Μ ( 3) = ( μb μs μs μs μc ) .

(19)

fS

7

The fourth automaton A( 4) includes B2, S4, S6 and S8. It is the only automaton in this SAN model, which does not contain a circuit breaker. Its infinitesimal generator is defined by

Λ( 4) =

(f

λ B2 b

(20)

Similarly, infinitesimal generators of A( 5) and A( 6) may be defined by

) λ );

λ T1 t

fS λ s fC λ c ; 9 4

6

λ T2 t

fS λ s fC 10

(5)

(i ) (1)

∀ i = 1, 2, 5, 6;

3

(27)

5

Functional transition rates associated with A( 4) are

( ) ( ) ( ) = ( s ( ) = 0 ) ∩ ( s ( ) ≠ 5 ) ∩ ( s ( ) ≠ 3) ; = ( s ( ) = 0 ) ∩ ( s ( ) ≠ 3 ) ∩ ( s ( ) ≠ 3) ; = ( s ( ) ≠ 3) ∩ ( s ( ) ≠ 5) ∩ ( s ( ) ≠ 3) .

fS

4

fS

6

fS

8

2

3

6

3

2

6

2

3

6

(28)

Functional transition rates associated with A( 5) are

(

) (

)

fT = s ( 3) = 0 ∩ s ( 1) ≠ 3 ; fS = fC = fT .

5 c

1

T Μ ( i) = ( μt μs μc ) ; i = 5, 6.

(21)

SAN descriptor of the system may be written as

9

4

1

(29)

And finally, functional transition rates associated with automaton A( 6) are as follows:

(

) (

) (

)

fT = s ( 4) = 0 ∩ s ( 3) ≠ 5 ∩ s ( 2) ≠ 3 ;

6

2

Q = ⊕ Q ( i). i=1 g

(22)

State space, generated by (22) is a set of 6-tuples

s=

(4)

(5)

fC = fS .

2

T Μ ( 4) = ( μb μs μs μs ) .

(f Λ( ) = ( f

(1)

fB = s ( 2) ≠ 3 ∩ s ( 3) ≠ 5 ∩ s ( 6) ≠ 3 ;

)

fS λ s fS λ s fS λ s ; 4 6 8

Λ( 5) =

) ≠ 3) ; ≠ 3) ; ≠ 3) ;

fB = s (1) ≠ 3 ∩ s (5) ≠ 3 ; 3

;

) ( = 0 ) ∩ (s = 0) ∩ (s = 0 ) ∩ (s

17

{ ( s )}; i = 1, 6; s = 0, 3; s = 0, 5; s ( i)

j

3

4

= 0, 4.

fS

10

= fC = fT . 5

(30)

2

5.3. SAN of independent automata

(23)

E.g., 6-tuple (0;0;0;0;1;0) denotes a state which means failed transformer T1, while other items are operating. It is easy to see from (23), that the size of global system state space is equal to 42⋅6⋅5⋅42 = 7680.

5.2. Preventative failure rates

The implementation of SAN with independent automata is straightforward – it suffices to change all functional transition rates into constants. In this case it can be achieved by replacing all functions in (18)–(21) by 1. Assuming the independence of automata would also allow representing the entire branch as a single component. For example, it is easy to verify that the matrix of type (9) can be replaced by that of type (10), if respective failure and repair rates are chosen as follows: ni

In this chapter we define functions used in (18)–(21). These functions model preventative circuit breaker operation under failure. All functions have values 0 or 1; therefore they can be defined by the use of logical predicates. At first we define functional transition rates of the first automaton A( 1) . Let us assume that L is disconnected and cannot fail, if

λ=

∑ λi, i=1

μ=

λ n

λ

∑i =i 1 μi

. (31)

i

This transformation allows for a significant reduction of statespace and could be useful if Markov chain reliability model had a large number (millions or more) of states.

1

any of S3, S5, S7, B1, C3 (all belong to A( 3)) or C4 (denoted by state 3 of A( 5)) is under repair. These conditions can be expressed by defining fL as 1

(

) (

fL = s ( 3) = 0 ∩ s ( 5) ≠ 3 1

)

(24)

The same conditions apply to S1 and C1, therefore:

fS = fC = fL . 1

1

(25)

1

Similarly, functional transition rates associated with the second automaton A( 2) can be defined as

(

) (

) (

)

fL = s ( 4) = 0 ∩ s ( 3) ≠ 5 ∩ s ( 6) ≠ 3 ; fS = fC = fL . 2

2

2

2

(26)

Functional transition rates associated with A( 3) can be written as

5.4. Including shared load and repair capacity SAN model described in the previous chapter can be easily modified in order to include shared load and repair capacity. Suppose that failure rate of transformer increases 50%, if another transformer is out (we must emphasize that it is a theoretical example and we do not have statistical data to support these numbers). This may be modeled by multiplying λt in (21) by the function

⎧ ⎪ 1, if s ( j) = 0 ; fT s = ⎨ ( i; j) = ( 1; 6), ( 2; 5). i ⎪ ⎩ 1.5, otherwise ;

(

)

(32)

Suppose that the repair capacity is also limited – repair rate of an item decreases 50% if any other item is under repair. This condition can be included in SAN model by multiplying each repair rate of infinitesimal generators (18)–(21) by the following function:

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

18

⎧ if ∩ s ( j) = 0 = 1 ; ⎪ 1, j≠i gi = ⎨ i = 1, 6. ⎪ ⎩ 0.5, otherwise ;

(

)

(33)

Table 3 Preventative disconnections due to the failure of the items of ring bus configuration. Failed L1, S2, S3, B2 L2, S4, S8, B3 T1, S1, S5, B1 T2, S6, S7, B4

6. Reliability modeling of ring bus substation configuration

Disconnected

Failed

Disconnected

C1-S1, C2-S3, C1-S2, C3-S5,

C1 C2 C3 C4

C2-S4, C1-S1, C1-S2, C2-S3,

C2-S4 C4-S7 C3-S6 C4-S8

C3-S6 C4-S7 C4-S8 C3-S5

In this chapter we consider a reliability model of the ring bus configuration (see Fig. 2). We use the same approach to create SAN reliability model – system division into individual automata is performed according to failure modes and effect analysis, while each automaton is described by an arrowhead matrix with functional transition rates. The circuit breaker operation under the repair is presented in Table 3.

System state space S consists of 4-tuples, and is defined as follows

6.1. SAN descriptor of the system

S = s ( i) ; i = 1, 4; si = 0, 5.

Table 3 suggests the use a SAN of 4 stochastic automata to model reliability of ring bus configuration. We choose a first automaton A( 1) to model items L , B , S , S and C . Its local in-

E.g., a 4-tuple (3;0;0;0) denotes a system state with failed disconnector S2, while other items of the substation are operating.

1

2

2

3

given by

(34)

)

fB λb fS λ s fS λ s fC λ c ; 3 4 8 2 (35)

The third automaton A( 3) models items T1, B1, S1, S5 and C3. Its local infinitesimal generator Q( 3) is defined by vectors

(36)

( 4)

Finally, automaton A which models the part consisting of T2, B4, S6, S7 and C4, is defined by:

λ T2 t

)

3

2

2

fS

3

L1 C1

S2

( ) = f = ( s ( ) = 0 ) ∩ ( s ( ) ≠ 5) ; = ( s ( ) = 0 ) ∩ ( s ( ) ≠ 5) .

(40)

( 2)

A

and

2

3

1

fS

4

4

C2

4

fS

8

3

(41)

S3

C2

S4

(

)

fT = fB = s ( 1) ≠ 5 ; 1

1

(

) (

)

fS = s ( 1) = 0 ∩ s ( 2) ≠ 5 ; 1

(

fS = fC = s 3

( 4)

) (

=0 ∩ s

( 1)

)

≠5 .

(42)

Finally, the functions presented in (37) are defined as follows:

(

)

fT = fB = s ( 3) ≠ 5 ; 2

L2 B2

4

fL = fB = s ( 4) ≠ 5 ;

(37)

S1

2

C1

5

fB λb fS λ s fS λ s fC λ c ; 4 6 7 4

Μ ( 4) = Μ ( 3) .

B1

2

Functions presented in (36) are

Μ ( 3) = ( μt μb μs μs μc ) .

(f

fS

)

T

Λ( 4) =

1

fB λb fS λ s fS λ s fC λ c ; 1 1 5 3

λ T1 t

( ) = f = ( s ( ) = 0 ) ∩ ( s ( ) ≠ 5) ; = ( s ( ) = 0 ) ∩ ( s ( ) ≠ 5) .

fL = fB = s ( 2) ≠ 5 ;

Functions associated with the second automaton presented in (35) may be written as

T Μ ( 2) = ( μl μb μs μs μc ) .

(f

(39)

Functional transition rates of ring bus substation model are defined in the same way, as described in the previous chapter. At first, functional rates associated with the first automaton A( 1) are

models L2, B2, S4, S8 and C8, is described by vectors

Λ( 3) =

(38)

6.2. Preventative failure rates

Similarly to (34) the rest of automata are described by arrowhead infinitesimal generators, therefore it is possible to write them in block matrix form (9). E.g., infinitesimal generator Q( 2) , which

λ L2 l

i=1 g

1

⎛* f λ f λ f λ f λ f λ ⎞ L1 l B2 b S2 s S3 s C1 c ⎟ ⎜ ⎜μ * 0 0 0 0 ⎟ ⎜ l ⎟ ⎜ μb 0 * 0 0 0 ⎟ 1) ( Q =⎜ ⎟. * 0 0 0 ⎟ ⎜ μs 0 ⎜ ⎟ * 0 0 0 ⎟ ⎜ μs 0 ⎜μ * ⎟⎠ 0 0 0 ⎝ c 0

(f

4

Q = ⊕ Q ( i).

( )

finitesimal generator is given by

Λ( 2) =

SAN descriptor of the system may be written as

B3

fS

6

4

(

) (

(

fS = fC = s 7

)

= s ( 3) = 0 ∩ s ( 1) ≠ 5 ; 4

( 2)

) (

=0 ∩ s

( 3)

)

≠5 .

(43)

B4 S5 T1

C3

S6 T2

S7

C4

S8

Fig. 2. Ring bus configuration of the substation.

6.3. SAN of independent automata If independence of different automata is assumed, it can be modeled in the same way as in previous chapter – i.e., all functional transition rates in (40)–(43) must be constants equal to 1.

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

19

In this case probability of the event A1 is given by

6.4. Including shared load and repair capacity

8

We assume the same effect of shared load of the transformer and shared repair capacity of the items at fault as in the previous chapter. In case of ring bus substation configuration, we multiply λt in (36) and (37) by the following function:

⎧ ⎪ 1, if s ( j) = 0 ; fT s = ⎨ ( i; j) = ( 1; 4), ( 2; 3). i ⎪ ⎩ 1.5, otherwise ;

(

)

(44)

Shared repair capacity can be modeled by using the functions defined in (33), if i = 1, 6 is replaced by i = 1, 4 . E.g., vector Μ( 3) may be written as T

Μ ( 3) = ( g3μt g3μb g3μs g3μs g3μc ) .

(45)

7. Estimating system measures of SAN reliability models Stochastic automata networks approach can be easily implemented together with classical reliability methods such as minimal set paths and structure functions.

Pr( A1) = ∪ Pr( Pi ).

Exact calculation of (48) by using general addition theorem requires estimating 28 − 1 = 63 terms. Calculation of Pr( A2 ) would be even more cumbersome task, because it consists of 216 − 1 = 255 different terms. Estimation of these components requires a thorough analysis and it is not a trivial task. 7.2. SAN representation of minimal path sets It is easy to associate a structure function to a minimal path set by using state variables s( i) of SAN descriptor. E.g., P1 of sectionalized bus configuration is successful if automata A( 1), A( 3) and A( 5) are in the operating state. Thus, the following binary function, whose arguments are system state vector s (23) of SAN model, can represent minimal path set P1 of sectionalized bus:

ϕ1( s) =

∏ ( s ( j ) = 0) ;

j = 1, 3, 5. (49)

j

The rest of minimal path sets of sectionalized bus configuration can be represented by the following binary functions:

φ2( s) =

7.1. Minimal path sets

(48)

i=1

∏ ( s ( i) = 0) ∩ ( s ( 3) ≠ 5);

i = 1, 3, 5;

∏ ( s ( i ) = 0 ) ∩ ( s ( j ) ≠ 3) ;

i ≠ 2, 5; j = 2, 5;

∏ ( s ( i ) = 0 ) ∩ ( s ( j ) ≠ 3) ;

i ≠ 1, 6; j = 1, 6.

i

We assume that transformers are connected to power lines. It is easy to see from Fig. 1, that 4 different minimal path sets exists for sectionalized bus configuration:

i, j

φ4( s) =

(50)

i, j

P1 ¼{L1,S1,C1,S3,B1,S7,C4,S9,T1} ∪{S5,C3}; P2 ¼{L2,S2,C2,S4,B2,S8,C5,S10,T2} ∪{C3,S6}; P3 ¼{L1,S1,C1,S3,B1,S5,C3,S6,B2,S8,C5,S10,T2} ∪{C4,S7,S4,C2}; P4 ¼{L2,S2,C2,S4,B2,S6,C3,S5,B1,S7,C4,S9,T1} ∪{C1,S3,S5,S8}.

Similarly we define the structure functions φi i = 1, 8 associated to the minimal path sets Pi of the ring bus configuration. The argument of these functions is state space vector (39).

(

If a fault occurs on S5 or C3, we assume that breakers C1 and C3 will operate thus breaking the continuity of P1. Therefore P1 contains S5, C3. For similar reasons P2 contains S5 and C3; P3 contains C4, S7, S4 and C2; P4 contains C1, S3, C5 and S8. We denote as Ai an event that at least i ( i = 1, 2) transmission lines are connected to a transformer. Thus, probability of the event A1 may be written as

⎛ 4 ⎞ Pr( A1) = Pr⎜ ∪ Pi⎟. ⎝ i=1 ⎠

φ3( s) =

φ1( s) =

∏ ( s ( i) = 0) ∩ ( s ( 2) ≠ 5);

)

i = 1, 3;

i

φ2( s) =

∏ ( s ( i ) = 0) ;

i = 1, 2, 3, 4;

i

φ3( s) =

∏ ( s ( i) = 0) ∩ ( s ( 4) ≠ 5);

i = 1, 2, 3;

∏ ( s ( i) = 0) ∩ ( s ( 1) ≠ 5);

i = 2, 3, 4;

∏ ( s ( i) = 0) ∩ ( s ( 2) ≠ 5);

i = 1, 3, 4;

∏ ( s ( i) = 0) ∩ ( s ( 3) ≠ 5);

i = 1, 2, 4;

i

φ4( s) =

i

(46)

φ5( s) =

i

An exact calculation of (46) by using the general addition theorem requires estimating the probability of 24 − 1 = 15 different terms. Probability of the event A2 can be expressed as

Pr( A2 ) = Pr( P1 ∩ P2 ∪ P1 ∩ P3 ∪ P2 ∩ P4 ∪ P3 ∩ P4 ).

(47)

Similarly, it is easy to see from Fig. 2, that 8 different minimal path sets exists in the ring bus configuration. They can be written as follows: P1 ¼{L1,B2,S2,C1,S1,B1,T1} ∪{S3,C2,S5,C3}; P2 ¼{L1,B2,S3,C2,S4,B3,S8,C4,S7,B4,S6,C3,S5,B1,T1} ∪{L2,S1,C1,S2,T2}; P3 ¼{L2,B3,S4,C2,S3,B2,S2,C1,S1,B1,T1} ∪{ L1,S5,C3}; P4 ¼{L2,B3,S8,C4,S7,B4,S6,C3,S5,B1,T1} ∪{S4,C2,C1,S1,T2}; P5 ¼{L1,B2,S2,C1,S1,B1,S5,C3,S6,B4,T2} ∪{S3,C2, S7,C4}; P6 ¼{L1,B2,S3,C2,S4,B3,S8,C4,S7,B4,T2} ∪{L2,S2,C1, S6,C3}; P7 ¼{L2,B3,S4,C2,S3,B2,S2,C1,S1,B1,S5,C3,S6,B4,T2} ∪{L1,S7,C4 S8, T1}; P8 ¼ {L2,B3,S8,C4,S7,B4,T2} ∪{S4,C2,S6,C3}.

φ6( s) =

i

φ7( s) =

∏ ( s ( j ) = 0) ;

i = 1, 2, 3, 4;

i

φ8( s) =

∏ ( s ( j) = 0) ∩ ( s ( 3) ≠ 5);

i = 2, 4.

i

(51)

7.3. Implementation of system reliability measures estimation It is easy to estimate the values of structure functions, associated with each minimal path set, once steady-state probabilities are calculated. Below we present a pseudocode for estimation of Pr(A1) for a sectionalized bus configuration, assuming the independence of automata. for i1 ¼0:3 for i2 ¼ 0:3

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

20

for i3 ¼0:5 for i4 ¼0:4 for i5 ¼0:3 for i6 ¼0:3 f1: ¼(i1 ¼0)*(i3 ¼0)*(i5 ¼0); f2: ¼(i2 ¼0)*(i4 ¼0)*(i6 ¼0)*(i3‡5); f3: ¼(i1 ¼0)*(i3 ¼0)*(i4 ¼0)*(i6 ¼ 0)*(i3‡3)*(i3‡3); f4: ¼(i2 ¼0)*(i3 ¼0)*(i4 ¼0)*(i5 ¼ 0)*(i1‡3)*(i6‡3); if (f1) or (f2) or (f3) or (f4) pr¼ psþ p(ind) end if end for i6 end for i5 end for i4 end for i3 end for i2 end for i1

Table 5 Sectionalized bus configuration reliability modeling results. 1 Set of parameters Reliability model

Event

Minimal path sets

A1 A2 A1 A2 A1 A2 A1 A2 A1 A2

SAN (independent) SAN (f, prev.) SAN (f, shar.) SAN (f, prev.þ shar.) 2 Set of parameters Reliability model Minimal path sets SAN (independent)

The meaning of variables used in pseudocode: ij denotes the state of automaton A( j), j = 1, 6 ; fi denotes the value of structure function, associated to path Pi is successful; pr denotes the probability Pr(A1); p(ind) denotes a steady-state probability of a system. Number ind denotes an array index, which can be calculated from ij and the number of states of each automaton ni. In this case, ind could be estimated from

(

)

⎛ 6 ⎞ ind = ∑ ⎜⎜ i j ∏ nk ⎟⎟. ⎠ j=1 ⎝ k=j+1 6

(52)

The main advantage over the use of general addition theorem is that it checks all system states, using only standard functions which can easily be implemented in most of programming languages. It does not require any special algorithms based on graph theory [37].

8. Modeling results We calculated probabilities Pr( A1) and Pr( A2 ) for sectionalized bus and ring bus substation configurations with two different sets of failure and repair rates. The first set of estimated model parameters (see Table 4) is statistical data collected by Lithuanian Energy Institute. The second set of parameters was chosen from [38]. In this case we chose higher failure rates and average repair times, in order to get more visible difference between different reliability modeling methods. Probabilities of events A1 and A2 for the sectionalized bus configuration was estimated by using two different modeling approaches: probability calculation of standard minimal set paths method and Markov chain modeling. First, we used minimal set paths and the general addition theorem for exact calculation of Pr( A1) and Pr( A2 ). The probabilities that an item j is operating was calculated from pj = λj⋅8760/t j . Table 4 Failure and repair rates. Set 1

λl λs λc λb λt

Set 2 0.0816 0.0002 0.0067 0.0002 0.0114

tl ts tc tb tt

70 8 48 2 48

λl λs λc λb λt

0.15 0.001 0.08 0.01 0.02

tl ts tc tb tt

100 15 160 6 200

SAN (f, prev.) SAN (f, shar.) SAN (f, prev.þ shar.)

Event A1 A2 A1 A2 A1 A2 A1 A2 A1 A2

Average time per year h min 0 19 3 11 0 19 3 11 0 19 3 11 0 20 3 12 0 19 3 11

s 30 4 30 2 29 1 35 20 30 1

Average time per year h min 13 0 72 1 12 59 71 54 12 55 71 44 13 37 72 46 12 55 71 44

s 49 0 38 54 9 21 15 41 53 37

The rest of reliability models are SANs, whose descriptors are built with different model assumptions. The first model is the SAN of independent automata, while the rest of three models include functional transition rates, defined in previous chapters. Two separate models include: 1) preventative circuit breaker operation; 2) shared load and repair capacity. The last SAN model was built by putting all functional transition rates together (both preventative failures and shared load/ repair capacity). Modeling results are presented in Table 5 (probabilities are interpreted as an average time per year). Two conclusions can be derived from Table 5. First, the difference between the modeling results of minimal set path method and SAN of independent automata seems to be very small. The difference of about 1 and 5 min per year for events A1 and A2 respectively (relative error is 0.0015 in both cases) can be seen only in case of the modes with less reliable equipment. The difference among the SAN models is also not very big. As expected, the probabilities of the model with preventative failure rates are lower than those of SAN with independent automata, but the difference seems to be almost invisible, especially with the first set of model parameters. The effect of shared load and repair capacity is slightly more visible. A relative difference between this and independent automata model is higher than 5% modeling event A1, though absolute difference is only about 1 min per year. Assuming both preventative failure rates and shared load/capacity seems to cancel each other out. Modeling results of the ring bus configuration are presented in Table 6. We did not estimate the probabilities by using general addition theorem for this substation configuration. Modeling results of ring bus substation brings similar conclusions as the sectionalized bus configuration. The effect of functional transition rates in SAN model basically can be seen only in case of higher failure rates and repair time.

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22 Table 6 Ring bus configuration reliability modeling results. 1 Set of parameters Reliability model

Event

SAN (independent)

A1 A2 A1 A2 A1 A2 A1 A2

SAN (f, prev.) SAN (f, shar.) SAN (f,prev.þ shar.) 2 Set of parameters Reliability model SAN (independent) SAN (f, prev.) SAN (f, shar.) SAN (f,prev.þ shar.)

Event A1 A2 A1 A2 A1 A2 A1 A2

Average time per year h min 0 39 2 51 0 39 2 51 0 39 2 51 0 39 2 51

s 5 28 4 28 17 43 16 41

Average time per year h min 25 39 59 8 25 33 59 0 25 57 59 28 25 45 59 13

s 42 24 34 39 34 49 11 4

9. Evaluating independence of automata Modeling results showed a relatively small difference between SANs of independent automata and SANs with functional transition rates. It is obvious that modeling system as a SAN of independent automata is a much easier task than SAN with functional transition rates. Assuming automata independence means not only an easier model solution by (7). System specification is also much easier, because in this case the definition of functional transition rates can be discarded. Therefore it would be beneficial to know a priori if the difference between two model assumptions leads to the significantly different modeling results. It is easy to see that functional transition rate has an effect only when failure occurs. Therefore, probabilities that 2 items are failed might indicate the difference between the models. Suppose that Pri( 2) denotes a probability that two items are failed if the system is modeled by SAN of independent automata, while Prf ( 2) denotes the same probability of SAN with functional transition rates. If the only functional transition rates are those of preventative failures, it is easy to see that Pri( 2) ≥ Prf ( 2). Therefore, if Pri( 2) is too small to have an impact on system reliability measures, there would be no point to consider the functional transition rates. In that case one could just assume independence of automata, which would simplify system specification and solution significantly. An exact estimation of Pri( 2) requires to solve SAN model of independent automata, however, a simpler criterion might be used. We propose a crude estimate of an upper bound for Pri( 2) which may be written as

Pri( 2) ≤ nf 2 ⋅ρ2 ,

(53)

where nf2 denotes a number of states with 2 failed items, ρ = λ max /μmin – i.e., it is a ratio of the highest failure rate and the lowest repair rate. Number nf2 depends on the SAN model. If SAN consists of automata, those infinitesimal generators are given by (16), it can be found from k−1 k−2

nf2 =

∑ ∑ ni nj ; i=1 j=i

(54)

where k denotes the number of automata; ni and nj denotes the number of items which are being modeled by automata A( i) and A( j).

21

Let us consider SAN model of ring bus configuration presented in Section 6. The model is a SAN of 4 automata whose sizes are 6. Thus from (54) we get that nf 2 = 150. If we use the first set of parameters presented in Table 4, from (53) follows Pri( 2) ≤ 0.0000638, which means about 33 min per year on average. However, if one assumes 10 times shorter incoming line (thus 10 lower value of λl ), from (53) to (54) we get Pri( 2) ≤ 0.0000012, which means only about 37 s per year on average. The average time would be even lower if more reliable equipment is used. In that case functional transition rates of preventative disconnection might be discarded with higher certainty. In general Pri( 2) ≥ Prf ( 2) does not hold if shared load and repair capacity is considered. However it is still possible to apply (53) if the values of λ max and μmin are changed accordingly to the functional transition rates. I.e. one must take the highest possible value of λ max and the lowest possible values of μmin from the definition of functions, which models the shared load and repair capacity.

10. Discussion The paper presents the use of SAN formalism in system reliability modeling. The proposed strategy of system division into individual automata allows for simplification of model specification. The use of decomposition techniques in Markov chain reliability modeling is not new and is probably unavoidable for a large system. The applications differ in types of modeled systems (e.g., nuclear power plants [23] or electrical bus networks [22]) and used methodology (Markov models can be specified using together with faulttree analysis [39], Petri nets formalism [40], etc.). The decomposition is natural if subsystems are independent, but it is much more complicated if they interact in any way. SAN formalism allows to describe the behavior of a system in a systematic way as we have demonstrated by the proposed technique in this paper. In general, it is intended that this research serves as an encouragement for a more widespread use of Markov chains and SAN formalism in reliability modeling. However, we have to note that the proposed methodology has some drawbacks. Specification of functional dependencies among different automata can be cumbersome even for a relatively small object, such as standard substation configurations presented in this paper. In the future research, it might be very advantageous to automate the process of system description based on graph theory algorithms, as it was done in [22,39]. Graph theory methods can also be applied for generation of minimal path sets, which, as has been shown in this study, can be merged with SAN formalism and structure functions for estimation of system measures. This could be achieved by adding some labeling schemes of each item with respect to different automata. We also believe that SAN application together with minimal cut sets could be implemented in a similar way. Modeling results showed that the difference between the model of independent automata and the models with functional transition rates is relatively small, especially in models with very reliable equipment. This means that probably the most cumbersome part, i.e. specification of functional transition rates, could be avoided altogether. In this case, assumption of independence is justified quantitatively, and modeling by SAN becomes a relatively easy task. One of the advantages of the proposed methodology is that it leads to a distinct structure, that is, the infinitesimal generator of each individual automaton is an arrowhead matrix. Although this property was not utilized in this study, it might be beneficial for application of efficient numerical methods. This property would be very important for lowering computation time for large models, when the storage of global system matrix and efficient calculation of steady-state probabilities requires special attention.

22

M. Šnipas et al. / Reliability Engineering and System Safety 157 (2017) 13–22

At the moment, we cannot state that the proposed methodology is the most efficient way for the system description. Other techniques might yield a lower number of system states or, alternatively, to simplify specification of functional dependencies between different automata. This study was opted for a trade-off between the size of state space and ability to describe interaction between different automata. In general, two main factors restrict the applicability of any Markov chain based methodology in reliability modeling. First, it is a rapid growth of the number of states in Markov model, which is often referred to as state space explosion. Although SAN formalism is suitable for simplifying this problem, it is not a complete solution. Depending on computational resources even the compact representation of transition matrix is not sufficient when the size of probability vector exceeds RAM capacity. To the best of our knowledge, a general compact form for a probability vector is not yet known; therefore, it remains the bottleneck in Markov chain modeling for large systems. In this case, state space reduction could be applied as it is common in reliability problems [41]. This study did not concentrate on these problems, and the presented examples were intended merely to demonstrate the proposed methodology. However, we believe that the advanced techniques, such as SAN formalism, allow for specification and solution of much larger Markov chain models in reliability problems than is commonly assumed. Second limitation to the application of Markov processes is its inapplicability to non-Markovian systems. In Markovian system, the time after which the system leaves any particular state does not depend on the time already spent in that state. This memoryless property is satisfied only when the time-spans between different events have an exponential distribution. Non-exponential distribution can be approximated by phase type distribution, which is also applicable to SAN formalism [42]. In such a way a non-Markovian system is replaced by a Markovian one, although it would further increase the size of state-space. In addition, this, most likely, would not be suitable in some situations, e.g., to model planned repairs in regular time intervals. In such cases, other approaches, such as the Monte Carlo method, should be applied.

[11]

[12] [13] [14]

[15]

[16]

[17]

[18]

[19] [20] [21]

[22]

[23]

[24] [25] [26]

[27]

[28]

Acknowledgments [29]

This work was supported by Grant (ATE-no. 04/2012) from the Research Council of Lithuania.

[30]

[31]

References [1] deMercado JB. Reliability prediction studies of complex systems having many failed states. IEEE Trans Reliab 1971 R  20:  223-30. [2] Shooman ML, Trivedi AK. A many-state Markov model for computer software performance parameters. IEEE Trans Reliab 1976:66–8. [3] Trivedi KS, Selvamuthu D. Markov modeling in reliability. Wiley Stats Ref: Statistics Reference Online. John Wiley & Sons, Ltd.; 2014. [4] Ferdous R, Khan FI, Veitch B, Amyotte PR. Methodology for computer-aided fault tree analysis. Process Saf Environ Prot 2007;85:70–80. [5] Malinowski J. A new efficient algorithm for generating all minimal Tie-sets connecting selected nodes in a Mesh-structured network. IEEE Trans Reliab 2010;59:203–11. [6] Yong L, Singh C. Reliability evaluation of composite power systems using Markov cut-set method. IEEE Trans Power Syst 2010;25:777–85. [7] Bouissou M, Bon J-L. A new formalism that combines advantages of fault-trees and Markov models: Boolean logic driven Markov processes. Reliab Eng Syst Saf 2003;82:149–63. [8] Billinton R, Lian G. Monte-Carlo approach to substation reliability evaluation. IEE Proc C: Gener Transm Distrib 1993;140:147–52. [9] Moazzami M, Hemmati R, Haghighatdar Fesharaki F, Rafiee Rad S. Reliability evaluation for different power plant busbar layouts by using sequential Monte Carlo simulation. Int J Electr Power Energy Syst 2013;53:987–93. [10] Billinton R, Gao Y, Karki R. Composite system adequacy assessment

[32] [33] [34] [35] [36] [37] [38] [39]

[40] [41]

[42]

incorporating large-scale wind energy conversion systems considering wind speed correlation. IEEE Trans Power Syst 2009;24:1375–82. Pindoriya NM, Jirutitijaroen P, Srinivasan D, Singh C. Composite reliability evaluation using Monte Carlo simulation and least squares support vector classifier. IEEE Trans Power Syst 2011;26:2483–90. Perninge M, Lindskog F, Soder L. Importance sampling of injected powers for electric power system security analysis. IEEE Trans Power Syst 2012;27:3–11. Shortle JF. Efficient simulation of blackout probabilities using splitting. Int J Electr Power Energy Syst 2013;44:743–51. Hou K, Jia H, Xu X, Liu Z, Jiang Y. A continuous time Markov chain based sequential analytical approach for composite power system reliability assessment. IEEE Trans Power Syst 2016;31:738–48. Leite da Silva AM, de Carvalho Costa JG, Chowdhury AA. Probabilistic methodologies for determining the optimal number of substation spare transformers. IEEE Trans Power Syst 2010;25:68–77. Soleymani S, Mosayebian ME, Mohammadi S. A combination method for modeling wind power plants in power systems reliability evaluation. Comput Electr Eng 2015;41:28–39. Karami-Horestani A, Hamedani Golshan ME, Hajian-Hoseinabadi H. Reliability modeling of TCR-FC type SVC using Markov process. Int J Electr Power Energy Syst 2014;55:305–11. Gilvanejad M, Abyaneh HA, Mazlumi K. Fuse cutout allocation in radial distribution system considering the effect of hidden failures. Int J Electr Power Energy Syst 2012;42:575–82. Haghifam MR, Manbachi M. Reliability and availability modelling of combined heat and power (CHP) systems. Int J Electr Power Energy Syst 2011;33:385–93. Guo H, Yang X. Automatic creation of Markov models for reliability assessment of safety instrumented systems. Reliab Eng Syst Saf 2008;93:829–37. Hirel C, Sahner R, Zang X, Trivedi K. Reliability and performability modeling using SHARPE 2000. In: Haverkort BR, Bohnenkamp HC, Smith CU, editors. Proceedings of the 11th international conference on computer performance evaluation modelling techniques and tools. TOOLS 2000. March 27–31, Schaumburg, IL, USA. Berlin, Heidelberg: Springer; 2000. p. 345–349. Habiballah IO, Ghosh-Roy R, Irving MR. Markov chains for multipartitioning large power system state estimation networks. Electr Power Syst Res 1998;45:135–40. Son KS, Kim DH, Kim CH, Kang HG. Study on the systematic approach of Markov modeling for dependability analysis of complex fault-tolerant features with voting logics. Reliab Eng Syst Saf 2016;150:44–57. Plateau B, Atif K. Stochastic automata network for modeling parallel systems. IEEE Trans Softw Eng 1991;17:1093–108. Brenner L, Fernandes P, Fourneau J-M, Plateau B. Modelling Grid5000 point availability with SAN. Electron Notes Theor Comput Sci 2009;232:165–78. Richoux WJ, Verghese GC. A generalized influence model for networked stochastic automata. IEEE Trans Syst Man Cybern Part A: Syst Hum 2011;41:10– 23. DeRemigio H, LaMar MD, Kemper P, Smith GD. Markov chain models of coupled calcium channels: Kronecker representations and iterative solution methods. Phys Biol 2008;5:14. El-Damcese MA. Markovian-model for systems with independent units. Microelectron Reliab 1995;35:1059–61. Tombuyses B. Automatic construction of Markov transition matrices under state and functional dependences. Reliab Eng Syst Saf 1999;64:49–58. Khatab A, Nahas N, Nourelfath M. Availability of K-out-of-N: G systems with non-identical components subject to repair priorities. Reliab Eng Syst Saf 2009;94:142–51. Valakevicius E, Snipas M, Radziukynas V. Markov chain reliability model of cogeneration power plant substation. Elektron Elektrotech 2013;19:61–6. Langville AN, Stewart WJ. The Kronecker product and stochastic automata networks. J Comput Appl Math 2004;167:429–47. Benoit A, Fernandes P, Plateau B, Stewart WJ. On the benefits of using functional transitions and Kronecker algebra. Perform Eval 2004;58:367–90. Hillston J, Kloul L. Formal techniques for performance analysis: blending SAN and PEPA. Form Asp Comput 2007 19:3-33. Stewart WJ, Atif K, Plateau B. The numerical-solution of stochastic automata networks. Eur J Oper Res 1995;86:503–25. Vanhuffel S, Park HS. Efficient reduction algorithms for bordered band matrices. Numer Linear Algebra Appl 1995;2:95–113. Vega M, Sarmiento HG. Algorithm to evaluate substation reliability with cut and path sets. IEEE Trans Ind Appl 2008;44:1851–8. Xi-Fan W, Xiu-Li W. Generalized outage table approach to reliability study of power station's bus schemes. IEEE Trans Power Deliv 1997;12:1480–5. Huang C-Y, Chang Y-R. An improved decomposition scheme for assessing the reliability of embedded systems by using dynamic fault trees. Reliab Eng Syst Saf 2007;92:1403–12. Koriem SM. Fast and simple decomposition techniques for the reliability analysis of interconnection networks. J Syst Softw 1999;45:155–71. Xiangyu P, Letian J, Guozhi X. Reliability modeling method based on state truncation and error bound analysis. In: Proceedings of the 2006 international conference on communications, circuits and systems; 2006. p. 2791–95. Sbeity I, Brenner L, Plateau B, Stewart WJ. Phase-type distributions in stochastic automata networks. Eur J Oper Res 2008;186:1008–28.