Transient stability probability assessment and statistical estimation

Electric Power Systems Research 67 (2003) 21 /33 www.elsevier.com/locate/epsr

Transient stability probability assessment and statistical estimation Flavio Allella, Elio Chiodo *, Davide Lauria Department of Electrical Engineering, University of Napoli ‘‘Federico II’’, Via Claudio 21, 80125 Napoli, Italy Received 14 November 2002; received in revised form 7 January 2003; accepted 15 February 2003

Abstract In the paper, a general analytical method for the probabilistic evaluation of power system transient stability is discussed and a new statistical inference approach for this evaluation is proposed. In particular, the transient stability probability (TSP) is defined and evaluated by taking into account the random nature of both the system loads and the fault clearing times (FCT). The paper is focused upon the aspect of statistical estimation of the TSP */a topic generally neglected in literature */on the basis of the obvious consideration that the parameters affecting the TSP (e.g. mean value and variance of loads, FCTs, etc.) are not known, but must be estimated. New properties of point and interval estimations of the TSP are derived and, in particular, an efficient ‘‘lower confidence bound’’ for the TSP estimation is proposed, based upon a suitable Beta probability distribution. In order to show the feasibility of the proposed approach, a numerical application to the Cigre test network is illustrated. Moreover, extensive Monte Carlo simulations to evaluate the estimator efficiency are performed. In the final part of the paper, also a practical example of possible application to the optimization of system design is illustrated. The application of the method is illustrated and performed by using the potential energy boundary surface method, but the estimation results hold their validity irrespective of the method adopted for the transient stability problem formulation and resolution. # 2003 Elsevier B.V. All rights reserved. Keywords: Power systems; Transient stability; Transient energy function; Probability theory; Statistical inference; Lognormal distribution; Normal distribution; Beta distribution; Confidence bounds; Coverage probability

1. Introduction The relevant transformation affecting the electric power systems environment stimulates the researchers to treat with more and more attention the classical analysis methodologies. In particular, the evolution from the actual power system operation to the deregulated environment involves more restrict stability margins [1 /3]; moreover, transmission lines are more difficult to be built nowadays, thus resulting in a need for transmission system operation with the highest possible levels of efficiency and near their maximum capacity. Transient stability is normally analyzed considering the effect on the system of large disturbances such as faults, loss of load, loss of generation, line switching, etc. [4 /6], and the application of probabilistic techniques for * Corresponding author. Tel.: /39-08-1768-3226; fax: /39-081239-6897. E-mail address: [email protected] (E. Chiodo). 0378-7796/03/$ - see front matter # 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0378-7796(03)00090-7

transient stability analysis was introduced in a series of papers by Billinton and Kuraganty, such as [7]. A complete account of the topic and of the analytical problems */still encountered nowadays*/challenging the researchers may be found in [4], while recent papers dealing with uncertainties in transient stability analysis are [8,9]. It is widely recognized that the stability probability is a fundamental index characterizing the dynamic behavior of the system and so its reliability, which is defined probabilistically [4], but often the assumptions about the assumed probability distributions (e.g. the abused Gaussian distribution) are not fully clarified or discussed. Moreover, the important practical aspect of the estimation of the relevant quantities, such as the transient stability probability (TSP) or other measures of ‘‘stability margin’’, is generally neglected, or dealt with approximate ‘‘sensitivity analyses’’ [4], probably because the numerical estimation of the various stability indexes requires to manage methodologies, which are not very familiar in the area of electric power systems. In the paper, an effort

22

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

is made in order to develop a TSP estimator which allows to obtain the required characteristics of robustness and with better performance with respect to the other ones proposed until now. The estimation problem is dealt with in the framework of the general methods for the TSP evaluation [5,6,10], whereas potential energy boundary surface (PEBS) method [11 /15] is employed for the determination of Vcr due to its inner feasibility. In fact, since this method avoids computing controlling u.e.p., it allows obtaining in a simple way the proposed estimation technique. It has to be highlighted that PEBS method could be not very adequate in some cases, but this does not affect the generality of the proposed estimation technique.

/ Cq a credible disturbance in the time interval under consideration (q/1, . . ., z); / S the event of power system transient stability in the time interval under consideration; / P(S) the probability of transient stability; / Tx the fault clearing time (FCT); / Ty the critical clearing time (CCT). Supposing that the credible disturbances form a set of z exhaustive and mutually exclusive random events, the probability of transient stability, in the fixed time interval, can be expressed by means of the total probability theorem [4,18]: r P(S) 

z X

P(Cq )P(Tx B Ty ½Cq )

(1)

q1

2. Probabilistic approach to transient stability assessment The need for the application of probabilistic techniques for transient stability is motivated by the random nature of: / / / /

the the the the

steady-state operating conditions of the system; fault time of occurrence; fault type and location; fault clearing phenomenon.

In fact, the steady-state operating condition*/which heavily affects stability */strongly depends on the load, which is a random process. This is especially evident in planning studies, as the one considered in the present analysis. Also the time interval by which the fault is cleared (fault clearing time (FCT)) should be regarded as a random variable (rv) [7,9,16,17]. In two previous papers [16,17], some theoretical results from probability theory have been utilized in order to develop an analytical approach to transient stability evaluation of electrical power systems. In these papers, the stability of the system is described by means of the critical clearing time (CCT); the dependence of this parameter upon the load demand and then the instability probability, with respect to the system parameters, are analytically computed, so providing a quantitative tool for devising an adequate planning management. However, these and other papers only hinted at, or lacked, the basic argument of the probabilistic parameters’ estimation, which cannot be avoided in current planning and operation, and will be dealt with in this paper. In this section, the probabilistic tools for TSP evaluation are briefly introduced. At the planning stage, the system operating conditions must be necessarily forecasted: in particular, system loads cannot be exactly predicted and have to be properly considered as rvs. As a general method for the TSP evaluation, it is possible to introduce the following set of basic quantities:

Of course, the probability distributions of the rvs CCT and FCT depend on the occurred disturbance. The introduction of the probability of disturbance occurrence is useful to reduce the number of contingencies to be taken into account, since the most unlikely ones can be omitted thus limiting the computational burden [3,19]. Without losing generality, in the following, a single disturbance C is supposed to occur in the time interval under consideration. By this assumption, relation (1) can be written, omitting*/for sake of notation simplicity */the conditioning event C, as follows: r P(Tx BTy )

(2)

being (Tx, Ty), respectively, the FCT and CCT corresponding to the disturbance C. The model in (2) is formally equivalent to the expression of reliability (hence the symbol r, while R is adopted for its estimator) in a ‘‘stress /strength model’’ [20,21], where Tx is the stress acting on a system and Ty its strength, so that the system does not fail until Tx B/Ty. Stress /strength theory models and estimation properties will in fact be used in the sequel.

3. Critical clearing time evaluation The CCT depends on the electrical power system configuration and on the load values at the instant of fault occurrence, this functional link has to be taken into account in the procedure. To reach this aim, an approach based upon the transient energy function (TEF) method [5,6,15] is presented. This approach permits to appraise the security margins as regards the possibility to lose synchronism. For the application of the above said method, it is necessary: / to single out the TEF function for the post-fault stable system;

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

/ to estimate the stability region of the post-fault stable equilibrium point; / to verify that, at the FCT, the point representing the system operating conditions in the state space belongs to this region, (i.e. Tx B/CCT). Being x/[q, 6]T the state vector of the system, the TEF function V(x) is [5]: V(x)

Ng Ng X 1X Mi 62i  Pi (qi qSi ) 2 i1 i1



N g 1 X

Ng X

i1

ji1

[Cij (cos qij cos qSij )

qi qj



g

Dij cos qij d(qi qj )]Vk Vp

(3)

qSi qSj

Mi MT

PCOA

(6)

The application of the PEBS approach consists of the following aspects: the dot product f(q)T ×/ (q/qS) is negative inside the PEBS and positive outside, and the potential energy Vp(q) attains its maximum just upon this surface. This value represents a good approximation of the critical energy Vcr for the considered disturbance [10]. The critical energy Vcr can be calculated by solving the following equation:  Ng dVp X  6i E2i Gpf ii Pmi dt i1  Ng X pf  (Cpf sin q D cos q ) ij ij ij ij j1;j"i

It is defined as a sum of the system kinetic energy, Vk, independent of the network and associated with the relative motion of generator rotors, and the system potential energy, Vp, associated to the changes in rotor potential energy relative to center of angles (COA), changes in magnetic stored energy of branches and changes in energy dissipated in branches. The stability region V of the post-fault stable equilibrium point is estimated by determining a critical value Vcr of the energy function, so that the region S / {x:V(x)5/Vcr} is placed inside the exact stability region V. In the paper, the determination of Vcr is obtained by applying the PEBS approach, permitting a fast calculation of Vcr. The PEBS is the boundary surface of the estimated stability region S:

(7)

0

with a suitable Taylor’s series expansions around the pre-fault stable equilibrium point [3,10]. The instant of time tPEBS, at which Eq. (7) is satisfied, is the time of PEBS crossing: Vp(tPEBS) /Vcr. Then the CCT can be identified as the instant of time at which the sum of kinetic and potential energies is equal to the critical energy: Vcr Vk (CCT)Vp (CCT)

(8)

that is: Ng 1X

2 

˙ @Vp (q) q @q

(4)

Mi 62i (CCT)

i1

Ng X

Pi (qi (CCT)qSi )

i1

N g 1 X

Ng X

[Cij (cos qij (CCT)cos qSij )

i1 ji1 qi (CCT)qj (CCT)

The PEBS can be constructed drawing, from the postfault stable equilibrium point, a number of rays in every direction in the angle space with COA as reference; searching, along each ray, for the first point where the potential part Vp(q) of the TEF function attains its relative maximum. Hence, the boundary surface of interest can be characterized by joining the points qi obtained on these rays [11,12]. In mathematical terms, the description of the PEBS can be obtained by setting, along the rays emanating from the post-fault stable equilibrium point, the directional derivative of Vp(q) to zero. This leads to the following mathematical description of the PEBS: T

f i (q)Pi Pei 

23

S

f(q) × (qq )0

(5)

where f(q)/(f1,. . .,fNg)T is the vector of the accelerating powers of rotors, during the system fault on conditions, that is:

g



Dij cos qij (CCT)d(qi qj )]

qSi qSj



Ng X

Pi (qi (tPEBS )qSi )

i1



N g 1 X

Ng X

[Cij (cos qij (tPEBS )cos qSij )

i1 ji1 qi (tPEBS )qj (tPEBS )



g

Dij cos qij d(qi qj )]

(9)

qSi qSj

Also in Eq. (9) the computation of the state-space vectors q(t) and 6(t) is obtained by using a limited Taylor’s series expansion around the pre-fault stable equilibrium point. Hence, the problem of the power system transient stability can be formalized by two eq. (7) and (9), providing the critical energy and the CCT.

24

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

4. FCT and CCT modeling and stability probability evaluation by means of Lognormal distributions In this section, the FCT pdf and the consequent stability probability evaluation are discussed. In the relevant literature, the Gaussian assumption for the FCT is widely adopted [9] mainly as a consequence of the lack of experimental data. The Gaussian distribution, however, is not a very adequate and flexible choice, as discussed in [16], and more flexible distributions, such as the Lognormal (LN) one, should be considered in order to describe positive rvs [22]. In this paper, as in [17], the FCT Tx is considered as a LN rv, with the following pdf [18,23]:   1 [ln(t)  a]2 p ffiffiffiffiffi ffi exp t]0 (10) f(t; a; b) bt 2p 2b2 The above pdf is zero for tB/0; the two parameters (a,b; b /0) are, respectively, the scale and the shape parameters of the LN pdf: they are, respectively, equal to the mean value and the standard deviation of the Normal rv: X /ln(Tx). The mean value m and the standard deviation s of the LN distribution, as functions of the shape and scale parameters, can be expressed by the following equations:   1 b2 m exp a ; s mfexp(b2 )1g2 (11) 2 In the following, this kind of rv will be indicated as ‘‘LN(a, b)’’. The LN assumption is not too restrictive since its pdf is very flexible, as it can assume a large variety of shapes. Furthermore, if the shape coefficient b is very small, the LN pdf tends to be symmetrical and can satisfactorily approximate also a Gaussian model [23]. The presence of a right tail in the LN pdf accounts for the possibility of relatively large clearing times compared with the expected value; thus, the LN assumption corresponds to a conservative approach resulting appropriate when the exact distribution, as it often occurs, is not precisely known [16]. The parameters a and b may be determined as functions of the values of m and s, about which some information is generally known in practice: denoting by v /s/m the coefficient of variation of Tx, the relations expressing a and b as functions of m and s can be easily deduced from Eq. (11): 1

b fln(1v2 )g2 ;

aln(m)

b2 2

(12)

Moreover, it has been found that, due to its great flexibility, the LN distribution can often well approximate also the CCT, Ty, in accordance to previous results obtained for the single machine case in [16] and for the multi-machine case in [17], where the (extended)

equal area criterion has been adopted. A confirmation of this result by means of the PEBS method is provided in Appendix B. The analytical description of both the critical and the FCT by means of LN probability distributions permits a fast and efficient evaluation of the system TSP r/ P(Tx B/Ty). In fact, assuming that: / Tx is a LN(a, b) rv with parameters: a /ax, b /bx, / Ty is a LN(a, b) rv with parameters: a /ay, b /by, / Tx and Ty are */as reasonable on physical grounds*/ statistically independent,and introducing the auxiliary Gaussian rvs: / X /ln(Tx), Y /ln(Ty), it is well known [4,18,23,24] that the rv: W YX is still normally distributed, and the mean value and standard deviation of W can be, respectively, computed as follows: mW E[W] E[Y]E[X]  ay ax ; sW (Var[X]Var[Y ])1=2 (b2x b2y )1=2 ; Consequently, denoting: u mW =sW (ay ax )=(b2x b2y )1=2

(13)

by elementary properties of the Normal cumulative distribution function (cdf), the TSP is easily evaluated as: r P(Tx BTy )P(X BY) P(W 0)  1F(mW =sW ) 1F(u) and, finally: r F(u)

(14)

where F(w) is the standard normal cdf, defined over the real line as: w

F(w) 

g



1 pffiffiffiffiffiffi exp(z2 =2)dz 2p

(15)

As a numerical example, from the computations in the Appendix B the following estimated parameters (mean value and S.D.) of the FCT, expressed in seconds, are obtained: my /E[Ty]/0.5281 s; sy /(Var[Ty])1/2 / 0.0112 s, to which correspond, through application of Eq. (12): ay E[Y] 0:6387; b2y Var[Y]4:494104 :

(16a)

Assuming that the CCT follows a LN distribution with mean value and std.: mx /E[Tx]/0.4850 s; sx / (Var[Tx])1/2 /0.10mx, (i.e. a cv of 10%), the corresponding parameters of X are:

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

ax  E[X] 0:7286; b2y  Var[X] 9:960103

(16b)

so: u (ay ax )=(b2x b2y )1=2 0:8811 and then the TSP is evaluated as: rF(u) 0:8109:

25

So, let: X (X1 ; . . . ; Xn ); ¯ where Xk ln(Txk );

k 1; . . . ; n

be a random sample of n elements constituted by the natural logarithms of a FCT sample, which can be obtained from field or laboratory data on the system protection components, relevant to the assumed kind of fault (three-phase fault in the application of the Appendix B, to which the FCT estimates are referred here, but the same approach can be used for every given fault), and let:

This is only a point estimate of the ‘‘true’’ TSP, as discussed in the following section: its value */here reported only to illustrate the computations */may be considered rather small, but it may be realistic in long time planning horizon as considered here1. However, the general analytical relation of the TSP as a function of system parameters allows providing means for its improvement, as it will be exemplified in the last section. By examination of the analytical relation:

where

rF[(ay ax )=(b2x b2y )1=2 ]

Yk ln(Tyk );

it is easy to show that, since F is an increasing function of its argument, then, as intuitive, r increases as the mean CCT my (and so ay, see Eq. (12)) increases, and/or the mean FCT mx decreases (and so ax), for given values of their variances. It may be also verified that the TSP also decreases with the FCT and CCT variances as long as ay /ax, (condition which is always satisfied in practice, of course, recalling that F(w) /0.5 only if w /0). This consideration has the practical implication that, if uncertainty in load values (which entails uncertainty in the CCT) and/or in FCT is neglected, the TSP may be undesirably overestimated.

be a random sample of m elements constituted by the natural logarithms of the a CCT sample, which can be obtained (as in the present study) from simulation of system loads */according to their assumed distributions */and resolution of PEBS model, still considering the given fault. Referring, for easier notation, to estimated quantities with capital letters, the most widely adopted estimators (Ax, Bx, Ay, By) of the above four parameters are */for the well known properties of the Maximum Likelihood (ML) estimation [25 /27] */are given, for the LN variables under study, by [23]:

Y (Y1 ; . . . ; Ym ); ¯

Ax (1=n)

k 1; . . . ; m;

n X

Xk ;

k1

5. Statistical inference for system transient stability probability



Bx  (1=n)

5.1. Introduction

1 Moreover, it has to be taken into account that such conditional TSP is referred (see Appendix B) to a three-phase fault, which is rather unlikely, the overall TSP being evaluated by Eq. (1).

(Xk Ax )2

1=2

(17)

k1

Ay (1=m) In this section, the problem of assessing the properties of statistical estimates of the system TSP in the given time horizon is addressed. The aim is to establish both a point and an interval estimate of the unknown probability: r /F(u), where u is defined by Eq. (13), given that the parameters (ax, bx, ay, by) of the two LN distributions must be estimated on the basis of the available random samples of the CCT (Txk: k /1, . . ., n) and of the FCT (Tyk: k /1, . . ., m).

n X

m X

Yk ;

k1



By  (1=m)

m X

2

(Yk Ay )

1=2

(18)

k1

In practice, such estimators coincide with the sample estimators of the mean values (Ax and Ay) and standard deviations (Bx and By) of the Normal rv: X /ln(Tx) and Y /ln(Ty), and possess some desirable properties, such as the one of consistency [25]. Moreover, the log /mean estimators Ax and Ay are also unbiased estimators of ax and ay, respectively. However, although the properties and also the pdf of the above estimators are known and extensively studied in the statistical literature, the same is not true for the ML estimator of r, which is obviously given [20] by:

26

R F[(Ay Ax )=(B2x B2y )]1=2

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

(19)

For a comprehensive examination of such estimator’s theoretical properties, [20,24] provide excellent summaries, also from a bibliographical point of view. We are here interested only in some of the properties of R useful for the TSP, and also in inferring a possible ‘‘lower confidence bound’’ (LCB) for r, i.e. a rv RL such that: P(RL Br) h

(20)

where h is a given, sufficiently large, probability value, assuring that the ‘‘true’’ TSP; r, is larger than the ‘‘lower bound’’ RL value with confidence level of 100h% [25 / 27]. In the application presented in this paper, the value h /0.90 is chosen (only for sake of example). In practice, the establishment of such a lower limit is very useful in order to verify if the designed system satisfies prescribed reliability requirements, which will be desirably issued in the new energy market by adequate guidelines. The classical statistical approach, in which r is a deterministic */but unknown */quantity, is adopted here. A different paradigm, the well known Bayesian one [26], considers the unknown r as a rv; in [20] also the Bayesian estimation approach is examined and the same topic is discussed, in a hypothesis-test framework, in [28]. The determination of a LCB is not an easy matter in general, and in the case under study no analytical solution exists. In this paper, a value of RL numerically computed by means of a Beta distribution */denoted RLB */is proposed. 5.2. Approximating the point and interval estimate of the TSP by the Beta pdf Given that the estimators Ax, Ay, Bx, By are rv, also the TSP ML estimator, R, given by Eq. (19) is a rv, but its exact distribution cannot be obtained analytically. However, its good properties have been confirmed in the application: for instance, in Table 2 (to be discussed in detail below) the estimation results, obtained from N/ 1e4 Monte Carlo samples of different sample size n and m (for sake of simplicity, only equal sample sizes: m /n are chosen for the present illustration, their common values being reported as n in the tables), and different values of the true TSP2 are reported and it is apparent that the bias of the estimator, defined as: bias (R) / E[R]/r, estimated from the simulation results (by means of the sample mean of the R values for each simulation) is always very small, in practice negligible for large sample sizes (n]/100) as expected on the basis 2 The three different reported values of r correspond to the same numerical values of the mean and std. of the FCT as reported in (16a), and three different mean values of the FCT, i.e. 0.485, 0.450, 0.400 s, with a S.D. of 10% of the mean in each of the three cases. Of course, in the simulation procedures they were considered as unknown.

of ML theory, but also for moderate sample sizes (n/ 30). In [20], large sample properties (valid asymptotically) of the ML estimator, U, of the argument u in Eq. (13): U [(Ay Ax )=(B2x B2y )]1=2

(21)

are derived, which allow to use the normal pdf for the large sample approximation for the pdf of U, and */as a consequence */the following approximate, large sample, 100h% LCB is obtained still by means of the standard Normal cdf F already introduced as follows (see [20] for details): RLN  F[(Ay Ax )=(B2x B2y )1=2 zh sR =(mn)1=2 ] (22) where: s2R K[B2x =nB2y =m0:5H(Ay Ax )2 =(B2x B2y )2 ] H (n1)B4x =n2 (m1)B4y =m2 K (mn)=(B2x B2y ) and zh is the h-quantile of the standard normal pdf, i.e. the (unique) solution of: F(zh )h: For instance, if h /0.90, then zh /1.2816. This will be referred as ‘‘Normal’’ (hence the suffix ‘‘N’’ in the above equation for RL) LCB for r. In the present research work, a different approach is taken, as suggested in the context of reliability estimation in [29], which consists in the approximate characterization of the point estimate R by means of a Beta distribution; here the same approach is extended also for obtaining the LCB estimate RL, with satisfactory results, as shown in the following developments also with comparisons with the above estimate RLN, for different sample sizes and TSP true values. First, it is considered the point estimate R: being this a rv in (0, 1), a reasonable choice for its characterization is the approximation of its true pdf with a suitable Beta distribution, which is very flexible for describing rv in (0, 1), being capable */as well known */of producing a large variety of shapes. The Beta is in fact the most commonly used distribution for describing random reliability values [26,29]. The analytical expression of the Beta pdf, as a function of the values r assumed by a rv R in (0, 1), [18] is: 0 rp1 (1  r)q1 G(pq) (0BrB1) (23) f(r; p; q) @ G(p)G(q) 0 elsewhere where p and q are positive shape parameters. In particular, if p /q /1, a uniform distribution over (0, 1) is obtained. Mean value and variance of the Beta

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

distribution are given by:  p q mB  ; s2B m2B pq p(p  q  1)

p? M?(M?M?2 V?)=V?; (24)

In order to choose the approximating Beta pdf for R, an adequate choice of the two parameters (p,q) must be made, for instance */as proposed in [29] */equating the above Beta statistical parameters (mB, s2B) to the mean value (M) and variance (V) of the estimator R, thus obtaining the following equations: p M(MM2 V)=V;

q p(1M)=M

(25)

The mean value M and variance V of the rv R, however, cannot be obtained from its (unknown) distribution, but an excellent approximation for them is obtained */still following [29] */by expanding the function R /F(U) in Eq. (19) in a Taylor series about the mean up to second-order terms, thus obtaining: M?F(U)0:5f(U)[(Ay Ax )B3 ]L V? f(U)2  [B2x =nB2 B2y =mB2 0:5(Ay Ax )2 B4x =n1 B6 0:5(Ay Ax )2 B4y =m1 B6 ]

(26)

being: m1 m1;

n1 n1;

B (B2x B2y )1=2 ; L B2x =nB2y =m0:5(Ay Ax )2 B4x =n1 B4 1:5B4x =n1 B2 0:5(Ay Ax )2 B4y =m1 B4 1:5B4y =m1 B2 ; F(x) and f(x), respectively, the standard Normal cdf and pdf; U the estimator of (13). In order to verify the validity of the proposed Beta approximation, a Monte Carlo simulation has been performed with N/104 independent random samples, each of the same sizes: n /m /104. The goodness of fit tests [25] gave excellent statistical results for different assumed values of the true TSP r. In Table 1, the values of the Kolmogorov /Smirnov Statistics [25], D, for three values of r are reported, being: Dsup(½Fs (r)FB (r; p?; q?)½; 0BrB 1)

q? p?(1M?)M?

27

(28)

Being the critical value of the test statistics for N/104 equal to D* /0.0107 at significance level 0.20, the Beta hypothesis for R is quite acceptable, with p-value larger than 0.20 [25,27]. This is graphically confirmed by histograms such as the one in Fig. 1, relative to r/ 0.8109, and the plot of the sample cdf: y/Fs(R) against the Beta cdf: x/FB(R), evaluated for the same observed (ordered) values of R, corresponding to the same value of r, reported in Fig. 2: in case of good fitting, the plot should appear as the straight line y /x, as indeed happens; same results were obtained for a wide range of r values, with the Beta approximation becoming better and better (which implies decreasing values of D) as r approaches 1, as also apparent from Table 1. It has been also verified that the Normal approximation */ though also acceptable for a wide range of r values */ does not performs well and in particular the Normal hypothesis has been generally rejected at 20% significance level for high values of r (r /0.99). In fact, the pdf of R becomes more asymmetric, or ‘‘skewed’’ */as intuitive */as r approaches 1, and the Beta pdf is capable of assuming a wide range of ‘‘skewness’’ values (see [22] for a precise definition of skewness), while the Normal pdf skewness is always equal to 0. Based on these results and on properties of interval estimation [25 /27], the same Beta distribution has been used for inference on the LCB, estimated as the (1/h) quantile of the above Beta pdf, i.e.: RLB F1 B (1h; p?; q?)

(29)

i.e. the inverse function of the above Beta cdf FB(x; p?, q?) evaluated in (1/h), or the solution of:

(27)

where Fs(r) is the sample cdf of R [25]; FB(x; p, q), the Beta cdf with parameters p and q, evaluated in x (i.e. the integral of the pdf in Eq. (23) over the interval (0, x));

Table 1 Some Kolmogorov /Smirnov statistics values, D, for the approximation of the sample distribution of R (ML estimate of the TSP) with the proposed Beta distribution for three values of r (true TSP) r D

0.8109 0.0074

0.9469 0.0045

0.9972 0.0037

Fig. 1. Histogram of N /1e4 simulated values of the TSP estimate, R, and its approximation with a Beta pdf with parameters: p?/120.35, q?/28.23.

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

28

Fig. 2. Plot of the sample cdf: y/Fs(R) against the Beta cdf: x/ FB(R;p?,q?) for the same sample as in Fig. 1.

Table 2 Synthesis of the estimation results */obtained from N /1e4 Monte Carlo samples of different sample size n and m (n/m /30, 100, 1000, 1e4) */in terms of bias of the TSP estimator R, estimated values of the 90% LCB according the proposed Beta quantile RLB, and the usual normal approximation RLN n

Bias of R

RLB

RLN

CPB

CPN

30 100 1e3 1e4 r/0.8109 30 100 1e3 1e4 r/0.9469 30 100 1e3 1e4 r/0.9972

7.690e/4 /1.960e/4 /7.630e/5 /5.055e/5

0.7219 0.7654 0.7975 0.8067

0.7318 0.7680 0.7978 0.8068

0.9287 0.9154 0.9058 0.9012

0.8871 0.8952 0.9017 0.9000

/0.0022 /8.267e/4 /5.057e/5 6.142e/6

0.8868 0.9189 0.9396 0.9448

0.8966 0.9224 0.9400 0.9448

0.9268 0.9277 0.9084 0.9019

0.8889 0.8929 0.8973 0.8967

/0.0012 /3.602e/4 /3.704e/5 /5.454e/6

0.9853 0.9927 0.9962 0.9969

0.9864 0.9932 0.9963 0.9969

0.8798 0.9098 0.9139 0.9022

0.8785 0.8902 0.8984 0.8969

Coverage probabilities (CP) are also shown for the Beta (CPB) and the normal (CPN) approximation. The three tables are referred to different assumed values of the true TSP (r /0.8109; r/0.9469; r/ 0.9972) as indicated below each table.

1h FB (RLB ; p?; q?) This LCB estimation performed very well, as shown in Table 2, since the estimated coverage probabilities [30], i.e. the relative frequency with which the property (RLB B/r) is verified in the repeated sampling, are always very close to the desired confidence level 0.90. In particular, they are generally larger than this value, implying conservative estimates, which is a desirable property [30]. In the same table, for sake of comparison, also coverage probabilities for the normal approxima-

tion (RLN) are reported, which also performs well for large sample size (as expected), but is less efficient than the Beta approximation for moderate or medium sample size. For instance, from the first part of Table 2 (case r /0.8109), in correspondence of a sample size n /30, the coverage probability CPB, namely P(RLB B/r), is 0.9287, instead the value obtained for RLN in the same case is CPN /P(RLN B/r) /0.8871 B/0.90. It has to be remarked that, while the sample size n of the FCT sample is limited by practical aspects (limitedness of field data, cost and time requirements of laboratory tests, etc.), in principle the sample size m of the CCT, being its values obtained from simulation of the loads, could be*/apparently */chosen arbitrarily large, being limited of course only by computational time requirements, so to approximate the true CCT distribution to any desired degree of precision. In fact this is not true, as also m has practical and theoretical limitations, being the input data for load simulation (e.g. mean value and variance of each load) of course always known with a certain degree of statistical uncertainty, as they must be estimated from data in a load forecasting procedure, which is in practice the first step of the statistical estimation procedure. In practice, the optimal value of m should be the highest compatible with the statistical errors of load parameters estimates. To summarize the proposed statistical inference method, it is opportune to highlight that it is based only upon the estimation of four parameters (ax,bx,ay,by) of FCT and CCT distributions based on two samples of such rv. The estimation is performed by means of the ML estimators of (17) and (18), which are the inputs for the TSP point estimate, R, of (19) and the LCB estimate, RLB, of (29). The parameters (p,q) of the estimated Beta pdf are easily computable from the above four LN parameters, through the parameters (M?, V?). The whole procedure is very simple and fast, requiring only some elementary algebraic computations and two numerical computations: the inversion of the beta cdf and the evaluation of the standard normal cdf (appearing in the R estimate computation and in those involving the parameters M?, V?): both computations are largely available in all statistical or mathematical software packages. Moreover, the procedure allows also the easy estimation of the whole distribution of R (i.e. the Beta pdf above illustrated), which is a powerful tool for evaluating the efficiency of the estimator (by means, for instance, of the distribution variance and quantiles) and establishing any desired confidence interval and/or hypothesis testing procedure for the TSP. The method is useful also in the design stage, since it allows an easy implementation of a rational procedure for the optimal choice or improvement of system parameters in order to satisfy prescribed constraints to

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

the TSP, as shown in the next and final section of the paper.

that, based on a sample of m/40 sample FCT, the ML estimates of the new assumed LN FCT distribution are4: A?x 0:7896;

Finally, a simple numerical example is presented here for illustrating a practical application of the above methodology to the evaluation of system stability and to its improvement by means*/in the design stage */of an appropriate setting of the protection system parameters, here summarized by the log /mean, Ax, and logS.D., Bx, of the FCT. It is supposed that the same values used in the previous numerical example (see Eqs. (16a) and (16b)) are known, being estimated from the available protection system data for the FCT parameters: B2x 9:960e3

(30)

corresponding to a mean FCT of 0.485 s with a S.D. of 10%. It is also assumed that the following estimates of the CCT parameters, from the simulation of m /1000 load values for each bus and successive PEBS resolution, are obtained3: Ay 0:6388;

By 0:02118

(31)

On the basis of these values, the estimated TSP is given, as easily computed from Eq. (19), by: R 0:8107

(32)

It is supposed that a probabilistic stability constraint is prescribed to the system true TSP, r, summarized by the following statement in terms of LCB: r0:85 with at least 90% confidence level. In other words, it is required that the 90% LCB of r is larger than 0.85. The estimated result in (32) clearly shows that the above constraint cannot be satisfied with the current system parameters’ values: hence, an intervention on the protection system components is performed in order to reduce the FCT values; in particular a new designed (mean) value of the FCT of 0.450 s is considered to be feasible. For obtaining the new point estimate and LCB of the TSP, a random sample of the new protection system FCT is needed, since a random variability naturally exists around the designed value: it is supposed 3

B?x 0:0980

while the values of the CCT parameters, Ay and By, remain of course unchanged as in (31). So, the new estimated TSP value is given by:

6. A practical numerical example

Ax 0:7286;

29

In fact these values were obtained for this example from a Monte Carlo simulation of 1000 random numbers from the LN distribution of the CCT as deduced in the Appendix B, i.e. with mean value 0.5281 s, and S.D. 0.00112 s. 4 In fact these values were obtained from a Monte Carlo simulation of 40 random numbers from a new LN distribution of the FCT, with mean value 0.450 s and S.D. of 10%.

R?0:9338 which shows the improvement of the TSP in terms of point estimate. More important is the LCB estimate, obtained by evaluating the estimated mean value M? and variance V? of the ML estimate. Applying the Eq. (28) as functions of the FCT parameters (A ? x, B ? x) and of the CCT parameters (Ay, By), with the sample size values: m /1000, n /40, the following values are obtained: M? 0:9178;

V? 0:0008

From these values, the following estimated parameters (p?,q?) of the Beta pdf, approximating the pdf of the estimator R, are computed through application of (28): p? 82:00;

q?7:340

The Beta pdf of (23) with such parameters values is shown in Fig. 1. It constitutes the best summary of the estimator’s statistics. In particular, the 90% LCB, RLB, of (29) is numerically computed as: RLB F1 B (0:10; p?; q?)0:8792 This LCB is graphically represented in Fig. 3 by the R-abscissa point, x, corresponding to an area under the pdf curve of 0.10 as measured from the origin to the point x. The result implies that, with the available data and estimation procedure, it can be assessed that r/0.8792 with confidence 0.90, so that the reliability constraint is now satisfied. It is highlighted that the estimation procedure does not require any new resolution of the PEBS equations. If the constraint was not satisfied and no improvement of the FCT was feasible, other more complex design actions could be put under study, e.g. network topological changes in order to obtain larger FCT: also in this case the above said estimation method allows an easy quantitative evaluation of the TSP improvement.

7. Conclusion In the paper an effort has been made to perform a probabilistic transient stability assessment, with reference to planning studies. In particular, the statistical estimation of the TSP is examined starting from the point that the TSP parameters are not known a priori, but have to be estimated from available data. The proposed statistical inference method is based only upon

30

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

Fig. 3. Estimated Beta pdf of the TSP estimate, R, for the numerical example.

the estimation of four parameters of CCT and FCT distributions based on two samples of such rv. The estimation is performed by means of the ML estimators of the parameters, which are the input for the TSP point estimate, R, and the LCB estimate, RLB. The whole procedure has been proven very simple and fast, requiring only some elementary algebraic and numerical computations. Moreover, the procedure allows also the easy estimation of the whole distribution of R (i.e. the Beta pdf above illustrated), which is a powerful tool for evaluating the efficiency of the estimator and establishing any desired confidence interval and/or hypothesis testing procedure for the TSP. As it is illustrated in the last part of the paper, the method could be addressed to the design stage, since it allows an easy implementation of a rational procedure for the optimal choice or improvement of system parameters in order to satisfy prescribed constraints to the TSP, which should become a standard in the new energy market. It is remarked that the proposed estimation approach */illustrated in the paper with reference to the PEBS method */can be supported by every other method for the transient stability problem. Possible extensions of the method include the evaluation of the risk consequent to the event of instability [31,32]. More kinds of fault or other random quantities (such as fault location) can be introduced in the model of Eq. (1), following a general methodology such as the one developed in [32].

Appendix A: Main symbols and acronyms CCT cdf cv FCT Ng f pdf LN ML PEBS rv S.D. r R TEF m, s di Mi v0 vi Pmi Ei Bij, Gij Vk

critical clearing time cumulative probability distribution function coefficient of variation fault clearing time number of generation buses Frequency probability density function Lognormal (distribution) maximum likelihood potential energy boundary surface random variable standard deviation theoretical value of TSP maximum likelihood estimate of the TSP, r transient energy function mean value and standard deviation of a rv rotor angle of the i-th synchronous machine inertia of the i-th synchronous machine synchronous speed of rotation relative angular velocity of the rotor of the i-th generator relative to the synchronous velocity mechanical power of the i-th generator electromotive force behind the transient reactance of the i-th generator components of i-th elements of the reduced short circuit admittance matrix for the network kinetic energy

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

Vp Mr

/

potential energy Ng ai1 Mi total inertia N

d0

qi 6i qij Cij Dij Pi PCOA Pei

/

g aj1 Mi dj Mt

di/d0 vi/v0 qi/qj EiEjBij EiEjGij Pmi/E2j Gij Ng /a i1 (Pmi Pei ) N g /a i1;j"1 (Cij sin qij Dij cos q ij )/ /electric active power of the i-th generator immediately after the occurrence of a fault

Meaning of superscripts s refers to the final post-fault stable equilibrium pf refers to the post fault network configuration

Appendix B: Numerical application for the evaluation of the CCT distribution In order to show the method for the evaluation of the CCT distribution, a numerical application is presented in this section, with reference to the Cigre test system [33], whose single line diagram is reported in Fig. 4. Data are expressed in p.u., with reference power 100 MV A. Series compensation is applied, for sake of computation, to the line 8 /9, with a degree of compensation K /0.4. Assuming that the contingency C is a three-phase to ground fault of the line 6 /8 near the bus 8, the CCT, Ty has been computed by random simulation of the system peak loads, which are properly considered as predicted with a certain degree of uncertainty: this is often accepted when the analysis time horizon extends over many years. More precisely, active peak loads are

31

considered as statistically independent Gaussian rvs; the mean values of these random loads are chosen equal to the constants assigned to the Cigre network loads, and their standard deviations are chosen, for sake of illustration, equal to 10% of the mean values (i.e. the coefficients of variation are all equal to 10%). Reactive powers are attributed to the load buses considering constant power factors (holding those attributed to the test network). This implies that the load reactive powers are Gaussian rvs as well, whose mean values and standard deviations depend on the above said power factors. As a consequence of these assumptions, the CCT cannot be considered as a constant, but it is a rv whose statistical parameters have been evaluated by means of Monte Carlo simulations. In particular, for each of the system active peak loads, N/104 independent random samples have been generated according to the Gaussian distribution with parameters m and cv/s/ m of Table 3. The random Gaussian reactive loads are generated coherently assuming, as previously said, constant power factors: they are reported in Table 3 as well. For each simulation cycle, a random value of the CCT, Ty has been computed. The estimated values, by the simulated N-dimensional sample, of the CCT parameters are the following (the superscript ‘‘e’’ denotes the estimate of the generic parameter): me E[Ty ]0:5281 s;

se 0:0112 s

They constitute the ML estimates [4] of the corresponding theoretical parameters of the assumed LN distribution. It has been found that a LN distribution can well approximate the simulated CCT, Ty: considering indeed a LN distribution with mean value equal to 0.5281 s and standard deviation equal to 0.0112 s, the frequency histogram of the simulated values of the CCT, Ty and the corresponding theoretical LN frequency distribution are superimposed in Fig. 5 to give a visual evidence of the goodness of fitting. A quantitative evaluation of the satisfactory approximation with the LN distribution is shown in the following by the comparison between theoretical and sample-estimated values of the quantiles [25]. For instance, the 0.50-quantile of the CCT, Ty, computed by the theoretical LN distribution, is given by: Table 3 Statistical parameters of the loads at load buses of the test network

Fig. 4. Single line diagram of Cigre test system.

Busbar

m/E[P] (MW)

cv

cos 8

8 9 10

100 230 90

0.10 0.10 0.10

0.8944 0.8461 0.8944

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33

32

Fig. 5. Monte Carlo histogram of the CCT, Ty sample, with Gaussian peak loads, and fitting theoretical Lognormal frequency distribution. Table 4 Theoretical and sample-estimated quantiles from a sample of n/1e4 simulated values of the CCT, Ty Theoretical parameters

N-sample estimates

Abs(Error) (%)

Ty Ty Ty Ty Ty Ty

Tey Tey Tey Tey Tey Tey

0.076 0.000 0.037 0.019 0.039 0.240

Ty

0.50 /0.5280 0.95 /0.5467 0.90 /0.5425 0.10 /0.5138 0.05 /0.5099 0.01 /0.5026

0:50 0:5280;

0.50 /0.5284 0.95 /0.5467 0.90 /0.5423 0.10 /0.5139 0.05 /0.5101 0.01 /0.5038

namely P(Ty B 0:5280) 0:50;

while the population estimate of the same quantile, obtained by the simulated N-dimensional sample, is: Tey

0:50 0:5284;

with a corresponding relative error of 7.6e/004, or 0.076%. Results of this numerical application are summarized in Table 4, whose last column clearly shows the excellent fitting of the LN distribution, with errors less than 0.25% also for the extreme-tail quantiles (i.e. Ty,p with p near to 0 or 1) which are the most sensitive to estimation errors [25]. It has to be highlighted that the PEBS method, chosen for its property of rapidity [6,34,35], in some cases [36,37], mainly when the FCT is much shorter than the actual CCT, can exhibit over-estimation problems, due to drastic system trajectory change after fault clearance. However, as said, this does not affect the generality of the proposed estimation technique, whose further devel-

opment could be, for instance, the application to methodologies such as the one presented in [38].

References [1] K.B. Kilani, R.A. Schlueter, Trends in model development for stability studies in power systems, Electric Power Systems Research 53 (1) (2000) 207 /215. [2] D. Chattopadhyay, D. Gan, Market dispatch incorporating stability constraints, Electric Power and Energy Systems 23 (1) (2001) 459 /469. [3] F. Allella, D. Lauria, A fast optimal dispatch with global transient stability constraint, IEE Proceedings-Generation, Transmission and Distribution 148 (5) (2001) 471 /476. [4] G.J. Anders, Probability Concepts in Electric Power Systems, Wiley, New York, 1990. [5] M.A. Pai, Power System Stability Analysis by the Direct Method of Lyapunov, North Holland Publishing Company, 1981. [6] M. Pavella, P.G. Murthy, Transient Stability of Power Systems. Theory and Practice, Wiley, 1994. [7] R. Billinton, P.R.S. Kuruganty, Probabilistic assessment of transient stability, IEEE Transactions on PAS 100 (5) (1981) 2163 /2170. [8] N.G. Bretas, L.F.C. Alberto, Transient stability of power systems: robustness with respect to parameter uncertainties, IEEE Power Engineering Society Winter Meeting 2 (2002) 1105 /1112. [9] G.M. Huang, Y. Li, Power systems reliability indices to measure impacts caused by transient stability crises, IEEE Power Engineering Society Winter Meeting 2 (2002) 766 /771. [10] A.A. Fouad, V. Vittal, S. Rajagopal, V.F. Carvalho, M.A. ElKady, J.V. Mitsche, M.V. Pereira, Direct transient stability analysis using energy functions: application to large power networks, IEEE Transactions on Power Systems PWRS-2 (1) (1987) 37 /44. [11] N. Kakimoto, Y. Ohnogi, H. Matsuda, W. Shibuya, Transient stability analysis of large scale power systems by Lyapunov’s direct method, IEEE Transaction on PAS PAS-103 (1) (1984) 160 /167.

F. Allella et al. / Electric Power Systems Research 67 (2003) 21 /33 [12] N. Kakimoto, Y. Ohsawa, M. Hayashi, Transient stability analysis of electric power systems via Lure-type Lyapunov function, Part I and II, Transaction on IEE of Japan 98(5/6) (1978) pp. 63 /79. [13] P.P. Varaiya, F.F. Wu, Direct methods for transient stability analysis of power systems: recent results, Proceedings of the IEEE 73 (12) (1985) 1703 /1715. [14] H.D. Chiang, F.F. Wu, P.P. Varaiya, Foundations of the potential energy boundary surface method for power system transient stability analysis, IEEE Transactions on Circuits and Systems 35 (6) (1988) 712 /728. [15] T.S. Chung, F. Da-Zhong, A fast approach to transient stability estimation using an improved potential energy boundary surface method, Electric Power Systems Research 34 (1) (1995) 47 /55. [16] E. Chiodo, F. Gagliardi, D. Lauria, Probabilistic approach to transient stability evaluation, IEE Proceedings-Generation, Transmission and Distribution 141 (5) (1994) 537 /544. [17] E. Chiodo, D. Lauria, Transient stability evaluation of multimachine power systems: a probabilistic approach based upon the extended equal area criterion, IEE Proceedings-Generation, Transmission and Distribution 141 (6) (1994) 545 /553. [18] A. Papoulis, Probability, Random Variables, Stochastic Processes, third ed, McGraw Hill, New York, 1991. [19] E. Chiodo, F. Gagliardi, M. La Scala, D. Lauria, Probabilistic on-line transient stability analysis, Proceedings of the IEE-C 146 (2) (1999) 176 /180. [20] R.A. Johnson, Stress /strength models for reliability, in: P.R. Krishnaiah, C.R. Rao (Eds.), Handbook of Statistics, vol. 7, North-Holland, Amsterdam, 1988, pp. 27 /53. [21] I.A. Ushakov, R.A. Harrison, Handbook of Reliability Engineering, Wiley, New York, 1994. [22] D.R. Cox, D. Oakes, Analysis of Survival Data, Chapman & Hall, London, 1984. [23] E.L. Crow, K. Shimizu, Lognormal Distributions, Marcel Dekker, New York, 1988. [24] S. Weerahandi, R.A. Johnson, Testing reliability in a stress / strength model when X and Y are normally distributed, Technometrics 34 (1) (1988) 83 /91. [25] V.K. Rohatgi, Statistical Inference, Wiley, New York, 1984. [26] H.F. Martz, R.A. Waller, Bayesian Reliability Analysis, Wiley, New York, 1982. [27] D.R. Cox, D.V. Hinkley, Theoretical Statistics, Chapman & Hall, 1974. [28] S.B. Nandi, A.B. Aich, Hypothesis-test for reliability in a stress / strength model with prior information, Microelectronics and Reliability 37 (4) (1997) 694 /695. [29] H.F. Martz, L.L. Hamm, W.H. Reed, P.Y. Pan, Combining mechanistic best-estimate analysis and Level 1 probabilistic risk assessment, Reliability Engineering and System Safety 39 (1993) 89 /108. [30] R.R. Sitter, C. Wu, A note on Woodruff confidence intervals for quantiles, Statistics and Probability Letters 52 (4) (2001) 353 /358. [31] J.D. Mc Calley, V. Vittal, N. Abi-Sambra, An overview of riskbased security assessment, Proceedings of the 1999 IEEE/PES Summer Meeting (1999) 173 /178. [32] V. Vittal, J.D. Mc Calley, V. Van Acker, F. Fu, N. Abi-Samra, Transient instability risk assessment, Proceedings of the 1999 IEEE/PES Summer Meeting (1999) 206 /211. [33] V. Miranda, J.N. Fidalgo, J.A. Pecas Lopes, L.B. Almeida, Real time preventive actions for transient stability enhancement with a hybrid neural network-optimization approach, IEEE Transactions on Power Systems 10 (2) (1995) 1029 /1035.

33

[34] M.B. Djukanovic, D.J. Sobajic, Y.H. Pao, An approach to the direct analysis of transient stability, European Transactions on Electrical Power (ETEP) 3 (3) (1993) 241 /247. [35] F. Shuti, H. Jiaxiao, N. Yixin, C. Jianlin, X. Youfang, Z. Zunlin, A practical on-line transient stability program, IEE Second International Conference on Advances in Power System Control, Operation and Management, Hong-Kong (1993) 884 /886. [36] P. Omahen, Fast transient stability assessment using corrective PEBS method, Proceedings of Sixth Mediterranean Electrotechnical Conference (IEEE sponsored), Ljubljana (Slovenia) (1991) 1408 /1411. [37] L.F.C. Alberto, F.H.J.R. Silva, N.G. Bretas, Direct methods for transient stability analysis in power systems: state of art and future perspectives, Proceedings of IEEE Porto Power Tech Conference, 10 /13th September 2001, Porto (Portugal). [38] L. Wehenkel, Machine learning approaches to power system security assessment, IEEE Expert-Intelligent System and Their Application 12(5) (1997) pp. 60 /72.

Biographies Dr Flavio Allella received the Electrical Engineering degree with honours from the University of Napoli (Italy) in 1998. He is a Ph.D. student in Electrical Engineering at the University of Napoli. His research areas include electrical power systems stability and reliability. Address: Dipartimento di Ingegneria Elettrica, Universita` di Napoli, via Claudio 21, 80125 Naples, Italy. Tel.: /39-08-1768-3227; fax: /39-081239-6897; E-mail: [email protected]. Professor Elio Chiodo, Ph.D., received the Electronics Engineering degree with honours from the University of Napoli (Italy), in 1985, and the Ph.D. degree in Statistics at the Department of Mathematical Statistics of the same University in 1992. He is Associate Professor at the Department of Electrical Engineering of the University of Napoli, Italy. His researches are devoted to engineering applications of probability theory and reliability. Address: Dipartimento di Ingegneria Elettrica, Universita` di Napoli, via Claudio 21, 80125 Naples, Italy. Tel.: /39-08-1768-3226; fax: /39-081239-6897; E-mail: [email protected]. Professor Davide Lauria received the Electrical Engineering degree with honours from the University of Napoli (Italy) in 1987. In 1998 he joined the Electrical Engineering Department of the University of Napoli where he is Associate Professor. His research areas include electrical power systems analysis, reliability and power electronics application in power systems. Address: Dipartimento di Ingegneria Elettrica, Universita` di Napoli, via Claudio 21, 80125 Naples, Italy. Tel.: / 39-08-1768-3212; fax: /39-08-1239-6897; E-mail: [email protected].