Printer Friendly

Markov chain reliability model of cogeneration power plant substation.

I. Introduction

Markov chain is an effective statistical modeling technique which can describe complex behavior of various stochastic systems and has a well-developed mathematical apparatus. Examples of Markov chain models can be found in computer and telecommunication networks [1], engineering [2] or biological systems [3]. Markov chains can also be used in reliability modeling.

However, some examples of Markov chains in reliability modeling deals with relatively small (less than 100 states) systems [4], assumes total independence of model units [5] or does not address model solution [6]. Markov chain models of industrial power systems can have thousands or even millions of states. This means that the use of efficient model specification techniques and fast computation algorithms is very important in reliability modeling.

In this paper, Markov chain is used to model the reliability of cogeneration power plant substation. In order to specify system behavior, stochastic automata networks (SAN) formalism [7] was used. Stochastic automata networks were successfully applied to model the availability of large computer networks [8]. We think that SAN formalism is suitable for specifying reliability models of power systems and estimating performance measures of a system under investigation.

One of the advantages of using Markov chain model is that it allows computing steady state probabilities of all system states, which helps to estimate probabilities of rare events and failure scenarios. This would be a difficult task in performing simulation, and would require a lot of CPU time or implementation of special modeling techniques [9].

The lack of statistical data is an important issue in reliability modeling, since parameter uncertainty can lead to the misestimation of system measures [10]. In this paper, Markov chain model is used to solve a reliability problem with different sets of parameters, which allows performing uncertainty and sensitivity analysis [11]. Necessary computations can be executed efficiently using a Markov chain model and iterative solution algorithms.

II. Markov Chain Models and SAN

In this paper stationary analysis of irreducible and homogenous continuous-time Markov chain is performed. Markov chain describes a system as a discrete set of states with possible transitions among them. In realistic models the size of state space can be large (thousands or millions of states), thus numerical modeling techniques must be applied. Numerical analysis can be divided in three main stages.

1) Generation of system states and transition matrix (an infinitesimal generator matrix). In reliability modeling it means specification of possible failure scenarios, failure and repair rates etc. Failure and repair rates can be evaluated statistically and stored in an infinitesimal generator matrix Q.

2) Computation of a steady-state probability vector n from the system of linear equations

[pi] x Q = 0. (1)

Equation (1) can also be interpreted as a left eigenvector problem. Since the infinitesimal generator matrix Q is singular, an additional condition is used, i.e. the sum of all probabilities must be equal to 1. It can be expressed in a matrix form as

[pi] x e = 1, (2)

where e denotes a column vector consisting of units.

Investigating the reliability model of a complex system the size of the infinitesimal generator matrix Q can be quite large. It means that numerical methods should be applied to find a steady-state solution. The standard methods to solve such type of problems are the following: direct algorithms, iterative and projective methods [12].

3) Computation of reliability measures using steady-state probabilities. Usual measures in reliability modeling include system availability, the mean time between system failures, frequency of system failure, etc. For example, system availability (i.e., probability that system is available) can be calculated according to the formula

P avail = [summation over (i [member of] A)] [[pi].sub.i], (3)

where A denotes the set of states in which the system is available.

The third stage requires a thorough reselection and analysis of all system states, though it is not as time consuming as the computation of steady-state solution.

One of the main problems in Markov chain modeling is the rapid growth of system states. For example, Markov chain model of the system, consisting of 10 parallel items (each can be in 2 possible states: failed or operating) has [2.sup.10] states. Moreover, each additional item doubles the size of state space--the phenomenon called state space explosion. This problem must be addressed in system specification and solution.

One of the methods to mitigate the state space explosion is the use of stochastic automata networks (SANs) formalism. SAN method is based on specifying the system by its division into smaller interacting subsystems, called automata. Each automaton [A.sup.(i)] is associated with its own state space and the transition among the states can depend on other automata. Interactions among different automata can be modelled by synchronising events and functional transition rates. The infinitesimal generator matrix Q of the entire network, i.e. SAN descriptor, can be represented as a sum of Kronecker products of infinitesimal generators [Q.sup.(i)] (each describes the behaviour of an individual automaton [13]). The Kronecker product of matrices A [member of] [[Real part][sup.m x n]] and B [member of][[Real part][sup.pxq]] is defined as the matrix


The Kronecker sum C = A [direct sum] B of two squared matrices A [member of] [[Real part][sup.m x m] and B [member of] [[Real part][sup.n x n] is defined as the matrix A [direct sum] B = A [cross product] [I.sub.n] + [I.sub.m] [cross product] B [part of] [[Real part][ x mn]] ([I.sub.n] and [I.sub.m] are identity matrices of size n and m respectively).

Since the Kronecker product is associative [13], i.e. (A [cross product] B)[[cross product] C = A [cross product] (B [cross product] C), the Kronecker sum of k square matrices [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be defined as follows


III. SAN Descriptor of Power Plant Substation Reliability Model

We assume a co-generative power plant substation (0.4/10.5 kV), consisting of two independent blocks. The first block has a single transformer (T1), and the second block has two transformers (T2-T3) connected to a busbar section (B2). Each block consists of transformers (T1-T3), switches (S1-S6), circuit breakers (C1-C3) and busbar sections (B1-B2). These items are connected with lines, which can also fail. Fig. 1. shows the diagram of co-generative power plant substation.

The reliability model of the power plant substation is created with the following assumptions:

1) Each item can be in one of two possible states, i.e. operating or failed;

2) The failed items are detected immediately and repair is initiated;

3) There is no limit on repair capacity;

4) The repaired item is as good as new;

5) An item can not fail if the power is disconnected;

6) The duration of failure and repair times are distributed according to exponential law.

The reliability model of the power plant substation is created with the following assumptions:

1) Each item can be in one of two possible states, i.e. operating or failed;

2) The failed items are detected immediately and repair is initiated;

3) There is no limit on repair capacity;

4) The repaired item is as good as new;

5) An item can not fail if the power is disconnected;

6) The duration of failure and repair times are distributed according to exponential law.

There are 14 different items connected with 12 line segments, it means that a Markov chain reliability model has [10.sup.26], i.e. more than 67 million states. However, most of those states are unreachable under the model assumptions. In this paper, we propose a model description which allows reducing the state space significantly.

Two methods to decrease the number of system states of the Markov chain are applied. Under the 5-th assumption, multiple failures in a consecutive branch are impossible, because electric current is off due to reparation. This allows us to define each cons consecutive branch as a single automaton, which can be described as a Markov chain with arrow shaped transition matrix.

The second method of state space reduction is lumping of the states with identical items. For example, if there are two switches in a consecutive branch, they can be represented by one state instead of two. In this case, identical failure rates add up and repair rate remains the same.

Failure rates of transformers, busbars, switches, circuit breakers and line segments are denoted as [[lambda].sub.t], [[lambda].sub.b], [[lambda].sub.s], [[lambda].sub.c] and [[lambda].sub.1] respectively. Similarly, the repair rates of those items are denoted as [[mu].sub.t], [[mu].sub.b], [[mu].sub.s], [[mu].sub.c] and [[mu].l]

The proposed model description leads to a network consisting of 4 automata. The first automaton [A.sup.(1)] describes the first block (T1-S1-C1-S1-B1). An infinitesimal generator of the first automaton is represented as follows


The first row in (6) signifies an operating state, while the others mean that one item has failed. In the first row of (6) failure rates [2[lambda].sub.s] and [4[lambda].sub.l] denote the lumping of the states. Diagonal elements (denoted as *) are negative sums of all row elements of the matrix (6).

In the model, the second block is described by three automata according to the resulting circuit breaker actions when items are being repaired. The second and the third automata [A.sup.(2)] and [A.sup.(30)] describe two identical parts of consecutive branches in the second block of the substation. [A.sup.(2)] and [A.sup.(3)] consist of (T2-S3-C2) and (T3-S5-C3) respectively. Infinitesimal generator of the second automaton is represented as


Since the second and third automata consists of identical items, so [Q.sup.(2)] = [Q.sup.(3)].

The fourth automaton represents the (S4-B2-S6) part of the power plant substation. The infinitesimal generator matrix of [A.sup.(4)] is represented as

In matrices (7) and (8) f refers to functional transition rates, because failure rates of these automata depend on each other. Each function can be expressed as


and [sA.sup.(i)] denotes the state of the i-th automaton--i.e., each function indicates, if the respective automaton is in the operating state. For example, if the transformer T2 is under repair, the branch (T2-S3-C2) is disconnected, while (T3-S5-C3) can still operate.

Global infinitesimal generator matrix Q of the whole reliability model can be expressed as a Kronecker sum of infinitesimal generators of automata


The subscript g in (10) denotes the generalization of Kronecker sum to matrices with functional transition rates.

The space of states of the system can be described as a 4tuple

([n.sub.1], [n.sub.2], [n.sub.3], [n.sub.4]); [n.sub.1] = [bar.0,5]; = [n.sub.i] = [bar.0,4], i = 2,3,4. (11)

For example, the state (0;1;0;0) means that the second automaton [A.sup.(2)] is in a failed state number 1 (which means that transformer T2 has failed), while every other automata are in the operating state.

IV. Modeling Results

After reducing the state space, the Markov chain reliability model of the power plant has only 750 states. In this case it is possible to store the infinitesimal generator matrix Q in RAM and to solve the model by the direct methods. However, since matrix Q is very sparse (most of its elements are zeros), more efficient approach is the use of sparse storage and iterative solution methods.

One of the advantages of Markov model is that it allows generating all possible system states and calculating steady states probabilities of the rarest failure scenarios. This task would be more difficult using simulation approach since it would require a lot of CPU time [9].

Statistical data collected by the Lithuanian Energy institute and from [14] were used to obtain modeling results. The model parameters (failure and repair rates per year) are presented in Table I.

Steady state solution was calculated with [10.sup.-15] precision, using Gauss-Seidel algorithm. In Table II we present some probable failure scenarios and their steady-state probabilities.

A special property of reliability modelling is the lack of statistical data to evaluate failure and repair rates precisely. Parameter uncertainty can lead to unreasonable conclusions and significant misestimation of system measures. Uncertainty analysis can mitigate the problem, but it requires repetitive model solution in order to estimate system measures with different sets of model parameters. Markov chain models have a certain advantage over some other techniques (e.g. simulation) if iterative solution algorithms are applied to find steady state probabilities.

Time requirements to perform the uncertainty analysis by simulation approach and Markov model is shown in Table III. N denotes the number of different sets of model parameters.

The Markov model has an advantage since steady-state solution vectors [[pi].sup.(i)] are relatively close to each other. Setting a vector [[pi].sup.(1)] as the first iteration step ensures fast convergence for the rest calculations. This can be described using pseudo-code:

1. ([p.sup.(1)], ..., [p.sup.(N)])[left

2. it [left arrow] Firstlteration ();

3.[[pi].sup(1)] [left arrow]SteadyStateProbabilities([p.sup.(1)], it);

4.[s.sup.(i)] [left arrow]SystemMeasures([[pi].sup.(1)]);

5.for i = 2 to N do

6. [[pi].sup.(i)] [left arrow] SteadyStateSolution([p.sup.(i)],

7. [s.sup.(i)] [left arrow] SystemMeasures([[pi].sup.(i)],

8. end for.

In this paper, uncertainty analysis was performed by the use of Markov chain model with 2 sets of parameters (size of each set is 1000) with different distribution laws. The first set of failure and repair rates have uniform distribution (Table IV) with probability density function


Random generation of model parameters and steady- state calculation were performed in C++ Builder, using PC with AMD Athlon 64 X2 dual core processor 4000+ 2.10 GHz and 896 MB of RAM physical address extension. It took 1.514 seconds of CPU time to perform the entire calculation.

System availability was estimated with different sets of model parameters, which allows evaluating the distribution of the system availability (Fig. 2).

It was assumed that system is available, if at least two transformers are operating. In that case, the set of states in which system is available consists of 4-tuples (11), satisfying the following condition

[3.summation over (i=1)] if ([n.sub.i] > 0) + 2 x if ([n.sub.4] > 0) [less than or equal to] 1. (13)

MATLAB statistics toolbox was used to analyze the modeling results. The experiments showed that beta distribution, with probability density function

f(x, a, b) = [GAMMA](a + b)/[GAMMA](a)[GAMMA](b) [x.sup.a - 1] [(1 - x).sup.b - 1], (14)

was the best fit the modeling data (Table V), since it provides the highest log likelihood value ([GAMMA] in (9) denotes the gamma function).

Beta distribution is also more suitable than left bounded distributions (e.g., lognormal or Weibull), because its domain (between 0 and 1) matches the range of probability.

The estimated average system availability is about 0.999081, which means that the power plant is shut down on average for 8 hours and 3 minutes per year.

Chi-square goodness-of-fit test affirmed our distribution fitting results. The estimated p-value was 0.1733, which indicates that null hypothesis, i.e. system availability has the beta distribution with parameters a = 11272.3 and b = 10.3635, can not be rejected at the standard 0.05 significance level.

The estimated distribution of system availability can be used as a part of a larger simulation model, because system availability can now be rapidly simulated by generating beta-distributed random variables.

Other numerical experiments were conducted under the assumption that model parameters are distributed normally. However, since failure and repair rates can not be negative numbers, one sided truncation (left tail) of normal distribution was used, which leads to the probability density function


In order to generate normally distributed model parameters with mean values and variations close to those of uniformly distributed parameters (Table IV), the following formulas were used:

[mu] = a + b/2, (16)

[sigma] = b - a/[square root of (12)], (17)

[x.sub.l] = a/10.

Parameters a and b in (12) mean minimum and maximum values of uniform distribution U(a; b). According to (16)-(18) and Table IV, transformer failure rates have truncated normal distribution with parameters [mu] = 0.015, [sigma] = 0.0058 and [x.sub.l] = 0.0005 .

Truncated normal random numbers were generated using acceptance-rejection method. In this case it requires similar amount of CPU time as the generation of the standard normal random numbers, since the left truncation [x.sub.l] is far from the mean value [mu].

As for the uniformly distributed model parameters, we computed steady state probabilities with 1000 different sets of normally distributed model parameters and estimated the distribution of system availability (Fig. 3).

Similarly as with uniformly distributed model parameters, beta distribution was the best fit (Table VI).

The estimated average probability of system availability is about 0.999091, which means that the power plant is shut down on average for 7 hours and 58 minutes per year: slightly less than in the case of uniformly distributed model parameters case.

Chi-square goodness-of-fit test confirmed our null hypothesis, i.e. system availability has beta distribution with parameters a=11933.8 and b=10.858. In this case estimated p-value 0.3247 is almost two times higher than in the case of uniformly distributed parameters case and it is significantly higher than the standard 0.05 level of significance. This proves that beta distribution provides a good fit to model the system availability.

Sensitivity analysis was performed in order to evaluate the model parameters which significantly contribute to the system availability. For this purpose we measured the correlation between randomly generated failure and repair rates and computed system availability. The following correlation coefficient were estimated: Pearson's r, Spearman's [rho] and Kendall's [tau]. Failure and repair rates which are significantly correlated to the system availability (with significance level 0.05) are presented in Table VII:

The results of correlation sensitivity analysis confirm the intuitive assumptions about the reliability model. The explanation of the negative correlation between system availability and failure rates (or its negative correlation between repair rates) is straightforward. The transformer has the most significant effect, since it has the highest failure rate and the lowest repair rate (which leads to longer average repair time). The fact that the number of line segments in reliability model exceeds the number of any other items could explain its significance.

Correlation sensitivity analysis of normally distributed model parameters and system availability showed similar results (Table VIII).

The same failure and repair rates are statistically significant in both sets of experiments. In both cases Pearson's and Spearman's correlation coefficients have higher values than Kendall's [tau].

V. Conclusions

A Markov chain reliability model of a cogeneration power plant substation was presented. Stochastic automata networks formalism with functional transition rates is suitable to specify the reliability behaviour of a system. The size of state space can be lowered significantly if suitable SAN descriptor is chosen.

The Markov chain reliability model and the iterative algorithms allow estimating rare failure scenarios with high precision. Repetitive model solution with different sets of model parameters can be performed efficiently by the use of iterative methods and Markov chain models. This property was used to perform uncertainty and sensitivity analysis of system availability of the power plant.

The obtained results showed that beta distribution is the best fit to model the system availability. Failure and repair rates of transformers and line segments have the most significant effect on system availability. 10.5755/j01.eee.19.5.1214

Manuscript received February 15, 2012; accepted January 15, 2013. This work was supported by grant (ATE--No. 04/2012) from the Research Council of Lithuania.


[1] R. P. Liu, G. J. Sutton, I. B. Collings, "A New Queueing Model for QoS Analysis of IEEE 802.11 DCF with Finite Buffer and Load", IEEE Trans. Wireless Communications, vol. 9, no. 8, pp. 2664-2675, Aug. 2010. [Online]. Available: TWC.2010.061010.091803

[2] M. K. Chang, S. Y. Lee, "Modeling of Single-Step Power Control Scheme in Finite-State Markov Channel and Its Impact on Queuing Performance", IEEE Trans. Vehicular Technology, vol. 58, no. 4, pp. 1711-1721, May 2009. [Online]. Available: 10.1109/TVT.2008.2004802

[3] A. Sakalauskaite, H. Pranevicius, M. Pranevicius and F. Bukauskas, "Markovian Model of the Voltage Gating of Connexin-based Gap Junction Channels", Elektronika ir Elektrotechnika (Electronics and Electrical Engineering), no. 5, pp. 103-106, 2011.

[4] A. M. Bazzi, A. Dominguez-Garcia, P. T. Krein, "Markov Reliability Modeling for Induction Motor Drives Under Field-Oriented Control", IEEE Transactions on Power Electronics, vol. 27, no. 2, pp. 534-546, 2012. [Online]. Available:

[5] M. A. El-Damcese, "Markovian-model for systems with independent units", Microelectronics and Reliability, vol. 35, no. 7, pp. 1059-1061, 1995. [Online]. Available: 2714(94)00062-S

[6] H. Guo, X. Yang, "Automatic creation of Markov models for reliability assessment of safety instrumented systems", Reliability Engineering & System Safety, vol. 93, no. 6, pp. 829-837, 2008. [Online]. Available:

[7] B. Plateau, K. Atif, "Stochastic Automata Network for Modeling Parallel Systems", IEEE Trans. Software Engineering, vol. 17, no. 10, pp. 1093-1108, Oct. 1991. [Online]. Available:

[8] L. Brenner, P. Fernandes, J. M. Fourneau, B. Plateau, "Modelling Grid5000 point availability with SAN", Electron. Notes Theor. Comput. Sci., vol. 232, pp. 165-178, 2009. [Online]. Available:

[9] J. K. Townsend, Z. Haraszti, J. A. Freebersyser, M. Devetsikiotis, "Simulation of rare events in communications networks", Communications Magazine, IEEE, vol. 36, no. 8, pp. 36-41, 1998. [Online]. Available:

[10] Y. Liang, M. A. J. Smith, K. S. Trivedi, "Uncertainty analysis in reliability modeling", in Proc. of the Annu.Reliability and Maintainability Symposium, 2001, pp. 229-234.

[11] O. Yong, J. B. Dugan, "Approximate sensitivity analysis for acyclic Markov reliability models", IEEE Trans. Reliability, vol. 52, no. 2, pp. 220-230, 2003.

[12] W. J. Stewart, Introduction to the Numerical Solution of Markov Chains. Princeton: Princeton University Press, 1994.

[13] A. N. Langville, W. J. Stewart, "The Kronecker product and stochastic automata networks", Journal of Computational and Applied Mathematics, vol. 167, no. 2, pp. 429-447, Jun. 2004. [Online]. Available:

[14] F. Roos, S. Lindah, "Distribution system component failure rates and repair times--an overview", in Conf. Rec. the Nordic Distribution and Asset Management Conference, Helsinki, Finland, 2004.

E. Valakevicius (1), M. Snipas (1, 2), V. Radziukynas (2)

(1) Department of Mathematical Research in Systems, Kaunas University of Technology, Studentu St. 50, LT-51368 Kaunas, Lithuania

(2) Laboratory of Systems Control and Automation, Lithuanian Energy Institute, Breslaujos St. 3, LT-44403 Kaunas, Lithuania


[[lambda].sub.t]   [[lambda].sub.s]   [[lambda].sub.c]
0.015                   0.0002             0.006
[[mu].sub.t]         [[mu].sub.s]       [[mu].sub.c]
146                      1095               182

[[lambda].sub.b]   [[lambda].sub.l]
0.002                   0.002
[[mu].sub.b]         [[mu].sub.l]
2190                     162


System state   Probability

(0; 1; 0; 0)    0.00010268
(5; 0; 0; 0)    0.00004936
(0; 4; 0; 0)    0.00003702
(0; 0; 3; 0)    0.00003295
(0; 0; 0; 3)    0.00002468



Model parameters                                         time

[p.sup.(1)] = ([[lambda].sup.(1).sub.1], ...,     [t.sup.(1).sub.s]
[[lambda].sup.(1).sub.1, ...,

[p.sup.(i) = ([[lambda].sup.(i).sub.l], ...,      [t.sup.(i).sub.s]
[[lambda].sup.(i).sub.n],                        [approximately equal
[[lambda].sup.(i).sub.l], ...,                   to][t.sup.(1).sub.s]
[[mu].sup.(i).sub.n]), I [member of] [bar 2, N]

Total time: [T.sub.s] [approximately equal to]
N x [t.sup.(l).sub.s]

Markov model

Model parameters                                         time

[p.sup.(1)] = ([[lambda].sup.(1).sub.1],          [t.sup.(1).sub.m]
[[lambda].sup.(1).sub.1], ...,

[p.sup.(i)] = ([[lambda].sup.(i).sub.l], ...,     [t.sup.(i).sub.m <
[[lambda].sup.(i).sub.n],                         [t.sup.(l).sub.m]
[[lambda].sup.(i).sub.1], ...,
[[mu].sup.(i).sub.n]), I [member of] [bar 2,N]

Total time: [T.sub.m] < N x [t.sup.(1).sub.m]


Failure rates        Distribution      Repair rates    Distribution

[[lambda].sub.t]    U(0.005; 0.025)    [[mu].sub.t]     U(50,100)
[[lambda].sub.s]   U(0.0001; 0.0003)   [[mu].sub.s]    U(548; 1644)
[[lambda].sub.c]    U(0.003; 0.009)    [[mu].sub.c]     U(84; 252)
[[lambda].sub.b]    U(0.001; 0.003)    [[mu].sub.b]    U(1095; 3285)
[[lambda].sub.l]    U(0.001; 0.003)    [[mu].sub.l]     U(81; 243)


Distribution fit
Distribution:                  Beta
Log likelihood:               6776.31
Domain:                      0 < y < 1
Mean:                        0.999081
Variance:                  8.13287e-008

Parameter estimation

Parameter                    Estimate
a                             11272.3
b                             10.3635


Distribution fit

Distribution:              Beta

Log likelihood:          6808.47
Domain:                 0 < y < 1
Mean:                    0.999091
Variance:              7.60278e-008

Parameter estimation

Parameter                Estimate
a                        11933.8
b                         10.858


Parameter             r       [rho]     [tau]

[[lambda].sub.t]   -0.8179   -0.8401   -0.6394
[[lambda].sub.l]   -0.1644   -0.1540   -0.1036
[[mu].sub.t]        0.4040    0.3614    0.2490
[[mu].sub.c]        0.1363    0.1450    0.0957
[[mu].sub.l]        0.1940    0.1918    0.1274


Parameter             r       [rho]      [tau]

[[lambda].sub.t]   -0.7991   -0.8100    -0.6176
[[lambda].sub.l]   -0.1605   -0.1608    -0.1080
[[mu].sub.t]        0.3813    0.3713     0.2544
[[mu].sub.c]        0.1175    0.1161     0.0777
[[mu].sub.l]        0.1815    0.1535     0.1038
COPYRIGHT 2013 Kaunas University of Technology, Faculty of Telecommunications and Electronics
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2013 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Valakevicius, E.; Snipas, M.; Radziukynas, V.
Publication:Elektronika ir Elektrotechnika
Article Type:Report
Geographic Code:4EXLT
Date:May 1, 2013
Previous Article:Constructed polynomial windows with high attenuation of sidelobes.
Next Article:The impact of feature extraction and selection on SMS spam filtering.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters