IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 20, NO.

2, MAY 2005 773

Cumulant-Based Probabilistic Optimal Power Flow

(P-OPF) With Gaussian and Gamma Distributions
Antony Schellenberg, William Rosehart, and José Aguado

Abstract—This paper introduces the cumulant method for the of the method in [5] to the P-OPF problem using a logarithmic
probabilistic optimal power flow (P-OPF) problem. By noting that barrier interior point method (LBIPM) [7]-type solution.
the inverse of the Hessian used in the logarithmic barrier interior
This paper is structured in the following manner. Section II
point can be used as a linear mapping, cumulants can be computed
for unknown system variables. presents information related to the Edgeworth form of the
Results using the proposed cumulant method are compared Gram–Charlier A series. It includes the A series itself, in addi-
against results from Monte Carlo simulations (MCSs) based on a tion to some background information on Tchebycheff–Hermite
small test system. The Numerical Results section is broken into polynomials and computation of A series coefficients. Next, in
two sections: The first uses Gaussian distributions to model system
Section III, an overview of the pure Newton step in the LBIPM
loading levels, and cumulant method results are compared against
four MCSs. Three of the MCSs use 1500 samples, while the fourth for numerical programming is provided. In Sections IV and V,
uses 20 000 samples. The second section models the loads with a the cumulant method is presented, in addition to the proposed
Gamma distribution. Results from the proposed technique are application to the P-OPF problem. Numerical results from a
compared against a 1000-point MCS. system based on the Matpower 9-bus system [8] using normally
The cumulant method agrees very closely with the MCS results
(Gaussian) and Gamma distributed independent random loads
when the mean value for variables is considered. In addition, the
proposed method has significantly reduced computational expense with the proposed cumulant method for P-OPF are detailed in
while maintaining accuracy. Section VI. Finally, conclusions are presented in Section VII.
Index Terms—Cumulants, optimal power flow (OPF), proba- Two appendixes are included to provide background infor-
bilistic optimization. mation in probability and statistics, focusing on moments and
cumulants, as well as information on Gaussian and Gamma dis-
tributions.
I. INTRODUCTION

O PTIMAL POWER FLOW (OPF) is a tool that has been

commonly used within the power systems industry for
many years [1] and has generally been addressed as a deter- II. GRAM–CHARLIER A SERIES
ministic optimization problem. However, it is becoming increas-
ingly important that solution methods to the optimal power flow The Gram–Charlier A Series allows many PDFs, including
problem be developed to address probabilistic quantities and, Gaussian and Gamma distributions, to be expressed as a series
thus, transform the optimal power flow problem into the proba- composed of a standard normal distribution and its derivatives.
bilistic optimal power flow (P-OPF) problem [2]. As a part of the proposed P-OPF method, distributions are re-
Probabilistic programming, or probabilistic optimization, is constructed with the use of the Gram–Charlier A Series. Addi-
concerned with the introduction of probabilistic randomness tional information can be found in [9]. The series can be stated
or uncertainty into conventional linear and nonlinear programs as follows:
[3]. However, the randomness introduced tends to have some
structure to it, and this structure is generally represented with
(1)
a probability density function (PDF) [4]. The goal of the
P-OPF problem is to determine the PDFs for all variables in
the problem. These PDFs are the distributions of the optimal where is the PDF for the random variable is the th
solutions. A typical example of an uncertain or probabilistic series coefficient, is the th Tchebycheff–Hermite, or
parameter in a P-OPF problem is bus loading. Hermite, polynomial, and is the standard normal distribu-
The cumulant method for probabilistic power flow was dis- tion function (see Appendix I-A).
cussed in [5] and [6]. The present paper briefly outlines the fun-
The Gram–Charlier form uses moments to compute series co-
damentals of the cumulant method and presents the adaptation
efficients, while the Edgeworth form uses cumulants. Since the
work presented in this paper is based on cumulants, only Edge-
Manuscript received June 10, 2004; revised November 17, 2004. Paper no. worth’s form of the A series is discussed.
TPWRS-00300-2004.
A. Schellenberg and W. Rosehart are with the University of Calgary, Calgary, Throughout this section, the operator is defined as the
AB T2N 1N4, Canada (e-mail: schellen@[Link]; rosehart@[Link]- derivative with respect to to simplify notation.
[Link]). The remainder of this section is devoted to discussion of the
J. Aguado is with the University of Malaga, 29071 Malaga, Spain (e-mail:
jaguado@[Link]). Hermite polynomials and the computations of A series coeffi-
Digital Object Identifier 10.1109/TPWRS.2005.846184 cients in (1).
0885-8950/$20.00 © 2005 IEEE
774 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 20, NO. 2, MAY 2005

A. Tchebycheff–Hermite Polynomials of the PDF is broken into its series representation and equated
There are two different forms for Hermite polynomials [10]. with the Gram–Charlier A series (1).
The first is based on and the second on . The The PDF, as an exponential, is written in the following form
second form is the same as the PDF for a standard normal dis- using cumulants [9]:
tribution and is more convenient for this application. To avoid
confusion, the notation is used in this paper to denote the (14)
use of the second type of Hermite polynomial.
Since the PDF for a normal distribution is an exponential where is the th derivative of the unit normal distribution,
term, taking derivatives successively returns the original func- is the th cumulant, and is the standard unit normal
tion with a polynomial coefficient multiplier. These coefficients PDF. A complete derivation of (14) can be found in [9].
are referred to as Tchebycheff–Hermite, or Hermite, polyno- Expanding (14) as an exponential series yields
mials.
To illustrate how the Hermite polynomials are generated, the
first four derivatives of the standard unit normal distribution are
taken as follows:

(2)
(3)
(4)
(5) (15)

(6) If each of the terms are expanded individually and grouped

based on powers of , the following result is obtained:
where is the th derivative.
The Tchebycheff–Hermite polynomials are the polynomial
coefficients in the derivatives. Using the results of the first four
derivatives in (2)–(6), the first five Tchebycheff–Hermite poly-
nomials are written as follows:

(7) (16)
(8)
(9) Using the relationship for Hermite polynomials in (12) to re-
(10) place the powers of in (16) gives
(11)

Because of the structure of (2)–(6), the highest power coef-

ficient of the odd derivatives, i.e., the third, fifth, seventh, etc.,
are negative. Equations (7)–(11) have been formed following
the convention that the equations relating to the odd derivatives (17)
are multiplied by negative one, such that the coefficient of the
highest power is positive [9]. Returning to the definition for the Gram–Charlier A series in
Therefore, the th Tchebycheff–Hermite polynomial can be (1) and expanding the summation yields
symbolically written as
(18)
(12)
Comparing (17) and (18), the values for the coefficients can
In addition, a recursive relationship is available to determine
be determined. Based on the equations presented, the first seven
third-order and higher polynomials
terms of the Edgeworth form of the A series are presented in
(13) Table I.

III. INTRODUCTION TO THE LBIPM

B. Edgeworth A-Series Coefficients The LBIPM has been used in a wide variety of applications.
Given the cumulants for a distribution in standard form, i.e., It is a method for the solution of constrained optimization prob-
zero mean and unit variance, the coefficients for the Edgeworth lems and makes use of a variety of techniques from many con-
form of the A series can be computed. In order to find the equa- tributors. Detailed information for the LBIPM can be found in
tions for the A series coefficients, an exponential representation [7], [11], and [12].
SCHELLENBERG et al.: CUMULANT-BASED PROBABILISTIC OPTIMAL POWER FLOW 775

TABLE I the moment generating function for the random variable

A SERIES COEFFICIENT EQUATIONS can be written as follows:

(22a)
(22b)

(22c)

Since are independent

(23a)
(23b)

The cumulants for the variable can be computed using the

cumulant generating function, as defined in Appendix II, in
The algorithm proposed in this paper uses some matrices and
terms of the component variables as follows:
components computed as part of the Newton step in the LBIPM.
Therefore, a very brief outline of the Newton step and its com-
putation in an LBIPM framework is included. (24a)
(24b)
A. Gradient, Hessian, and Newton Direction
Any optimal solution must satisfy the condition that the gra- (24c)
dient of the Lagrangian is zero. This condition is known as (24d)
the Karush–Kuhn–Tucker (KKT) first-order condition for op-
timality. To compute the second-order cumulant, the zero-, first-, and
A Newton–Raphson-based solver can be used to solve second-order derivatives of the cumulant generating function for
this system of nonlinear equations. In this case, the the random variable are computed as
Newton–Raphson method attempts to solve the system using
the following relationship: (25)

(19)
(26)
where is a vector of primal and dual variables in the optimiza-
tion problem, and and are the Hessian and the gra- (27)
dient of the Lagrangian, respectively, evaluated for the variable
values at the current iteration. The update step is computed Evaluating (27) at gives
from the Newton–Raphson method and is known as the Newton
step or, alternatively, the Newton direction if magnitude is nor- (28)
malized. System variables are updated via the following equa-
tion: Third- and higher order cumulants can be computed following
the same procedure. In general, the th-order cumulant for ,
(20) a linear combination of independent random variables, can be
determined with the following equation:
Step lengths are chosen to ensure that the resulting point remains
within the feasible space. (29a)
(29b)
IV. CUMULANT METHOD
The cumulant method relies on the behavior of random vari- where the exponent denotes the th derivative with respect
ables and their associated cumulants when they are combined in to .
a linear fashion. This section discusses the formation of random
variables from a linear combination of others and the role cu- V. ADAPTATION OF THE CUMULANT METHOD
mulants play in this combination. TO P-OPF PROBLEM
Given a new random variable , which is the linear combina- The cumulant method is adapted from the basic derivation
tion of independent random variables, above to accommodate the P-OPF problem when an LBIPM-
type solution is used. The Hessian of the Lagrangian is neces-
(21) sary for the computation of the Newton step in the LBIPM. The
776 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 20, NO. 2, MAY 2005

inverse of the Hessian, however, can be used as the coefficients random variables, and PDFs are reconstructed using the
for the linear combination of random bus loading variables. Gram–Charlier/Edgeworth Expansion theory [5].

A. Inverse Hessian as Linear MAP D. Computation and Use of the Statistical Step
The pure Newton step is computed at iteration of the The proposed algorithm makes use of statistical information
LBIPM using the following equation: during every iteration of the LBIPM.
The statistical step is computed at each iteration in the fol-
(30) lowing fashion.
where is the vector of variables, and are Hessian 1) Distributions for system variables are reconstructed
and gradient of the Lagrangian, respectively, evaluated at , based on computed cumulants and the Edgeworth
and is the pure Newton step. Replacing with form of the A series.
in (30) and rearranging gives 2) The difference between the variable values for the cur-
rent iteration and the peak values of the distributions is
(31) computed.
A general linear equation can be written as 3) A step, known as the statistical step, is set equal to this
difference.
(32) The statistical step is combined with the pure Newton step in
a linear fashion to produce a step that is applied to system vari-
where is the slope, is the variable, and is the -intercept. ables. The emphasis on the Newton step increases as the LBIPM
Noting the similarities in the form of (32) and (31), the matrix progresses toward a solution to ensure good convergence be-
, the inverse Hessian, in (31) contains the multipliers havior. The original variable update (20) can be rewritten in the
for a linear combination of PDFs for random bus loads. Alterna- following fashion to include the linear combination of the pure
tively stated, the negative inverse Hessian is a linear map from Newton and statistical steps:
one variable to another.
(35)
B. Including Random Loads
where and are the scalar weighting for the linear combina-
It is necessary to introduce the cumulants related to the
tion used to vary the emphasis between the pure Newton step
random loads into the system in such a way that the cumu-
and the statistical step.
lants for all other system variables can be computed. Some
It is particularly noteworthy that when Gaussian distributions
characteristics of the gradient of the Lagrangian are used to
are used, there is no obvious statistical step available since the
accomplish this.
mean of the distribution corresponds to the peak in all cases.
When the gradient of the Lagrangian is taken, the power flow
Therefore, the statistical step in this case is always zero.
equations appear unmodified in this vector. Therefore, cumulant
In the case of non-Gaussian distributions, including Gamma
models in the bus loads map directly into the gradient of the La-
distributions, the peak of the PDF does not generally correspond
grangian. For the purposes of mapping, the mismatch vector,
to the mean. Consequently, a statistical step is available for use
in (31), is replaced by a new vector containing the cu-
in the procedure described.
mulants of the random loads in the rows corresponding to their
The statistical step is introduced to place greater emphasis
associated power flow equations.
on optimizing around parameter settings that are more likely to
C. Generalized Results occur. In this paper, the weighting factors for the statistical step
have been arbitrarily chosen, and a range of different weightings
The linear mapping information contained in the inverse Hes- have been tested. However, convergence problems developed
sian can be used to determine cumulants for other variables when the statistical step was heavily weighted compared to the
when bus loading is treated as a random variable. If Newton step. It was found, in the simulations, that when the ini-
is written in the following form tial weighting of the statistical step was greater than ,
where , the solution tended to fail. This is expected
since, when using a primal-dual interior point approach, the
.. .. .. .. .. (33) Newton step is required to converge to the optimal solution.
. . . . . For all systems tested in this paper, the program converged to
the same point regardless of values of and , subject to them
then the th cumulant for the th variable in is computed using being below the threshold discussed above. Therefore, final re-
the following equation: sults using the statistical step were identical to the results if the
distribution mapping was applied only once after the optimiza-
(34) tion was completed.
where is the th element in , and is the th cumulant
for the th component variable. VI. NUMERICAL RESULTS
In the proposed cumulant method for P-OPF, the cumulants The proposed cumulant method was tested using random
for unknown random variables are computed from known bus loads with the mean value set at the nominal bus loading
SCHELLENBERG et al.: CUMULANT-BASED PROBABILISTIC OPTIMAL POWER FLOW 777

TABLE II
MEAN VALUE COMPARISON TABLE

level. Problems based on the Matpower nine- and 118-bus 1) Mean Values: Table II contains all results related to the
systems [8] are used to show general trends and characteristics mean values for system variables. Values for the system vari-
under random loading conditions using the proposed cumulant ables are included, in per unit, as well as a comparison between
method. Of particular interest are the optimal distributions for the cumulant method and each of the four MCSs presented as
the decision variables. an absolute percent difference. Columns labeled “value” are the
Two different and independent sets of results are included. actual value of the variable in p.u., while columns titled “differ-
In the first set, loads are modeled in the nine-bus and 118-bus ence” are the absolute percent difference between the cumulant
systems using Gaussian distributions with variances such that method and the MCS results.
the 99% confidence interval is equal to % of the nominal The results for the mean value of the distributions using the
loading value. The second set models loads using Gamma dis- cumulant method are, in general, well within 1% of the values
tributions such that the variance is 15% of the nominal loading found using MCS. With the exception of the reactive power
value, and only results for the nine-bus system are presented. For generation at bus 2, the maximum percent difference between
all problems, the problem converges when the barrier parameter the mean from any of the MCSs and the cumulant method is
is less than as computed using the complementary gap. 0.4437% and occurred for the angle at bus 8 in the third 1500
sample MCS. The results for reactive power generation at bus
A. Gaussian Distributions 2, however, had a very high percent difference between MCS
The results for the Gaussian distributions are divided into and cumulant method. It is noteworthy that the minimum abso-
three sections. The first section presents the results for the mean lute percent difference of 62.3% for reactive power generation
value of the distributions and includes discussion about these re- at this bus resulted from a difference of 0.0005 p.u. Similarly,
sults. The second section presents and discusses the results for the maximum difference of 79.17% was from a difference of
the variance of the distributions. The first and second sections only 0.0012 p.u. Although the percent difference is for reactive
use the nine-bus problem, while the third section presents results power generation at this bus is high, the actual error in per unit
using the 118-bus system to illustrate the cumulant method’s is small. Therefore, the percent error is a somewhat misleading
performance as the system size increases. measure for these situations.
In all cases, the mean value for bus loading was taken at the In general, the cumulant method approximation matches very
nominal loading value from the Matpower problems, and the well with the MCS results with respect to the mean values.
variance is such that the 99% confidence interval is 10% of the 2) Variance Values: The covariances are calculated numer-
nominal value. ically from the results in the MCSs. The results from the MCSs
For the nine-bus system, a total of four Monte Carlo simula- are presented in Table III along with the results using the pro-
tions (MCSs) are included. Three were run using 1500 samples posed cumulant method to allow for a direct comparison. Again,
and one with 20 000 samples. The raw results are included in columns labeled “value” are the actual value of the variable in
both sections in addition to the comparison between the cumu- p.u., while columns titled “difference” are the absolute percent
lant method results and the MCS results. difference between the cumulant method and the MCS results.
778 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 20, NO. 2, MAY 2005

TABLE III
VARIANCE VALUE COMPARISON TABLE

Of particular interest is the fact that the angle at bus 1 has TABLE IV
118-BUS SYSTEM GAUSSIAN DISTRIBUTION RESULTS SUMMARY:
zero variance for all cases. This phenomenon results from the MPE—MEAN PERCENT ERROR, MAPE—MEAN ABSOLUTE PERCENT ERROR
fact that bus 1 is used as the angle reference bus and is fixed at
precisely zero.
In most cases, the percent difference between the MCSs and
the cumulant method results is less than 6%. The notable excep-
tion to this statement is the percent difference in variances for
voltage variables, which is substantially higher, in general, than
other system variables. In particular, the variance in bus volt-
ages at buses 1, 4, 6, and 8 is much higher than 6%. Although
the percent difference is much higher for voltage, the difference general, the difference in the absolute mean values for system
in the actual variance value is not. The worst absolute percent variables are less than 1.5% compared to the MCSs. The differ-
difference for all variance results occurred for the voltage mag- ence in the variances is between 2% and 8.01%. Fig. 1 shows
nitude at bus 6. However, the actual value of the variance for this the PDF for the objective function in the 118-bus problem. In
bus voltage is extremely small, i.e., less than . Therefore, this problem, the objective is a linear active power generation
this variable is almost deterministic, compared to others in the cost function.
problem, and can be treated as such without any significant loss
in statistical information. B. Gamma Distributions
3) 118-Bus System: For the 118-bus system, the results are
aggregated since the number of buses and variables is too high Results included in this section are based on Gamma-dis-
to present them individually. The results presented in Table IV tributed random loads with the mean at the nominal load and
provide an illustration of how the proposed cumulant performs the variance 15% of the nominal value. MCSs are performed
when applied to larger systems. Results have been tabulated in with 1000 samples, and these results are taken as the reference
terms of mean and variance values since the number of variables solution.
in the system exceeds 300. The cumulant method is compared Reconstructions using the Gram–Charlier A series can be per-
against a MCS consisting of 1500 samples. formed with any number of cumulants. Results are presented
The column in Table IV labeled MPE denotes the mean per- here for solutions using up to fifth- and ninth-order cumulants.
cent error, that is, the average error with the sign considered. As the cumulant order increases, the computational expense for
In contrast, the column labeled MAPE (mean absolute percent the reconstruction also increases.
error) takes the absolute value of the individual percent errors One of the benefits of the proposed algorithm is a substantial
prior to computing the average. As in the nine-bus system and reduction in computational expense while maintaining a high
discussed in Section VI-A2, several variables are effectively de- level of accuracy. Table V demonstrates the difference in com-
terministic and are, therefore, not analyzed as probabilistic. In putation time for several different algorithms. Based on the re-
SCHELLENBERG et al.: CUMULANT-BASED PROBABILISTIC OPTIMAL POWER FLOW 779

Fig. 2. Bus 2 Active Power Generation—Fifth-Order Reconstruction.

Fig. 1. PDF for 118-bus system. Histogram—Monte Carlo, Line—Cumulant Method.

TABLE V
COMPUTATIONAL EXPENSE WITH GAMMA DISTRIBUTIONS

sults presented, the proposed algorithm is several times faster

than an MCS.
MCS run times are highly dependent on the number of points
used. In this case, only 1000 samples are taken, and in general,
a higher number of points results in better solutions. Compar-
atively, the proposed methodology depends only on the cumu-
lant order used. Generally, a higher number of cumulants should
result in better approximations. However, additional complica-
tions can occur as order increases including negative probabili- Fig. 3. Bus 5 Voltage Magnitude—Fifth-Order Reconstruction.
ties in the tail sections and multimodal distributions [9]. Histogram—Monte Carlo, Line—Cumulant Method.
Figs. 2 and 3 show the reconstructed distribution compared
with the Monte Carlo results. In these figures, the solid black that more simulations will be required to address uncertainty in
line is the result using the proposed cumulant method, and the the sampling process.
histogram is the raw MCS results. Non-Gaussian behavior is ev-
ident in both of the figures, and alternative methodologies that
only consider Gaussian distributions would be unable to accu- VII. CONCLUSION
rately model the lopsided shape of these distributions. This paper introduces the cumulant method for the P-OPF
Of note is the fact that the approximation in Fig. 3 extends far problem. By noting that the inverse of the Hessian used in the
to the right of the MCS results. This occurs because the proposed logarithmic barrier interior point can be used to perform linear
cumulant method returns a continuous result that does not con- mapping, cumulants can be computed for unknown system vari-
sider operational limitations in the distributions, although oper- ables.
ational limitations are enforced during the optimization process. Results using the cumulant method are compared against the
In this case, the PDF should be truncated on the right-hand side results from MCS for nine- and 118-bus systems. The proposed
to limit the maximum value for the voltage magnitude. algorithm has a substantial reduction in computational expense
The assumptions made in the proposed system do not change while maintaining a high level of accuracy.
based on system size. Therefore, larger systems can be analyzed
in the same fashion. It is expected that accuracy would remain
APPENDIX I
high, and there will be increasing computational performance
RANDOM VARIABLE DISTRIBUTIONS
gain of the proposed algorithm compared with Monte Carlo. The
computational performance will increase because the individual The material presented in this paper makes use of two dif-
optimizations in the MCS will take longer, in addition to the fact ferent types of PDFs: Gaussian and Gamma. This appendix
780 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 20, NO. 2, MAY 2005

briefly introduces these distributions and some important prop- and then uses these results to develop the necessary cumulant
erties of each. Further information is available in [4] and [13]. relationships.
The expected value of a random variable is defined as
A. Gaussian Distributions
The Gaussian Distribution, also known as the Normal Distri- (40)
bution [13], is commonly used in a variety of different areas. It
is a simple distribution and is characterized by its mean and where is the PDF of .
variance , according to the following relationship: The th-order raw moment is defined in the following
manner:
(36)
(41)
where is the random variable.
The standard unit normal distribution is defined as the It is possible to compute the raw moments through the use of the
Gaussian distribution with zero mean and unit variance. In a moment generating function . Mathematically, this func-
Gaussian distribution, the peak of the PDF always occurs at the tion is stated as [4]
mean value.
(42)
B. Gamma Distributions The th raw moment is computed from the moment generating
Another frequently used distribution is known as the Gamma function by taking the th derivative with respect to and eval-
distribution. This distribution is characterized by three vari- uating at . For example, the third raw moment can be
ables: the random variable in addition to two non-negative computed as follows:
shape parameters.
The general formula for the Gamma distribution is as follows (43a)
[4]:
(43b)
(37)
otherwise (43c)
where is the random variable, and and are the shape pa- (43d)
rameters. The notation stands for the complete Gamma
(43e)
function, which can be written in the following manner [4]:
The cumulant generating function, denoted by , is often
(38) written in terms of the moment generating function , as
follows [4]:
In the case that is an integer, the complete Gamma function
can be written in the following simplified form [4]: (44)

(39) The cumulant generating function is employed in the same

manner as the moment generating function; successive deriva-
The Gamma distribution has several important properties. tives are taken with respect to and evaluated at . The
First, the shape is controlled by two independent non-negative th cumulant is denoted as .
parameters. The values for these two parameters strongly affect
the shape of the resulting PDF. Second, the PDF can only be B. Illustrative Example
nonzero for positive values of the random variable. In other Consider the standard normal distribution based on (36) with
words, the PDF is bounded on one side. Third, the mean value zero mean and unit variance. The moment generating function
and the peak value of the PDF are generally different. in (42) can be developed using the definition of expected value
from (40) as follows:
APPENDIX II
PROBABILITY AND STATISTICS BACKGROUND (45a)

This appendix includes an overview of information related to (45b)

moments and cumulants for PDFs. More detailed information is
available in [4]. (45c)

A. Moments and Cumulants

(45d)
Moments and cumulants are both measures of a PDF. Since
equations related to cumulants are often developed in terms of
(45e)
moments, this section begins with an introduction to moments
SCHELLENBERG et al.: CUMULANT-BASED PROBABILISTIC OPTIMAL POWER FLOW 781

[4] A. Papoulis and S. Pillai, Probability, Random Variables, and Stochastic

(45f) Processes, 4th ed. New York: McGraw-Hill, 2002.
[5] P. Zhang and S. T. Lee, “Probabilistic load flow computation using the
method of combined cumulants and Gram–Charlier expansion,” IEEE
(45g) Trans. Power Syst., vol. 19, no. 1, pp. 676–682, Feb. 2004.
[6] W. Tian, D. Sutanto, Y. Lee, and H. Outhred, “Cumulant based prob-
(45h) abilistic power system simulation using Laguerre polynomials,” IEEE
Trans. Energy Convers., vol. 4, no. 4, pp. 567–574, Dec. 1989.
[7] G. Torres and V. Quintana, “An interior-point method for nonlinear op-
Taking the result in (45h) and computing the cumulant gener- timal power flow using voltage rectangular coordinates,” IEEE Trans.
ating function based on (44) yields Power Syst., vol. 13, no. 4, pp. 1211–1218, Nov. 1998.
[8] R. Zimmerman and D. Gan. (1997) Matpower. Power Systems Engi-
neering Research Center (PSERC). [Online]. Available: [Link]
(46) [Link]/matpower/
[9] M. G. Kendall and A. Stuart, The Advanced Thoery of Statistics, 4th
ed. New York: Macmillan, 1977.
To compute the th cumulant, successive derivatives of the [10] J. B. Calvert, “Hermite polynomials,”,
cumulant generating function are taken and evaluated at . [Link] 2004.
Suppose the third cumulant for a standard normal distribution is [11] G. Astfalk, I. Lustig, R. Marsten, and D. Shanno, “The interior-point
method for linear programming,” IEEE Softw., vol. 9, no. 4, pp. 61–68,
of interest and needs to be computed. The first three derivatives Jul. 1992.
of (46) with respect to are taken as follows: [12] J. D. Weber, “Implementation of a Newton-based optimal power flow
into a power system simulation environment,” Univ. Illinois at Urbana-
(47) Champaign, [Online]. Available: [Link]
ment%20Library/Newton-based%20Optimal%20Power%[Link],
(48) 1997.
[13] E. W. Weisstein. Eric Wesstein’s World of Mathematics [Online]. Avail-
(49) able: [Link]

The third cumulant can be found by evaluating (49) at

(50)
Antony Schellenberg received the [Link]. degree in electrical engineering from
For this example, the third cumulant is zero. The same process the University of Calgary, Calgary, AB, Canada, in 2002. He is currently a grad-
uate student pursuing the Ph.D. degree at the University of Calgary.
can be repeated for any cumulant of interest.
All cumulants of order three and higher are zero in this ex-
ample since the derivative of zero is zero. In fact, this result is
true for any general Gaussian distribution, and all higher order
cumulants, third or greater, are zero. Consequently, cumulants William Rosehart received the Ph.D. degree in electrical engineering from the
can, in some sense, be considered as a measure of the departure University of Waterloo, Waterloo, ON, Canada, in 2001.
He is currently an Assistant Professor in the Department of Electrical and
from normality. Computer Engineering, University of Calgary, Calgary, AB, Canada.

REFERENCES
[1] M. Huneault and F. Galiana, “A survey of the optimal power flow liter-
ature,” IEEE Trans. Power Syst., vol. 6, no. 2, pp. 762–770, May 1991.
[2] G. Viviani and G. Heydt, “Stocahstic optimal energy dispatch,” IEEE José Aguado received the Ph.D. degree from the University of Malaga, Malaga,
Trans. Power App. Syst., vol. PAS-100, no. 7, pp. 3221–3228, Jul. 1981. Spain, in 2001.
[3] J. R. Birge and F. Louveaux, Introduction to Stochastic Program- Currently, he is an Associate Professor in the Department of Electrical Engi-
ming. New York: Springer-Verlag, 1997. neering at the University of Malaga.