IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 19, NO. 3, AUGUST 2004 1483
Evaluation of Uncertainty in Dynamic
Simulations of Power System Models:
The Probabilistic Collocation Method
James R. Hockenberry, Member, IEEE, and Bernard C. Lesieutre, Member, IEEE
Abstract—This paper explores the use of a new technique, the
probabilistic collocation method (PCM), to enable the evaluation
of uncertainty in power system simulations. The PCM allows the
uncertainty in transient behavior of power systems to be studied
using only a handful of simulations. The relevant theory is outlined
here and simple examples are used to illustrate the application of
PCM in a power systems setting. In addition, an index for identi-
fication of key uncertain parameters, as well as an example with a
more realistic power system, are presented.
Index Terms—Power system modeling, power system simulation,
uncertainty.
I. INTRODUCTION
T
IME-STEP simulation techniques form an important class
of tools for power system analysis. They are employed
whenever dynamic phenomena are to be studied, from the stand-
alone analysis of individual components or the study of fast elec-
tromagnetic transients, to the behavior of large-scale systems
over many time scales. (We use the words “time-step simula-
tion” to refer to methods that emulate the dynamic response of
a system typically represented by ordinary differential and al-
gebraic equations, and to distinguish from the usage of “simu-
lation” to refer to techniques to evaluate probabilistic problems
using Monte Carlo and related methods.) The results of such
studies are used to make decisions concerning the structure,
tuning, and operation of the system. The constant use and impor-
tance of these tools motivates the equally constant research on
improving models and simulation algorithms. The contribution
of this paper is to demonstrate practical evaluation of uncertain-
ties in power system models using presently available simula-
tion tools. This method will enable detailed uncertainty studies
that have been infeasible due to computational limitations. By
directly considering uncertainties in the models, the additional
information gained will lead to less need for conservative oper-
ation of the power grid.
The literature concerning uncertainty analyses in power sys-
tems is vast (see the bibliographies in [1]–[7]). The referenced
Manuscript received January 9, 2004.
J. R. Hockenberry was with the Massachusetts Institute of Technology, Cam-
bridge, MA 02139 USA. He is now with Robert Bosch, GmbH, Stuttgart 70049,
Germany (e-mail: james.hockenberry@de.bosch.com).
B. C. Lesieutre was with the Massachusetts Institute of Technology, Cam-
bridge, MA 02139 USA. He is now with the Lawerence Berkeley National Lab-
oratory, Berkeley, CA 94720 USA (e-mail: BCLesieutre@lbl.gov).
Digital Object Identifier 10.1109/TPWRS.2004.831689
books and papers cover a range of topics, many related to statis-
tical reliability studies. With regard to uncertainty studies ad-
dressing transient stability, there are fewer papers, and these
consider uncertainties related to disturbances and operating con-
ditions [8]–[14]. These do not consider uncertainties in param-
eter values, and the approaches combine deterministic simula-
tion techniques with stochastic analyses. In some cases [11],
[10] the dynamic simulation step is significantly enhanced by
the use of Lyapunov-like energy functions. However, these en-
ergy functions are based on simplified models of a power system
that lack significant detail. While these studies are important in
their own right, they do not address the problem we consider and
do not design the time-step simulations around the uncertainty
characterization.
The specific area of parameter uncertainties in detailed
time-domain simulation studies is nearly vacant. This is
understandable. Historically, a time-step simulation of a large
power system model required lengthy computer runs (
hours).
With faster computers these same simulations now can be
completed in minutes. (However, there is now new interest in
simulating larger and more detailed models which increases
the simulation time.) Still, practical techniques for the analysis
of uncertainty, such as Monte Carlo and its derivatives, require
many sample data obtained from repeated simulations. For
example, suppose that 1000 points are necessary, obtained from
1000 time-step simulations, then hours to days of computer
time will be required. Furthermore, each such study only
represents the analysis of a single event (line outage, generator
loss, etc.). Evaluating many contingencies while accounting
for uncertainties is prohibitively time intensive. Consequently,
such studies are not typically done. The results from nominal
simulations are incorporated into design and operation in a
conservative manner.
The recent literature on this topic is not completely empty.
In addition to our initial limited study [15], Hiskens
et al. have
used their “trajectory sensitivity” approach to approximate the
effect of uncertain parameter values on the outcome of time-step
simulations [16]. The method essentially employs an augmented
model that includes additional variables to represent the sensi-
tivity of specified state variables to select parameters and initial
conditions. The result of the time-step simulation of this larger
model yields the nominal trajectory and the sensitivities of the
trajectory to the aforementioned parameters. These linear sensi-
tivities then can be used to form a linear model with which to ap-
proximate the effect of uncertain parameters in a region about a
0885-8950/04$20.00 © 2004 IEEE
1484 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 19, NO. 3, AUGUST 2004
nominal trajectory without the need for further time-consuming
(time-domain) simulation.
We would like to advocate in this work a technique developed
by researchers who study global climate change. Similar to ours,
their time-domain simulation studies require significant com-
puter time and their models include parameters with uncertain
values. The so-called probabilistic collocation method (PCM)
[17], [18] (later renamed the “deterministic equivalent mod-
eling method) involves a few time-step simulations, the results
of which are used to develop a polynomial model that directly
maps the uncertain parameter to outcomes of interest. This
very
reduced-order model can be evaluated quickly to determine the
effect of the uncertainty. In their studies, the uncertain parameter
may be carbon dioxide emissions and the outcome the average
temperature in 2100. In our studies, the uncertainty may be re-
lated to load models or fault time and the outcome may be max-
imum frequency, angle, and voltage variations over the course
of the simulation (among many other possibilities) that can be
compared against operational and stability limits imposed on
these variables.
This approach is particularly appealing because it allows
the use of nonlinear models and the evaluation of complicated
output functions (maximum deviations). And, importantly,
these studies can be performed using existing commer-
cially-available power system simulation software.
The techniques we discuss are applicable to all studies that
rely on dynamic simulations: from fast time-scale EMTP sim-
ulations to slower time-scale large-scale system studies. In this
paper, we focus on system studies for the purpose of accurately
assessing operational limits while considering conditions during
which disturbances may lead to voltage and frequency viola-
tions. Instead of relying on operational safety margins, which
we hope are sufficiently conservative, a planning/operational
decision can be based on the probability of exceeding a set of
specified limits.
Our work on this topic is summarized in this paper. Here
we introduce and evaluate the PCM technique applied to small
power system models. We also address the additional challenges
associated with the analysis of very large system models con-
taining numerous uncertainties as well as presenting a large
system example. A more detailed treatment of some of the topics
presented here can be found in [19].
II. T
HEORY
PCM essentially creates polynomial models relating the un-
certain parameters of the system to the outputs of interest. The
power of the PCM method lies in its ability to select appropriate
simulation points to create a polynomial model which has the
same moments as a higher order model. Judiciously selected
simulations are carried out initially in order to determine the
coefficients of this polynomial model. Orthogonal polynomials
and Gaussian quadrature integration [20] are used to justify the
selection of suitable simulations, so these methods will be dis-
cussed before the actual development of PCM is presented. One
important property of PCM to note in the development is the
independence of the simulations necessary to determine the co-
efficients of the polynomial model from the output(s) of interest.
The same set of simulations is performed to determine the co-
efficients of the polynomial models for all of the outputs of in-
terest.
A. Orthogonal Polynomials
The following function is an inner product on the space of
polynomials:
(1)
where
is any nonnegative weighting function defined ev-
erywhere in a connected
. In the context of PCM, is a
probability density function (pdf) describing the uncertainty in
a system parameter. This particular inner product is the one used
for Gaussian quadrature integration and the probabilistic collo-
cation method.
Given this inner product, a pair of polynomials is said to be
orthogonal if their inner product is zero. Further, a set of poly-
nomials
is said to be orthonormal if and only if the following
relationship holds for all
in
(2)
Given a weighting function
, we are interested in a par-
ticular orthonormal set,
, of polynomials of increasing order
in which
is a polynomial of order . These polynomials are
unique and form a basis for all polynomials. Efficient recursive
methods exist for obtaining these polynomials [20].
We omit the proofs, but it can be shown that each
has ex-
actly
roots and all of the roots are contained in [20]. These
roots play a pivotal role in the probabilistic collocation method.
B. Gaussian Quadrature Integration
Gaussian quadrature integration is a numerical integration
technique for integrals of the form:
(3)
where
is a polynomial and is a nonnegative weighting
function. If
is a pdf, this integral is the expected value of
. The main result of Gaussian quadrature integration is the
following exact formula for calculating this integral:
(4)
where the
are constants which only depend on the weighting
function
and the are constants in the region of integra-
tion. The formula is exact for all polynomials
of order less
than or equal to
. This result is somewhat surprising. The
polynomial
itself could be determined using samples,
but only
samples of are needed to compute the integral.
The
are computed using the orthonormal polynomial set
, where is used as the weighting function. The constants
are the roots of , which exist and are contained in . Fur-
thermore, the
only depend on , since only depends on
. We show the independence of the from and the
correctness of (4) with these
constructively.
HOCKENBERRY AND LESIEUTRE: EVALUATION OF UNCERTAINTY IN DYNAMIC SIMULATIONS OF POWER SYSTEM MODELS 1485
The polynomials in up to and including order constitute
an orthonormal basis for the space of all polynomials of degree
less than or equal to
. Therefore, a polynomial of order
can be expressed in terms of these orthonormal polynomials
using constant coefficients
and
(5)
If this expression is expanded by multiplying through by
,
the result is a sum of 2n linearly independent polynomials (not
necessarily orthogonal), which proves that such a representation
is always feasible. Since
is a constant, the integral is easily
determined by orthogonality
(6)
Finally, we build a linear set of equations by evaluating
in (5) at the roots of
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
(7)
By inverting this matrix, we derive an expression for
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
(8)
If we define
to be 1 (which is customary) and assume
is a pdf (which integrates to 1), then is the desired result
(9)
The weights
are given by the last row of the matrix in (8).
Since
is determined solely by , both and are inde-
pendent of
.
The independence of the
from the particular polynomial
for which we are calculating the integral is analogous to the
desired independence between the set of parameter values for
which we perform simulations and the particular output variable
of interest. One set of
suffices for all polynomials of order less
than or equal to
and one set of simulations is sufficient
for all outputs of interest.
C. Monte Carlo Method
A brief description of the simple, “brute force” Monte Carlo
method for a single uncertain parameter is helpful as a basis of
comparison before we discuss PCM.
In its simplest form, the Monte Carlo method is a method of
repeated trials. As an example, let
be a pdf describing some
uncertain parameter
(10)
(11)
For some desired output,
, which is a function of that pa-
rameter the expected value of the output is computed using the
following integral:
(12)
A simplistic Monte Carlo method approximates the integral by
generating random numbers using
and then performing a
simple average of the resulting answers
(13)
where the
are generated according to using a pseu-
dorandom number generator. Given sufficiently large
, this
method is guaranteed to converge to the actual expected value as
a direct result of the law of large numbers [21]. Unfortunately,
the required
may be quite large. Since all of the moments
of
(e.g., the variance) are expected values, the law of large
numbers applies to these statistics as well under the Monte Carlo
method.
The size of
is dependent on the variance of the output vari-
able. A larger variance means a larger number of simulations
are required. Variance reduction techniques reduce the size of
the required
for a particular output by modifying the system
and pseudorandom numbers in such a way that the expected
value of the desired output is accurately generated but with a sig-
nificantly reduced variance. However, these methods have the
major drawback that a new set of simulations are then neces-
sary for each output variable and each moment, since only the
expected value of one output is accurately computed.
D. PCM—Single Uncertain Parameter
The probabilistic collocation method is a polynomial mod-
eling technique; the desired output is described as a polynomial
in the uncertain parameter of the system. After this model is
identified, a standard, simplistic Monte Carlo method can be
applied to the polynomial model. The problem with (13) is that
each
is computationally expensive, since it represents a
separate simulation of the power system. If this function
can be modeled reasonably accurately by a polynomial ,
an essentially unlimited number of samples
can be com-
puted because no simulations are involved once the polynomial
has been identified. The only simulations necessary are those to
identify the polynomial. If the same simulations can be used to
identify the polynomials for all outputs of interest, the number
of simulations is very limited in contrast to variance reduction
techniques. Therein lies the power of the probabilistic colloca-
tion method.
Concretely, let
be the uncertain parameter, be the pdf
describing this parameter and
be the output of interest.
PCM creates a polynomial model of the form
(14)
where the
are constants. The model parameters could be iden-
tified by performing many simulations using various
and ap-
plying a least squares algorithm. However, we are trying to per-
1486 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 19, NO. 3, AUGUST 2004
form an absolute minimal number of simulations, since each
simulation is expensive in terms of time.
Gaussian quadrature integration says that if this polynomial is
fit using exactly the right simulations, the expected value of
is
identical to the expected value of
as long as is a polynomial
of order
or less. PCM fits the polynomial model with sim-
ulations of the system at the
indicated by Gaussian quadrature
integration. If a high-order polynomial is a reasonable model
for the relationship between uncertain parameter and output of
interest, PCM yields extremely good results for the expected
value. Not surprisingly, the higher order moments are also well
approximated though a firm theoretical basis is lacking.
Following this line of logic, let us represent
using the
polynomials in
without loss of generality, where is of
order
(15)
The following linear system of equations can be solved to de-
termine the
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
(16)
where the
are chosen to be the roots of [to exploit the
benefits of Guassian quadrature integration explored in (5)–(9)].
If
, the expected value is simply and no further
calculations are necessary. Similar relationships exist for higher
order moments of
. For example
(17)
In general, this approach is heuristic. We cannot usually guar-
antee that the actual relationship between the uncertain param-
eter and the output of interest is exactly a polynomial. In ad-
dition, we may be interested in further statistical information
about the output, not just its moments. As previously noted,
higher order moments are also expected values, so the argu-
ments for the expected value hold just as well for the variance,
although a higher order polynomial model may be required, but
further statistical information, such as the exact probability that
the output lies in some interval, may not be accurate. PCM is
designed to calculate the coefficients of the polynomial model
with the bare minimal number of simulations while also trying
to reproduce the moments of the output with high fidelity by
modeling the polynomial particularly well in the regions that are
more probable. Given that, one could also expect that the PCM
model would perform well when used to compute the proba-
bility of events in the high probability region of the output.
III. S
IMPLE ILLUSTRATIVE EXAMPLES
A. RC Circuit
To illustrate the probabilistic collocation method, a simple
first-order circuit, such as the one shown in Fig. 1, may be
helpful because of its familiarity and its tractability to analytical
Fig. 1. Simple RC circuit with voltage step input.
Fig. 2. Pdf for resistance
R
.
techniques. The response voltage is described by a first-order
differential equation
(18)
Given zero initial conditions and the step input shown, the
output is a familiar exponential response
(19)
For this example, the resistance is treated as uncertain and the
capacitance is treated as if it were known exactly. The voltage
at time
is chosen as the output of interest
(20)
Before any uncertainty analysis can be performed, a descrip-
tion of the parameter uncertainty must be available. We use a
somewhat unusual probability density function to describe the
uncertainty in
to demonstrate that the probabilistic colloca-
tion method can accommodate any desired pdf
(21)
This pdf is shown graphically in Fig. 2. For completeness, we
compute the expected value and variance of
(22)
(23)
Since an algebraic expression for
as a function of is avail-
able, the pdf for
can be derived directly, along with its as-
sociated moments. The expected value, variance and standard
deviation of
are summarized in Table I.
HOCKENBERRY AND LESIEUTRE: EVALUATION OF UNCERTAINTY IN DYNAMIC SIMULATIONS OF POWER SYSTEM MODELS 1487
TABLE I
E
XPECTED
VALUE,VARIANCE, AND
STANDARD
DEVIATION FOR
RC-CIRCUIT
EXAMPLE
At this point, we would like to demonstrate how to apply
PCM to this problem. In normal practice, one would not use
PCM in such a situation, since the relationship between the un-
certain parameter and the output of interest is known analyti-
cally, but having the exact answer for comparison purposes is
helpful pedagogically. The first step is to obtain the first few or-
thogonal polynomials based on
(see [20] for a straight-
forward recursive algorithm for the construction of orthogonal
polynomials)
(24)
(25)
(26)
(27)
The roots of
yield two values of , which we can use
to develop a linear approximation of the relationship between
and . The roots of and the resulting values of are as
follows:
(28)
(29)
These two points are sufficient to uniquely define a line, which
is the linear PCM approximate model
(30)
As noted earlier, the expected value, variance and standard devi-
ation are directly available from the coefficients and the results
are shown in Table I. The accuracy of these results is startling
considering that they are based on a linear approximation and
that approximation was created using two points selected based
on the pdf for
and not based on the relationship between
and .
As a basis of comparison, we also examine another standard
way to obtain a linear approximation of
, namely, a Taylor
series approach. Instead of creating a line using two widely sep-
arated values of
, we select a nominal value for and an-
other value of
close to the nominal value, which results in
a line tangential to the actual relationship between
and at
that point; in other words, the approximation is local. One might
select
for the nominal value, since the pdf for is max-
Fig. 3. Comparison of linear approximations using PCM and Taylor series
approximation.
imal there. Another good choice is , since that is the
expected value of
. The PCM linear model as well as these two
linear approximations are shown in Fig. 3.
As Table I shows, the PCM results compare favorably with
those based on the Taylor series approximations. The results
using
as the nominal point are uniformly poor since this
approximation diverges substantially from the actual
for larger
values of
. The approximation based on a nominal value of
yields a comparable variance and standard deviation
to that using PCM since their slopes are nearly identical. But,
the expected value as computed using PCM is more accurate
since the Taylor series approximation in this example always
lies below the actual graph of
; the PCM model is not limited
to being tangential to
.
One can also create higher order polynomial approximations
using PCM. To create a quadratic model, the roots of
are used to find the coefficients. The roots of and the corre-
sponding values of
are as follows:
(31)
(32)
(33)
The quadratic PCM model is as follows:
(34)
The moments can be computed directly from the coefficients,
just as for the linear model, and the results are in Table I. The
quadratic PCM model is able to reproduce the actual moments
almost exactly.
B. Multiple Machine Power System
A power system with multiple machines has richer behavior
than the previously presented example and we close this sec-
tion with such a system. The power system model comprises a
two-area system adapted from [22] and shown in Fig. 4. We ex-
plore the amount of power which can be transferred between the
