University of Southern California Dissertations and Theses

Uncertainty management for complex systems of systems: applications to the future smart grid

(USC Thesis Other)

Uncertainty management for complex systems of systems: applications to the future smart grid

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content UNCERTAINTY MANAGEMENT FOR COMPLEX SYSTEMS OF SYSTEMS:
APPLICATIONS TO THE FUTURE SMART GRID
by
Hadi Meidani
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CIVIL ENGINEERING)
December 2012
Copyright 2012 Hadi Meidani
Dedication
Dedicated to my unboundedly caring mother, Zohreh, and my insightful father, Mehdi...
and to my lovely wife, Marmar.
ii
Acknowledgements
I would like to express my utmost gratitude and most sincere appreciation to my PhD advisor,
Professor Roger Ghanem, for his unbounded inspiration and motivation towards excellence in re-
search endeavors. Words cannot do justice to the extent of academic supports he has passionately
provided in preparing me to tackle scientific challenges. I feel truly indebted to him for investing
his trust in me and introducing me to a plethora of intriguing technical subjects.
During the course of my PhD study at the University of Southern California, I had the un-
paralleled opportunity to conduct research with other prominent researchers in the field of civil
engineering. I would like to acknowledge the academic and financial support provided for me
by Professor Maria Todorovska during the first years at USC. Through my research with her and
Professor Trifunac, I was introduced to the essential components of first-rate technical research,
a major factor for my success for which I will be forever grateful.
I also acknowledge the members of my dissertation and qualifying exam committee, Prof.
Sami Masri, Prof. Amy Rechenmacher, Prof. Burcin Becerik-Gerber and Prof. Giuseppe Caire,
for their valuable comments and suggestions. I especially thank Prof. Sami Masri for his insight-
ful general guidance during my PhD study.
Studying at USC was fortunate for me, in that I met wonderful friends and colleagues from
which I have learned a great deal. Among them, I would like to particularly acknowledge Prof.
Maarten Arnst, for patiently and passionately introducing me to new ideas and also for sharing of
his computer codes.
iii
Studying abroad necessitates a multitude of sacrifices, the greatest of which is having to live
distantly from my loved ones; especially from my parents, Zohreh Karkhanehchi and Mehdi
Meidani, and my two lovely sisters, Elham and Sarah. I owe my everything to my wonderful
parents, and a verbal gratitude would be in despair considering the great deal of sacrifice on their
side, especially during my PhD years.
Finally, I would like to extend my special thanks to my passionate and loving companion,
Marmar Ghadiri, whose unending emotional supports made this seemingly challenging journey
absolutely enjoyable.
University of Southern California, Los Angeles, CA Hadi Meidani
August 2012
iv
Table of Contents
Dedication ii
Acknowledgements iii
List of Tables vii
List of Figures viii
Abstract xiii
Chapter 1 Introduction 1
1.1 Outline and Summary of Contribution . . . . . . . . . . . . . . . . . . . . . . . 4
Chapter 2 Dynamics of the Smart Grid 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Power Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Dynamics of the Power Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4 Small Signal Stability of the Smart Grid . . . . . . . . . . . . . . . . . . . . . . 11
Chapter 3 Uncertainty Quantification for Markov Chain Models 17
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.3 Construction of Uncertain Markov chains . . . . . . . . . . . . . . . . . . . . . 23
3.3.1 Time-invariant RTM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3.2 Time-variant RTM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 MaxEnt Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5 Numerical Illustration I - A 3-state Markov Chain . . . . . . . . . . . . . . . . . 30
3.6 Numerical Illustration II - Markov Chain Deterioration Model . . . . . . . . . . 39
3.6.1 Data-driven MaxEnt Estimation of Deterioration Markov Chain Models . 40
3.6.2 Bayesian Update of the RTM . . . . . . . . . . . . . . . . . . . . . . . . 41
Chapter 4 Stochastic Markov Chain Model for Power Demand 44
4.1 Introduction and Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2 Electric Load Models in Power System Equations . . . . . . . . . . . . . . . . . 47
4.3 Markov Chain Model for Power Demand Based on Acitivity Data . . . . . . . . 50
4.4 Uncertainty Quantification for Markov Chain Model of Power Demand . . . . . 58
v
Chapter 5 Stochastic Model for Optimally Scheduled Generation 60
5.1 Introduction and Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 Graph Theory Notations and Definitions . . . . . . . . . . . . . . . . . . . . . . 62
5.3 Discrete-time Consensus Problem . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.4 Consensus over Random Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.4.1 Random Consensus Dynamics . . . . . . . . . . . . . . . . . . . . . . . 65
5.4.2 Probabilistic Characterization ofP

(!) . . . . . . . . . . . . . . . . . . 67
5.5 Distributed Optimization of Power Generation . . . . . . . . . . . . . . . . . . . 72
Chapter 6 Stability of the Smart Grid under Uncertainty 79
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.2 Stochastic Eigenvalue Problem - Symmetric Matrices . . . . . . . . . . . . . . . 81
6.2.1 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6.2.2 Solution Algorithms for Dominant Eigenpairs . . . . . . . . . . . . . . . 86
6.2.3 Analysis of Stochastic Power Iteration . . . . . . . . . . . . . . . . . . . 89
6.2.4 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.2.5 Model Reduction Using Stochastic Eigenpairs . . . . . . . . . . . . . . . 108
6.3 Stochastic Eigenvalue Problem - Non-symmetric Matrices . . . . . . . . . . . . 116
6.4 Polynomial Chaos Representation for Random Transition Matrices . . . . . . . . 119
6.5 Impact of Demand Uncertainty on the Small Signal Stability . . . . . . . . . . . 121
6.6 Impact of Supply Uncertainty on the Small Signal Stability . . . . . . . . . . . . 123
Chapter 7 Conclusion 132
Bibliography 134
Appendices 140
vi
List of Tables
4.1 The discrepency metric
(P
1
;P
2
)
between the zones for the first set of respondents 54
4.2 The discrepency metric
(P
1
;P
2
)
between the zones and between the first and sec-
ond sets of respondents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.1 Initial conditions and parameters of the 4 generation units . . . . . . . . . . . . . 76
6.1 Values of relative L
2
-error and KL-divergence for the four dominant random
eigenpairs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.2 Initial conditions and parameters of the 4 generation units . . . . . . . . . . . . . 124
vii
List of Figures
1.1 Schematic of the hierarchy of the controllable electrical load in future Smart Grid. 3
2.1 Schematic of a single-machine infinite bus system . . . . . . . . . . . . . . . . 10
3.1 The standard deviation in the 200 estimates of asymptotic distributions (taken to
be the QoI) versus the length of the observation window - in each time window,
200 QoI’s are estimated based on the 200 sampled trajectories. . . . . . . . . . 23
3.2 The 3-state Markov chain modeling the HIV progression; states are identified
by the CD4-cell count ranges. Shown on the edges are the observed transition
frequencies between the corresponding states, which are taken to be the mean
transition probabilities in the MaxEnt estimation . . . . . . . . . . . . . . . . . 31
3.3 The joint pdf forp
22
andp
23
computed by MaxEnt RTM estimation with mean
values. The mean values are p
22
= 0:24 and p
23
= 0:18. . . . . . . . . . . . . 32
3.4 The joint pdf forp
22
andp
23
computed by MaxEnt RTM estimation with mean
values and standard deviations. The same mean values with 30% coefficient of
variation are considered. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.5 Stationary distribution for the three cell count states,
1
1
;
1
2
;
1
3
in time-
invariant Markov chain using MaxEnt RTM estimation with mean values. . . . 33
3.6 Stationary distribution for the three cell count states,
1
1
;
1
2
;
1
3
in time-
invariant Markov chain using MaxEnt RTM estimation with mean values and
standard deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.7 The scatter in the random progression path in a time-invariant Markov chain
using MaxEnt RTM estimation with mean values and standard deviation. . . . . 34
3.8 The scatter in the random progression path in a time-variant Markov chain using
MaxEnt RTM estimation with mean values and standard deviation. . . . . . . . 35
viii
3.9 The pdf of the probability that the system starting in state 1 will be at state
2 in exactly 10 steps. Circle marker refers to the average value of this random
quantity of interest. Cross marker refers to the same quantity of interest obtained
by the deterministic transition matrix. Time-invariant RTM implementation is
used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.10 The pdf of the expected first passage time from state lowest cell count state
to the highest one. Circle marker refers to the average value of this random
quantity. Cross marker refers to the same quantity of interest obtained by the
deterministic transition matrix. Time-invariant RTM implementation is used. . . 37
3.11 The histogram of the random optimal allocation of service
1
obtained by 3000
samples of the MaxEnt RTM. Circle marker refers to the average value of this
random quantity. Cross marker refers to the optimal allocation obtained by
using the deterministic transition model. Time-invariant RTM implementation
is used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.12 The pdf’s and joint pdf’s corresponding to the independent random transition
probabilities outgoing states 1 to 5, plotted in figures (a) to (e), respectively. . . 41
3.13 The posterior joint pdf’s corresponding to the independent random transition
probabilities outgoing state 4 to states 5 and 6. Fig (a) refers to the case where
three transitions from state 4 to state 5 have been observed. Fig (b) corresponds
to the case where one transition from state 4 to state 5 and two transitions from
state 4 to state 6 are observed. . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.1 Activity profiles of 6 representative respondents during the time period 4 am to
12 am. The activity indices are defined as follows (1) away, (2) sleeping, (3)
home with high consumption (washing, cooking, etc.) and (4) home with low
consumption. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2 The dynamics of the activity pattern of the consumers. Y axes correspond to the
ratio of the respondents in each activity status. The solid lines correspond to the
ratios calculated from the full data. The dashed lines refer to the Markov chain
prediction where its transition matrix is estimated based on half of the data. . . 56
4.3 The histograms of the stochastic QoI for the activity prediction problem. The
QoI are the ratios of the population in each activity state. Circle marker refers
to the average value of this random quantity obtained by the stochastic model.
Cross marker refers to the QoI obtained by using the deterministic transition
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
ix
4.4 Comparison between the pdf’s of QoI (taken to be the three asymptotic dis-
tributions). The solid line corresponds to MaxEnt RTM estimation with mean
values and standard deviation (based on the observation of only 30 samples).
The dashed line corresponds to the 900 sample trajectories drawn from a per-
turbed transition matrix. The perturbation is assumed to be uniformly dis-
tributed around the mean values. . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.1 The random trajectories show the dynamics of the incremental costs in a random
communication topology governed by a MaxEnt probability measure. . . . . . 77
5.2 The probabity density function of the random power output of unit 2 sched-
uled under the uncertainty in the communication topology. As can be seen, the
scatter around the mean values progressively decreases, as the adjustment nears
convergence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.1 The three angles between the vector and coordinates, i.e.;; and
. . . . . 96
6.2 The 2DOF mass-spring system (after Ghosh and Ghanem (2008)). . . . . . . . 98
6.3 Comparison between the exact and approximate pdf’s of the phase angles of the
first and second eigenvector with respect to the first and second axes. . . . . . . 99
6.4 Comparison between the functional forms of the exact and approximate phase
angles of the first and second eigenvector relative to first and second axes. . . . 100
6.5 Comparison between the exact and approximate pdf’s of the first and second
eigenvalues of the 2DOF system. . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.6 Comparison between the functional forms of the exact and approximate solu-
tions for the first and second eigenvalues of the 2DOF system . . . . . . . . . . 101
6.7 Discretization of the cantilever beam with 10 elements and random elasticity
E(x;!) and random density(x;!) in the most general case . . . . . . . . . . 102
6.8 Comparisons between the distribution the first random eigenvector, obtained
from Monte Carlo and Subspace Iteration; shown are PDF of phase angles with
respect to (a) first, (b) fifth, (c) sixth, (d) 11
th
, (e) 14
th
, and (f) 16
th
coordinate
axes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
6.9 Comparisons between the distribution the fourth random eigenvector, obtained
from Monte Carlo and Subspace Iteration; shown are PDF of phase angles with
respect to (a) first, (b) fifth, (c) seventh, (d) ninth (e) 13
th
, and (f) 16
th
coordinate
axes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
x
6.10 Figures on each row correspond to a particular mode and show the norm of the
corresponding random eigenvector versus three uncertainty sources,
1
,
2
and

3
, which are independent standard normal random variables. For each value of
an uncertainty source, the error bar corresponds to the variation induced by the
other two independent uncertainty sources. . . . . . . . . . . . . . . . . . . . . 107
6.11 The pdf’s of the four dominant random eigenvalues. The pdf’s obtained by
the subspace iteration algorithm overlap with those obtained by Monte Carlo
simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.12 The random displacement response of the first degree of freedom in the 2DOF
oscillator. The gray lines refer to responses resulting from Monte Carlo samples
of the stiffness matrix. The blue lines refer to the response calculated at the
5 quadrature points. The blue lines are weighted based on the weight of the
corresponding quadrature points. . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.13 The pdf of the random displacement response of the first degree of freedom in
the 2DOF oscillator at time 2 sec. The solid lines refer to responses resulting
from Monte Carlo samples of the stiffness matrix. The dashed lines refer to the
collocation based solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.14 The random displacement of the first degree of freedom in the 2DOF system
generated from the stochastic modal decomposition compared with that ob-
tained by direct MC simulations of the equation of motion. Both random eigen-
vectors are included. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.15 The random displacement of the first degree of freedom in the 2DOF system
generated from the reduced model based on the first stochastic mode compared
with that obtained by direct MC simulations of the equation of motion. . . . . . 114
6.16 The random displacement of the point at 60 cm from the fixed end in the can-
tilever beam with random elasticity. The response is calculated based on a model
constructed based on the first random eigenmode and is compared with that ob-
tained by direct MC simulations of the equation of motion. . . . . . . . . . . . 115
6.17 The random displacement of the point at 60 cm from the fixed end in the can-
tilever beam with random elasticity. The response is calculated based on a model
constructed based on the first three random eigenmodes and is compared with
that obtained by direct MC simulations of the equation of motion. . . . . . . . 115
6.18 The random displacement of the point at 60 cm from the fixed end in the can-
tilever beam with random elasticity and uncertain external force. The response
is calculated based on a model constructed based on the first three random eigen-
modes and is compared with that obtained by direct MC simulations of the equa-
tion of motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
xi
6.19 The relationship between the random transition probabilitiesx
1
(!);x
2
(!);x
3
(!)
and the input uncertainties (random uniform variables)u
1
andu
2
. Figures on
the left column correspond to the PC representation of the transition probabili-
ties, whereas those on the right one refer to the Monte Carlo samples. . . . . . 125
6.20 The comparison between the distribution of the random transition probabilities.
Dashed line refers to the PC representation of the transition probabilities. Solid
line refers to the Monte Carlo samples. . . . . . . . . . . . . . . . . . . . . . . 126
6.21 The comparison between the distributions of the random electricity demand:
solid line refers to the Monte Carlo samples of the MaxEnt joint pdf, dashed
line refers to the case where the transition matrix is represented by a 3rd order
3 dimensional PCE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6.22 The Western System Coordinating Council (WSCC) 3-machine 9-bus grid used
in the small signal stability analysis (Milano (2008)) . . . . . . . . . . . . . . . 127
6.23 The probability density function of the real part of the random eigenvalues;
comparison between results obtained from Monte Carlo simulations and those
obtained from the 3-point quadrature based spectral projection. . . . . . . . . . 128
6.24 The probability density function of the real part of the random eigenvalues;
comparison between results obtained from Monte Carlo simulations and those
obtained from the 3-point quadrature based spectral projection. . . . . . . . . . 129
6.25 The probability density function of the real part of the random eigenvalues ob-
tained from Monte Carlo simulations. . . . . . . . . . . . . . . . . . . . . . . 130
6.26 The probabity density function of the real part of the random eigenvaluesob-
tained from Monte Carlo simulations. . . . . . . . . . . . . . . . . . . . . . . 131
xii
Abstract
Today, many of the infrastructures are composed of several coupled sub-systems, many of
which are by themselves complex, and further couplings introduce yet more complexities which
pose system operators with several challenging issues. In many of these infrastructure systems,
uncertainty is an intrinsic features, whereas in several others it is merely a mathematical ab-
straction referring to our ignorance due either to the unknown governing physics or to missing
information. A major challenge is the predictability of such complex systems of systems under
these uncertainties and also the quantification of confidence in the prediction. A predictive sci-
ence for these systems should (a) assess the uncertainties; that is, to identify and characterize
them in the input variables or model parameters; and (b) predict their impacts on the performance
of the system, which will eventually assist decision and policy making processes. This study
aims to contribute to these two pillars of a successful predictive model for complex infrastruc-
ture systems. We have focused on the uncertainty management of the future Smart Grid, as an
excellent example of a complex system of systems. The future Smart Grid will be a composition
of interacting systems, including the power grid, the weather system, the social networks, the
communication network, etc. Over the course of this study, the conceptual, mathematical and
algorithmic challenges relevant to the establishment of a predictive science for this system have
been examined and novel methodologies were developed to address these challenges. The scope
of this dissertation includes the challenges induced by the uncertainties in the demand and supply
sides of the future Smart Grid. Models for the characterization of these uncertainties have been
proposed and their impacts on performance metrics of the overall system have been quantified.
xiii
Chapter 1
Introduction
Complex Systems of Systems (CSoS) are critically significant eco-socio-economic-technical sys-
tems which greatly influence our lives at different scales (Glass et. al. (2011)). As a first step
before articulating the characteristics of these systems, we define the following terms
System: A set of interacting entities that serve a common goal
System of Systems: A set of (interacting) systems that cannot be reduced to entities and may
or may not serve a common objective
Complexity: Used for various characteristics such as emergent behavior from individual
interections, nonlinearities, chaotic behavior and may also refer to uncertainties.
CSoS are usually large and irreducible systems spanning a wide range of physics with dynam-
ics spanning a wide range of time scales. Examples of these systems are cities and megacities,
infrastructures, healthcare systems, energy systems, etc. The characterization of these systems
are challenging, as their behaviors are often hysteretic and/or recursive with initial conditions that
are not repeatable. Furthermore, they involve people as a significant component which hinders
1
the experimentation on the system due to the associated risks and costs. It is, however, the aspira-
tions of the engineers to yet control and design these systems that have motivated efforts towards
developing or adopting predictive models for these systems.
Future infrastructure systems will be complex systems of interacting systems with wide ex-
change of information between people and systems. The future data intake for only a single
electricity utility company is projected to annually exceed 800 TB. There is need for computa-
tional algorithms that efficiently process this “tsunami” of data and therefore assist people and
system managers with their decisions. There is also need for new scientific tools to address the
unprecedented challenges in the analysis of these multiscale systems of interacting systems.
The power grid, amongst other infrastructures, is becoming increasingly congested. For in-
stance, despite the 25% increase in electricity demand since 1990, construction of transmission
facilities has decreased by about 30%. According to Energy Information Administration, global
electricity consumption will double in the next 25 years. Therefore, there is a vital need for opti-
mized energy consumption. An area for energy efficiency improvement is the demand side man-
agement, where the electricity load is controlled by encouraging consumers to schedule their con-
sumption towards increased efficiency. In order for such management to materialize, the condition
of the connectivity between the components of the power grid should be sufficient, a necessity that
currently is within reach, hence providing unprecedented opportunity for the establishment of a
revolutionized architecture of the power grid referred to as the Smart Grid architecture.
The future Smart Grid will be designed to entail capabilities for the incorporation of the
renewable energy resources, dynamic pricing, and many distributed control and command proto-
cols. For instance, the hierarchy of the scales for the controllable electricity demand is depicted
in Fig. 1.1, where at each scale different types of factors drive the electricity amount consumed
2
Figure 1.1: Schematic of the hierarchy of the controllable electrical load in future Smart Grid.
at the associated scale. With the advent of such revolutions, the Smart Grid will get imposed
into new sources uncertainties. First, renewable energy resources will impose weather-induced
fluctuations into the grid operation, a major factor to be accounted for in the generation schedul-
ing for load balance. Second, the demand side management relies on the consumer reaction to
the electricity price, a behavior that is complex, uncertain and therefore quite elusive. Third, the
uncertainties in the topology of the dynamic communication network due to the introduction of
the “Plug-and-Play” devices introduce new challenges in the distributed control of the grid.
In this study, we aim to address challenges in the management of two outstanding sources
of the uncertainties; (1) uncertainties in the prediction of the electricity consumption prediction
of a population of consumers, and (2) uncertainties in the distributed economic dispatch of the
generation due to the dynamic nature of the connectivities. The uncertainty management that we
pursued in this study involves three major steps. First, after the identification of the significant
3
uncertainty sources in the system, one needs to characterize and represent these uncertainties base
first principles, accepted physical laws, and available data from the experimentation or observa-
tion of the component of interest. This uncertainty characterization and representation involves
determining the shape of the probability distribution functions or determining the coefficients in
a Polynomial Chaos expansion.
Second, given the characterized uncertainties in the system parameters and/or initial or bound-
ary conditions, the second step, uncertainty propagation, involves in quantifying the impact of the
input uncertainties on the quantities of interest and other decision making variables. Finally, such
quantified uncertainties can significantly assist the decision and policy makers with their robust
decisions in the face of inevitable uncertainties.
1.1 Outline and Summary of Contribution
This dissertation is structured in the following way. Chapter 2 includes the basics of the power
system theory used in the formulation of the stability analysis of these systems. In particular, in
this study, the small signal stability condition of the grid is chosen as the performance measure
of interest for the grid. The formulation of this stability problem based on the linearization of the
power system governing equations is elaborated in Section 2.1.
In this study, we investigate the uncertainties in the demand and supply. The electricity de-
mand is modeled as a Markov chain. In Chapter 3, we propose a new probabilistic model for the
parameters of the Markov chain. First, we describe the mathematical settings based on which a
Markov chain is constructed with a random transition matrix (RTM). Next, we present the Maxi-
mum Entropy (MaxEnt) formulation used to construct the probability measure of the RTM based
4
on the observed values for the mean and standard deviation of transition rates. As illustration,
the implementation of the MaxEnt RTM is demonstrated on a simplified model of infection dy-
namics, where the implications of the uncertainty representation of the Markov chain model on
ensuing decision processes are discussed. The application of the proposed probabilistic model to
the deterioration modeling of infrastructure systems is also illustrated.
Chapter 4 concerns models for the random electricity demand incorporating human factors.
Electricity consumption can be modeled by a Markov chain. Due to the elusive and random
nature of the human behavior in electricity consumption, it will be too simplistic to prescribe a
single deterministic Markov chain to represent the whole system. Variations induced due to the
inevitable human-induced uncertainties will remain unexplored. Also, the measurement noise or
insufficient data make the case for the uncertainty treatment of these models. In this chapter, the
MaxEnt formalism of Chapter 3 is adopted for the characterization of these uncertainties. The
application of the proposed model to the random electricity demand prediction is illustrated.
In Chapter 5, we study the distributed optimization for the economic dispatch of generation
units under the uncertainties in the communication networks. We make use of the reformulation
of optimization problem as a consensus problem. A quantitative measure is proposed for the
characterization of the random graph representing the condition of the connecitivities. Finally, in a
4-generator example, the application of the proposed random consensus problem is demonstrated.
In Chapter 6, we investigate the impact of the characterized uncertainties in Chapters 4 and
5 on the performance of the grid and the quantities on interest. Specifically, we are interested in
studying the small signal rotor angle stability of the power systems. Two sections are included on
the eigenvalue analysis of symmetric and non-symmetric matrices, the latter one corresponding to
the case of the small signal stability problem. Using the proposed stochastic eigenvalue analysis,
5
the uncertainties in the quantities of interest, which are the real parts of the eigenvalues, are
quantified in a representative 9-bus power grid, providing opportunities for the robust decision
making based on the stochastic stability assessment.
6
Chapter 2
Dynamics of the Smart Grid
2.1 Introduction
The electric power grid is a hierarchical structure which transmits the electricity generator at
power plants through transmission and distribution networks to residential or industrial con-
sumers. In what follows, we describe the basics of the power theory, which serves as the founda-
tion for the formulation of the small signal stability problem of the power systems.
2.2 Power Flow
We start with the formulation of power flow in an AC power grid, with the purpose of calculating
the node voltage phasors when the injected complex power is given. In a linear network of n
nodes, the injected current at node i, denoted by
_
I
i
, is related to the injected voltages at other
nodes
_
V
j
, through the admittance matrixY , as follows
_
I
i
=
n
X
j=1
Y
ij
_
V
j
i = 1; 2; ;n: (2.1)
7
The relation of node current and power reads
_
I
i
=
P
i
jQ
i
V

i
; (2.2)
whereP
i
andQ
i
are respectively the active and reactive components of the power injected at node
i. The relation between node power and voltage will then be given by
P
i
+jQ
i
_
V
i
=
n
X
j=1
Y

ij
V

j
i = 1; 2; ;n; (2.3)
where superscript denotes complex conjugate. After manipulation, the governing power flow
equation will be given by
P
i
=V
i
X
j2i
V
j
(G
ij
cos
ij
+B
ij
sin
ij
)
Q
i
=V
i
X
j2i
V
j
(G
ij
sin
ij
B
ij
cos
ij
)
(2.4)
whereY
ij
= G
ij
+jB
ij
, and
_
V
i
= V
i
e
j
i
, and
ij
=
i

j
. The Newton-Raphson method is
used to numerically solve the power flow equations. In these equations, for each node, we have
four variables, V;;P andQ, and two equations. Thus, two variables for each node should be
given as input. Based on the type of input variables, we deal with two different types of nodes:
PQ nodes: Usually load buses are PQ nodes, since the load powers are assumed fixed.
PV nodes: Usually generator buses are PV nodes, since there is controllable reactive power
to keep the voltage magnitude fixed.
8
2.3 Dynamics of the Power Systems
The dynamics model of synchronous machines is similar to the equation of (linear) motion. That
is, as we relate the acceleration to the force in the linear motion, for synchronous machines, we
derive the relation between speed deviation and total torque. The resulting equation is called the
swing equation.
Consider a generator running at the synchronous speed !
s
with respect to other generators
produces an electromagnetic torqueT
e
. In a no-loss case, the driving mechanical torqueT
m
under
steady-state conditions equals the electromagnetic torque. Any disturbance to the synchronous
behavior will depart the system from the steady-state condition and create an inertia in the system
J
d
2

dt
2
=T
m
T
e
; (2.5)
where is the angular displacement of the rotor with respect to the stationary reference axis on the
stator, andJ is the combined moment of inertia of the prime mover and the generator. Consider a
single-machine infinite bus system (Fig. 2.1), where the generator voltage angle leads that of the
infinite bus by. If we represent this angular displacement relative to the synchronously rotating
reference frame rotating at a constant angular velocity!
s
, we get
= !
s
t + (2.6)
d
2

dt
2
=
d
2

dt
2
; (2.7)
9
g
Figure 2.1: Schematic of a single-machine infinite bus system
where we have represented the rotor acceleration in terms of the rotor position at time t = 0.
Then Eq. 2.5 can be rewritten as
J
d
2

dt
2
=T
m
T
e
; (2.8)
and in terms of power as
M _ ! =P
m
P
e
; (2.9)
whereM = J! is called the inertia constant, and! =
_
. AlthoughM is not constant because
of the deviations of!, near the stability the deviations can be assumed to be small, and thusM
is considered to be constant and evaluated at the synchronous speed!
s
. The swing equation can
also be derived in the following similar form
2H
!
s
_ ! =P
m
P
e
; (2.10)
whereH is a given machine property. If the damping is to be considered, the swing equation will
be given by
M _ ! =P
m
P
e
D
0
!; (2.11)
whereD
0
is the damping coefficient.
10
Now, we describe the calculation of output electric powerP
e
. The power transferred between
the machine and the bus is given by the following highly nonlinear equation called power-angle
relation
P
e
=
E
G
E
B
X
sin (2.12)
whereE
G
andE
B
denote the voltage at the generator and the bus, andX is the inductive reac-
tance. Based on Eq. 2.12, in a multimachine system ofn generators, the electrical power output
of machinei is given by
P
ei
=E
2
i
G
ii
+
X
j2i
E
i
E
j
(B
ij
sin
ij
G
ij
cos
ij
) i = 1; 2; ;n: (2.13)
where
ij
=
i

j
.
2.4 Small Signal Stability of the Smart Grid
For the power system to remain operational, it should satisfy certain stability conditions. An
IEEE/CIGRE task force has been assigned to unify the definitions and classifications of the sta-
bility of power system (Kundur et al. (2004)). According to this task force the stability of the
power system is defined as:
“Power system stability is the ability of an electric power system, for a given initial operating
condition, to regain a state of operating equilibrium after being subjected to a physical distur-
bance, with most system variables bounded so that practically the entire system remains intact.”
Also in this study, the stability is classified based on the main system variable in which insta-
bility can occur, the size of the physical disturbances, and the time span that should be considered
11
for the stability study. Three main classes of stability are then defined as (1) rotor angle stability,
(2) voltage stability, and (3) frequency stability.
Rotor angle stability means the ability of the synchronous machines to restore synchronism
when faced with disturbance, and is a problem of electromechanical oscillation of synchronous
machines. Small-signal rotor angle stability is the ability the synchronous machines to remain
synchronized or ”in-step” when faced with a small disturbance. The small disturbance assumption
permits the use of linearized equation of motion. Large-disturbance or transient rotor angle
stability, on the other hand, is the ability of the synchronous machines to maintain synchronism
when a large transient disturbance is introduced to the system. In this section, we explain the
formulatation of the small-signal stability problem, which will later be used as the performance
metric of the grid and will be studied under the impacts of demand and supply uncertainties.
In order to study the small signal stability of the power systems, we need to first construct
the governing equation for the dynamics of the rotor angles. These governing equations can be
represented in a general differential-algebraic form (Eq. 2.14). The state equation consists of
differential equations for the generator and excitation system, and algebraic equations are derived
from the power flow equation and include the stator and network variables.
_ x =f(x; y; u; p)
0 =g(x; y; p)
(2.14)
12
where x and y are vectors of state (differential) variables and algebraic variables, respectively; u
is the input vector and p is the parameter vector. The linearized form of Eqs. 2.14 will be given
by
_ x =Ax +By
0 =Cx +Dy
(2.15)
In a system ofn generators, we have x = [x
1
; ; x
n
]
T
, where x
i
2R
s
i
is the state variables
of machinei, ands
i
denotes the number of total state variables considered for the dynamics of
machinei. Depending on the type of the assumptions, different sets of differential variables may
become relevant. One class of assumptions (see, e.g. Sauer and Pai (1990)) leads to the widely
used formulation called ‘the classic model’, in which for each generatori we only consider two
state variables, i.e. x
i
= [! ]
T
2 R
2
. There are more detailed models, however, where other
state variables, such asE
0
d
andE
0
q
, can also be considered. In general, we can write x = [!]
T
where is the vector of additional state variables.
The set of algebraic variables in y can also differ from one model to another, but will always
include the variables of the network power flow equations, i.e. the magnitude and angle of the
(complex) voltage values at the network nodes, denoted byV and, respectively. In more detailed
models, it could also involve the currentsI
d
andI
q
.
We start with the description of the formulation of the small-signal rotor angle stability prob-
lem for a single-machine infinite-bus. The linearization of Eq. 2.12 and Eq. 2.5 yields
M
d
2
(
0
+

)
dt
2
=P
m

E
G
E
B
X
sin(
0
+ ); (2.16)
13
which with the assumption that is small can be rewritten as
M
d
2

dt
2
+P
max
cos
0
() = 0: (2.17)
P
max
cos
0
is the slope of the power-angle curve at
0
and is called the synchronizing torque
coefficient, denoted byK
s
which plays an important role in the stability.
In a multimachine system ofn generators, as mentioned earlier, the electrical power output of
machinei is given by
P
ei
=E
2
i
G
ii
+
X
j2i
E
i
E
j
(B
ij
sin
ij
G
ij
cos
ij
) i = 1; 2; ;n: (2.18)
Where
ij
=
i

j
. Having the increments
i
=
i0
+
i
and
ij
=
ij0
+
ij
, we can
linearize these equations about the operating condition
ij0
and arrive at
M
i
d
2
(
i
)
dt
2
+
X
j2i
E
i
E
j
(B
ij
cos
ij0
G
ij
sin
ij0
)
ij
= 0 i = 1; ;n (2.19)
where E
i
is the constant voltage of machine i, and Y
ij
= G
ij
+jB
ij
is the elements of the
admittance matrix given by a reduction scheme so that only the generator nodes are present in the
network.
Eq. 2.19 is a set ofn1 independent second order equations, since all the angles are calculated
relative to a single reference generator. Let us take then-th generator to be the reference generator.
If we subtract then-th equation from thei-th equation, we will get
14
d
2
(
in
)
dt
2
+
!
s
2H
X
j2i
P
ij

ij

!
s
2H
X
j2n
P
nj

nj
= 0 i = 1; 2; ;n 1; (2.20)
whereP
ij
=E
i
E
j
(B
ij
cos
ij0
G
ij
sin
ij0
).
Let x
1
2R
n1
denote the vector of relative generator angle changes, x
1
= [
1n
; ;
(n1)n
]
T
.
The time derivatives of these angle changes are therefore given by x
2
= _ x
1
. The block matrix
representation of Eq. 2.19 is as follows
2
6
6
4
_ x
1
_ x
2
3
7
7
5
=
2
6
6
4
0 I
A 0
3
7
7
5
2
6
6
4
x
1
x
2
3
7
7
5
; (2.21)
where A
ij
is constructed from Eq. 2.20. The eigenvalue analysis of the (2n 2) (2n 2)
matrix in Eq. 2.21 is carried out for the small signal stability study. As mentioned earlier, this
representation of the generators is called the classical model representation, where the only two
state variables are and!.
Network Equations
As mentioned earlier, the algebraic part of the governing equation (Eq. 2.14) is derived from the
power flow equations. It essentially establishes the interdependencies between algebraic vari-
ables, i.e. the voltage magnitudes and angles at the network buses. In other words, it relates the
dynamics of individual generators to the rest of the network. As mentioned earlier, the set of
algebraic variables will always contain the voltage values at the generator or load buses. These
15
variables are in fact time-varying. Based on the modeling assumption, the dynamics of these
variables can be ignored or retained (Kundur (1994))
16
Chapter 3
Uncertainty Quantification for Markov Chain
Models
3.1 Introduction
Markov chains are ubiquitous in modeling phenomena that exhibit random fluctuations while
evolving in time. Their value lies in their simplicity, both of interpretation and implementation.
Indeed, the mere justification of finite memory with constant structure, allows the postulation of
a Markovian model from which consistent sample paths for futures of the system can be real-
ized. Decision processes associated with this behavior can also be readily formalized yielding
elegant analytical insight. These Markovian models are typically parameterized by their initial
state and the probabilities of transition between states. These so-called transition probabilities are
usually estimated from the evolution of the Markov chain observed along a number of sample
paths. Given the dichotomy between the simplicity of Markovian models, and the highly com-
plex natural, man-made, and social phenomena to which they have been adapted, a framework for
assessing the uncertainty in associated mode-based inferences is essential. One can contemplate
17
two constructions of Markov chains that provide different perspectives on these uncertainties. The
first construction maintains the belief that a Markov model is a suitable description of the under-
lying motives such as physical or social processes. In this case uncertainty is attributed to lack of
knowledge (due perhaps to insufficient observations) of the true values of transition probabilities,
which are themselves treated as random. An ensemble of Markov chains is thus obtained, that
is in one-to-one correspondence with the ensemble of transition matrices. The second construc-
tion recognizes that a Markov model is merely a convenient representation of the processes of
interest, and develops adaptation strategies as information is acquired. Accordingly, the transition
matrix of the Markov chain is permitted to vary in time, effectively, refitting a new Markov chain.
Clearly, in this construction, the amount of data available to re-fit a Markov chain is more limited
than in the first construction, and uncertainty in the estimated transition probability can be signif-
icant. In this dissertation, we will refer to these two constructions as the time-invariant and to the
time-variant models, respectively. We present a methodology for ascertaining the uncertainty in
transition probabilities of Markov chains, and demonstrate its value in quantifying the uncertainty
in the predicted steady-state dynamics and associated decisions.
Without loss of generality, we will focus our attention on finite-state discrete-time Markov
chains. The time evolution of the state probabilities from a time instance to the next is determined
by the entries of the transition matrix which describe the transition rates between two states of the
chain. These transition rates are usually estimated as deterministic model parameters using least
squares, Maximum Likelihood, Most a-Posteriori, or Maximum Entropy (MaxEnt) arguments.
Within the confines of probability theory, our objective, then, is to characterize a probability
measure for the transition matrix that is consistent, in some sense, with mathematical constraints
and available evidence.
18
Cooman et al. (2009) constructed a Markov chain based on the concept of imprecise prob-
abilities. They formed the credal sets for the probability mass function and the transition proba-
bilities, whose extreme points were used in estimating upper and lower bounds on the expected
behavior of the Markov chain, yielding results similar to the classical Perron-Frobenius theorem
(Luenberger (1979)) on the asymptotic behavior of the associated Markov chains. No probabilis-
tic characterization of the state distributions of the Markov chain was introduced, a factor that
distinguishes it from our formulation. Alternatively, Li and Si (2008) and Abbad and Filar (1992)
investigated the dynamics of Markov chains under parametric perturbations of transition matrices.
This approach is, however, restrictive in that the mathematical form of the parametric perturbation
together with the number and range of the parameters should be manually specified. In our for-
mulation, the variability in the transition matrices follows from imposing the Maximum Entropy
principle where only summary statistics based on observation are exploited and other subjective
assumption may not pollute the estimation.
A different approach proposed in the literature is based on characterizing a support domain
for random transition matrices (RTM). In Nilim and El Ghanoui (2005), the authors characterized
this domain by considering first the deterministic Maximum Likelihood estimate of the transition
matrix. They then formed a set of transition matrices whose members have likelihoods that are
within a prescribed range with respect to the maximum likelihood. They also proposed a simi-
lar construction by considering those transition probabilities that are within a Kullback-Leibler
distance with respect to a given prior transition probability. These proposed support regions are
subsequently used in the control of the Markov decision processes. This approach, however, does
not permit the probabilistic characterization of the fluctuations in Markov Chain’s behavior, hence
the limitation in its inferential value.
19
Our approach involves characterizing a probability measure over the set of all mathematically
admissible transition matrices, constraining them to be valid stochastic matrices. We estimate
this measure using the MaxEnt principle given the information on the mean values and possibly
higher moments. We make use of no other assumptions in the estimation procedure. With a
MaxEnt RTM available, a Markov chain can be constructed in two different ways, describing
two distinct behaviors. In the first one, referred to as Markov chain with time-invariant RTM,
a single sample of the random transition matrix is assumed to govern the Markov chain’s time
evolution. In the second implementation, referred to as Markov chain with time-variant RTM, at
each time step, a new independent sample is drawn from the RTM based on which the Markov
chain makes a one-step transition. The second variation is also referred to as Markov chains in
random environment (Cogburn (1984, 1991), Lasserre (1997), Lu and Mukherjea (1997), Orey
(1991), Stenflo (2001)).
This section is structured as follows. After motivating our research using a simple numerical
example, we describe the mathematical settings based on which a Markov chain is constructed
with an RTM. Two different settings are discussed for the time-invariant and time-variant imple-
mentations. Next, we present the MaxEnt formulation used to construct the probability measure
of the RTM based on the observed values for the mean and standard deviation of transition rates.
As illustration, the implementation of the MaxEnt RTM is demonstrated on a simplified model of
infection dynamics, where the implications of the uncertainty representation of the Markov chain
model on ensuing decision processes are discussed.
20
3.2 Motivation
The stochastic behavior described through Markovian dynamics is typically attributed to white
noise disturbances which, when closely viewed within a multiscale context, could be further ex-
plained by subscale interactions. It is thus the effect of these unmodeled physics, or more gener-
ally causal relationships, that is captured by the random evolution of a Markov chain. As indicated
above, additional uncertainties accrue in the process of calibrating the transitions of these chains
from finite samples. While standard approaches to analyzing Markovian behavior aggregates the
effect of these uncertainties, their segregation provides a path towards addressing each of them
appropriately.
First, and through a simple numerical example, we seek to emphasize the significance of
representing and propagating uncertainty in transition matrices. Consider a single dynamical
system with three possible states and the corresponding transition matrix known to be
P =
2
6
6
6
6
6
6
4
0:92 0:04 0:04
0:58 0:24 0:18
0:23 0:23 0:54
3
7
7
7
7
7
7
5
(3.1)
We argue that one cannot deterministically recover the true transition matrix by observing a
single realization or even a few realizations from a Markov chain. To do so, we generate 200
sample trajectories, each including 1000 time-steps, based on the known transition matrix. We
21
then deterministically estimate the transition matrix entries for each trajectory, using transition
frequencies. That is, the deterministic transition rate of thek-th trajectory,p
k
ij
, is calculated by
p
k
ij
=
c
k
ij
P
j
c
k
ij
; (3.2)
wherec
k
ij
denotes the number of transitions from statei to statej in thek-th trajectory recorded
over a given time window. We expect, as the size of the time window increases, the bias in
200 estimates of the transition matrix to decrease. However, limited resources often impede the
long time observation of a chain. This will result in considerable errors in predicting the “exact”
transition matrix, even in the rare cases that such underlying matrix uniquely exists.
In order to investigate the impacts of limited observation on the predicted Quantities of Interest
(QoI), we consider a set of time windows with different size, ranging from 20 to 1000. For each
fixed time window size, we estimate the asymptotic distributions for each of the 200 trajectories.
This asymptotic distribution, indicated by ^ , is our QoI, which is known to satistfy, for a given
transition matrixP ,
^ = ^ P: (3.3)
Fig. 3.1 shows the standard deviation in the 200 estimated QoI versus the size of observation win-
dow, where significant scatter in the estimated QoI can be observed. One can, therefore, conclude
that the transition matrix estimation based on observation data involves inevitable variability and
error, hence the motivation for treating these matrices as random.
22
0 200 400 600 800 1000
0
0.05
0.1
0.15
0.2
0.25
observation window size
standard deviation

π
1
π
2
π
3
Figure 3.1: The standard deviation in the 200 estimates of asymptotic distributions (taken to be the
QoI) versus the length of the observation window - in each time window, 200 QoI’s are estimated
based on the 200 sampled trajectories.
3.3 Construction of Uncertain Markov chains
As indicated above, probabilistic transition matrices can be motivated by either time-invariant
or time-variant systems. The mathematical setting for these two implementations is discussed
next, in the context of finite-state discrete-time Markov chains. Consider the probability triple
(S;F
S
;), whereS =f1; ;ng is the sample space (finite set of states),F
S
is a-algebra on
S and = [
1
; ;
n
] is the probability measure (probability mass function) on (S;F
S
), which
is referred to as the state distribution. LetX =fX
1
;X
2
;g denote the states that the Markov
chain resides in at different times, with X
t
2S; t = 1; 2; . Let P
t
denote the transition
matrix at timet of this chain with itsij-th component,p
t
ij
, defined by
p
t
ij
= Pr(X
t+1
=jjX
t
=i): (3.4)
23
Transition matrices, also referred to as stochastic matrices, are square matrix with non-negative
components satisfying
P
n
j=1
p
t
ij
= 1; 8i2S. The dynamics of the state distributionsf
t
g is
assumed to be governed by the following finite-memory equation

t+1
=
t
P (3.5)
The set of admissible probability measures on the setS, denoted by and the set of admissible
transition matrices, denoted byP, can be formed as follows
=f2R
1n
j
n
X
i=1

i
= 1;
i
2 [0; 1]g
P =fP2R
nn
j
n
X
j=1
[P ]
ij
= 1 8i2S;p
ij
2 [0; 1]g
(3.6)
3.3.1 Time-invariant RTM
Consider (

0
;F

;) with

0
=P. The RTM, denoted byP (!), is aP-valued random matrix
whoseij-th component determines the transition probabilities from statei to statej for!2

0
. A
realization of time-invariant Markov chain is constructed with a single transition matrix sampled
at!
s
from the probability measure ofP (!) and is denoted byX
s
. The dynamics of this chain is
governed by the transition matrix whoseij-th component is defined as
[P (!
s
)]
ij
= Pr(X
t+1
s
=jjX
t
s
=i;! =!
s
) 8t (3.7)
24
Associated with the sample Markov chains constructed based onP (!
s
), there exists a stationary
distribution ^
TI
s
which is its left eigenvector, i.e.
^
TI
s
= ^
TI
s
P (!
s
) (3.8)
Using sufficient samples of the RTM, one can estimate the pdf of the stationary distribution
^
TI
(!) which will predict the asymptotic behavior of the Markov chain by
^
TI
(!) = ^
TI
(!)P (!) a.s. (3.9)
3.3.2 Time-variant RTM
In this section, we present the mathematical setting for the time-variant RTM implementation
(Cogburn (1984), Takahashi (1969)). Consider the probability triple (

0
;F

;) defined in the
previous section. We define the product probability triple (
;B;) where its elements are given
by

=f [!
1
;!
2
; ] : !
t
2

0
; 8t2N
+
g
B =F

F

=
(3.10)
Let the mappingP
t
:
!P be thet-th transition matrix, defined asP
t
(!) = !
t
. In other
words, one deals with a set of independent random transition matricesfP
1
(!);P
2
(!);g each
with identical probability measure. Theij-th component of theP-valued random matrixP
t
(!)
is denoted by [P
t
(!)]
ij
and governs the transition probabilities from statei to statej in one step
25
in the random environment!
t
. Therefore, if the chain starts off at positionx
1
the probability that
it will visitx
1
;x
2
; ;x
T
at times 1; 2; ;T is given by
P
!
x
1
[X
1
=x
1
;X
2
=x
2
; ;X
T
=x
T
] =
1
x
1
(x)[P
1
(!)]
x
1
x
2
[P
2
(!)]
x
2
x
3
[P
T1
(!)]
x
T1
x
T
(3.11)
where 1
x
1
is the indicator function for statex
1
.
A Markov chain in random environment is then defined as a sequence of probability measures
f
t
2 :t2N
+
g that are generated by the RTMP
t
(!) according to the following definition

t+1
(!) =
t
(!)P
t
(!) -a.s. (3.12)
3.4 MaxEnt Formulation
We seek to characterize the probability measure on the support of random matrices,P, using
MaxEnt principle. We assume the rows of the RTM to be statistically independent from each other.
This assumption implies that the change in the transition pattern outgoing from node i has no
impact on the transition pattern outgoing from nodej. Therefore, in this exposition, for the sake of
notation brevity, we formulate the uncertainty representation for a given row of a transition matrix,
called a transition row hereinafter. To this end, let the random vector r(!) := [r
1
(!); ;r
n
(!)]
denote a particular transition row of the RTM.
Our MaxEnt formalism for the estimation of probability measure for the multivariate r(!)
relies on two classes of available information; (1) the support for the probability measure, and
26
(2) the statistical moments of r(!) obtained from data. With regards to availability of the sta-
tistical moments, we consider two different cases in this work. The first case corresponds to
situations where only the mean values of the entries ofr(!) are available from the observation
data, whereas the second case corresponds to the situations where we also include the second
statistical moments. In what follows, we will elaborate on the implications of each case on the
corresponding MaxEnt probability measure and, henceforth, on the Markov chain behavior.
Let [ m
1
; ; m
n
] denote the mean values forr(!). These mean values can be calculated
using the transition records (i.e. using Eq. 3.2). The sole use of mean values may fail to provide
enough information about the variability of the transition probabilities. Using two conceptual
cases, we will justify the need for including second moments, in addition to mean values, in the
MaxEnt formulation.
First, consider a Markov chain which is observed for a long time period. Let us further assume
that this observation time period consists of a number of time windows of smaller size close to
the characteristic time scale of the system. In general, the mean transition probabilities estimated
based on one time window may differ from those calculated in another window. By incorporating
the second order statistics of the observation data obtained from different time windows, one can
account for the level of the inter-window fluctuations around the mean value.
Second, consider a dynamical system that is composed of multiple heterogeneous constituents.
Let us assume that the overall behavior of the system is to be represented by a single Markov chain
and only certain types of the constituents are to be observed. In order to capture variability in the
transition frequencies from one type to another, a robust parameter estimator needs to accom-
modate the second order statistics, hence the justification for including standard deviation in the
MaxEnt estimation.
27
In summary, one can think of the total observation data to be decomposed into patches; tem-
poral patches in the first case, and categorical patches in the second one. To account for the
variabilities, we calculate the standard deviation, based on the temporal or categorical patches,
for the transition probabilities on the rowr and denote them by [s
1
; ;s
n
].
Having the information on the mean values and standard deviations, our objective is to charac-
terize the MaxEnt probability measure knowing that the random multivariater(!) takes values in
. Letf
r
(r
1
; ;r
n
) denote a joint probability density function for the multivariater(!). Given
the one-sum constraint, we can relate then-th component of the multivariate to other components
using
r
n
(!) = 1
n1
X
i=1
r
i
(!) a:s: (3.13)
In other words, in each random row, there are indeedn 1 independent components, and thus we
consider an (n 1)-dimensional multivariate denoted byr
0
(!) = [r
1
(!); ;r
n1
(!)] to rep-
resent the random row. Therefore, we can fully characterize the uncertainty in then-dimensional
random row by the lower dimensional joint pdf f
r
0(r
1
; ;r
n1
) together with the relation in
Eq. 5.14.
Definition [MaxEnt jpdf] Letf
r
0(r
1
; ;r
n1
) denote a joint pdf onn 1 components of
the transition rowr
0
(!). The MaxEnt joint density function, f

r
0
(r
0
), among the family of the
joint probability distribution, is the one given by
f

r
0(r
0
) = arg max
f
r
0

Z
V
0
f
r
0(r
0
) lnf
r
0(r
0
)dr
0
(3.14)
28
s:t:
Z
V
0
f
r
0(r
0
)dr
0
= 1
Z
V
0
r
i
f
r
0(r
0
)dr
0
= m
i
i = 1; ;n 1
Z
V
0
r
2
i
f
r
0(r
0
)dr
0
=s
2
i
+ m
2
i
i = 1; ;n 1
(3.15)
where the supportV
0
is characterized by
V
0
=
(
r
0
R
n1
j
n1
X
i=1
r
0
i
1; r
0
0
)
(3.16)
A solution to equations 3.14-5.22 is obtained by solving the following optimization problem,
sup
;;
inf
f
r
0
L(f
r
0;;;) (3.17)
whereL() is the associated Lagrangian,;2R
n1
, and2R
n1
are the Lagrange multipliers
associated, respectively, with the first, the second and the third of equations 5.22.
For the case where only knowledge of mean values are used, the last constraint in the Eqs. 5.22
is not used. The MaxEnt distributionf

in this case is of the form
f

r
0(r
0
) =e
(1)
exp
"
n1
X
i=1

i
r
i
#
1
V
0(r
0
) (3.18)
29
where 1
V
0 is the indicator function forV
0
. If the standard deviations are also included the MaxEnt
solution will instead have the following form
f

r
0(r
0
) =e
(1)
exp
"
n1
X
i=1

i
r
i
+
n1
X
i=1

i
r
2
i
#
1
V
0(r
0
) (3.19)
The solution of the dual Lagrangian problem involves solving 2n 1 equations for 2n
1 Lagrange multipliers. We resort to numerical integration in order to calculate the multiple
integrals due to the difficulty imposed by the supportV
0
.
3.5 Numerical Illustration I - A 3-state Markov Chain
We now present a numerical example that illustrates the behavior of Markov chains with a RTM.
We consider the Swiss HIV cohort study (Ledergerber et al. (1994)) where the progression of HIV-
infected subjects at the greatest risk of developing the Mycobacterium Avium Complex (MAC)
infection is studied. The progression is considered to involve transitions between three states
categorized by CD4-cell count ranges (with and without AIDS). The transitions between these
three states are observed between 1993 and 1995 in six-month intervals. Based on the observed
transition records, the average transition matrix is found to be
P =
2
6
6
6
6
6
6
4
0:92 0:04 0:04
0:58 0:24 0:18
0:23 0:23 0:54
3
7
7
7
7
7
7
5
(3.20)
The diagram representing the Markov chain showing the mean values is depicted in Fig. 3.2.
In the MaxEnt estimation based on mean values and standard deviations, we artificially induce
30
Figure 3.2: The 3-state Markov chain modeling the HIV progression; states are identified by
the CD4-cell count ranges. Shown on the edges are the observed transition frequencies between
the corresponding states, which are taken to be the mean transition probabilities in the MaxEnt
estimation
standard deviations for these transition rates by assuming a constant coefficient of variation of
30% for all nine transition rates.
Three MaxEnt joint pdf’s were then independently obtained for the rows of the RTM. Figs. 3.3
and 3.4 depict the MaxEnt joint pdf for the two independent random components on the second
row, i.e. for the transition rates outgoing state 2, and respectively correspond to the MaxEnt
formulation with only mean values and that with mean values and standard deviations. It is of
worth noting that the inclusion of the standard deviation have resulted in a joint density function
that is more concentrated around the associated mean value.
We next investigate the impact of the uncertainty representation of transition rates on the
behavior of the Markov chain. We first focus on the time-invariant implementation and quantify
the uncertainty in the stationary or asymptotic state distribution as a major QoI. Figs. 3.5 and 3.6
show the variability in these QoI where, as expected, less scatter is observed when second order
statistics is included in the MaxEnt density estimation.
31
0 0.5 1
0
0.5
1

p
22

p
23
0.5
1
1.5
2
2.5
3
3.5
Figure 3.3: The joint pdf forp
22
andp
23
computed by MaxEnt RTM estimation with mean values.
The mean values are p
22
= 0:24 and p
23
= 0:18.
0 0.5 1
0
0.5
1

p
22

p
23
0
5
10
15
20
25
30
35
Figure 3.4: The joint pdf forp
22
andp
23
computed by MaxEnt RTM estimation with mean values
and standard deviations. The same mean values with 30% coefficient of variation are considered.
32
0 0.2 0.4 0.6 0.8 1
0
1
2
3
4
5
6
7
8
9
10
π
pdf

π
1
π
2
π
3
Figure 3.5: Stationary distribution for the three cell count states,
1
1
;
1
2
;
1
3
in time-invariant
Markov chain using MaxEnt RTM estimation with mean values.
0 0.2 0.4 0.6 0.8 1
0
5
10
15
20
25
π
pdf

π
1
π
2
π
3
Figure 3.6: Stationary distribution for the three cell count states,
1
1
;
1
2
;
1
3
in time-invariant
Markov chain using MaxEnt RTM estimation with mean values and standard deviation.
33
0 2 4 6 8 10
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
π
2
Time step [six months]

deterministic
MC: mean
MC: mean+std
MC: mean−std
Figure 3.7: The scatter in the random progression path in a time-invariant Markov chain using
MaxEnt RTM estimation with mean values and standard deviation.
The uncertainty in the time evolution of the state distributions is studied next. Specificaly, we
focus on the dynamics of the state distribution corresponding to state 2, denoted byf
t
2
g, when
the initial state distribution is assumed to
0
= [0:7 0:2 0:1]. We study this dynamics in a
time-invariant and time-variant construction of Markov chain using the MaxEnt RTM. Fig. 3.7
compares the resulting random trajectory in a time-invariant construction with the deterministic
trajectory obtained based on the average transition matrix. On the other hand, Fig. 3.8 refers to
the time-variant chain where at each time step an independent realization of the random transition
matrix is used. An ensemble of 400 Markov chains with time-variant RTM is considered based
on which the scatter in the state distribution
2
shown in the figure is calculated. As the chain
marches through time and consequently traces of initial conditions fade away, one can expect the
variability in these state distributions to depend on a finite-length stretch in the pseudo-stationary
regime of the probabilistic trajectory. Thus, the probabilistic descriptors, such as pdf’s, in this
pseudo-stationary region are expected not to differ significantly from one time step to another; a
fact that is evident in Fig. 3.8.
34
0 5 10 15 20
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
π
2
Time step [six months]

deterministic
MC: mean
MC: mean+std
MC: mean−std
Figure 3.8: The scatter in the random progression path in a time-variant Markov chain using
MaxEnt RTM estimation with mean values and standard deviation.
Besides the stationary state distribution, another important decision-supporting variable is the
probability that a system transits from a specific state to another in a given number of steps. This
variable for the transition from statei to statej in exactlym steps is given byp
(m)
ij
= [P
m
]
ij
, i.e.
theij-th component of matrixP
m
. In the time-invariant RTM implementation, we have
p
(m)(TI)
ij
(!) = [P (!)
m
]
ij
a:s:; (3.21)
whereas for the time-variant case we have
p
(m)(TV )
ij
(!
1m
) = [P
1
(!)P
2
(!)P
m
(!)]
ij
a:s:; (3.22)
where!
1m
= [!
1
; ;!
m
]. Fig. 3.9 depicts the pdf of transition probability from a lowest cell
count state (state 1) to next cell count state (state 2) in 10 steps in a time-invariant implementation.
Since, we are dealing with the product of a number of independent identically distributed random
matrices with their common mean value equal to the average matrix of Eq. 3.20, the average of
35
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16
0
5
10
15
20
25
P
12
(10)
(ω)
pdf

pdf ofP
(10)
12
(ω)
E[P
(10)
12
(ω)]
¯
P
(10)
12
Figure 3.9: The pdf of the probability that the system starting in state 1 will be at state 2 in exactly
10 steps. Circle marker refers to the average value of this random quantity of interest. Cross
marker refers to the same quantity of interest obtained by the deterministic transition matrix.
Time-invariant RTM implementation is used.
the estimated random QoI is expected to coincide with the QoI obtained from the deterministic
model with the average matrix. This is clearly verified in the figure.
Next, we consider the uncertainty quantification of another important decision variable called
the first passage time, which is the time it takes the dynamical system to reach a particular state
starting from different initial states. In our example, lett
3
= [t
13
t
23
]
T
denote the expected first
passage times to state 3, from states 1 and 2, which is given by
t
3
= (I +Q +Q
2
+Q
3
+ )1
= (IQ)
1
1;
(3.23)
where 1 = [1 1]
T
and
Q =
2
6
6
4
p
11
p
12
p
21
p
22
3
7
7
5
: (3.24)
36
0 20 40 60 80 100 120 140
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
t
13
TI
(ω) [steps]
pdf

pdf oft
TI
13
(ω)
E[t
TI
13
(ω)]
t
TI
13
(deterministic)
Figure 3.10: The pdf of the expected first passage time from state lowest cell count state to the
highest one. Circle marker refers to the average value of this random quantity. Cross marker refers
to the same quantity of interest obtained by the deterministic transition matrix. Time-invariant
RTM implementation is used.
For time-invariant and time-variant RTM implementation we have
t
TI
3
(!) = (IQ(!))
1
1 a:s:;
t
TV
3
(!
1m
) = [I +Q
1
(!) +Q
1
(!)Q
2
(!) +
+Q
1
(!)Q
m
(!)
| {z }
m times
] 1 asm!1 a:s:;
(3.25)
where Q
t
(!) is constructed from the random transition matrix P
t
(!) according to Eq. 3.24.
Fig. 3.10 shows the pdf of the expected first passage time from the lowest cell count to the highest
cell count for the time-invariant case. The nonlinear form of this QoI induces the discrepancy
between the average value of the random QoI and the QoI obtained from the deterministic model
that is apparent in the figure.
In order to further emphasize the significant impact of the uncertainty treatment of the Markov
chains on policy making, a resource allocation problem for this disease progression example
is considered. Let us assume that two types of medical services, denoted by
1
and
2
, are
37
available to enhance the conditions of the populations in states 1 and 2. The policy maker is
then concerned with the optimal allocation of each medical service according to the following
optimization problem
[

1
;

2
] = arg max
[
1
;
2
]
2
X
i=1
G
i
(
1
;
2
)
2
X
i=1
C
i
(
i
)
s:t:
i

i

i
i = 1; 2
12
1
+
2
=B
(3.26)
where
1
= 0:5,
1
= 1,
2
= 5,
2
= 12 are the bounds on the two types of medical
service. The total budget to be invested on the services is assumed to beB = 12. The objective
is to maximize the utility minus the implementation costs. The utility, which is a measure of
acceptable health conditions, is assumed to be the sum of mean first time passage to the third
state, i.e.G = [G
1
G
2
]
T
=t
3
. The implementation cost for each service is assumed to be given
byC
1
(
1
) = exp(
1
) andC
2
(
2
) = 0:0001 exp(
2
). Let us assume that the implementation of
these medical services will impact the disease progression by changing the transition rates in the
Q matrix (Eq. 3.24) according to
Q =
2
6
6
4
p
11
(1 +
1
=100) p
12
(1
1
=50)
p
21
(1 +
2
=100) p
22
(1 +
2
=100)
3
7
7
5
: (3.27)
We seek to investigate the implication of the random treatment of the transition rates on the
optimal resource allocation for a time-invariant disease progression. To this end, we solve the
optimal resource allocation problem (Eqs. 3.26) for the average transition matrix and also for
3000 samples obtained from the MaxEnt probability measure of the RTM. The histogram of
38
0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0
50
100
150
200
250
300
350
400
450
α
1
Figure 3.11: The histogram of the random optimal allocation of service
1
obtained by 3000
samples of the MaxEnt RTM. Circle marker refers to the average value of this random quantity.
Cross marker refers to the optimal allocation obtained by using the deterministic transition model.
Time-invariant RTM implementation is used.
Fig. 3.11 shows the variability in the optimal allocation of medical service 1 resulting from the
3000 samples. As can be seen, there is a wide scatter in the optimal allocation. More importantly,
the random model provides the policy makers with a point estimate (i.e. a mean value) that is
profoundly different from the one obtained by the deterministic model. This further accentuates
the significance of uncertainty quantification in decision making procedures.
3.6 Numerical Illustration II - Markov Chain Deterioration Model
In order to optimize the maintenance policy for infrastructure systems and structures, it is essential
to construct sufficiently accurate models that can assist the decision making procedures concern-
ing the life cycle condition of the systems of interest. The deterioration and degradation behavior
of infrastructure systems plays a significant role in directing the optimal maintenance strategies.
To model this deterioration processes, accurate numerical models that can quantitatively assess
and predict the future conditions of these systems should be developed. The constructing of these
39
deterioration models rely on the experts’ judgments as well as measurement data, both of which
induce great uncertainties into the formalism. Thus, it would be too simplistic to conjecture that
a deterministic model would suffice to capture the true deterioration process.
To this end, probabilistic models have been developed to model the deterioration process of
infrastructure and building components (Marcous et al. (2002), Zhang and Vidakovic (2005)).
These models are in the form of Markov chains, which are attractive models widely used for
dynamical systems whose transition between different states through time is considered to be
random. In the model construction, a finite set of health conditions can be considered for the
system of interest, and the deterioration is parameterized by the transition probabilities between
the health states.
3.6.1 Data-driven MaxEnt Estimation of Deterioration Markov Chain Models
As the example, let us consider the random deterioration for concrete bridge decks protected
with rigid overlays (Marcous et al. (2002)). The subjective assessment of transitions have been
acquired by experts resulting in the following transition matrix between six different damage
conditions
P =
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
0:71 0:27 0:02 0 0 0
0 0:93 0:07 0 0 0
0 0 0:78 0:13 0:09 0
0 0 0 0:92 0:06 0:02
0 0 0 0 0:58 0:42
0 0 0 0 0 1
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
(3.28)
We treat these transition rates as the mean values. Since no information has been provided on
the second order statistics or the scatter in the experts’ opinions, we assumed a 30% coefficient
40
of variation for these values. Figs. 3.12 show the estimated MaxEnt pdf’s and joint pdf’s for the
independent transition rates.
p
11
p
12
0.5 0.6 0.7 0.8 0.9
0
0.1
0.2
0.3
0.4
0.5
(a)
0 0.2 0.4 0.6 0.8 1
0
5
10
15
20
p
23
pdf
(b)
p
33
p
34
0.5 0.6 0.7 0.8 0.9
0
0.1
0.2
0.3
0.4
0.5
(c)
p
45
p
46
0 0.05 0.1
0
0.01
0.02
0.03
0.04
0.05
0.06
(d)
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
3
p
56
pdf
(e)
Figure 3.12: The pdf’s and joint pdf’s corresponding to the independent random transition prob-
abilities outgoing states 1 to 5, plotted in figures (a) to (e), respectively.
3.6.2 Bayesian Update of the RTM
The proposed MaxEnt probability density function for the random transition matrices can be
incorporated in a Bayesian framework as the prior density function, based on which a posterior
density is estimated using the observation on the damage conditions of the infrastructures. Let us
consider a case where a Markov chain with the MaxEnt RTM has been calibrated for an ensemble
41
of infrastructure systems of a particular type. Letd denote the array of observations that has been
made available. According to Bayes’ theorem one can then write
^
f
P
(Pjd)/f
d
(djP )f

P
(P ) (3.29)
where f

P
(P ) is the MaxEnt probability density function for transition matrix P and Pr(djP )
is the likelihood function for observationd. Let us assume that three transitions were observed
during the experimentations, which ares
0
1
tos
00
1
,s
0
2
tos
00
2
ands
0
3
tos
00
3
. Then, one can write
^
f
P
(Pjd)/ [P ]
s
0
1
s
00
1
[P ]
s
0
2
s
00
2
[P ]
s
0
3
s
00
3
f

P
(P ) (3.30)
In order to investigate the impact of new observation data on these bridges on the posterior
update of the probability density functions of this RTM, let us consider two different hypothet-
ical cases. In the first case, let us consider that through the experimentations, three transitions
from state 4 to state 5 have been observed, whereas in the second case, observation include one
transition from state 4 to state 5 and two transitions from state 4 to state 6. As can be seen in
Fig. 3.13, the observation of the first case has resulted in a shifted concentration of probability
density for higher values ofp
45
which is consisted with what one expect based on the observation
data. However, no shift is apparent in thep
46
direction. In other words, the probability density of
the random variablep
46
was not informed given such an observation. In the second case, however,
one can see that both transitionsp
45
andp
46
have been informed.
42
p
45
p
46
0 0.05 0.1
0
0.01
0.02
0.03
0.04
0.05
0.06
(a)
p
45
p
46
0 0.05 0.1
0
0.01
0.02
0.03
0.04
0.05
0.06
(b)
Figure 3.13: The posterior joint pdf’s corresponding to the independent random transition prob-
abilities outgoing state 4 to states 5 and 6. Fig (a) refers to the case where three transitions from
state 4 to state 5 have been observed. Fig (b) corresponds to the case where one transition from
state 4 to state 5 and two transitions from state 4 to state 6 are observed.
43
Chapter 4
Stochastic Markov Chain Model for Power
Demand
4.1 Introduction and Literature Review
One of the capacities of the future Smart Grid to be realized in the demand side management. It
consists in managing the grid and improving its condition by modifying the electricity demand,
as opposed to the way conventional protocols were defined. Before the advent of the Smart
Grid technologies, in the face of an instability problem, the sole option for grid operators was to
increase the generation capacity. Such response needs to take effect fast enough before further
implications in the grid operation. Thus, a reserved generation is scheduled to spin at a minimal
speed so that it can be brought up to the desired generation level quite fast. In the future grid, a
large portion comes from dispatchable generation units, and in order to compensate for the load
variation, the operation should conventionally be more dependent on reserve units. The operation
of this reserved generation is costly and the fact that it may not be needed at all makes this option
uneconomical.
44
A more economical protocol aims at the modification of the electrical load, instead of increas-
ing the generation capacity. A major incentive in the modification of the loads is the electricity
price. If consumers realize the increase in the price during the peak hours and actually are forced
to pay the difference, they may reconsider the operation of certain appliances, such as dishwash-
ers, dryers, etc. Therefore, dynamic pricing can be an ideal framework for the demand-side
management. A major challenge, however, is to predict the ‘demand response’, i.e. the way con-
sumers react to the dynamic pricing. As a prerequisite to the demand response model, an accurate
high resolution model of the consumption behavior of the population under study should be con-
structed. Aside from the importance of such predictive models for demand side management, they
also serve as a vital component for other power system managements. The demand forecast at an
intraday resolution can greatly assist in maintenance, operations, and capital planning, as well as
in potential load shedding planning (Cottet and Smith (2003)).
There are various modeling approaches for the demand forecast (— (1980), Bunn (1982),
Gross and Galiana (1987), Mahmoud et al. (1981)). One of the most commonly used approach
is based on the time series analysis or system identification, including autoregressive (AR) and
autoregressive moving average (ARMA) processes (Hagan and Behr (1987),Moghram and Rah-
man (1989)). In order to account for the nonlinearities in the load parameters, support vector
machines has been proposed to solve the nonlinear regression problem of load forecast (Pai and
Hong (2005)). Neural networks, as another artificial intelligence approach, have also been applied
to the load forecast problems (Hippert et al. (2001, 2005)) Kalman filter has also been proposed
to predict electrical load which involves the state estimation of the load dynamics (Moghram and
45
Rahman (1989)). In all these approaches, the model fails to provide insights and produce quan-
titative measure as to how the total demand depends on influential underlying factors such as the
activity patterns of the users.
Our objective in this study is to form a bottom-up predictive model that relies on parameters
that have physically meaningful interpretation. Specifically, we aim to incorporate the activ-
ity patterns of the consumers as the influential factor for our model construction. The motiva-
tion for such formalism is two-fold. First, the activities of the consumers are indeed the actual
driver for their consumption behavior. Second, with the aid of new smart phones and gadgets,
the activity patterns of individuals have become increasingly more observable, providing an un-
precedented opportunity for the validation and parameter estimation for activity-based demand
prediction models.
In this chapter, after the discussion on different approaches to include the electricity load
in the power system equations, in Section 4.3, we explain the Markov chain model to be used
for demand prediction. This model has been extended to represent a heterogeneous population
of consumers. The theoretical development for such extension is included in Section3, where
the behavior of Markov chains associated with the uncertainties in its parameters, i.e. transition
probabilities, is described using a probabilistic model for the transition matrices. The principle
of Maximum Entropy (MaxEnt) is used to characterize the probability measure of the transition
rates. Section3 also includes the application of the proposed model to a disease progression
problem and deterioration of infrastructure systems. Finally, in Section 4.4, we illustrate the
application of this model to the residential electricity demand prediction problem based on an
actual time-use dataset obtained in 2011.
46
4.2 Electric Load Models in Power System Equations
Constant-impedance Load Model, Internal Generator Bus A widely used approach is to
consider loads to be constant impedance and therefore represent the loads by equivalent admit-
tance values. So, for a system of n generators and m loads, an enlarged (n +m) (n +m)
admittance matrix, denoted by

Y is constructed. This global admittance matrix is then reduced so
that the resulting admittance matrix corresponds to the generator nodes only. The reduction tech-
nique is called the Kron reduction. We denote byI
G
2C
n
andV
G
2C
n
the vector of currents
and voltages at the generator nodes and letV
L
2C
m
denote the vector voltage values at the load
buses. Since the current injection at load buses are zero, then we can represent the current-voltage
relation in the following block-matrix form
2
6
6
4
I
G
0
3
7
7
5
=
2
6
6
4
Y
GG
Y
GL
Y
LG
Y
LL
3
7
7
5
2
6
6
4
V
G
V
L
3
7
7
5
: (4.1)
We can eliminateV
L
and derive the reduced admittance matrix (annn matrix) as

Y
r
=Y
GG
Y
GL
Y
1
LL
Y
LG
: (4.2)
This equivalent reduced matrix is then used in the state equation (Eq. 2.19). This way, we have
captured the effect of voltage interdependence through the use of equivalent admittance matrix.
The final set governing equations will only involve differential equations with no algebraic part.
Constant-impedance Load Model, Jacobian of Load Flow Another approach is to preserve
the power flow model (i.e. Eqs. 2.4). Let x denote the vector of perturbed state variables, z
47
the vector of perturbed algebraic variables of the generators (voltage magnitudes and currents),
and w the vector of perturbed algebraic variables of the load buses, i.e. voltage magnitudes and
angles. Consistent with the notation in Eq. 2.14, y = [z; w]
T
. The block matrix representation of
the linearization of Eqs. 2.14 is as follows
2
6
6
6
6
6
6
4
_ x
0
0
3
7
7
7
7
7
7
5
=
2
6
6
6
6
6
6
4
A B
C
D
11
D
12
D
21
D
22
3
7
7
7
7
7
7
5
2
6
6
6
6
6
6
4
x
y
w
3
7
7
7
7
7
7
5
(4.3)
whereD
22
is in fact the Jacobian of the network power flow model, denoted byJ
fl
. The Jacobian
of the system can then be given by
J
sys
=ABD
1
C (4.4)
whereD is given by
D =
2
6
6
4
D
11
D
12
D
21
J
fl
3
7
7
5
(4.5)
Another alternative is to take the output power at the generator as the coupling variable, in-
stead of the use of rotor angles in the above formulation (Cardell and Ilic (2004)). The use of
output power as the coupling variable is more interesting to the participants who want to track the
dynamics of the interconnected generators. LetP
G
2 R
n
denote the vector of output powers at
48
the generators. The state equation of this variable in terms of load disturbances,
_
P
L
2 R
m
and
generator frequencies!
G
2R
n
by
_
P
G
=K
P
!
G
+D
P
_
P
L
; (4.6)
where thenn matrixK
P
and thenm matrixD
P
are derived from linearization of the load
flow equations, based on
K
P
= J
GG
+J
GL
C
!
(4.7)
D
P
= J
GL
J
1
LL
(4.8)
C
!
= J
1
LL
J
LG
; (4.9)
whereJ
ij
=
@P
N
i
@
j
, and subscriptsG andL denote generation and load nodes, respectively.P
N
i
is
the load flow constraint written for nodei. With the matricesK
P
,D
p
andC
!
a network is fully
characterized.
Other Load Models In general, the real power drawn by a load bus is a nonlinear function of
the voltage and frequency. There are load models that are based on the static voltage characteristic
of the load. They essentially relate the power at a load bus to its voltage magnitude. Polynomial
or exponential relations can be considered (see Wang et al. (2009) - page 501). Also, there are
models that assuming constant voltage at the loads relate the real power to the derivative of the
voltage angle
i
at load busi. In this case the voltage angles are the new state variables (Bergen
and Hill (1981)).
49
4.3 Markov Chain Model for Power Demand Based on Acitivity Data
In order to model domestic electricity consumption, several bottom-up high resolution models
have been proposed (Capasso et al. (1994), Paatero and Lund (2006), Swan and Ugursal (2009)).
A major challenge in using such models is that their mathematical forms should be simple enough
to be computationally tractable. Another drawback is that they rely on a large number of assump-
tions and further depend heavily on the availability of detailed data, which is usually not the case.
Among the high resolution models for residential electricity consumption, there are models
that incorporate the behavioral information of the consumers by using the time use (TU) data (see,
e.g. Richardson et al. (2008, 2010), Widen and Wackelgard (2010), Widen et al. (2009)). TU
data can serve as a highly realistic estimator of the behavioral factor in the consumption pattern.
In these models, the TU data is used for the estimation of the parameters of a Markov chain
representing the activity dynamics of a group of surveyed individuals. These parameters are the
transition probabilities which are deterministically estimated based on the frequency of observed
transitions (more details in the Chapter 3). As will be discussed in next section, due to inherent
variability of the population or insufficient data, deterministic estimation of these parameters may
result in the infidelity of the Markov chain model. In next section, we propose a Markov chain
model with random transition probabilities, which will account for the uncertainties induced by
the missing information or the heterogeneity and variability of the consumer populations.
The data used in our study is the 2011 American Time Use Survey (ATUS) data provided
by the Bureau of Labor Statistics of the US Department of Labor. It is based on the survey of
12480 individuals about their daily activities. Specifically, the time and duration of each activity
are recorded. A detailed description of the activity is declared by the respondents based on a
50
hierarchical indexing of the activities. The time resolution of the activity reports is 1 minute. The
raw data is obtained for a time period between 4 am and 12 am. The ATUS data has been available
since 2003.
The first step in analyzing the data is the choice of the resolution in the activity space. It
usually involves clustering the large number of detailed activity identifiers into a small number of
activity groups. In our study, since the electricity consumption of the individuals are concerned,
the status of the consumer is divided into 4 states with different consumption levels; (1) away,
(2) sleeping, (3) home with high consumption (washing, cooking, etc.) and (4) home with low
consumption. Figs. 4.1 show sample activity trajectories for 6 different respondents.
51
4 am 9 am 2 pm 7 pm 12 am
0
1
2
3
4
5
Activity index
Time
(a)
4 am 9 am 2 pm 7 pm 12 am
0
1
2
3
4
5
Activity index
Time
(b)
4 am 9 am 2 pm 7 pm 12 am
0
1
2
3
4
5
Activity index
Time
(c)
4 am 9 am 2 pm 7 pm 12 am
0
1
2
3
4
5
Activity index
Time
(d)
4 am 9 am 2 pm 7 pm 12 am
0
1
2
3
4
5
Activity index
Time
(e)
4 am 9 am 2 pm 7 pm 12 am
0
1
2
3
4
5
Activity index
Time
(f)
Figure 4.1: Activity profiles of 6 representative respondents during the time period 4 am to 12 am.
The activity indices are defined as follows (1) away, (2) sleeping, (3) home with high consumption
(washing, cooking, etc.) and (4) home with low consumption.
52
The transition probabilities are the obtained based on the observed transition frequencies,
as described in the previous section, over a time period between 6 a.m. and midnight. The
observation period should be large enough to reveal realistic transitions. In order to also account
for the diurnal changes of the transition probabilities, four different time zones are considered;
(1) 6 a.m. to 11 a.m., (2) 11 a.m. to 16 p.m., (3) 16 p.m. to 20 p.m. and (4) 20 p.m. to 12 a.m.
Our objective is to obtain transition matrices over 10 min intervals. Since the time resolution of
the input data is 1 minute, the estimated transition matrix should be raised to the power of 10 to
give the 10-min transition matrix. The resulting average transition matrices based on the first half
of the respondents are as follows
P
(1)
611
=
2
6
6
6
6
6
6
6
6
6
6
4
0:8719 0:0083 0:0896 0:0303
0:0173 0:8795 0:0928 0:0104
0:1377 0:0091 0:7821 0:0711
0:1211 0:0070 0:0945 0:7774
3
7
7
7
7
7
7
7
7
7
7
5
; P
(1)
1116
=
2
6
6
6
6
6
6
6
6
6
6
4
0:8548 0:0057 0:0949 0:0446
0:0276 0:8887 0:0626 0:0211
0:1544 0:0115 0:7637 0:0705
0:1032 0:0094 0:0718 0:8156
3
7
7
7
7
7
7
7
7
7
7
5
P
(1)
1620
=
2
6
6
6
6
6
6
6
6
6
6
4
0:7794 0:0044 0:1444 0:0718
0:0172 0:8867 0:0626 0:0334
0:0851 0:0115 0:7878 0:1156
0:0463 0:0061 0:0659 0:8817
3
7
7
7
7
7
7
7
7
7
7
5
; P
(1)
2024
=
2
6
6
6
6
6
6
6
6
6
6
4
0:6519 0:0760 0:1421 0:1300
0:0019 0:9956 0:0017 0:0009
0:0496 0:1797 0:6159 0:1549
0:0205 0:1010 0:0488 0:8298
3
7
7
7
7
7
7
7
7
7
7
5
As is evident from the values of the transition probabilities, the activity pattern in the four time
zones are distinguished from one another, as expected. Once should also expect that for the same
time zone, the discrepancies in the estimates due to different choice respondent groups should be
53
notably smaller than the inter-zone differences. In order to verify this, a second set of transition
matrix estimates based on the second half of the respondents were calculated and are as follows.
P
(2)
611
=
2
6
6
6
6
6
6
6
6
6
6
4
0:8892 0:0054 0:0823 0:0232
0:0185 0:8807 0:0895 0:0113
0:1332 0:0087 0:7966 0:0615
0:1325 0:0075 0:0934 0:7666
3
7
7
7
7
7
7
7
7
7
7
5
; P
(2)
1116
=
2
6
6
6
6
6
6
6
6
6
6
4
0:8816 0:0042 0:0815 0:0326
0:0265 0:8899 0:0590 0:0246
0:1504 0:0121 0:7765 0:0610
0:0902 0:0088 0:0694 0:8315
3
7
7
7
7
7
7
7
7
7
7
5
P
(2)
1620
=
2
6
6
6
6
6
6
6
6
6
6
4
0:7955 0:0046 0:1320 0:0679
0:0207 0:8769 0:0630 0:0394
0:0922 0:0123 0:7940 0:1015
0:0543 0:0073 0:0680 0:8704
3
7
7
7
7
7
7
7
7
7
7
5
; P
(2)
2024
=
2
6
6
6
6
6
6
6
6
6
6
4
0:6709 0:0720 0:1363 0:1208
0:0007 0:9965 0:0019 0:0009
0:0486 0:1633 0:6264 0:1617
0:0216 0:0945 0:0530 0:8310
3
7
7
7
7
7
7
7
7
7
7
5
In order to quantify the discrepancy, for the two matricesP
1
;P
2
2R
nn
, let
(P
1
;P
2
)
denote
the R
n
2
-norm of the difference between the two matrices. Tables 4.1 and 4.2 which tabulate
the values of this metric between different pairs of transition matrices verify the aforementioned
argument on the distinct characters of the activity patterns in different time zones.
Table 4.1: The discrepency metric
(P
1
;P
2
)
between the zones for the first set of respondents
P
(1)
611
P
(1)
1116
P
(1)
1620
P
(1)
2024
P
(1)
611
0.000 0.0548 0.1736 0.3658
P
(1)
1116
0.000 0.1475 0.3450
P
(1)
1620
sym 0.000 0.2981
P
(1)
2024
0.000
54
Table 4.2: The discrepency metric
(P
1
;P
2
)
between the zones and between the first and second
sets of respondents
P
(2)
611
P
(2)
1116
P
(2)
1620
P
(2)
2024
P
(1)
611
0.0253 0.0716 0.1500 0.3486
P
(1)
1116
0.0366 0.1235 0.3273
P
(1)
1620
sym 0.0251 0.2802
P
(1)
2024
0.0237
Given the estimated transition matrix based on half of the data, the Markov chain model can
be validated against the full array recorded activity statuses. Specifically, the activity status at 9
am is considered as the initial condition. The activities at 10 minute intervals in the next hour are
estimated by the Markov chain model. Figs. 4.2 compare these predicted values versus the actual
activity records of all 12480 respondents.
55
9:10 9:20 9:30 9:40 9:50 10:00
0.2
0.25
0.3
0.35
π
1
Time
(a)
9:10 9:20 9:30 9:40 9:50 10:00
0.2
0.3
0.4
0.5
0.6
0.7
0.8
π
2
Time
(b)
9:10 9:20 9:30 9:40 9:50 10:00
0.2
0.25
0.3
0.35
0.4
π
3
Time
(c)
9:10 9:20 9:30 9:40 9:50 10:00
0
0.05
0.1
0.15
0.2
π
4
Time
(d)
Figure 4.2: The dynamics of the activity pattern of the consumers. Y axes correspond to the ratio
of the respondents in each activity status. The solid lines correspond to the ratios calculated from
the full data. The dashed lines refer to the Markov chain prediction where its transition matrix is
estimated based on half of the data.
Given the above mean values, a MaxEnt probability measure for the random transition ma-
trices is obtained by incorporating the mean values observed from half of the respondents. Here,
we illustrate the implication of such uncertainty representation on the activity condition of the
consumers. Specifically, in the first time zone, we start with the deterministically known status of
the population at 9 am, and we investigate the uncertainty in the prediction for the activities at 11
am. The QoI is the ratio of the population in each one of the four activity states. Figs 4.3 depict
56
the histograms of the random QoI. As can be seen, due to the nonlinearity in the Markov model,
the average of the stochastic QoI does not coincide with the QoI obtained from the deterministic
model, a point that was elaborated on in Chapter 3. The impact of the probabilistic behavior of
the consumers given such uncertainty representation on the grid performance metrics is discussed
in Chapter 6.
0 0.2 0.4 0.6 0.8 1
0
200
400
600
800
1000
1200
π
1
(ω)
(a)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
0
100
200
300
400
500
600
700
800
π
2
(ω)
(b)
0 0.2 0.4 0.6 0.8 1
0
200
400
600
800
1000
1200
π
3
(ω)
(c)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0
200
400
600
800
1000
1200
1400
1600
1800
π
4
(ω)
(d)
Figure 4.3: The histograms of the stochastic QoI for the activity prediction problem. The QoI are
the ratios of the population in each activity state. Circle marker refers to the average value of this
random quantity obtained by the stochastic model. Cross marker refers to the QoI obtained by
using the deterministic transition model.
57
4.4 Uncertainty Quantification for Markov Chain Model of Power
Demand
As mentioned in the previous section, in order to account for the uncertainties due to the missing
information, and more importantly, in order to construct a global activity-based demand model
which represents an ensemble of heterogeneous users, we need to consider variability in the pa-
rameters of the Markov chain model described in Section 4.3.
In order to illustrate the application of the Markov chain with MaxEnt random transition
matrix, a heterogeneous multi-agent system is considered. Using a synthetic example, we seek to
show that with the aid of only few samples from the agents dynamics we can obtain sufficiently
accurate results, in terms of predicting QoI’s, with respect to the result obtained from the full
system. To this end, let us consider a system that consists of 900 heterogeneous agents with
different transition properties. Let us assume that the mean values of the set of 900 transition
matrices is given by
P =
2
6
6
6
6
6
6
4
0:92 0:04 0:04
0:58 0:24 0:18
0:23 0:23 0:54
3
7
7
7
7
7
7
5
: (4.10)
As discussed earlier, one would obtain different sampled trajectories even from a fixed (de-
terministic) transition matrix. In order to induce heterogeneity among the agents, we consider
a transition matrix for each agent slightly perturbed around the nomial values given above. The
perturbations are assumed to follow a uniform distribution in a range centered at the mean value
with a size equal 10% of the mean value. For each agent, then, we sample a transition matrix
58
0 0.2 0.4 0.6 0.8 1
0
5
10
15
20
25
π
pdf

π
1
π
2
π
3
Figure 4.4: Comparison between the pdf’s of QoI (taken to be the three asymptotic distributions).
The solid line corresponds to MaxEnt RTM estimation with mean values and standard deviation
(based on the observation of only 30 samples). The dashed line corresponds to the 900 sample
trajectories drawn from a perturbed transition matrix. The perturbation is assumed to be uniformly
distributed around the mean values.
and then synthesize a trajectory of length 300. Based on these large array of trajectories, we can
calculate the “exact” QoI’s, such as asymptotic distribution, first passage times, etc.
In order to evaluate the MaxEnt prediction, we observe the trajectory of only 30 agents based
on which we calculate the mean and standard deviation values for transition frequencies. A Max-
Ent RTM is then estimated based on this information. In order to validate the MaxEnt estimation,
the asymptotic distribution corresponding to the three states (i.e.
1
1
,
1
2
and
1
3
) obtained from
the MaxEnt RTM is compared with those obtainedd from the “exact” one. Fig. 4.4 shows that
the pdf of the random QoI estimated by the MaxEnt distribution is in good agreement with that
obtained from large-size sampling of the underlying system.
59
Chapter 5
Stochastic Model for Optimally Scheduled
Generation
5.1 Introduction and Literature Review
The advent of the Smart Grid technologies will provide an unparalleled opportunity for energy
efficiency and sustainability by integrating more renewable energy resources and providing the
capacity for distributed generation. Such improvements, however, will materialize at the expense
of exposing the electric grid to major sources of uncertainties and complexities on the genera-
tion side. The scope of this chapter is the study of the uncertainties induced by the distributed
management and optimization of the generation.
The future Smart Grid will integrate communication and control network, including compo-
nents such as sensors, distributed controllers, and actuators, introducing a large-scale networked
control system. There will be a shift in paradigm from central control architecture to the dis-
tributed sensor design. The reason for this shift is two-fold. First, in future Smart Grid, one will
60
deal with a variety of physical grid topology as well as the communication topology. The intro-
duction of the “Plug-and-Play” technology will further complicate the identification of the true
topology. Second, in such large systems, the central control unit will need a high level of connec-
tivity with the distributed actuators; hence the computational burden and the escalated severity of
the consequences of a loss of link or a cyber attack. Distributed control, on the other hand, proves
more robust to the connectivity losses.
Several central-unit optimization methods have been proposed for the economic dispatch
problem. Among them, the gradient search method is one the conventional techniques adopted to
this end when convex cost functions are concerned (Wood and Wollenberg (1996)). Non-convex
cost functions have required development of techniques, such as Genetic Algorithms (Bakirtzis
et al. (1994)), the particle swarm (Gaing (2003)). As mentioned earlier, limitations due to the
centrally controlled architecture of these frameworks have motivated research efforts towards de-
velopment of the distributed optimization techniques, such as consensus framework.
Consensus problems have been widely studied and used in computer science and distributed
computing. Its origin dates back to 1960s, where they first appear in management science and
statistics (DeGroot (1974)). Its wide use in systems and control theory and its application to
networked systems was due to the pioneering work of (Borkar and Varaiya (1982), Tsitsiklis et al.
(1986)), which was followed by applications to the fields ranging from animals science Conradt
and Roper (2005), Olfati-Saber (2006) to vehicle formation problem (Olfati-Saber et al. (2007))
and sensor fusion (Gupta et al. (2005), Olfati-Saber and Shamma (2005)).
Solving a consensus problem, or implementing a consensus protocol, for a network of indi-
viduals means finding an agreement state based on a certain quantity of interest for the population.
61
It involves the information exchange between the individuals and their local neighbors in the net-
work, and alleviates the need for a central processing and control unit. The theoretical framework
for the consensus problems has been studied in Saber and Murray (2003) and Olfati-Saber and
Murray (2004).
The application of the consensus problem for the economic dispatch of the generators in a
Smart Grid was first introduced in Zhang and Chow (2011). It involves the reformulation of the
optimization problem that is solved at the central control unit. The objective is to find the genera-
tion levels at different units that result in the minimum generator costs, subject to the constraint of
demand satisfaction. In the central control paradigm, control signals are communicated between
the generators and the central unit. Such optimization problem can be recast as a distributed prob-
lem by the implementation of the consensus protocol, as will be described in Section 5.5. This
chapter is organized as follows. Section 5.2 includes basics of the graph theory, which sets the
foundation for the introduction of the consensus over networks in Section 5.3. A stochastic formu-
lation of the consensus problem is introduced in Section 5.4, which can address the uncertainties
induced by the uncertain and switching nature of the communication networks.
5.2 Graph Theory Notations and Definitions
We denote a directed graph (digraph) byG(V;E), where V =f1; 2; ;ng is a finite set of
vertices, and E2 V V is a set of edges, which consists of ordered pairs of vertices. In an
undirected graph,E is a set of unordered pairs of vertices. A directed (undirected) graph is called
strongly connected if there exists a directed path (a path) from every vertex to every other one.
62
The adjacency matrix of a weighted graph is the matrixA = [a
ij
], witha
ij
> 0, if (i;j)2 E,
anda
ij
= 0, otherwise. In an unweighted grapha
ij
= 1, if (i;j)2E, anda
ij
= 0, otherwise.
The Laplacian matrix of a weighted directed graph is the matrixL = [l
ij
], withl
ij
given by
l
ij
=
8
>
>
<
>
>
:
P
n
k=1;k6=i
a
ik
j =i
a
ij
j6=i
(5.1)
According to Greshgorin theoretm, the eigenvalues of the Laplacian matrix falls within a close
disk centered at in the complex plane with a radius equal to the graph’s maximum degre. The
Laplacian of undirected graphs is symmetric positive semi-definite, and we have
0 =
1

2

2n
2: (5.2)
The second smallest eigenvalue of Laplacian,
2
is called the algebraic connectivity of the asso-
ciated graph and determines the convergence rate of the consensus algorithms.
The digraph of an nn real matrix M = [m
ij
], denoted by (M) is the directed graph
formed on then vertices inV =f1; ;ng, in which there exists a directed edge the verticesi
andj if and only ifm
ij
> 0. We call a matrix (and a vector) nonnegative if all its components
are nonnegative. According to the Perron-Frobenius Theorem (Luenberger (1979)), if the digraph
of a nonnegative matrixM is strongly connected (i.e. M is nonnegative and irreducible), then
its spectral radius,(M), which is the Perron-Frobenius eigenvalue, is a simple eigenvalue ofM
associated with its only positive eigenvector.
63
A matrix whose row sum to 1 is called a (row) stochastic matrix. If a stochastic matrix P
satisfiesP
m
> 0 for some positivem, i.e.P is primitive, then we have the asymptotic result that
lim
k!1
P
k
= 1v
T
, wherev is a positive column vector that also satisfies 1
T
v = 1.
5.3 Discrete-time Consensus Problem
Consider the graphG with LaplacianL. Letx = [x
1
; ;x
n
]
T
denote the information encoded
on each vertex on the graph. The consesus problem consists in the dynamics of this informa-
tion via consecutive local interactions between the vertices until they reach a convergence. The
continuous-time consensus problem is governed by the following first order dynamical system
_ x
i
(t) =
n
X
j=1;j6=i
a
ij
(x
j
(t)x
i
(t)) i = 1; ;n (5.3)
or, in compact form, by
_ x(t) =Lx(t) (5.4)
The discrete-time form of the consensus problem is govern by the following iteration written for
each vertex (Olfati-Saber et al. (2007))
x
i
(k + 1) =x
i
(k) +
n
X
j=1;j6=i
a
ij
(x
j
(k)x
i
(k)) i = 1; ;n (5.5)
which can be rewritten as
x(k + 1) =P

x(k); (5.6)
64
with
P

=IL (5.7)
whereI is the identity matrix, is the time step. The matrixP

is refered to as Perron matrix
of graphG with parameter. It can be shown (Olfati-Saber et al. (2007)) that if 0 < < 1=,
where is the maximum degree of graphG, thenP

is a nonnegative row-stochastic matrix, i.e.
n
X
j=i
[P

]
ij
= 1; i = 1; ;n: (5.8)
Furthermore, ifG is alos strongly connected, thenP is primitive, and it can be shown (Olfati-Saber
et al. (2007)) that a consensus can be asymptotically reached given any initial statex(0).
5.4 Consensus over Random Graphs
We are interested in the consensus behavior in a dynamic (switching) communication network,
where the local interagent interactions are subject to change. The proposed formualtion could
also represent cases where the communication links can become unavailable or may be subject to
noise or delays. We make use of probability theory in representing these variabilities.
5.4.1 Random Consensus Dynamics
The consensus dynamics is governed by the graph Laplacian,L, which determines the state of the
communication network. Therefore, our objective is to characterize a probability measure for the
65
graph Laplacian in order to incorporate the uncertainties in the communication network. Based
upon this probabilistic representation for matrixL, the random Perron matrix is defined by
P

(!) =IL(!); (5.9)
and consequently the probabilistic description of the consensus dynamics can be estimated.
The mathematical setting for our formulation follows the setting used in the theory of Markov
chains in random environment first introduced in (Cogburn (1984)). LetL denote the set of all
admissible Laplacian matrices for a directed graph, i.e.
L =fL2R
nn
: L1 = 0g (5.10)
whereR
nn
is the set of allnn matrices with real components. In the case of an undirected
graph, one deals withL
sym
defined as follows
L
sym
=fL2R
nn
sym
: L1 = 0g (5.11)
whereR
nn
sym
is the set of allnn symmetrix matrices with real components. Let (L;B
0
) denote
theB
0
-measurable set of Laplacian matrices at time 0, and let
=L
Z
denote the event space.
LetL
t
:
!L be thet-th realization of random graph Laplacian. We can alternatively write
L(!
t
) :=L
t
(!). LetB =B
1
1
denote the-field generated byL
t
between times1 and +1.
We can then construct the probability space (
;B;), with being the probability measure.
Following these notations, the random graph Laplacian, denoted byL(!
t
) is aL-valued ran-
dom matrix whoseij-th component determines the status of the interconnections between nodes
66
i and j, realized based on the sample !
t
drawn from the sample space
at time t. Given the
instance of communication network at timet, the corresponding Perron matrix will be
P

(!
t
) =IL(!
t
): (5.12)
Therefore, using the notation!
ij
:= [!
i
; ;!
j
], the random consensus variable at timek +s
evolving from the consensus state at timek can be estimated using the following equation
x
k+s
(!
1s
) = [P

(!
k+s
)P

(!
k+2
)P

(!
k+1
)]x
k
(!
1k
) (5.13)
5.4.2 Probabilistic Characterization ofP

(!)
Undirected graphs with unconstrained communication
We use the Maximum Entropy (MaxEnt) principle in order to estimate the probability measure
for the random Laplacian
L(!) =
2
6
6
6
6
6
6
6
6
6
6
4
l
11
(!) l
12
(!) l
1n
(!)
l
22
(!) l
2n
(!)
sym
.
.
.
.
.
.
l
nn
(!)
3
7
7
7
7
7
7
7
7
7
7
5
By definition (Eq. 5.1), the diagonal components of Laplacian matrices depend on the non-
diagonal at their respective row as follows
l
ii
(!) =
n1
X
j=1;j6=i
l
ij
(!) a:s: (5.14)
67
The non-diagonal components on our definition of graph Laplacian take values in [0; 1], unlike
the case for conventional Laplacians where non-diagonal components are constrained to be either
zero or one. In fact, the componentl
ij
(!) is a measure of the communication strength between
the agentsi andj.
Note that L(!) takes value in the space of symmetric matrices, and the diagonal terms are
dependent random variables according to Eq. 5.14. We further assume that non-diagonal compo-
nents are independent from each other, that is the existance/failure of a communication between a
given pair is independent of that between any other pair. Therefore, there are onlyn(n 1)=2 in-
dependent random variables in annn random matrix. If we characterize the probability measure
for these independent random components, with the use of symmetery and zero-sum property, we
have practically characterized the probability measure of the whole matrix. The resulting random
matrix will then have symmetric and zero-sum properties almost surely.
Let us assume a number of summary statistics are available from observing the communica-
tion topology over a given time period. Let us assume, without loss of generality, that at each time
instance, the status of the interagent communications is represented by a non-weighted digraph,
with the adjacency matrixA
t
= [a
t
ij
]. During an observation time window of lengthT , we can
therefore determine the first and second statistical moments forl
ij
(!) according to

l
ij
=
1
T
T
X
t=1
a
t
ij
(5.15)
68
Now, let us divide the total observation period into aW observation windows, and let

l
w
ij
denote
the mean values obtained in time windoww. One can then compute the second order statistics as
a measure of variability in the data. Specifically, the standard deviation forl
ij
(!) can be given by

2
ij
=
1
W
W
X
w=1
(

l
w
ij

l
ij
)
2
(5.16)
We seek to estimate the probability measure for the componentl
ij
(!) given these statistics,
using MaxEnt principle. MaxEnt estimation is the optimal estimator for probability measures in
that it is least biased with respect to the missing information. In fact, it only relies on the available
information, such as statistical moments, and remains uncomittal to any other prior assumption.
For a random variable x(!) taking value in [0; 1], given its first and second statistical mo-
ments, x and, respectively, the MaxEnt pdf, denoted byf

x
(x), is one that solves the following
optimization problem
f

x
(x) = arg max
f

Z
1
0
f
x
(x) lnf
x
(x)dx
s.t.
Z
1
0
f
x
(x)dx = 1
Z
1
0
xf
x
(x)dx = x
Z
1
0
x
2
f
x
(x)dx =
2
+ x
2
(5.17)
69
Instead of using standard numerical algorithms to solve the above optimization problem, we
construct and consequently solve its dual. The Lagrangian of the MaxEnt optimization problem
can be constructed as
L(f
x
;;;) =
Z
1
0
f
x
(x) lnf
x
(x)dx
+ (
Z
1
0
f
x
(x)dx 1)
+ (
Z
1
0
x
i
f
x
(x)dx x)
+ (
Z
1
0
x
2
f
x
(x)dxs
2
x
2
)
(5.18)
The dual problem is then given by
sup
;;
D(;;) = sup
;;
inf
fx
L(f
x
;;;)
For the case where only knowledge of mean values are used, we won’t have the last constraint
in the Eqs. 5.22 and its associated Lagrange multipliers. The MaxEnt distributionf

in this case
belongs to the exponential family of the following form
f

x
(x) = exp( 1 +x) 1
[0;1]
(x)
where 1
[0;1]
is the indicator function for the set [0; 1]. If the standard deviations are also included,
the MaxEnt solution will have the following form
f

x
(x) = exp( 1 +x +x
2
) 1
[0;1]
(x):
70
Directed graphs with limited communication resources
The random Laplacian of a directed graph is not symmetric (i.e.l
ij
(!)6=l
ji
(!)) necassarily. Let
us assume that at each time instance, the number of communication links any agenti can afford
to establish is a known constant, denoted byU
i
. Therefore, we can write
n
X
j=1;j6=i
l
ij
(!) =U
i
a:s: i = 1; ;n: (5.19)
Let l
i
(!) taking values in [0; 1]
n
denote the array of ith row on the random Laplacian. The
admissible set for then 1 non-diagonal random entries on rowi can be formed as follows
V =
8
<
:
[l
i1
; ;l
i(i1)
;l
i(i+1)
; ;l
in
] [0; 1]
n1
:
n
X
j=1;j6=i
l
ij
=U
i
9
=
;
: (5.20)
On each row, there are only n 2 independent random variables and one can express the k
i
-
th component, l
ik
i
, as a function of the other n 2 components, and thus deal with a lower
dimensional admissible set defined by
V
0
=
8
<
:
l
0
i
:= [l
i1
; ;l
i(i1)
;l
i(i+1)
; ;l
i(k
i
1)
;l
i(k
i
+1)
; ;l
in
] [0; 1]
n2
:
n
X
j=1;j6=i;j6=k
i
l
ij
U
i
9
=
;
:
(5.21)
71
The MaxEnt probability density function for the random rowl
0
i
, denoted byf

l
0
i
, is then given by
f

l
0
i
= arg max
f
l
0
i

Z
V
0
f
l
0
i
(l
0
i
) lnf
l
0
i
(l
0
i
)dl
0
i
s:t:
Z
V
0
f
l
0
i
(l
0
i
)dl
0
i
= 1
Z
V
0
l
ij
f
l
0
i
(l
0
i
)dl
0
i
=

l
ij
j = 1; ;n;j6=i;j6=k
i
Z
V
0
l
2
ij
f
l
0
i
(l
0
i
)dl
0
i
=
2
ij
+

l
2
ij
j = 1; ;n;j6=i;j6=k
i
(5.22)
For the case where only knowledge of mean values are used, the MaxEnt distribution f
l
0
i
belongs to the exponential family of the following form
f

l
0
i
(l
0
i
) = exp
2
4

i
1 +
n
X
j=1;j6=i;j6=k
i

ij
l
ij
3
5
1
V
0(l
0
i
) i = 1; ;n;
where 1
V
0 is the indicator function. If the standard deviations are also included the MaxEnt
solution will have the following form
f

l
0
i
(l
0
i
) = exp
2
4

i
1 +
n
X
j=1;j6=i;j6=k
i

ij
l
ij
+
n
X
j=1;j6=i;j6=k
i

ij
l
2
ij
3
5
1
V
0(l
0
i
) i = 1; ;n:
5.5 Distributed Optimization of Power Generation
The economic dispatch of the generation in a multi-unit grid under the uncertainty in the commu-
nication links between the controllers is considered. The objective is to minimize to total cost of
gerenation subject to the constraint that the total demand should be met. The optimal condition
can be achieved through a decentralized communication protocol. In this approach, instead of a
central control units, the the local controllers at each generator through communication with one
72
another achieve the global optima. As a result, the schedule for each generation unit is gradually
adjusted based on their relative cost conditions with respect to the rest of the units, hence the
name Incremental Cost (IC) .
In what follows, the consensus formulation of the dynamics of this schedule adjustment pro-
posed in Zhang and Chow (2011) will be described. Let us consider a power grid withn genera-
tion units. LetC
i
(P
Gi
) denote the generation cost associated with theith unit, as a function of its
output power. The economic dispatch problem involves in finding the array of generation output
powersP

G
=fP

G1
; ;P

Gn
g that solves the following global optimization problem
P

G
= arg min
n
X
i=1
C
i
(P
Gi
)
s:t: P
D
=
n
X
i=1
P
Gi
;
(5.23)
whereP
D
is the power demand, andC
i
(P
Gi
) is the cost function for uniti assumed to be given
by
C
i
(P
Gi
) =
i
+
i
P
Gi
+
i
P
2
Gi
; (5.24)
where
i
;
i
and
i
are constant.
In the decentralized approach for the economic dispatch, one of the generation units is des-
ignated as the leader, which ensures that at each iteration the demand power is served. The
incremental cost at uniti, at timet, given its output powerP
t
Gi
is given by

t
i
=
dC
i
(P
t
Gi
)
dP
t
Gi
: (5.25)
73
The optimal solution is achieved if all the units arrive at the same incremental cost. The discrete-
time linear consensus formulation of this condition will therefore read

t+1
=P

t
; (5.26)
where P

is defined in Eq. 5.9. Based on these increments, new values for output powers are
obtained using Eq. 5.24. The adjusted power outputs may not meet the power demand, hence the
power imbalance given by
P
t+1
=P
d

n
X
i=1
P
t+1
Gi
: (5.27)
It is on the designated leader unit to compensate for this imbalance. To this end, the incremental
cost for the leader unit,
l
, is adjusted by the power imbalance as follows

k+1
l
=
k+1
l
+P
t+1
(5.28)
where is a time factor that determines the speed of convergence.
Let us now consider a dynamic communication topology, i.e. a case where the inter-unit links
are subject to unknown change over time. One can then adopt the random treatment of consensus
problem proposed earlier to arrive at a distributed optimization that accounts for the uncertainty in
74
the communication topology. The following reformulated equations will thus govern the random
dynamics
C
i
(P
Gi
(!)) =
i
+
i
P
Gi
(!) +
i
P
Gi
(!)
2

t+1
(!) =P

(!)
t
(!)
P
t+1
(!) =P
d

n
X
i=1
P
t+1
Gi
(!)

k+1
l
(!) =
k+1
l
(!) +P
t+1
(!)
(5.29)
The power system studied in this section is the 4-generation example used in Zhang and
Chow (2011). Initial conditions and parameters of the cost function corresponding to these units
are shown in Table 5.1. The total demand is assumed to be P
D
= 900MW . We assume the
communication links between the units are subject to unknown change. However, let us assume
a history of observation on the communication topology yeilds the following average adjacancy
matrix

A =
2
6
6
6
6
6
6
6
6
6
6
4
0 0:2 0:8 0:7
0:2 0 0:3 0:1
0:8 0:3 0 0:4
0:7 0:1 0:4 0
3
7
7
7
7
7
7
7
7
7
7
5
:
A maximum entropy probability density is obtained for the associatedP

(!). Given a time step
of size 0.1 sec, Figs. 5.1 show the random trajectories of the 4 incremental costs until conver-
gence. Figs. 5.2 depict the pdf’s of the random power dispatch incorporating the uncertainty in
the communication links at different times. These probability distributions provide quantitative
measure of the scatter in the status of the scheduled generation. Such quantified scatter can assist
75
Table 5.1: Initial conditions and parameters of the 4 generation units
Unit
i

i

i
P
Gi
(0)(MW )
1 561 7.92 0.001562 450
2 310 7.85 0.001940 300
3 78 7.80 0.004820 100
4 100 7.70 0.005120 50
grid operator and other entities account for the confidence in their prediction of the generation
and assist them with making robust decisions.
76
0 2 4 6 8 10
8.9
8.95
9
9.05
9.1
9.15
9.2
9.25
9.3
9.35
time (sec)
λ
1

mean
mean+std
mean−std
(a)
0 2 4 6 8 10
8.9
8.92
8.94
8.96
8.98
9
9.02
time (sec)
λ
2

mean
mean+std
mean−std
(b)
0 2 4 6 8 10
8.75
8.8
8.85
8.9
8.95
9
9.05
9.1
time (sec)
λ
3

mean
mean+std
mean−std
(c)
0 2 4 6 8 10
8.2
8.3
8.4
8.5
8.6
8.7
8.8
8.9
9
9.1
9.2
time (sec)
λ
4

mean
mean+std
mean−std
(d)
Figure 5.1: The random trajectories show the dynamics of the incremental costs in a random
communication topology governed by a MaxEnt probability measure.
77
250 260 270 280 290 300 310
0
0.01
0.02
0.03
0.04
0.05
0.06
P
G2
(MW) @t=0.5 sec
pdf
(a)
250 260 270 280 290 300 310
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
P
G2
(MW) @t=1 sec
pdf
(b)
250 260 270 280 290 300 310
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
P
G2
(MW) @t=2 sec
pdf
(c)
250 260 270 280 290 300 310
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
P
G2
(MW) @t=3 sec
pdf
(d)
250 260 270 280 290 300 310
0
0.2
0.4
0.6
0.8
1
1.2
1.4
P
G2
(MW) @t=4 sec
pdf
(e)
250 260 270 280 290 300 310
0
1
2
3
4
5
6
7
8
P
G2
(MW) @t=6 sec
pdf
(f)
Figure 5.2: The probabity density function of the random power output of unit 2 scheduled under
the uncertainty in the communication topology. As can be seen, the scatter around the mean
values progressively decreases, as the adjustment nears convergence.
78
Chapter 6
Stability of the Smart Grid under Uncertainty
6.1 Introduction
In this chapter, we investigate the impacts of the uncertainties modeled in Chapters 4 and 5 on the
performance metrics of the Smart Grid which can assist grid operators with their robust decision
making processes, and therefore close the loop for a comprehensive uncertainty management. To
this end, the small signal stability condition of the power grid is chosen as the quantitiy of interest.
The demand and supply uncertainties will impact the system behavior, and the quantification of
this impact is the subject of this chapter.
In what follows, based on the basics of the basics of power system theory discussed in Sec-
tion 2.1 , we formulate the small signal stability analysis of the power system. This analysis
involves solving a stochastic eigenvalue problem. Eigenvalue problems appear in a wide range of
applications and have been studied extensively in this dissertation. Specifically, we report on the
developed solution algorithms for the stochastic eigenvalue analysis for symmetric operators in
Section 6.2 and for non-symmetric operators in Section 6.3. A sample 9-bus power grid is chosen
in order to illustrate the uncertainty propagation proposed in this study.
79
If the impacts of the sources of uncertainty in the operation of the Smart Grid are not suit-
ably accounted for, aside from the failure in providing the decision makers with the quantitative
confidence information about the output analytics, the credibility of the analysis results can be
degraded due to the nonlinearities in the power system models. We seek to formulate the stochas-
tic small signal stability problem under the uncertainties in the load and generation, as described
in Chapters 4 and 5. As discussed earlier, small signal stability concerns the capability of the
system to remain in equilibrium when faced with small disturbance. It does not deal with the
characteristics of the disturbance, nor with the mechanical input of the generators. It is a function
of the operating condition and network properties.
The operating (or equilibrium) condition is indeed the solution of the power flow model.
Given the random active and reactive power at load buses, the random voltage magnitudes and
angles are computed which are then used to calculate the Jacobian of the state matrix. So, as op-
posed to the damping and inertia parameters,M andD that are present explicitly in the Jacobian,
the random load and network properties enter the problem non-explicitly through the rotor angles
and admittance matrixY .
In what follows, we describe the mathematical and algorithmic developements for the solution
of stochastic eigenvalue problems. We first describe the algorithms inspired by the power method
for the solution of the eigenvalues and eigenvectors of symmetric matrices. We will also discuss
a method for the eigen-analysis of non-symmetric matrices, which is the case for the stochastic
small signal stability problem which involves the non-symmetric Jacobian matrix of the DAE in
Eq. 2.14.
80
6.2 Stochastic Eigenvalue Problem - Symmetric Matrices
The robust predictive models for natural and engineered systems should incorporate the significant
variabilities in the behavior of these systems induced because of the inherent variability of system
components, inadequacy of the input-output models and inaccuracy of the numerical implemen-
tations (finite-dimensional approximation). Towards this end, efficient uncertainty quantification
techniques have been developed that represent the uncertainty in the input, propagate it through
the numerical models, and finally produce a probabilistic description for the system outputs and
Quantities of Interest (QoI). In addition to the classical sampling techniques, such as Monte Carlo
simulation and its improved derivatives, such as Markov chain Monte Carlo techniques, recently,
there has been growing interest in the non-sampling spectral technique, named Polynomial Chaos
Expansion (Debusschere et al. (2003), Ghanem and Spanos (2003), Maitre et al. (2001), Soize
and Ghanem (2004), Xiu and Karniadakis (2003)). In this stochastic projection technique, the
governing equation is projected onto truncated polynomial chaos coordinates, resulting in a finite-
dimensional approximation of the sample space, which is amenable to numerical calculation.
Our objective in this work is to characterize the eigenspace of linear systems with parametric
uncertainties. Such characterization can be used as the basis in various model reduction tech-
niques, and thus significantly improve the computational cost of the analysis of large uncertain
systems. Our particular interest is on random matrices that arise from finite-dimensional approx-
imation (discretization) of continuous operator, such as partial differential operators, involving
parametric uncertainty, as well as on those corresponding to a discrete random system, such as
multi-degree of freedom oscillators. Generally, there is no closed-form expression for these ran-
dom matrices.
81
Among the first solution approaches for the eigenvalue analysis of such matrices, the pertur-
bation technique was proposed in (Collins and Thomson (1969)), where the solution eigenpair
was calculated based on few low-order perturbations. This approach will give inaccurate results
in cases where higher-order perturbations are needed, or when the perturbations are not small
enough, disqualifying the Taylor series approximation. Statistical sampling is another approach
first proposed in (Shinozuka and Astill (1972)). It is computationally intensive, however, espe-
cially in cases where a prohibitively large number of samples of independent random variables
are required. A dimensional decomposition approach was proposed in (Rahman (2006)) for the
probabilistic characterization of the random eigenvalues using a lower-variate approximation of
these random variables which proves superior over the perturbation techniques.
Recently, Polynomial Chaos Expansion techniques have been applied to the solution of ran-
dom eigenvalue problem (Ghanem and Ghosh (2007), Ghosh and Ghanem (2008), Red-Horse and
Ghanem (1999)), where given a PC representation for the random matrix, the solution eigenvalue
and eigenpairs are assumed to have similar spectral representation in the same vector space. The
weak form of the eigenvalue problem is then obtained by Galerkin projection of the eigenvalue
equation on the PC orthogonal basis, spanning the Hilbert space of input random variables, lead-
ing to a system of nonlinear algebraic equations. The PC coefficient of the solution eigenpairs
are then obtained by the numerical solution of these systems. The resulting PC solution lends
itself readily to accurate uncertainty representation of QoI and sensitivity indices (Pellissetti and
Ghanem (2004)), and is greatly advantageous to other proposed methods.
In our approach, given a PC representation for the random matrix, instead of solving the
system of algebraic equations for the PC coefficients, we utilize the idea of power iteration scheme
in order to calculate the PC coefficients of dominant random eigenpair. Also, we propose an
82
algorithm that calculates a prescribed number of dominant random eigenpairs, and thus arrive at
an optimal approximation for a lower-dimensional subspace in which the random operator can be
represented.
In what follows, first in Section 6.2.1, we describe the mathematical setting and provide the
theoretical background for the proposed spectral approximation. Sections 6.2.2 includes the pro-
posed solution algorithms for the random eigenvalues and eigenvectors of a random matrix.
6.2.1 Problem Definition
Consider an experimental setting defined by the probability triplet (
;F;P). Let :
7! R
d
be anR
d
-valued random variable and letF

F be the-algebra induced by. Further, denote
the probability density of byd. We consider to refer to the source of uncertainty in an
underlying mathematical problem, such as random coefficients or boundary conditions in a partial
differential equation describing the physics of the problem. We will further assume that these
random variables are statistically independent, with the implication thatd =
Q
d
i=1
d
i
where

i
is the probability measure of the real-valued random variable
i
. LetA :
7! Mat
sym
nn
(R) be
a (F

;)-measurable matrix-valued random variable with values in Mat
sym
nn
(R), that is the set of
symmetric real-valued square matrices of sizen. We seek the solution of the following stochastic
eigenvalue problem
A()
s
() =
s
()
s
() s = 1; ;n a.s. (6.1)
83
where superscript s indexes the invariant subspace associated with the pair (
s
;
s
). Suppose
there exists a functional representation for the random matrix in the following form
A() =
P
0
X
i=0
A
i

i
(); (6.2)
whereA
i
2 Mat
nn
(R) is deterministic, andf g is an orthonormal basis in
def
= L
2
(
;F

;).
The analytical form in Eq. 6.2 could be synthesized from a Karhunen-Loeve (KL) expansion
of an underlying stochastic process, in which case the random variables would correspond to
the KL random variables, the deterministic coefficientA
i
would depend on theith eigenvector
and eigenvalue of the covariance matrix and
i
() would be a set of d-orthogonal random
variables (Das et al. (2009)). Equation (6.2) could also be obtained directly from a Polynomial
Chaos expansion (PCE) of the underlying stochastic process, with coefficients estimated from
available information (Desceliers et al. (2006), Ghanem and Doostan (2006)). We assume that
such representation is available to the analyst.
Stochastic fluctuations in entries of matrix A will induce corresponding fluctuations in its
eigenvalues and eigenvectors. Thus
s
:
7! R and
s
:
7! R
n
areF

-measurable random
variables. With mild conditions on matrix A, both
s
and
s
can be assumed to be square-
integrable and are thus in . As a basis in , we choose polynomials orthogonal with respect to
the measure ofd. The approximate representation for the eigenpairs is then given by

s
()
P
X
i=0

s
i

i
(); s = 1; ;n

s
()
P
X
i=0

s
i

i
(); s = 1; ;n;
(6.3)
84
where
s
i
and
s
i
are the chaos coefficients given by

s
i
=h
s
;
i
i; s = 1; ;n;

s
i
=h
s
;
i
i; s = 1; ;n;
(6.4)
with the following notational convention that will be used throughout this dissertation,
hf;
i
i
def
=
Z Z

Z
R
d
f()
i
()d(); (6.5)
for any integrable functionf defined onR
d
.
Therefore, the Galerkin projection on the stochastic coordinates is formulated as
h
P
0
X
i=0
A
i

i
P
X
i=0

i

i

P
X
i=0

i

i
P
X
i=0

i

i
;
k
i = 0; k = 0; ; 2P; (6.6)
where, for the sake of notational simplicity, we have dropped the superscripts denoting the eigen-
mode. This weak formulation leads to a nonlinear system of equations in the following form
P
0
X
i=0
P
X
j=0
A
i

j
c
ijk
=
P
X
i=0
P
X
j=0

i

j
c
ijk
k = 0; ; 2P; (6.7)
wherec
ijk
=h
i

j

k
i andhi is the expectation operator. The random eigen-solution is ob-
tained by solving Eq. 6.7 for the deterministic coefficients
i
and
i
. In Ghanem and Ghosh
(2007) and Ghosh and Ghanem (2008), two numerical schemes are proposed for solving this set
of equations; one based on the Newton-Raphson method, and one based on an optimization prob-
lem. Instead of solving for the chaos coefficients in Eq. 6.7, we adapt the power iteration method
to calculate the PCE of the random eigenpairs using efficient iterative algorithms. Two different
85
algorithms are proposed. First one produces the PCE of the dominant random eigenpair, and the
second one is capable of estimating a dominant subspace of arbitrary dimension.
6.2.2 Solution Algorithms for Dominant Eigenpairs
Among the solution algorithms for the deterministic eigenvalue problem, methods based on power
iteration are well-adapted to high-performance computing as they can be formulated efficiently as
matrix-vector (mat-vec) multiplication. We extend this iterative scheme to the stochastic eigen-
value problem.
The main idea behind the Power Iteration algorithm is that an initial vector under consecutive
application of the matrixA will converge to the direction of the dominant eigenvector, provided
proper normalization at each step. The only case where convergence is not achieved is when
the initial vector has zero component in the direction of the actual dominant eigenvector. The
deterministic Power Method is summarized in Algorithm 1.
Algorithm 1 Deterministic power iteration method
Step 1. Choose an initial normal vectoru
(0)
k = 0
while>
tol
do
Step 2.y
(k)
=Au
(k)
Step 3.u
(k+1)
=
y
(k)
ky
(k)
k
=ku
(k+1)
u
(k)
k
k =k + 1
end while
Step 4. = (u
(k1)
)
T
A(u
(k1)
)
86
In the stochastic version, we iteratively apply the random matrix to an initial random vector,
where both the random matrix and vector are represented in the PCE form. Next, we explain the
modification to each step in the algorithm to make it applicable for the stochastic case.
Initialization and mat-vec multiplication We start with a PC representation for an initial (nor-
mal) random vector
u
(0)
() =
P
X
i=0
u
(0)
i

i
() (6.8)
We apply the random matrixA() to the left and compute the product
y
(0)
() =A()u
(0)
(); (6.9)
using Galerkin projection - that is, we replace the left hand side with a PCE representation for
y():
P
X
i=0
y
(0)
i

i
() =A()u
(0)
(); (6.10)
and then calculate its deterministic chaos coefficient as
y
(0)
k
=
P
0
X
i=0
P
X
j=0
A
i
u
(0)
j
c
ijk
(6.11)
wherec
ijk
=h
i

j

k
i.
Normalization We seek to characterize the eigenpairs (
s
(!);
s
(!)) such that the eigenvec-
tors are normalized almost surely, i.e.
k
s
(!)k
R
n = 1; a.s.; s = 1; ;n; (6.12)
87
wherekk
R
n denotes the l
2
-norm inR
n
. Therefore, at the kth iteration in our algorithm, we
normalize the random vectory
(k)
() and obtain a new vectoru
(k+1)
() that is normalized in the
following sense
ku
(k+1)
()k
R
n = 1; a.s.; (6.13)
which is equivalent to
u
(k+1)
() =
y
(k)
()
ky
(k)
()k
R
n
a.s.
=
P
P
i=0
y
(k)
i

i
()
k
P
P
i=0
y
(k)
i

i
()k
R
n
a.s.
(6.14)
The PC coefficient ofu
(k+1)
() is approximated as follows
u
(k+1)
i
=
Q
X
q=1
y
(k)
(
q
)
ky
(k)
(
q
)k
R
n

i
(
q
)w
q
; (6.15)
wheref
1
; ;
Q
g andfw
1
; ;w
Q
g are quadrature points and associated weights, respec-
tively. Using the quadrature integration we will incur a small error as follows
u
(1)
i
=hu
(1)
;
i
i +
quad
i
;
(6.16)
which
quad
i
denotes the integration error which decreases asQ increases.
Convergence Criterion The algorithm stops when the following inequality in terms of the PCE
coefficients of the vectoru
(k)
() is satisfied
P
X
i=0
ku
(k)
i
u
(k1)
i
k
R
n
tol
(6.17)
88
Eigenvalue Calculation Once the PCE of the dominant eigenvector is calculated, Eq. 6.7 can
be solved for the PC coefficients of the corresponding random eigenvalue.
The steps in the proposed Stochastic Power Iteration are summarized in Algorithm 2:
Algorithm 2 Stochastic power method
Step 1. Choose a set of deterministic coefficientsfu
(0)
0
; ;u
(0)
P
g
k = 0
while>
tol
do
Step 2.y
(k)
r
=
P
P
0
i=0
P
P
j=0
A
i
u
(k)
j
c
ijr
r = 0; ;P
Step 3. Normalizefy
(k)
0
; ;y
(k)
P
g and obtainfu
(k+1)
0
; ;u
(k+1)
P
g
=
P
P
i=0
ku
(k+1)
i
u
(k)
i
k
k =k + 1
end while
Step 4. Solve Eq. 6.7 withf
i
g =fu
(k1)
i
g forf
0
; ;
P
g.
6.2.3 Analysis of Stochastic Power Iteration
In order to prove convergence for the stochastic power iteration algorithm, we will make use of
the following two Lemmas. The proofs of these Lemmas is included in the Appendix.
Lemma 1 The PC representation of the exact dominant eigenvector is a fixed point of iterations
defined by Algorithm 2.
Lemma 2 Starting with an initial vector composed of the first two dominant eigenvectors, i.e.
u
(0)
() =
1

1
()+
2

2
() =
1
P
i

1
i

i
()+
2
P
i

2
i

i
(), the solution of the algorithm
eventually converges, in direction, to the PCE of the dominant eigenvector.
89
Theorem Given the PC representation for the random matrixA(), one of the following two
statements is true for the solution obtained by Algorithm 2
lim
k!1
P
X
i=0
ku
(k)
i

1
i
k
R
n = 0 as P;Q!1
lim
k!1
P
X
i=0
ku
(k)
i
+
1
i
k
R
n = 0 as P;Q!1
(6.18)
wheref
1
i
g are the chaos coefficients of the true PCE of the dominant eigenvector (
1
() =
P

1
i

i
()), andfu
(k)
i
g are the chaos coefficient in a P -th order expansion, obtained using a
quadrature-based normalization atQ quadrature points
Proof In the light of Lemmas 1 and 2, the proof easily follows.
Algorithm for Dominant Subspaces
In many applications, the analyst is interested in calculating not only the dominant eigenpair, but
a dominant subspace. In this section, we propose an iterative algorithm which finds the subdomi-
nant eigenpairs. This algorithm is developed based on the Subspace Iteration algorithm.
In the Power Iteration technique, if instead of a single initial vector, one chooses a set of
m vectors and apply the matrix to them multiplicatively, all the vectors will converge to the
dominant eigenvector. However, if we orthogonalize thesem vectors at each step, the resulting
vectors will evolve in independent directions and will eventually converge, under this process, to
the eigenvectors spanning the dominantm-dimensional subspace. Algorithm 3 summarizes the
steps of this method for the calculation of the first m dominant eigenpairs of the matrix A2
Mat
nn
(R):
90
Algorithm 3 Deterministic subspace iteration
Step 1 Choose an initial set ofm normalized vectorsf
1
u
(0)
; ;
m
u
(0)
g
k = 0
while>
tol
do
Step 2 (
s
y
(k)
) =A (
s
u
(k)
); s = 1; ;m andY = [
1
y
(k)
; ;
m
y
(k)
]
Step 3 Compute the QR factorization,QR =Y , and set [
1
u
(k+1)
; ;
m
u
(k+1)
] =Q
=
P
m
s=1
k
s
u
(k+1)

s
u
(k)
k
k =k + 1
end while
Step 4
s
= (
s
u
(k1)
)
T
A(
s
u
(k1)
) s = 1; ;m
Here,V is the set ofm vectors inR
n
, which are distinguished by their left superscripts, and
i

c
is thei-th most dominant eigenvector given by a converged iteration.
In what follows, we explain the way each step is adapted to the random eigenvalue problem.
First, we initialize the set ofm random vectors by setting their chaos coefficients. The initial set
ofm random vectors resulting from these coefficients need not be orthogonal. For example, we
can choose the eigenvectors of the mean matrix as the initial guess. In Step 2, we solve the weak
form of the product similar to that in the Step 2 of the Power Iteration.
QR Factorization In order to decompose them random vectors we have developed a stochastic
variation of the Gram-Schmidt decomposition. We use a quadrature-based decomposition, where
given a set ofm random vectors, a new set ofm random vectors are generated that are orthogonal
at every point in an acceptably large subset of the stochastic domain.
91
Let us assume that we have a set of m random vectors denoted byf
1
y(); ;
m
y()g,
where
i
y() : R
d
! R
n
, and for each random vector a PCE is available. The objective is to
compute the PCE of a new set of m “orthogonal” random vectorsf
1
q(); ;
m
q()g where
i
q() : R
d
! R
n
. Ideally, these new orthogonal random vectors can be obtained by satisfying
the Gram-Schmidt equalities almost surely, as follows
k
q() =
k
y()
k1
X
j=1
j
q()
h
k
y();
j
q()i
R
n
h
j
q();
j
q()i
R
n
k = 1; ;m a.s. (6.19)
Let the random vector
jk
() denote the term under the summation,
jk
() =
h
k
y();
j
q()i
R
n
h
j
q();
j
q()i
R
n
(
j
q()) a.s.; (6.20)
Eq. 6.19 can then be rewritten as
k
q() =
k
y()
k1
X
j=1
jk
() k = 1; ;m a.s. (6.21)
We rely on the weak form of this equation for computational convenience. Thus, given the PCE
of random vector
jk
() in the form,
jk
() =
P
i
jk

i

i
(), it is easy to see that
k
q
i
=
k
y
i

k1
X
j=1
jk

i
k = 1; ;m; i = 0; ;P : (6.22)
The coefficients
jk

i
can be readily evaluated as,
jk

i
=h
jk
;
i
i =
1
N
N
X
r=1
jk
(
(r)
)
i
(
(r)
); (6.23)
92
where
(r)
refers to a realization of random vector.
Algorithm 4 summarizes the steps in the proposed stochastic subspace iteration scheme.
Algorithm 4 Stochastic subspace iteration
Step 1 Choose a set of chaos coefficients for m normalized random vectors
ff
1
u
(0)
i
g; ;f
m
u
(0)
i
gg
k = 0
while>
tol
do
Step 2
s
y
(k)
r
=
P
P
0
i=0
P
P
j=0
A
i
(
s
u
(k)
j
)c
ijk
s = 1; ;m
Step 3 Compute the stochastic QR factorization (Eq. 6.21, usingf
s
y
j
g =f
s
y
(k)
j
g) and set
f
s
u
(k+1)
j
g =f
s
q
j
g
=
P
m
s=1
P
P
i=0
k
s
u
(k+1)
i

s
u
(k)
j
k
k =k + 1
end while
Step 4 Solve Eq. 6.7 withf
i
g =f
s
u
(k1)
i
g forf
s

i
g s = 1; ;m.
Measure of discrepancy
We will consider two measures of discrepancy between representations of the solution to the
random eigenvalue problem. The first measure is consistent with the Hilbert space structure in
which we have formulated the problem. The second measure involves a comparison of some
relevant probability density functions.
Error in L
2
convergence We define the L
2
-error as the norm of the difference between two
PCE representations. This error definition is used to evaluate the performance of the algorithm.
93
Using superscriptsex to denote an exact quantity, andPC to denote a PCE approximation of that
quantity, we introduce the following two error measures for the eigenvectors and eigenvalues,

=ke

k
2
R
n

=
Z
k
ex
()
PC
()k
2
R
nd =
P
X
i=0
k
ex
i

PC
i
k
2
R
n; (6.24)
and

=je

j
2

=
Z
j
ex
()
PC
()j
2
d =
P
X
i=0
j
ex
i

PC
i
j
2
: (6.25)
Clearly, in most cases of interest,
ex
and
ex
are not available, and they are obtained by suitable
approximation of the integrals in the above definitions. We will rely on a Monte Carlo-based
approximation of these integrals, resulting in the following approximations to the error measures

1
N
N
X
i=1
k
MC
(
(i)
)
PC
(
(i)
)k
2
R
n; (6.26)
and

1
N
N
X
i=1
j
MC
(
(i)
)
PC
(
(i)
)j
2
; (6.27)
whereN refers to the number of samples in the Monte Carlo scheme, andf
(i)
g refers toN inde-
pendent samples of random variable. Here,
MC
(
(i)
) and
MC
(
(i)
) refer to the deterministic
eigenpairs associated with realization
(i)
.
Error in convergence in distribution Our second measure of discrepancy consists in compar-
ing probability density functions of the random eigenpairs obtained through our proposed ap-
proximation scheme and the probability density functions obtained through an extensive Monte
94
Carlo sampling scheme. We introduce this error to show convergence in distribution for the ran-
dom eigenpairs. Letf
PC

() andf
ex

() denote the probability density functions of
PC
() and

ex
(), respectively. The solution is said to converge in distribution to the exact random eigen-
value, if the following condition holds
f
PC

()f
ex

() = 0; (6.28)
at every 2 R
+
at which f
ex

is continuous. In order to quantify the error in satisfying this
condition, the Kullback-Leibler (KL) divergence is used, which forf
PC

andf
ex

is given by
D

KL
=
Z
1
0
f
ex

() ln
f
ex

()
f
PC

()
d (6.29)
In order to investigate the convergence in distribution for random eigenvectors, the distribu-
tion of the phase angles between the random vector and the Cartesian coordinate axes are chosen
for comparison (Fig. 6.1). Letf
PC
1
(); ;
PC
n
()g andf
ex
1
(); ;
ex
n
()g denote the co-
sine angles corresponding to random eigenvectors
PC
() and
MC
(), respectively, with their
associated pdf’s denoted byff
PC

1
(
1
); ;f
PC
n
(
n
)g andff
ex

1
(
1
); ;f
ex
n
(
n
)g. The KL di-
vergence for a random eigenvector is then defined as follows
D

KL
def
=
1
n
n
X
i=1
Z
2
0
f
ex

i
(
i
) ln
f
ex

i
(
i
)
f
PC

i
(
i
)
d
i
: (6.30)
95
Figure 6.1: The three angles between the vector and coordinates, i.e.;; and
.
Here, we rely on the test points of the kernel smoothing density estimates in order to approximate
the KL-divergence terms as follows
D

KL

1
N
N
X
k=1
f
ex

(
(k)
) ln
f
ex

(
(k)
)
f
PC

(
(k)
)
D

KL

1
n
n
X
i=1
1
N
N
X
k=1
f
ex

i
(
(k)
i
) ln
f
ex

i
(
(k)
i
)
f
PC

i
(
(k)
i
)
;
(6.31)
wheref
(k)
g andf
(k)
i
g refer to theN test points used in the kernel smoothing density estimators.
6.2.4 Numerical Examples
In what follows section, we investigate the performance of the proposed algorithms. First, a
low-dimensional system is chosen so that the eigenpairs can be derived analytically. Specifically,
a 2 Degree-Of-Freedom (DOF) mass-spring oscillator is chosen where the stiffness of one of
the springs is random. As a second example, a finite element model of an elastic cantilever
beam with random elasticity is considered. Since the random eigenpairs of this system cannot be
96
calculated analytically, the results of MC simulations are assumed to yield a sufficiently accurate
approximation to the exact PCE coefficients and are used to verify our procedure.
2 DOF Mass-Spring Oscillator
A 2 DOF oscillator example is considered (Fig. 6.2), where the stiffness of the first spring, denoted
byk
1
is the only random parameter in the system, given by the following equation
k
1
(!) =

k
1
+(
2
1) (6.32)
where is a standard normal random variable. This particular form generates a random variable
k
1
(!) that is positive almost surely (see Desceliers et al. (2005)). The numerical values used for
the model parameters are as follows:

k
1
=k
2
= 10N=m; k
3
= 4N=m; m
1
=m
2
= 1kg and
= 2N=m. The exact solution for this eigenvalue problem can be derived analytically as

1
=

k
1
+(
2
1) +
q

k
2
1
+ (
2
1)(2

k
1
20) 20

k
1
+
2

4
2
2

2
+
2
+ 164
2
+ 9

2
=

k
1
+(
2
1)
q

k
2
1
+ (
2
1)(2

k
1
20) 20

k
1
+
2

4
2
2

2
+
2
+ 164
2
+ 9
Given the aforementioned numerical values, the exact eigenvalues are

1
= 13 +
2
+
p

4
2
2
+ 17

2
= 13 +
2

p

4
2
2
+ 17
A second order PC expansion is considered for the random eigenpairs. Here, we report the
results obtained by the random subspace algorithm. Fig. 6.3 shows the comparison between the
97
Figure 6.2: The 2DOF mass-spring system (after Ghosh and Ghanem (2008)).
PDF’s of the phase angles of the exact and approximate random eigenvectors. Fig. 6.4 compares
the exact and approximate functional forms of the random phase angles versus the uncertainty
source. Similarly, Fig. 6.5 compares the pdf’s of the exact and approximate random eigenvalues,
whereas Fig. 6.6 compares the functional forms of these two. It can be seen that the solutions for
the random eigenvalues and eigenvectors are in close agreement with the exact values.
98
−10 0 10 20 30 40 50 60
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
1st eigenvector
θ
1
(degree)
pdf

exact
subspace
(a)
30 40 50 60 70 80 90 100
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
1st eigenvector
θ
2
(degree)
pdf

exact
subspace
(b)
30 40 50 60 70 80 90 100
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
2nd eigenvector
θ
1
(degree)
pdf

exact
subspace
(c)
−10 0 10 20 30 40 50 60
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
2nd eigenvector
θ
2
(degree)
pdf

exact
subspace
(d)
Figure 6.3: Comparison between the exact and approximate pdf’s of the phase angles of the first
and second eigenvector with respect to the first and second axes.
99
−3 −2 −1 0 1 2 3
0
10
20
30
40
50
60
1st eigenvector
ξ
θ
1
(degree)

exact
subspace
(a)
−3 −2 −1 0 1 2 3
30
40
50
60
70
80
90
1st eigenvector
ξ
θ
2
(degree)

exact
subspace
(b)
−3 −2 −1 0 1 2 3
30
40
50
60
70
80
90
2nd eigenvector
ξ
θ
1
(degree)

exact
subspace
(c)
−3 −2 −1 0 1 2 3
0
10
20
30
40
50
60
2nd eigenvector
ξ
θ
2
(degree)

exact
subspace
(d)
Figure 6.4: Comparison between the functional forms of the exact and approximate phase angles
of the first and second eigenvector relative to first and second axes.
100
15 20 25 30 35 40 45 50
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
λ
1
pdf

exact
subspace
(a)
5 10 15 20 25 30 35
0
0.2
0.4
0.6
0.8
1
1.2
1.4
λ
2
pdf

exact
subspace
(b)
Figure 6.5: Comparison between the exact and approximate pdf’s of the first and second eigen-
values of the 2DOF system.
−5 0 5
15
20
25
30
35
40
45
50
55
60
65
ξ
λ
1

exact
subspace
(a)
−5 0 5
5
10
15
20
25
30
35
ξ
λ
2

exact
subspace
(b)
Figure 6.6: Comparison between the functional forms of the exact and approximate solutions for
the first and second eigenvalues of the 2DOF system
Cantilever Beam with Random Elasticity
The finite element model of a one-dimensional cantilever beam with random elasticity,E, is next
considered (Fig. 6.7). The length of the beam is divided into 10 elements of equal size. Only
rotational and transverse degrees of freedom are retained and axial deformation is neglected.
101
Figure 6.7: Discretization of the cantilever beam with 10 elements and random elasticityE(x;!)
and random density(x;!) in the most general case
The spatial variation in the elastic modulus of the beam is represented by the random process,
E(x;!), that follows a lognormal distribution. Specifically,
E(x;!) = exp(g(x;!)) (6.33)
where g(x;!) is a Guassian random field with its mean equal to zero and its autocorrelation
function given by
r(x
1
;x
2
) = sinc

x
1
x
2
2L

2
(6.34)
whereL is the correlation length, which in this example is set equal to the length of the beam. A
coefficient of the variation of 10% is considered for the Gaussian process. Consider the following
KL expansion for the Gaussian random field
g(x;!) =
d
X
i=0
g
i
(x)
i
(!); (6.35)
whered is the order of the expansion andf
i
g are a set of independent Gaussian random variables.
Given this KL expansion, the lognormal random field can now be represented by the following
P th-orderd-dimensional PC representation (Roger and Ghanem (1999))
E(x;!) =
P
X
i
E
i
(x)
i
((!)); (6.36)
102
where
E
i
(x) =
h
i
()i
h
2
i
i
exp
2
4
g(x) +
1
2
d
X
j=1
g
i
(x)
2
3
5
;
i
=
i
g
i
: (6.37)
Using a PCE representation with Hermite polynomials of order 2 and dimension 3 for the elastic
modulus, the PCE of the random stiffness matrix of the finite element model, which takes values
in Mat
2020
(R) is formed. The random eigenvalue problem is then given by
K()() =()(): (6.38)
The 2nd order 3-dimensional PCE form of the four most dominant eigenpairs were calcu-
lated using subspace iteration algorithm. First, we evaluate the performance of this algorithm in
estimating the random eigenvectors. To this end, the pdf’s of the phase angles of the random
eigenvectors with respect to the Cartesian coordinate axes are plotted. Figs. 6.8 and 6.9 depict
a representative set of these pdf’s for the first and fourth random eigenvectors with respect to
different axes, which show agreement with results obtained from Monte Carlo simulations.
103
60 65 70 75 80 85 90 95
0
0.02
0.04
0.06
0.08
0.1
0.12
θ
1
(degree)
pdf

MC
Subspace method
(a)
40 50 60 70 80 90 100
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
θ
5
(degree)
pdf

MC
Subspace method
(b)
89.4 89.6 89.8 90 90.2
0
1
2
3
4
5
6
θ
6
(degree)
pdf

MC
Subspace method
(c)
60 65 70 75 80 85 90 95
0
0.05
0.1
0.15
0.2
0.25
θ
11
(degree)
pdf

MC
Subspace method
(d)
89.2 89.4 89.6 89.8 90 90.2
0
1
2
3
4
5
θ
14
(degree)
pdf

MC
Subspace method
(e)
89.4 89.6 89.8 90 90.2
0
1
2
3
4
5
6
7
8
θ
16
(degree)
pdf

MC
Subspace method
(f)
Figure 6.8: Comparisons between the distribution the first random eigenvector, obtained from
Monte Carlo and Subspace Iteration; shown are PDF of phase angles with respect to (a) first, (b)
fifth, (c) sixth, (d) 11
th
, (e) 14
th
, and (f) 16
th
coordinate axes
104
65.5 66 66.5 67 67.5 68 68.5
0
0.5
1
1.5
2
2.5
3
3.5
θ
1
(degree)
pdf

MC
Subspace method
(a)
65 70 75 80 85 90 95
0
0.02
0.04
0.06
0.08
0.1
0.12
θ
5
(degree)
pdf

MC
Subspace method
(b)
62 64 66 68 70 72 74 76
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
θ
7
(degree)
pdf

MC
Subspace method
(c)
60 65 70 75 80 85 90 95
0
0.02
0.04
0.06
0.08
0.1
θ
9
(degree)
pdf

MC
Subspace method
(d)
60 65 70 75 80 85 90
0
0.05
0.1
0.15
0.2
0.25
θ
13
(degree)
pdf

MC
Subspace method
(e)
87 87.5 88 88.5 89
0
1
2
3
4
5
θ
16
(degree)
pdf

MC
Subspace method
(f)
Figure 6.9: Comparisons between the distribution the fourth random eigenvector, obtained from
Monte Carlo and Subspace Iteration; shown are PDF of phase angles with respect to (a) first, (b)
fifth, (c) seventh, (d) ninth (e) 13
th
, and (f) 16
th
coordinate axes
105
In the proposed algorithms, the normalization of the random vectors was imposed at few
quadrature points. In order to investigate whether the random vectors have unit norm almost any-
where in the stochastic domain, the norms of the four solution eigenvectors are plotted. Figs. 6.10
show a close agreement with the almost sure unit norm constraint of Eq. 6.12. In order to evaluate
the performance of the subspace algorithm in estimating the random eigenvalues, the pdf’s of
the solution random eigenvalues are plotted in Fig. 6.11 which show perfect agreement with the
results obtained by Monte Carlo simulations.
106
−3 0 3
0.95
1
1.05
1.1
1.15
|| φ
1
||
ξ
1
(a)
−3 0 3
0.95
1
1.05
1.1
1.15
1.2
|| φ
1
||
ξ
2
(b)
−3 0 3
0.95
1
1.05
1.1
1.15
|| φ
1
||
ξ
3
(c)
−3 0 3
0.95
1
1.05
1.1
1.15
|| φ
2
||
ξ
1
(d)
−3 0 3
0.95
1
1.05
1.1
1.15
1.2
1.25
|| φ
2
||
ξ
2
(e)
−3 0 3
0.98
1
1.02
1.04
1.06
|| φ
2
||
ξ
3
(f)
−3 0 3
0.99
1
1.01
1.02
1.03
1.04
|| φ
3
||
ξ
1
(g)
−3 0 3
0.98
1
1.02
1.04
1.06
1.08
|| φ
3
||
ξ
2
(h)
−3 0 3
0.98
1
1.02
1.04
1.06
1.08
|| φ
3
||
ξ
3
(i)
−3 0 3
0.998
1
1.002
1.004
1.006
|| φ
4
||
ξ
1
(j)
−3 0 3
0.995
1
1.005
1.01
1.015
|| φ
4
||
ξ
2
(k)
−3 0 3
0.995
1
1.005
|| φ
4
||
ξ
3
(l)
Figure 6.10: Figures on each row correspond to a particular mode and show the norm of the
corresponding random eigenvector versus three uncertainty sources,
1
,
2
and
3
, which are
independent standard normal random variables. For each value of an uncertainty source, the error
bar corresponds to the variation induced by the other two independent uncertainty sources.
107
150 200 250 300 350 400 450 500 550
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
λ
pdf

1st mode
2nd mode
3rd mode
4th mode
Figure 6.11: The pdf’s of the four dominant random eigenvalues. The pdf’s obtained by the
subspace iteration algorithm overlap with those obtained by Monte Carlo simulations.
Furtheremore, Table 6.1 tabulates values for relative L
2
-error and errors in converging in
distribution for the four dominant random eigenpairs. The relativeL
2
-errors are taken to be the
errors given by Eqs. 6.24 and 6.25 normalized by the corresponding mean values given by the
MC simulations. These small errors are further evidences of the acceptable performance of the
subspace iteration algorithm.
6.2.5 Model Reduction Using Stochastic Eigenpairs
In this section, we present a framework for the analysis of the structural dynamics under the
uncertainties in the material properties and loading. The uncertainty formulation of the second
order differential equation of motion, when the mass and/or stiffness matrices are given by a PCE
is discussed. This formulation can be used in the dynamics analysis of wind turbine blades. Wind
turbine members can be modeled using the Blade Element Momentum method, where the wind
108
Table 6.1: Values of relativeL
2
-error and KL-divergence for the four dominant random eigenpairs
Mode
RelativeL
2
-error KL-divergence

1 1.1010
3
5.3910
2
4.2710
8
1.1210
2
2 8.1910
4
4.2810
2
1.4410
8
1.9310
2
3 1.5010
4
1.5910
2
1.9110
9
1.4710
2
4 1.1210
4
0.6110
2
1.1010
9
1.0010
2
load profile is given as functions of the drag and lift coefficients together with the angle of attack.
The blade is discretized into a finite number of elements, of which each member takes different
values of external loading. The fluctuation in the wind property can then be captured by the
random treatment of the drag and lift coefficients. In a subsequent generic cantilever example, we
show how the uncertainties in the external loading can be incorporated in the calculation of the
stochastic response.
LetM(
M
) =
P
P
M
i=0
M
i

i
(
M
),K(
K
) =
P
P
K
i=0
K
i

i
(
K
) andf(t;
f
) =
P
P
f
i=0
f
i
(t)
i
(
f
)
denote the PC expansions for the random mass matrix, random stiffness matrix and the random
external load, where
M
;
K
and
f
are vector valued random variables referring to the uncer-
tainty sources in the density, stiffness and external load properties, with the dimensionsd
M
,d
K
andd
f
respectively. The stochastic undamped equation of motion will then read
P
M
X
i=0
M
i

i
(
M
) x(t;!) +
P
K
X
i=0
K
i

i
(
K
)x(t;!) =
P
f
X
i=0
f
i
(t)
i
(
f
) (6.39)
109
The random displacement response, x(!) can be subsequently characterized by the following
form
x(t;!) =
P
X
i=0
x
i
(t)
i
() (6.40)
where = [
M
;
K
;
f
]
T
.
In order to reduce the model, the free vibration response of the system is considered leading
to the following generalized eigenvalue problem
M(
M
)
s
(
0
) =
s
(
0
)K(
K
)
s
(
0
) (6.41)
where
0
= [
M
;
K
]
T
with its corresponding measure given by

0. In order to transform this
equation into a regular eigenvalue problem we seem to characterize a random matrixA(
0
) =
P
P
A
i=0
A
i

i
(
0
) such that
A(
0
) =K(
K
)
1
M(
0
)

0 a:s: (6.42)
Assuming smooth function form between the random matrices and their corresponding uncer-
tainty sources, if a solution A
0
(
0
) is obtained such that the above equality holds true at few
quadrature points, one can expect that such approximate solution yields accurate statistical mo-
ments and moreover is close enough to the true solutionA(
0
) in a large region of the stochastic
domain. A quadrature-based matrix inversion is therefore performed in a similar way to the dis-
cussed orthogonalization scheme.
110
Having the random eigenpairsf
s
;
s
i
(
0
)g corresponding to random mass and stiffness ma-
trices, the equation of motion can be decomposed into the following sets of random modal equa-
tions

s
(
0
)
T
M(
M
)
s
(
0
) a
s
(t;) +
s
(
0
)
T
K(
K
)
s
(
0
)a
s
(t;) =
s
(
0
)
T
f(t;
f
)
s
(
0
)
(6.43)
wherex(t;) =
P
s
a
s
(t;)
s
(
0
) are the random variables are represented by their correspond-
ing PCE’s.
In order to illustrate the formulation, the DOF oscillator example discussed in previous sec-
tion is considered. The numerical values used for the model parameters are as follows:

k
1
= 6,
k
2
= 1 N=m ; k
3
= 5 N=m ; m
1
= m
2
= 1 kg and = 2 N=m. A second order PC
expansion is considered for the random eigenpairs. We now investigate the structural dynam-
ics of this 2DOF system. First, the response of the system is calculated based on the values of
the stiffness matrix calculated at 5 quadrature points. Fig. 6.12 shows the random displacement
response corresponding of the first degree of freedom captured by the 5 quadrature-based trajec-
tories. Figs. 6.13 depicts the pdf of this displacement att = 2 and shows that the random solution
synthesized based on the few quadrature-based responses is close in distribution to the precise
behavior obtained by sufficiently large number of MC simulations.
111
0 2 4 6 8 10
−1
0
1
2
3
4
5
x
1
t
Figure 6.12: The random displacement response of the first degree of freedom in the 2DOF
oscillator. The gray lines refer to responses resulting from Monte Carlo samples of the stiffness
matrix. The blue lines refer to the response calculated at the 5 quadrature points. The blue lines
are weighted based on the weight of the corresponding quadrature points.
0 0.5 1 1.5 2 2.5 3 3.5
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
pdf
y
1
@ t=2

MC
PC−quadrature
Figure 6.13: The pdf of the random displacement response of the first degree of freedom in the
2DOF oscillator at time 2 sec. The solid lines refer to responses resulting from Monte Carlo
samples of the stiffness matrix. The dashed lines refer to the collocation based solution.
We now focus on the model reduction using the random eigenvector bases. The dynamic
response of the 2DOF oscillator is calculated based on the modal decomposition of the equation of
112
motion. Fig. 6.14 show the dynamic response calculated by using both random eigenvector bases
can closely capture the precise random behavior obtained by the large number of MC simulations.
In order to investigate the behavior resulting from a reduced model, Fig. 6.14 evaluates the random
response obtained by including only the dominant random eigenmode.
0 2 4 6 8 10
−0.4
0
0.4
0.8
1.2
1.6
1.8
y
1
t

Exact response (68.2% confidence)
Mean of exact response
Mean of approx. response
Mean +/− std of approx. response
Figure 6.14: The random displacement of the first degree of freedom in the 2DOF system gener-
ated from the stochastic modal decomposition compared with that obtained by direct MC simula-
tions of the equation of motion. Both random eigenvectors are included.
113
0 2 4 6 8 10
−0.4
0
0.4
0.8
1.2
1.6
1.8
y
1
t

Exact response (68.2% confidence)
Mean of exact response
Mean of approx. response
Mean +/− std of approx. response
Figure 6.15: The random displacement of the first degree of freedom in the 2DOF system gener-
ated from the reduced model based on the first stochastic mode compared with that obtained by
direct MC simulations of the equation of motion.
Next, the cantiliver in the previous section with random elesticity is considered. The dominant
structural modes can be obtained by solving the following random eigenvalue problem
K()
1
M() =()(): (6.44)
The 2nd order 3-dimensional PCE form of the four most dominant eigenpairs were calcu-
lated using subspace iteration algorithm, based upon which the random dynamics response of the
cantilever to a constant external loading is calculated. Fig. 6.16 compares the random vertical
displacement at 60 cm from the fixed end obtained by Monte Carlo simulation with that obtained
by the reduced model using only the first random eigenpair. The bias in estimating the mean
values is further removed when the top three random eigenpairs are used as the reduction bases,
as is evident in Fig. 6.17
114
0 1 2 3 4 5 6
−0.002
0
0.004
0.008
0.012
0.016
0.02
0.022
y
1
t

Exact response (68.2% confidence)
Mean of exact response
Mean of approx. response
Mean +/− std of approx. response
Figure 6.16: The random displacement of the point at 60 cm from the fixed end in the cantilever
beam with random elasticity. The response is calculated based on a model constructed based on
the first random eigenmode and is compared with that obtained by direct MC simulations of the
equation of motion.
0 1 2 3 4 5 6
−0.002
0
0.004
0.008
0.012
0.016
0.02
0.022

Exact response (68.2% confidence)
Mean of exact response
Mean of approx. response
Mean +/− std of approx. response
Figure 6.17: The random displacement of the point at 60 cm from the fixed end in the cantilever
beam with random elasticity. The response is calculated based on a model constructed based on
the first three random eigenmodes and is compared with that obtained by direct MC simulations
of the equation of motion.
115
Next, for the same cantilever example, we consider the case where the external loading is
stochastic. We consider a time and spatial invariant constant vertical load at each node in the
cantilever. For the spatial varying case where a spatial correlation function is given, the procedure
discussed before for the random elastic property can be utilized. In this numerical example,
the external load at each point is assumed to be a standard normal random variable with unit
mean and 30% coefficient of variation. Fig. 6.18 shows the stochastic response of the same point
calculated in the previous example. The observed increase in the uncertainty of the structural
response is an indicator of the non-trivial impact of the load uncertainties. For lower levels of
load uncertainties, the variabilities due to the load uncertainties may fall within the uncertainties
due to the structural properties and does not increase the response uncertainty. As it can be seen,
with the use of only three dominant eigenmode, and with incorporating few quadrature points (5
Gauss-Hermite quadrature point for each stochastic dimension), we have been able to obtain an
accurate approximation of the stochastic response.
6.3 Stochastic Eigenvalue Problem - Non-symmetric Matrices
In the eigenvalue analysis of non-symmetric matrices, one cannot apply the method developed
in Section 6.2, as the eigenvalues may be negative or complex. In this section, we explain the
non-intrusive approach for the eigen-analysis of random matrices, when the PC representations of
its components are available. In what follows, we explain the methodology in a generic setting,
where the more generalized notion of a nonlinear map is used which encompasses the eigenvalue
analysis as a special case.
116
0 1 2 3 4 5 6
−0.005
0
0.005
0.01
0.015
0.02
0.025
y
1
t

Exact response (68.2% confidence)
Mean of exact response
Mean of approx. response
Mean +/− std of approx. response
Figure 6.18: The random displacement of the point at 60 cm from the fixed end in the cantilever
beam with random elasticity and uncertain external force. The response is calculated based on a
model constructed based on the first three random eigenmodes and is compared with that obtained
by direct MC simulations of the equation of motion.
Let and be two scalar-valued random variable, satisfying the nonlinear relationship
(!) =L((!)) a.s. (6.45)
Let us further assume that the PCE of the random variable is available and given by
P
i

i

i
().
The non-intrusive approximation of the nonlinear map, first, involves considering a PCE
P
i

i

i
()
for the random variable whose coefficients are to be calculated. A weaker condition to be im-
posed on the PC form of the two random variables can be thought of as
X
i

i

i
(
q
) =L(
X
i

i

i
(
q
)) q = 1; ;Q (6.46)
117
forQ quadrature points, chosen based on the probability measure of. The unknown coefficients
f
0
; ;
P
g can then be calculated as

j
= h
X
i

i

i
;
j
i (6.47)
=
Q
X
q=1
L(
X
i

i

i
(
q
))
j
(
q
)w
q
(6.48)
An alternative approach to the quadrature integration is the approximation of the stochastic inte-
grals with Monte Carlo samples. Specifically, givenM Monte Carlo samplesf
1
; ;
M
g, we
will have

j
= h
X
i

i

i
;
j
i (6.49)
=
1
M
M
X
m=1
L(
X
i

i

i
(
m
))
j
(
m
) (6.50)
with the following weak identity
X
i

i

i
(
m
) =L

X
i

i

i
(
m
)
!
m = 1; ;M (6.51)
The nonlinear mapL can be considered as the solution of the eigenvalue problem. In next sec-
tions we will report on the accuracy of this solution method for solving for the real parts of the
eigenvalues of non-symmetric matrices in the small signal stability analysis of power systems.
118
6.4 Polynomial Chaos Representation for Random Transition Ma-
trices
The uncertainty quantification in problems involving the MaxEnt probability measure for the
random transition matrices, discussed in Chapter 4, only sample-based techniques, i.e. Monte
Carlo simulations, can be used. This will leave the computational and analytical potential of
the non-sampling techniques, such as polynomial chaos expansion, unrealized. In order to enjoy
the rich analytical capabilities of the methods relying on the functional representation of random
variables, a PC representation of the MaxEnt RTM should be derived. In this section, we present
the way the joint pdf representation of the RTM can be used in order to obtain the associated PC
form.
As mentioned in Chapter 4, the rows of the RTM are assumed to be independent. Thus,
for simplicity purposes, we will present the formulation for a single row of the random matrix.
Conforming to the notation of Section 3, letx :
! R
n
denote the random vector whose
MaxEnt joint pdf is know and denoted byP

x
(x). As discussed earlier, the components of this
vector are not independent. Our objective is to find a PCE,
P
i
x
0
i

i
(!) such that
x(!) =
X
i
x
0
i

i
(!) P

x
(x) a.s. (6.52)
In order to explain the PC transformation of the random vector with dependent components,
we start with the description of the sampling from the joint pdfP

x
(x). The first step involves con-
structing the sequence of marginal CDFs. LetF
1
(x
1
) denote the marginal CDF associated with
the random variable x
1
. Subsequently, the series of conditional marginal CDFs, F
2j1
(x
2
jx
1
=
119
X
1
);F
3j1;2
(x
3
jx
1
=X
1
;x
2
=X
2
); ;F
nj1;(n1)
(x
n
jx
1
=X
1
; ;x
n1
=X
n1
) are calcu-
lated. Based on this series of marginal CDFs, a sequence of random variable,x
1
;x
2
;x
3
; are
generated based on the seeds (uniform random variables)u
1
;u
2
;u
3
; . Therefore, we can write
u
1
F
1
1
! x
1
u
2
F
1
2j1
! x
2
u
3
F
1
3j1;2
! x
3
.
.
.
For each sample ofx, we essentially transform a set ofn independent uniform random vari-
able through a nonlinear map to the set ofn dependent components of random vectorx, i.e.
u
L
!x (6.53)
whereu is the vector of sampled independent uniform variables. We can approximate the non-
linear map,L, by a polynomial chaos expansion based on Legendre polynomials. To this end, we
make use the non-intrusive approximation of the nonlinear map, described in Section 6.3.
In our case where the input to the nonlinear map was a set of independent uniform random
variables, we first consider aP th order PCE of the random vector in the RHS of Eq. 6.53, and
denote it by
P
i
x
i

i
(u), wheref
i
g are the set of orthonormal Legendre polynomials. The
unknown coefficientsfx
i
g can the be given by
x
i
=
1
M
M
X
m=1
L(u
m
)
i
(u
m
) i = 0; ;P (6.54)
120
In order to illustrate the procedure, and for the presentation efficiency, let us consider the
3-state Markov transition matrix associated with the disease progression example in Chapter 4.
Therefore, our objective is to find for each row of the transition matrix, a PC form of the random
vectorx which takes values inR
2
. The third component on each row can be readily given as a
function of the other two, almost surely. In what follows, we evaluate the PC approximation of
the MaxEnt probability measure forx. Figs. 6.19 show good agreement between the functional
relationships of thex
1
andx
2
versus the corresponding input uncertaintiesu
1
andu
2
for Monte
Carlo samples of the MaxEnt joint pdf and the PC approximation of the random vector. Such
agreement is also observed between the distributions of the transition probabilities in Figs. 6.20.
6.5 Impact of Demand Uncertainty on the Small Signal Stability
In this section, we consider the impact of the uncertainty in the demand model described in Chap-
ter 4 on the performance of the power system. Specifically, we investigate how the real parts of the
eigenvalues of the linearized power system equations are affected by the demand uncertainties. In
particular, let us consider the 16-20 pm time zone of the day. Let us further assume that the uncer-
tainties are limited to the transition behaviors outgoing the fourth activity state, i.e. home with low
consumption level. In other words, we consider the fourth row of the random transition matrix to
be random. On this row, theP
44
component can be readily calculated given the values for other
three components. Let us also assume that at a certain time in the third time zone a deterministic
estimate of the activity status is known to be
0
= [0:0242 ; 0:0012 ; 0:0061 ; 0:9685]. Our
objective is to find total random demand 10 minutes from the current time. In order to convert the
121
activities into the electricity consumption, the electricity demand (in per unit values) associated
with the four activity status is assumed to be [5; 12; 850; 650] 2:17e
6
.
In order to quantify the impact of the uncertainties on the stability condition of the power
system, a 3rd order PC form is calculated for the random transition matrix, and subsequently for
the random total demand, based on which, the 3rd order 3 dimensional PCE of the total demand at
the time of interest can be readily obtained. Fig. 6.21 compares the pdf’s of the random demand
obtained by 5000 Monte Carlo samples and that by the PC representation of the transition matrix.
The power grid chosen for this analysis is the Western System Coordinating Council (WSCC)
3-machine 9-bus grid (Fig. 6.22). In total, there are 24 state variables. Prior to performing the
eigenvalue analysis, a power flow calculation is carried out, based upon which the system of
nonlinear DAE (Eq. 2.14) is linearized. As the input to the power flow calculation, one needs
to input the the electricity loads at the load buses and assign the power level at the generation,
among other parameters, such as network or generator properties.
The calculated per unit random demand are assigned to the three load buses. The non-intrusive
(Stochastic collocation) scheme described in Section 6.3 is then applied in order to quantify the
impact of the demand uncertainties on the real part of the eigenvalues of the linearized power
system equations, as the quantities of interest. Specifically, a 3-point Gauss-Legendre quadra-
ture is used in order to sample the three dimensional PCE of the random demand, based upon
which deterministic small signal stability analyses are performed. An open source MATLAB
toolbox, called PSAT, was used in order to compute the network equilibrium condition as well as
the system eigenvalues (Milano (2005, 2008)). The computed QoI at the quadrature points are
then projected back on the chaos coordinates leading to the functional representation for these
quantities.
122
Figs. 6.23 and 6.24 compares the pdf’s of the QoI (real parts of the random eigevalues) ob-
tained by sampling from the corresponding estimated PCE and those obtained by 5000 MC sim-
ulations. Out of the total 24 physical modes, there are 10 deterministic modes, associated with
which the resulting PCE contained only the first term (the mean value). In these figures, a good
agreement between the respective results is evident, indicating the accuracy of the non-intrusive
PC-based eigenvalue analysis. In addition to the observed accuracy in the distributions, the unique
advantage of PC-based machinery lies in its functional treatment of random variables which estab-
lishes the foundation for efficient sensitivity analysis, design procedure and other post-processing
procedures.
6.6 Impact of Supply Uncertainty on the Small Signal Stability
A similar approach is pursued in order to quantify the impact of generation uncertainties induced
by the random communication network on the stability condition of the power system. The same
9-bus grid depicted in Fig. 6.22 where observed status of the connectivity network results in the
following adjacency matrix

A =
2
6
6
6
6
6
6
4
0 0:2 0:4
0:2 0 0:3
0:4 0:3 0
3
7
7
7
7
7
7
5
:
The initial condition as well as the cost function parameters are reported in Table 6.2 The
random generation level immediately after introducing a total demand ofPd = 700 , is computed
based on the consensus-based methodology discussed in Chapter 5.
123
Table 6.2: Initial conditions and parameters of the 4 generation units
Unit
i

i

i
P
Gi
(0)(MW )
1 78 7.80 0.04820 100
2 310 7.85 0.01940 300
3 561 7.92 0.01562 100
Based on the proposed stochastic consensus formulation of Chapter 5, the uncertainty in the
scheduled generation for the next time step is quantified and assigned to Bus number 1 and 3 in the
WSCC 9-bus grid. Monte Carlo samples of the resulting probability distributions for the sched-
uled generation are obtained based on which deterministic small signal stability analyses were
carried out. Figs. 6.25 and 6.26 depict the distribution of the real parts of the random eigenvalues.
Very small scatter is observed in the resulting eigenvalues. One can conclude that the uncertainty
induced by the random connectivities between the generators, as modeled in Chapter 5, does not
have significant impact on the stability condition of the multimachine grid.
124
−1 1
−1
1
x
1
(ω) − PC
u
1

u
2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(a)
−1 1
−1
1
x
1
(ω) − MC
u
1

u
2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(b)
−1 1
−1
1
x
2
(ω) − PC
u
1

u
2
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
(c)
−1 1
−1
1
x
2
(ω) − MC
u
1

u
2
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
(d)
−1 1
−1
1
x
3
(ω) − PC
u
1

u
2
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
(e)
−1 1
−1
1
x
3
(ω) − MC
u
1

u
2
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
(f)
Figure 6.19: The relationship between the random transition probabilities x
1
(!);x
2
(!);x
3
(!)
and the input uncertainties (random uniform variables) u
1
and u
2
. Figures on the left column
correspond to the PC representation of the transition probabilities, whereas those on the right one
refer to the Monte Carlo samples.
125
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0
2
4
6
8
10
12
x
1
(ω)
pdf

PC
MC
(a)
−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
0
1
2
3
4
5
6
x
2
(ω)
pdf

PC
MC
(b)
−0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
0
1
2
3
4
5
6
7
x
3
(ω)
pdf

PC
MC
(c)
Figure 6.20: The comparison between the distribution of the random transition probabilities
x
1
(!);x
2
(!);x
3
(!). Dashed line refers to the PC representation of the transition probabilities.
Solid line refers to the Monte Carlo samples.
126
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
0
1
2
3
4
5
6
7
Random demand
pdf

PC
MC
Figure 6.21: The comparison between the distributions of the random electricity demand: solid
line refers to the Monte Carlo samples of the MaxEnt joint pdf, dashed line refers to the case
where the transition matrix is represented by a 3rd order 3 dimensional PCE.
WSCC 3-machine, 9-bus system (Copyright 1977)
Bus 9
Bus 8
Bus 7
Bus 6 Bus 5
Bus 4
Bus 3 Bus 2
Bus 1
Figure 6.22: The Western System Coordinating Council (WSCC) 3-machine 9-bus grid used in
the small signal stability analysis (Milano (2008))
127
−0.74 −0.73 −0.72 −0.71 −0.7 −0.69 −0.68 −0.67
0
10
20
30
40
50
60
70
80
Eigenvalue − mode #4,5
pdf

PC
MC
(a)
−0.2 −0.195 −0.19 −0.185 −0.18 −0.175 −0.17
0
20
40
60
80
100
120
140
Eigenvalue − mode #6,7
pdf

PC
MC
(b)
−5.334 −5.332 −5.33 −5.328 −5.326 −5.324 −5.322
0
50
100
150
200
250
300
Eigenvalue − mode #12,13
pdf

PC
MC
(c)
−5.23 −5.22 −5.21 −5.2 −5.19 −5.18 −5.17 −5.16 −5.15 −5.14
0
5
10
15
20
25
30
35
40
45
50
Eigenvalue − mode #14
pdf

PC
MC
(d)
−3.42 −3.4 −3.38 −3.36 −3.34 −3.32 −3.3 −3.28 −3.26 −3.24
0
50
100
150
200
250
Eigenvalue − mode #15
pdf

PC
MC
(e)
−0.446 −0.445 −0.444 −0.443 −0.442 −0.441 −0.44 −0.439 −0.438
0
200
400
600
800
1000
1200
Eigenvalue − mode #16,17
pdf

PC
MC
(f)
Figure 6.23: The probability density function of the real part of the random eigenvalues; compar-
ison between results obtained from Monte Carlo simulations and those obtained from the 3-point
quadrature based spectral projection.
128
−0.441 −0.4405 −0.44 −0.4395 −0.439 −0.4385 −0.438 −0.4375 −0.437 −0.4365 −0.436
0
100
200
300
400
500
600
700
800
Eigenvalue − mode #18,19
pdf

PC
MC
(a)
−0.435 −0.43 −0.425 −0.42
0
50
100
150
200
250
Eigenvalue − mode #20,21
pdf

PC
MC
(b)
Figure 6.24: The probability density function of the real part of the random eigenvalues; compar-
ison between results obtained from Monte Carlo simulations and those obtained from the 3-point
quadrature based spectral projection.
129
−0.51 −0.505 −0.5 −0.495 −0.49 −0.485 −0.48 −0.475 −0.47
0
10
20
30
40
50
60
70
Real part of the eigenvalue − mode #4,5
pdf
(a)
−0.22 −0.215 −0.21 −0.205 −0.2 −0.195
0
20
40
60
80
100
120
Real part of the eigenvalue − mode #6,7
pdf
(b)
−5.213 −5.212 −5.211 −5.21 −5.209 −5.208 −5.207 −5.206
0
50
100
150
200
250
300
350
400
450
Real part of the eigenvalue − mode #8,9
pdf
(c)
−5.38 −5.36 −5.34 −5.32 −5.3 −5.28
0
50
100
150
200
250
300
Real part of the eigenvalue − mode #10,11
pdf
(d)
−5.38 −5.36 −5.34 −5.32 −5.3 −5.28
0
50
100
150
200
250
300
Real part of the eigenvalue − mode #12,13
pdf
(e)
−5.68 −5.66 −5.64 −5.62 −5.6 −5.58
0
5
10
15
20
25
30
35
Eigenvalue − mode #14
pdf
(f)
Figure 6.25: The probability density function of the real part of the random eigenvalues obtained
from Monte Carlo simulations.
130
−3.335 −3.33 −3.325 −3.32 −3.315 −3.31 −3.305
0
20
40
60
80
100
120
140
Eigenvalue − mode #15
pdf
(a)
−0.424 −0.422 −0.42 −0.418 −0.416 −0.414 −0.412 −0.41
0
50
100
150
200
250
300
Real part of the eigenvalue − mode #16,17
pdf
(b)
−0.408−0.4075−0.407−0.4065−0.406−0.4055−0.405−0.4045−0.404
0
100
200
300
400
500
600
700
800
900
Real part of the eigenvalue − mode #18,19
pdf
(c)
−0.364 −0.362 −0.36 −0.358 −0.356 −0.354
0
50
100
150
200
250
300
350
400
Real part of the eigenvalue − mode #20,21
pdf
(d)
Figure 6.26: The probabity density function of the real part of the random eigenvaluesobtained
from Monte Carlo simulations.
131
Chapter 7
Conclusion
In this dissertation, we presented solution methodologies for the future Smart Grid as a Complex
System of System. An uncertainty management framework for this CSoS was carried out and its
quantitative results were demonstrated. This framework consists of three essential elements; (a)
the uncertainty representation in the model inputs which in this study were the parameters in the
electricity demand model and the connectivities of the distributed generators, (b) the uncertainty
propagation of these uncertainties through a mathematical model which in this study was the
stochastic small signal stability problem of the smart grid and (3) quantification of the uncertainty
impacts on the model outputs and quantities of interests which was the real parts of eigenvalues
of the system.
Due to the large size of the CSoS, solution approaches developed for these systems should
enjoy computational efficiencies. In particular, in extending uncertainty quantification models for
these systems, sampling-based techniques, such as Monte Carlo simulations, will prove ineffi-
cient, due to the prohibitively large number of samples they require. Polynomial Chaos approach
is an ideal candidate for this purpose, as they usually require far less number of numerical solves.
132
The core of the developed methods of this study was based on this non-sampling approach, lend-
ing computational advantage to their potential application.
A major challenge in the management of the complex system is the fact that first-principle
physical models fail to capture many of the complex behaviors that are characteristics of these
system. Therefore, an ensemble of multi-physics multi-scale models that are data-driven is needed
to represent such systems. Uncertainties prevail in different components of these systems, making
it critical to develop models that address their impacts. The contributions in this dissertation can
be thought of as a representative approach and thus be generalized for other instances of CSoS.
Markov chains, in particular, due to their computational advantages and the rich theoretical
studies invested on them, are quite appealing among other approaches for modeling dynamical
components of CSoS. The proposed uncertainty treatment of these models in Chapter 3 can be
readily applied to these systems, or can with additional efforts be extended to other Markov chain
derivatives , such as Markov decision processes, hidden Markov models, or coupled Markov
models.
One of the major insights provided by this research is on how the estimation of the parameters
of Markov chains using insufficient data will involve error. We quantitatively illustrated how this
error results from in cases where the trajectory of the chain is observed over a short period of time.
The probabilistic model for the Markov transition matrices was the first effort to quantitatively
account for the scatter and error in these parameters, respecting the unity-sum constraint for the
rows of these matrices.
133
Bibliography
—. Load forecast bibliography phase i. Power Apparatus and Systems, IEEE Transactions on,
PAS-99(1):53 –58, jan. 1980. ISSN 0018-9510. doi: 10.1109/TPAS.1980.319608.
M. Abbad and J. Filar. Perturbation and stability theory for Markov control problems. Automatic
Control, IEEE Trans 37(9):1415-1420, 1992.
A. Bakirtzis, V . Petridis, and S. Kazarlis. Genetic algorithm solution to the economic dispatch
problem. Generation, Transmission and Distribution, IEE Proceedings-, 141(4):377 –382, jul
1994. ISSN 1350-2360. doi: 10.1049/ip-gtd:19941211.
A.R. Bergen and D.J. Hill. A structure preserving model for power system stability analysis.
Power Apparatus and Systems, IEEE Transactions on, PAS-100(1):25 –35, jan. 1981. ISSN
0018-9510. doi: 10.1109/TPAS.1981.316883.
V . Borkar and P. Varaiya. Asymptotic agreement in distributed estimation. Automatic
Control, IEEE Transactions on, 27(3):650 – 655, jun 1982. ISSN 0018-9286. doi:
10.1109/TAC.1982.1102982.
D.W. Bunn. Short-term forecasting: A review of procedures in the electricity supply industry.
The Journal of the Operational Research Society, 33(6):pp. 533–545, 1982. ISSN 01605682.
URLhttp://www.jstor.org/stable/2581037.
A. Capasso, W. Grattieri, R. Lamedica, and A. Prudenzi. A bottom-up approach to residential
load modeling. Power Systems, IEEE Transactions on, 9(2):957 –964, may 1994. ISSN 0885-
8950. doi: 10.1109/59.317650.
J. Cardell and M. Ilic. Maintaining stability with distributed generation in a restructured industry.
In Power Engineering Society General Meeting, 2004. IEEE, volume 2, pages 2142–2149, june
2004.
R. Cogburn. The Ergodic Theory of Markov Chains in Random Environments. Z. Wahrschein-
lichkeitstheorie verw. Gebiete 66:109-128, 1984.
R. Cogburn. On the Central Limit Theorem for Markov Chains in Random Environments Annals
of Probability 19(2):587-604, 1991.
J.D. Collins and W.T. Thomson. The eigenvalue problem for structural systems with statistical.
AIAA Journal, 7(4):642–648, 1969.
134
L. Conradt and T.J. Roper. Consensus decision making in animals. Trends in Ecology and
Evolution, 20(8):449 – 456, 2005. ISSN 0169-5347. doi: 10.1016/j.tree.2005.05.008. URL
http://www.sciencedirect.com/science/article/pii/S0169534705001564.
G. Cooman, F. Hermand and E. Quaeghebeur. Imprecise Markov chains and their limit behavior.
Prob Eng Inform Sci 23:597-635, 2009.
R. Cottet and M. Smith. Bayesian modeling and forecasting of intraday electricity load. Journal
of the American Statistical Association, 98(464):pp. 839–849, 2003. ISSN 01621459. URL
http://www.jstor.org/stable/30045335.
S. Das, R. Ghanem, and S. Finette. Polynomial chaos representation of spatio-temporal
random fields from experimental measurements. Journal of Computational Physics,
228(23):8726 – 8751, 2009. ISSN 0021-9991. doi: 10.1016/j.jcp.2009.08.025. URL
http://www.sciencedirect.com/science/article/pii/S0021999109004677.
B.J. Debusschere, H.N. Najm, A. Matta, O.M. Knio, R. Ghanem, and Olivier P. Le Maitre.
Protein labeling reactions in electrochemical microchannel flow: Numerical simulation and un-
certainty propagation. Physics of Fluids, 15(8):2238–2250, 2003. doi: 10.1063/1.1582857. URL
http://link.aip.org/link/?PHF/15/2238/1.
M.H. DeGroot. Reaching a consensus. Journal of the American Statistical Association, 69(345):
pp. 118–121, 1974. ISSN 01621459. URLhttp://www.jstor.org/stable/2285509.
C. Desceliers, R. Ghanem, and C. Soize. Polynomial chaos representation of a stochastic pre-
conditioner. International journal for numerical methods in engineering, 64(5):618–634, 2005.
C. Desceliers, R. Ghanem, and C. Soize. Maximum likelihood estimation of stochastic
chaos representations from experimental data. International Journal for Numerical Methods
in Engineering, 66(6):978–1001, 2006. ISSN 1097–0207. doi: 10.1002/nme.1576. URL
http://dx.doi.org/10.1002/nme.1576.
Z. Gaing. Particle swarm optimization to solving the economic dispatch considering the gener-
ator constraints. Power Systems, IEEE Transactions on, 18(3):1187 – 1195, aug. 2003. ISSN
0885-8950. doi: 10.1109/TPWRS.2003.814889.
R. Ghanem. Ingredients for a general purpose stochastic finite elements imple-
mentation. Computer Methods in Applied Mechanics and Engineering, 168(1-4):
19 – 34, 1999. ISSN 0045-7825. doi: 10.1016/S0045-7825(98)00106-6. URL
http://www.sciencedirect.com/science/article/pii/S0045782598001066.
R. Ghanem and P.D. Spanos. Stochastic finite elements: a spectral approach. Dover publica-
tions, 2003. ISBN 0486428184.
R. Ghanem and D. Ghosh. Efficient characterization of the random eigenvalue problem in a
polynomial chaos decomposition. International Journal for Numerical Methods in Engineering,
72(4):486–504, 2007.
135
R. Ghanem and A. Doostan. On the construction and analysis of stochastic models:
Characterization and propagation of the errors associated with limited data. Journal of Compu-
tational Physics, 217(1):63 – 81, 2006. ISSN 0021-9991. doi: 10.1016/j.jcp.2006.01.037. URL
http://www.sciencedirect.com/science/article/pii/S0021999106000386.
¡ce:title¿Uncertainty Quantification in Simulation Science¡/ce:title¿.
D. Ghosh and R. Ghanem. Stochastic convergence acceleration through basis enrichment of
polynomial chaos expansions. International Journal for Numerical Methods in Engineering, 73
(2):162–184, 2008.
R. Glass, T. Brown, A. Ames, J. Linebarger, W. Beyeler, S. Maffitt, N. Brodsky, P. Finley.
Complex Adaptive Systems of Systems (CASoS) Engineering: Mapping Aspirations to Prob-
lem Solutions. in 8th International Conference on Complex Systems and 6th IEEE Interational
Conference on Systems of Systems Engineering (SoSE), (SAND 2011-3354 C), 2011.
G. Gross and F.D. Galiana. Short-term load forecasting. Proceedings of the IEEE, 75(12):1558
– 1573, dec. 1987. ISSN 0018-9219. doi: 10.1109/PROC.1987.13927.
V . Gupta, B. Hassibi, and R.M. Murray. On sensor fusion in the presence of packet-dropping
communication channels. In Decision and Control, 2005 and 2005 European Control Con-
ference. CDC-ECC ’05. 44th IEEE Conference on, pages 3547 – 3552, dec. 2005. doi:
10.1109/CDC.2005.1582712.
M.T. Hagan and S.M. Behr. The time series approach to short term load forecasting. Power
Systems, IEEE Transactions on, 2(3):785 –791, aug. 1987. ISSN 0885-8950. doi: 10.1109/TP-
WRS.1987.4335210.
H.S. Hippert, C.E. Pedreira, and R.C. Souza. Neural networks for short-term load forecasting:
a review and evaluation. Power Systems, IEEE Transactions on, 16(1):44 –55, feb 2001. ISSN
0885-8950. doi: 10.1109/59.910780.
H.S. Hippert, D.W. Bunn, and R.C. Souza. Large neural networks for electricity
load forecasting: Are they overfitted? International Journal of Forecasting, 21(3):
425 – 434, 2005. ISSN 0169-2070. doi: 10.1016/j.ijforecast.2004.12.004. URL
http://www.sciencedirect.com/science/article/pii/S016920700400130X.
P. Kundur. Power System Stability and Control. McGraw-Hill Professional, 1994.
P. Kundur, J. Paserba, V . Ajjarapu, G. Andersson, A. Bose, C. Canizares, N. Hatziargyriou,
D. Hill, A. Stankovic, C. Taylor, T. Van Cutsem, and V . Vittal. Definition and classification of
power system stability ieee/cigre joint task force on stability terms and definitions. Power Sys-
tems, IEEE Transactions on, 19(3):1387 – 1401, aug. 2004. ISSN 0885-8950. doi: 10.1109/TP-
WRS.2004.825981.
J. Lasserre. Invariant Probabilities for Markov Chains on a Metric Space Stat. Probab. Lett.
34:259-265, 1997.
B. Ledergerber, J. von Overbeck, M. Egger, R. Luthy. The Swiss HIV Cohort Study: rationale,
organization and selected baseline characteristics. Soz Praventivmed. 39(6):387-94, 1994.
136
B. Li , J. Si. Robust Optimality for Discounted Infinite-Horizon Markov Decision Processes with
Uncertain Transition Matrices. Automatic Control, IEEE Trans 53(9):2112-2116, 2008.
G. Lu, A. Mukherjea. Invariant Measures and Markov Chains with Random Transition Proba-
bilities Probab. Math. Stat. 17(1):115-138, 1997.
D. Luenberger. Introduction to Dynamic Systems: Theory, Models, and Applications (Wiley,
New York), 1979.
A.A. Mahmoud, T.H. Ortmeyer, and R.E. Reardon. Load forecasting bibliography phase ii.
Power Apparatus and Systems, IEEE Transactions on, PAS-100(7):3217 –3220, july 1981. ISSN
0018-9510. doi: 10.1109/TPAS.1981.316650.
O.P. Le Maitre, O.M. Knio, H.N. Najm, and R. Ghanem. A stochastic projec-
tion method for fluid flow: I. basic formulation. Journal of Computational Physics,
173(2):481 – 511, 2001. ISSN 0021-9991. doi: 10.1006/jcph.2001.6889. URL
http://www.sciencedirect.com/science/article/pii/S0021999101968895.
G. Marcous, Z. Lounis, M.S. Mirza. Life-cycle Assessment of Highway Bridges in Proceedings
of the 2002 Taiwan-Canada Workshop on Bridges, 2002.
F. Milano. An open source power system analysis toolbox. Power Systems, IEEE Transactions
on, 20(3):1199 – 1206, aug. 2005. ISSN 0885-8950. doi: 10.1109/TPWRS.2005.851911.
F. Milano. Power system analysis toolbox: Documentation for psat version 2.0.0, 2008.
I. Moghram and S. Rahman. Analysis and evaluation of five short-term load forecasting tech-
niques. Power Systems, IEEE Transactions on, 4(4):1484 –1491, nov 1989. ISSN 0885-8950.
doi: 10.1109/59.41700.
A. Nilim and L. El Ghaoui. Robust Control of Markov Decision Processes with Uncertain Tran-
sition Matrices. Operations Reseach 53:780-798, 2005.
R. Olfati-Saber. Flocking for multi-agent dynamic systems: algorithms and theory. Auto-
matic Control, IEEE Transactions on, 51(3):401 – 420, march 2006. ISSN 0018-9286. doi:
10.1109/TAC.2005.864190.
R. Olfati-Saber and R.M. Murray. Consensus problems in networks of agents with switching
topology and time-delays. Automatic Control, IEEE Transactions on, 49(9):1520 – 1533, sept.
2004. ISSN 0018-9286. doi: 10.1109/TAC.2004.834113.
R. Olfati-Saber and J.S. Shamma. Consensus filters for sensor networks and distributed sensor
fusion. In Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC ’05.
44th IEEE Conference on, pages 6698 – 6703, dec. 2005. doi: 10.1109/CDC.2005.1583238.
R. Olfati-Saber, J.A. Fax, and R.M. Murray. Consensus and cooperation in networked multi-
agent systems. Proceedings of the IEEE, 95(1):215 –233, jan. 2007. ISSN 0018-9219. doi:
10.1109/JPROC.2006.887293.
S. Orey. Markov Chains with Stochastically Stationary Transition Probabilities Annals of Prob-
ability 19(3):907-928, 1991.
137
J.V . Paatero and P.D. Lund. A model for generating household electricity load profiles. In-
ternational Journal of Energy Research, 30(5):273–290, 2006. doi: 10.1002/er.1136. URL
http://dx.doi.org/10.1002/er.1136.
P. Pai and W. Hong. Support vector machines with simulated annealing algorithms
in electricity load forecasting. Energy Conversion and Management, 46(17):2669
– 2688, 2005. ISSN 0196-8904. doi: 10.1016/j.enconman.2005.02.004. URL
http://www.sciencedirect.com/science/article/pii/S0196890405000506.
M. Pellissetti and R. Ghanem. A method for the validation of predictive com-
putations using a stochastic approach. Journal of Offshore Mechanics and
Arctic Engineering, 126(3):227–234, 2004. doi: 10.1115/1.1782915. URL
http://link.aip.org/link/?JOM/126/227/1.
S. Rahman. A solution of the random eigenvalue problem by a dimensional de-
composition method. International Journal for Numerical Methods in Engineer-
ing, 67(9):1318–1340, 2006. ISSN 1097-0207. doi: 10.1002/nme.1665. URL
http://dx.doi.org/10.1002/nme.1665.
J. Red-Horse and R. Ghanem. Polynomial chaos representation of the random eigenvalue prob-
lem. In 40th Structures, Structural Dynamics, and Materials Conference, St. Louis, MS., April
12-15, 1999.
I. Richardson, M. Thomson, and D. Infield. A high-resolution domestic building
occupancy model for energy demand simulations. Energy and Buildings, 40(8):
1560 – 1566, 2008. ISSN 0378-7788. doi: 10.1016/j.enbuild.2008.02.006. URL
http://www.sciencedirect.com/science/article/pii/S0378778808000467.
I. Richardson, M. Thomson, D. Infield, and C. Clifford. Domestic electricity
use: A high-resolution energy demand model. Energy and Buildings, 42(10):
1878 – 1887, 2010. ISSN 0378-7788. doi: 10.1016/j.enbuild.2010.05.023. URL
http://www.sciencedirect.com/science/article/pii/S0378778810001854.
R.O. Saber and R.M. Murray. Consensus protocols for networks of dynamic agents. In American
Control Conference, 2003. Proceedings of the 2003, volume 2, pages 951 – 956, 4-6, 2003. doi:
10.1109/ACC.2003.1239709.
P.W. Sauer and M.A. Pai. Power system steady-state stability and the load-flow jacobian.
Power Systems, IEEE Transactions on, 5(4):1374 –1383, nov. 1990. ISSN 0885-8950. doi:
10.1109/59.99389.
M. Shinozuka and C.J. Astill. Random eigenvalue problems in structural analysis. AIAA Journal,
10(4):456–462, 1972.
C. Soize and R. Ghanem. Physical systems with random uncertainties: Chaos
representations with arbitrary probability measure. SIAM Journal on Scientific
Computing, 26(2):395–410, 2004. doi: 10.1137/S1064827503424505. URL
http://link.aip.org/link/?SCE/26/395/1.
138
O. Stenflo. Markov Chains in Random Environment and Random Iterated Function Systems
Trans. Amer. Math. Soc. 383(9):3548-3562, 2001.
L.G. Swan and V . Ismet Ugursal. Modeling of end-use energy consumption in the res-
idential sector: A review of modeling techniques. Renewable and Sustainable Energy
Reviews, 13(8):1819 – 1835, 2009. ISSN 1364-0321. doi: 10.1016/j.rser.2008.09.033. URL
http://www.sciencedirect.com/science/article/pii/S1364032108001949.
Y . Takahashi. Markov Chains with Random Transition Matrices. KODAI Math. Sem. Rep.
21:426-447, 1969.
J. Tsitsiklis, D. Bertsekas, and M. Athans. Distributed asynchronous deterministic and stochastic
gradient optimization algorithms. Automatic Control, IEEE Transactions on, 31(9):803 – 812,
sep 1986. ISSN 0018-9286. doi: 10.1109/TAC.1986.1104412.
X. Wang, Y . Song, and M. Irving. Modern Power Systems Analysis. Springer, 2009.
J. Widen and E. Wackelgard. A high-resolution stochastic model of domes-
tic activity patterns and electricity demand. Applied Energy, 87(6):1880 –
1892, 2010. ISSN 0306-2619. doi: 10.1016/j.apenergy.2009.11.006. URL
http://www.sciencedirect.com/science/article/pii/S0306261909004930.
J. Widen, M. Lundh, I. Vassileva, E. Dahlquist, K. Ellegard, and E. Wackel-
gard. Constructing load profiles for household electricity and hot water from
time-use datamodelling approach and validation. Energy and Buildings, 41(7):
753 – 768, 2009. ISSN 0378-7788. doi: 10.1016/j.enbuild.2009.02.013. URL
http://www.sciencedirect.com/science/article/pii/S0378778809000413.
A.J. Wood and B.F. Wollenberg. Power Generation, Operation and Control. J. Wiley and Sons,
NY ., 1996. ISBN 0486428184.
D. Xiu and G. Karniadakis. Modeling uncertainty in flow simulations via gen-
eralized polynomial chaos. Journal of Computational Physics, 187(1):137 –
167, 2003. ISSN 0021-9991. doi: 10.1016/S0021-9991(03)00092-5. URL
http://www.sciencedirect.com/science/article/pii/S0021999103000925.
Z. Zhang and M. Chow. Incremental cost consensus algorithm in a smart grid environment.
In Power and Energy Society General Meeting, 2011 IEEE, pages 1–6, july 2011. doi:
10.1109/PES.2011.6039422.
Y . Zhang and B. Vidakovic. Uncertainty Analysis in Using Markov Chain Model to Predict
Roof Life Cycle Performance in 10DBMC International Confrence On Durability of Building
Materials and Components, 2005.
139
Appendix A
In what follows, the detailed proofs for the two Lemmas are presented.
Lemma 1 The PC representation of the exact dominant eigenvector is a fixed point of iterations
defined by Algorithm 2.
Proof Let
1
() =
P

1
i

i
() denote the solution to the following equation
P
0
X
i=0
P
X
j=0
A
i

1
j
c
ijk
=
P
X
i=0
P
X
j=0

1
i

1
j
c
ijk
k = 0; ; 2P; (1)
where
1
is the largest random eigenvalue. We will show that if we start with the PCE of dominant
eigenvector, i.e.
u
(0)
i
=
1
i
k = 0; ;P (2)
the PCE coefficients of the output of the algorithm after one iteration will satisfy
P
X
i=0
ku
(1)
i

1
i
k
R
n! 0: (3)
Starting withu
(0)
i
=
1
i
, the mat-vec multiplication using Galerkin projection yields
y
(0)
k
=
P
0
X
i=0
P
X
j=0
A
i

1
j
c
ijk
k = 0; ;P
=
P
X
i=0
P
X
j=0

1
i

1
j
c
ijk
k = 0; ;P;
(4)
140
with the following pointwise equality
P
X
i=0
y
i
() =
P
X
i=0

1
i

i
()
P
X
i=0

1
i

i
() +
P
proj
() +
P
trunc
() a:s: (5)
where
P
proj
() represents theL
1
error incurred by the Galerkin projection ofP -th order PCE’s,
i.e.

P
proj
() =
2P
X
i=0
y
i

i
()
P
X
i=0

1
i

i
()
P
X
i=0

1
i

i
() a:s:; (6)
with the followingL
2
property (to be interpreted component-wise)
h
P
proj
;
i
i = 0 i = 0; ; 2P: (7)
Due to the analytic extension of the random eigenpairs in the parameter space (i.e. tensor product
of the supports off
i
g), we expect the following component-sizeL
1
property to hold true

P
proj
() = 0 a:s: as P!1: (8)
The other error term in Eq. 5,
P
trunc
(), refers to the order truncation, i.e.
P
X
i=0
y
i

i
() =
2P
X
i=0
y
i

i
() +
P
trunc
() a:s: (9)
with the following property

p
trunc
() = 0 a:s: as P!1: (10)
141
In the normalization step, using Eq. 6.15, the normalized vector at theqth quadrature point
can be written by
y
(0)
(
q
)
ky
(0)
(
q
)k
R
n
=
P
P
i=0
y
i

i
(
q
)
k
P
P
i=0
y
i

i
(
q
)k
R
n
=
P
P
i=0

1
i

i
(
q
)
P
P
i=0

1
i

i
(
q
) +
proj
(
q
) +
p
trunc
(
q
)
k
P
P
i=0

1
i

i
(
q
)
P
P
i=0

1
i

i
(
q
) +
proj
(
q
) +
p
trunc
(
q
)k
R
n
=
P
P
i=0

1
i

i
(
q
)
P
P
i=0

1
i

i
(
q
)
k
P
P
i=0

1
i

i
(
q
)
P
P
i=0

1
i

i
(
q
)k
R
n
as P!1
=
P
X
i=0

1
i

i
(
q
) as P!1;
(11)
where we used the fact that the assumption random vector that the dominant random eigenvector
has unit norm almost anywhere (Eq. 6.12).
Now, we will show that theith chaos coefficient of the normal random vectoru
(1)
() is equal
to that of
1
().
ku
(1)
i

i
k
R
n =k
Q
X
q=1
y
(0)
(
q
)
ky
(0)
(
q
)k
R
n

i
(
q
)w
q
h
1
;
i
ik
R
n
=k
Q
X
q=1

1
(
q
)
i
(
q
)w
q
h
1
;
i
ik
R
n as P!1
=kh
1
;
i
i +
Q
i
h
1
;
i
ik
R
n as P!1
kh
1
;
i
ih
1
;
i
ik
R
n +k
Q
i
k
R
n as P!1
k
Q
i
k
R
n as P!1
(12)
142
where
Q
i
is the small integration error incurred by the numerical quadrature integration in the
calculation of thei-th chaos coefficient usingQ quadrature points. Thus, we can show
ku
(1)
i

i
k
R
n! 0 as P;Q!1 (13)
Lemma 2 Starting with an initial vector composed of the first two dominant eigenvectors, i.e.
u
(0)
() =
1

1
()+
2

2
() =
1
P
i

1
i

i
()+
2
P
i

2
i

i
(), the solution of the algorithm
eventually converges, in direction, to the PCE of the dominant eigenvector.
Proof After the Galerkin mat-vec multiplication, we have
y
(0)
k
=
P
0
X
i=0
P
X
j=0
(
1
A
i

1
j
+
2
A
i

2
j
)c
ijk
k = 0; ;P
=
P
X
i=0
P
X
j=0
(
1

1
i

1
j
+
2

2
i

2
j
)c
ijk
k = 0; ;P;
(14)
Similarly to the proof of the previous lemma, the q-th normalized vector used in Eq. 6.15 will
then be given by
y
(0)
(
q
)
ky
(0)
(
q
)k
R
n
=

1
P
P
i=0

1
i

i
(
q
)
P
P
i=0

1
i

i
(
q
) +
2
P
P
i=0

2
i

i
(
q
)
P
P
i=0

2
i

i
(
q
)
k
1
P
P
i=0

1
i

i
(
q
)
P
P
i=0

1
i

i
(
q
) +
2
P
P
i=0

2
i

i
(
q
)
P
P
i=0

2
i

i
(
q
)k
R
n
as P!1
(15)
143
Theith coefficient of the PCE constructed based on these normalized samples, which is the output
of the algorithm at the end of first iteration will be given by
u
(1)
i
=
Q
X
q=1
y
(0)
(
q
)
ky
(0)
(
q
)k
R
n

i
(
q
)w
q
=
Q
X
q=1

1

1
(
q
)
p

2
1

1
(
q
)
2
+
2
2

2
(
q
)
2

1
(
q
)
i
(
q
)w
q
+
Q
X
q=1

2

2
(
q
)
p

2
1

1
(
q
)
2
+
2
2

2
(
q
)
2

2
(
q
)
i
(
q
)w
q
as P!1;
(16)
which gives the PCEu() =
P
u
i

i
() with the following pointwise property for someq
u
(1)
(
q
) =
P
X
i=0
u
(1)
i

i
(
q
)
=
P
X
i=0
2
4
Q
X
q=1
y
(0)
(
q
)
ky
(0)
(
q
)k
R
n

i
(
q
)w
q
3
5

i
(
q
)
=
P
X
i=0
"
h
y
(0)
ky
(0)
k
R
n
;
i
i +
Q
i
#

i
(
q
)
=
P
X
i=0
h
y
(0)
ky
(0)
k
R
n
;
i
i
i
(
q
) +
P
X
i=0

Q
i

i
(
q
)
=
y
(0)
(
q
)
ky
(0)
(
q
)k
R
n
as Q!1
(17)
At the second iteration, the mat-vec step will then involve calculating the new PCE fory
(1)
()
using the Galerkin projection according to
y
(1)
k
=
P
0
X
i=0
P
X
j=0
A
i
u
(1)
j
c
ijk
k = 0; ;P (18)
The approximate pointwise equivalent to the above mat-vec product at the quadrature point
q
reads
144
y
(1)
(
q
) =A(
q
)u
(1)
(
q
) +
P
proj
(
q
) +
P
trunc
(
q
)
=A(
q
)
y
(0)
(
q
)
ky
(0)
(
q
)k
R
n
+
P
proj
(
q
) +
P
trunc
(
q
) as Q!1
=A(
q
)
"

1

1
(
q
)
p

2
1

1
(
q
)
2
+
2
2

2
(
q
)
2

1
(
q
) +

2

2
(
q
)
p

2
1

1
(
q
)
2
+
2
2

2
(
q
)
2

2
(
q
)
#
as P;Q!1
=

1

1
(
q
)
2
p

2
1

1
(
q
)
2
+
2
2

2
(
q
)
2

1
(
q
) +

2

2
(
q
)
2
p

2
1

1
(
q
)
2
+
2
2

2
(
q
)
2

2
(
q
) as P;Q!1
(19)
and the normalized samples to be used in Eq. 6.15 will then be given by
y
(1)
(
q
)
ky
(1)
(
q
)k
R
n
=

1

1
(
q
)
2

1
(
q
) +
2

2
(
q
)
2

2
(
q
)
k
1

1
(
q
)
2

1
(
q
) +
2

2
(
q
)
2
(
q
)k
R
n
as P;Q!1; (20)
based on which the PC coefficients of the new iteration output is calculated. We can therefore
show that at the end of thek-th iteration we will have the following normalized sample at theq-th
quadrature point
y
(k)
(
q
)
ky
(k)
(
q
)k
R
n
=

1

1
(
q
)
k

1
(
q
) +
2

2
(
q
)
k

2
(
q
)
k
1

1
(
q
)
k

1
(
q
) +
2

2
(
q
)
k

2
(
q
)k
R
n
as P;Q!1; (21)
145
Thus, thei-th PC coefficient of thek-th iteration output will be given by
u
(k)
i
=
Q
X
q=1
y
q

i
(
q
)w
q
=
Q
X
q=1

1

1
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k

1
(
q
)
i
(
q
)w
q
+
Q
X
q=1

2

2
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k

2
(
q
)
i
(
q
)w
q
(22)
We can now write
P
X
i=0
ku
(k)
i

1
i
k
R
n =k
Q
X
q=1

1

1
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k

1
(
q
)
i
(
q
)w
q
+
Q
X
q=1

2

2
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k

2
(
q
)
i
(
q
)w
q

1
i
k
R
n as P;Q!1
k
Q
X
q=1

1

1
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k

1
(
q
)
i
(
q
)w
q

1
i
k
R
n+
k
Q
X
q=1

2

2
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k

2
(
q
)
i
(
q
)w
q
k
R
n as P;Q!1
k
Q
X
q=1

1

1
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k

1
(
q
)
i
(
q
)w
q

1
i
k
R
n+
Q
X
q=1
j
2
j
2
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k
j
i
(
q
)w
q
j as P;Q!1
(23)
146
But, we know
j
2
j
2
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k
<
j
2
j
2
(
q
)
k
p

2
1

1
(
q
)
2k
=j

2

1
j

2
(
q
)

1
(
q
)

k
j

2

1
j
k
(24)
where < 1 is
= inf
q

2
(
q
)

1
(
q
)
(25)
Therefore, we can write
lim
k!1
Q
X
q=1
j
2
j
2
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k
j
i
(
q
)w
q
j lim
k!1
j

2

1
j
k
= 0 (26)
Similary, we can write
lim
k!1

1

1
(
q
)
k
p

2
1

1
(
q
)
2k
+
2
2

2
(
q
)
2k
= lim
k!1

1

1
(
q
)
k
j
1
j
1
(
q
)
k
q
1 + (

2

1
)
2
(

2
(
q
)

1
(
q
)
)
2k
=1
(27)
Therefore, for
1
> 0, substituting Eqs. 24 and 27 in Eq. 23, we will get
lim
k!1
P
X
i=0
ku
(k)
i

1
i
k
R
nk
Q
X
q=1

1
(
q
)
i
(
q
)w
q

1
i
k
R
n as P;Q!1
=k
1
i
+
Q
i

1
i
k
R
n as P;Q!1
=k
Q
i
k
R
n as P;Q!1
= 0 as P;Q!1
(28)
147
For
1
< 0, we will have
lim
k!1
P
X
i=0
ku
(k)
i
(
1
i
)k
R
nk
Q
X
q=1

1
(
q
)
i
(
q
)w
q
+
1
i
k
R
n as P;Q!1
=k
1
i

Q
i
+
1
i
k
R
n as P;Q!1
=k
Q
i
k
R
n as P;Q!1
= 0 as P;Q!1
(29)
148

Abstract (if available)

Abstract Today, many of the infrastructures are composed of several coupled sub-systems, many of which are by themselves complex, and further couplings introduce yet more complexities which pose system operators with several challenging issues. In many of these infrastructure systems, uncertainty is an intrinsic features, whereas in several others it is merely a mathematical abstraction referring to our ignorance due either to the unknown governing physics or to missing information. A major challenge is the predictability of such complex systems of systems under these uncertainties and also the quantification of confidence in the prediction. A predictive science for these systems should (a) assess the uncertainties

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Risk assessment, intrinsic interpolation and computationally efficient models for systems under uncertainty

PDF

A polynomial chaos formalism for uncertainty budget assessment

PDF

Analytical and experimental studies in modeling and monitoring of uncertain nonlinear systems using data-driven reduced‐order models

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Comprehensive uncertainty quantification in composites manufacturing processes

PDF

User-centric smart sensing for non-intrusive electricity consumption disaggregation in buildings

PDF

The power of flexibility: autonomous agents that conserve energy in commercial buildings

PDF

Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Stochastic and multiscale models for urban and natural ecology

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Semantic modeling of outdoor scenes for the creation of virtual environments and simulations

PDF

Behavioral form finding using multi-agent systems: a computational methodology for combining generative design with environmental and structural analysis in architectural design

PDF

Studies into computational intelligence approaches for the identification of complex nonlinear systems

PDF

Studies into vibration-signature-based methods for system identification, damage detection and health monitoring of civil infrastructures

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

PDF

A framework for comprehensive assessment of resilience and other dimensions of asset management in metropolis-scale transport systems

PDF

The warehouse traveling salesman problem and its application

PDF

Studies into data-driven approaches for nonlinear system identification, condition assessment, and health monitoring

PDF

Probabilistic data-driven predictive models for energy applications

Asset Metadata

Creator Meidani, Hadi (author)

Core Title Uncertainty management for complex systems of systems: applications to the future smart grid

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Publication Date 03/04/2014

Defense Date 08/29/2012

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag consensus problem,Markov chains,maximum entropy,OAI-PMH Harvest,random eigenvalue problem,robust decision making,smart grid,Structural dynamics,systems of systems,uncertainty quantification

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ghanem, Roger G. (committee chair), Becerik-Gerber, Burcin (committee member), Caire, Giuseppe (committee member), Masri, Sami F. (committee member), Rechenmacher, Amy Lynn (committee member)

Creator Email meidani@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-94167

Unique identifier UC11289378

Identifier usctheses-c3-94167 (legacy record id)

Legacy Identifier etd-MeidaniHad-1183.pdf

Dmrecord 94167

Document Type Dissertation

Rights Meidani, Hadi

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA