1
Model Predictive Control based Real Time Power
System Protection Schemes
Licheng Jin, Member, IEEE, and Ratnesh Kumar, Fellow, IEEE, and Nicola Elia, Member, IEEE
Abstract—The objective of power system controls is to keep the
electrical flow as well as voltage magnitudes within acceptable
limits in spite of the load and network topology changes. The con-
trol of voltage level is accomplished by controlling the production,
absorption as well as flow of reactive power at various locations in
the system. This paper presents an approach to determine a real
time system protection scheme for maintaining voltage stability
following the occurrence of a contingency by means of reactive
power control. This approach is based on the Model Predictive
Control (MPC) theory. According to an economic criterion and
control effectiveness, a control switching strategy consisting of
a sequence of amounts of the shunt capacitors to switch is
identified for voltage restoration. The effect of the capacitive
control on voltage recovery is measured via trajectory sensitivity.
This approach is applied to the WECC system to enhance the
performance of voltage and to the 39 bus New England system
for preventing voltage collapse.
Index Terms—Model predictive control, trajectory sensitivity,
voltage stabilization, switching control, power system
I. INTRODUCTION
As a result of deregulation as well as increasing demands,
power systems operate close to their capacity. Although power
systems are designed with proper planning and with proper
stability margin, the instability can still occur under certain
severe disturbances. It is imperative that schemes for power
system protection be in place to mitigate their catastrophic
effects such as large scale shutdowns and collapses. The
objective of SPSs is to detect a potential instability or a
safety/security degradation of a power system and carry out
the necessary control actions to mitigate their effects (such as
a partial shutdown or a total collapse).
The traditional SPS is determined off-line and is rule based
[1], [2], [3]. A rule based system protection scheme relies on
voltage, or their rate of change levels, or line flow limits. For
example, if the measured voltage is lower than a specific value,
or the line flow exceeds the line rating limit, a predefined
SPS is triggered (such as adjustment of generator outputs or
load shedding). The limitation of the rule based SPSs lies in
the use of limited local information. In contrast, a real time
SPS computes and carries out control actions based on global
The work was supported in part by the National Science Foundation under
Grants NSF-ECCS-0424048, NSF-ECCS-0601570, NSF-ECCS-0801763, and
NSF-CCF-0811541.
L. Jin is with the Department of Electrical and Computer Engineering of
Iowa State University, Ames, IA 50011, USA and also with California ISO,
Folsom, CA, 95630 USA (email:ljin@caiso.com)
R. Kumar is with the the Department of Electrical and Com-
puter Engineering of Iowa State University, Ames, IA 50011, USA
(email:rkumar@iastate.edu)
N. Elia is with the the Department of Electrical and Computer Engineering
of Iowa State University, Ames, IA 50011, USA (email:nelia@iastate.edu)
state information in response to an impending contingency
detected by an online dynamic security assessment program.
Recent advances in monitoring, communication, and comput-
ing technologies have greatly facilitated the implementation of
real-time SPSs [4].
A real time system protection scheme for voltage stabiliza-
tion is studied in this work. The control of voltage level is
accomplished by controlling the production, absorption, and
flow of reactive power at various locations in the system. With
regard to a power system, sources and/or sinks of reactive
power, such as shunt capacitors, shunt reactors, synchronous
condensers, and static var compensators (SVCs) are used to
control voltage level. In literature, many algorithms [5], [6], [7]
have been developed to determine the amounts and locations
of shunt reactive power compensation devices needed for
maintaining a satisfactory voltage profile, while minimizing
their cost.
Most these work however are based on static analysis,
which means that the voltage performance criteria could be
met only if the system reaches a post-contingency stable
operating point. However, if the disturbances are sever, the
power system may lose stability. Under this situation, the
control strategy to restore the stable equilibrium point requires
a dynamic analysis.
Model predictive control has been applied in power system
voltage control based on dynamic analysis. [8] presents a
method of coordination of load shedding, capacitor switching
and tap changers using model preventive control. The predic-
tion of states is based on the numerical simulation of nonlinear
differential algebraic equations (DAEs) together with Euler
state prediction. A tree search method is adopted to solve the
optimization. [9] proposes a coordination of generator voltage
setting points, load shedding and ULTCs using a heuristic
search and the predictive control. The prediction of states is
based on the linearization of nonlinear DAEs. [10] presents
an optimal coordinated voltage control using model predictive
control. The controls used include: shunt capacitors, load
shedding, tap changers and generator voltage setting points.
The prediction of voltage trajectory is based on the Euler state
prediction. The optimization problem is solved by a pseudo
gradient evolutionary programming (PGEP) technique. In [11]
and [12], authors present a method to compute a voltage
emergency control strategy based on model predictive control.
The prediction of the output trajectories is based on trajectory
sensitivity. However, in these two papers, the authors employ a
simplified model predictive control, which computes the con-
trol actions only at the initial time and implements it over the
entire control horizon. A voltage stabilization control strategy
2
is also proposed in [13] based on load shedding, where the
objective function is to minimize the amount of load shedding
required to restore the voltages. It shows load shedding is
an effective voltage control under emergency condition. [14]
presents a MPC based voltage control design. The controls
are reference voltage of automatic voltage regulators and load
shedding.
In this paper, we propose computation of the optimal
strategies based on model predictive control (MPC). We utilize
shunt capacitors for control purposes as they are effective
means of voltage stabilization. The problem then becomes
one of determining capacitor switching sequence and amounts
given their locations and limits which are determined in a
prior planning stage (see for example [15]), together with the
requirements on the magnitudes of voltages. In this work, the
trajectory deviation and the cost of controls are simultaneously
minimized. Here, trajectory deviation refers to the deviation
of voltage trajectory from the nominal value. This is a multi-
objective optimization and a positively weighted convex sum
is chosen as the objective function. Trajectory sensitivities
are used to estimate the effect of controls on the voltage
behavior in a linear manner. Due to the use of model predictive
approach, the influence of each optimization is limited to one
step and the control gets recalculated and refined at each step,
the overall control strategy turns out to be sound and robust.
The features of our work, compared to the prior works dealing
with dynamic analysis, is summarized as follows:
• Trajectory sensitivity is used to compute a 1st-order
(linear) approximation of the effect of control without
having to linearize the system model. Further at each
control step, the trajectory sensitivity is updated (based
on a prediction of system trajectory starting from an
estimate of the current state under the control applied
in the past steps). This way of computing the effect
of control provides a better approximation as compared
to [8], [9], [10], where either system linearization or
numerical simulation of DAEs was used.
• Optimization minimizes costs of control as well as
voltage-deviations. In contrast [13], only considers the
amount of controls to restore the voltage.
• Optimization at each control step is a quadratic pro-
gramming problem, and hence can be efficiently solved.
In contrast [10] uses a pseudo gradient evolutionary
programming. [8], [9] use a tree search method.
• In contrast to [11] and [12], where the control action is
calculated only at the initial time, and remains the same
over the entire control horizon, a sequence of control
inputs is determined and only the first of them is applied
in our case.
• [14] presented a voltage control strategy based on ref-
erence voltage of automatic voltage regulators and load
shedding. [13] applied MPC for load-shedding computa-
tion. This paper presents an MPC-based shunt capacitor
control for voltage stabilization, a commonly used control
mechanism in North America power grid.
II. BACKGROUND
A. Model predictive control
Model Predictive Control (MPC) refers to a class of al-
gorithms that compute a sequence of manipulated variable
adjustments in order to optimize the future behavior of a
plant. An introduction to the basic concepts of MPC and a
formulation can be found in [16]. The principle of MPC is
graphically depicted in Fig. 1. Here x represents the state
variable that needs to be controlled to a specific range. The
available control is represented by variable u.
Predicted state
Manipulated input
State trajectory
until current time, x
k
Time
Magnitude
Input until
current time, u
k
t
k
t
k
+ T
s
t
k
+ T
c
t
k
+ T
p
Control horizon Tc
Prediction horizon Tp
Fig. 1. Principle of MPC
At a current time t
k
, the MPC solves an optimization
problem over a finite prediction horizon [t
k
, t
k
+ T
p
] with
respect to a predetermined objective function such that the
predicted state variable ˆx(t
k
+ T
p
) can optimally stay close to
a reference trajectory. The control is computed over a control
horizon [t
k
, t
k
+ T
c
], which is smaller than the prediction
horizon (T
c
≤ T
p
). If there were no disturbances, no model-
plant mismatch and the prediction horizon is infinite, one
could apply the control strategy found at current time t
k
for
all times t ≥ t
k
. However, due to the disturbances, model-
plant mismatch and finite prediction horizon, the true system
behavior is different from the predicted behavior. In order to
incorporate the feedback information about the true system
state, the computed optimal control is implemented only until
the next measurement instant (t
k
+ T
s
), at which point the
entire computation is repeated.
In a MPC, the optimization problem to be solved at time
t
k
can be formulated as follows:
min
ˆu
Z
t
k
+T
p
t
k
F (ˆx(τ), ˆu(τ))dτ (1)
subject to
˙
ˆx(τ) = f(ˆx(τ), ˆu(τ)), ˆx(t
k
) = x(t
k
) (2)
u
min
≤ ˆu(τ) ≤ u
max
, ∀τ ∈ [t
k
, t
k
+ T
c
] (3)
3
ˆu(τ) = ˆu(t
k
+ T
c
), ∀τ ∈ [t
k
+ T
c
, t
k
+ T
p
] (4)
x
min
(τ) ≤ ˆx(τ ) ≤ x
max
(τ), ∀τ ∈ [t
k
, t
k
+ T
p
] (5)
Here, T
c
and T
p
are the control and prediction horizon with
T
c
≤ T
p
. ˆx denotes the estimated state and ˆu represents
“estimated” control (The true state may be different and the
true control matches the estimated control only during the first
sampling period).
Equation (1) represents the cost function of the MPC opti-
mization. Equation (2) represents the dynamic system model
with initial state x(t
k
). Equations (3) and (4) represent the
constraints on the control input during the prediction horizon.
Equation (5) indicates the state operation requirement during
the prediction horizon.
B. Trajectory sensitivity
Consider a differential algebraic equation (DAE) of a sys-
tem,
˙x = f(x, y, u), x(0) = x
0
(6)
0 = g (x, y, u) (7)
where x is a vector of state variables, y is a vector of algebraic
variables, and u is a vector of control variables. Trajectory
sensitivity considers the influence of small variations in the
control u (and any other variable of interest) on the solution
of the state equations (6) and (7). Let u
0
be a nominal value
of u, and assume that the nominal system in (8) and (9) has
a unique solution x(t, x
0
, u
0
) over [t
0
, t
1
].
˙x = f(x, y, u
0
), x(0) = x
0
(8)
0 = g (x, y, u
0
) (9)
Then the system in Equations (6) and (7) has a unique solution
x(t, x
0
, u) over [t
0
, t
1
] that is related to x(t, x
0
, u
0
) as:
x(t, x
0
, u) = x(t, x
0
, u
0
) + x
u
(t)(u − u
0
) + H.O.T.(10)
y(t, x
0
, u) = y(t, x
0
, u
0
) + y
u
(t)(u − u
0
) + H.O.T.(11)
Here x
u
(t) =
∂x(t,x
0
,u)
∂u
is called the trajectory sensitivities
of state variables with respect to variable u and y
u
(t) =
∂y(t,x
0
,u)
∂u
is the trajectory sensitivities of algebraic variables
with respect to variable u.
The evolution of trajectory sensitivities can be obtained by
differentiating Equations (6) and (7) with respect to the control
variables u and is expressed as:
˙x
u
(t) = f
x
(t)x
u
(t) + f
y
(t)y
u
(t) + f
u
(t) (12)
0 = g
x
(t)x
u
(t) + g
y
(t)y
u
(t) + g
u
(t) (13)
Detailed information about trajectory sensitivity theory can
be found in [17]. The trajectory sensitivity can be solved
numerically. [18] provides a methodology for the computation
of trajectory sensitivity. When time domain simulation of a
power system is based on trapezoidal numerical integration,
the calculation of trajectory sensitivity requires solving a set
of linear equations, thus costing a little time. In our work, we
extended the Power System Analysis Tool [19] (a MATLAB
based tool) to do trajectory sensitivity calculation and the MPC
optimization.
Fig. 2 illustrates the application of trajectory sensitivity
in evaluating the effect of controls on system behavior. The
trajectory x
k
of the nominal system represents the behavior
under the control u
k
. When the control is increased by ∆u
k
1
at time t
k
, the change in predicted system behavior based
on sensitivity analysis at time t
l
, can be approximated as
∆x
kl
1
= x
l
u
k
1
∆u
k
1
. Here x
l
u
k
1
is the trajectory sensitivity of
the state variable at time t
l
with respect to the control at
time t
k
. Similarly if we increase the control by ∆u
k
n
at time
t
k
+ (n − 1)T
s
, the change in the state variable at time t
l
is
represented by ∆x
kl
n
= x
l
u
k
n
∆u
k
n
. Here, x
l
u
k
n
is the trajectory
sensitivity of the state variable at time t
l
with respect to the
control at time t
k
+ (n − 1)T
s
.
Magnitude
Step 1
Time
Step n
k
u
1
∆
k
n
u∆
kl
u
kl
uxx
k
11
1
∆=∆
k
n
l
u
kl
n
uxx
k
n
∆=∆
k
x
k
u
k
t
sk
Tt +
sk
Tnt )1( −+
l
t
Nominal system
trajectory
Fig. 2. Application of trajectory sensitivity in system behavior prediction
III. PROBLEM FORMULATION AND SOLUTION
The purpose of this work is to find an effective and eco-
nomic control strategy for controlling the shunt capacitors so
as to eliminate voltage instability following any pre-identified
contingency. For analyzing voltage performance following dis-
turbances, we model generator and automatic voltage regulator
(AVR) as well as aggregated exponential dynamic load models
[20], [21]. The overall power system is represented by a set
of differential algebraic equations (DAE) as in Equations (6)
and (7). Here x is a vector of states including state variables
in generator dynamic models, AVR models and dynamic load
models such as, rotor angles and angular speeds of generators,
outputs of AVRs, and active power recovery and reactive
power recovery of dynamic load models. y is a vector of
algebraic variables such as bus voltage magnitudes and phase
angles. The vector u indicates the output of shunt capacitors.
The computation is iterative over a finite control horizon,
where in each step a quadratic programming problem is solved
to compute the amounts of shunt capacitors to be added in that
step. The quadratic programming formulation is valid when
the capacitor control is continuous as in SVC. Even in the
case where capacitor control is discrete, we can still proceed
by assuming continuous control so as to compute an optimal
control by solving a quadratic programming relaxation. Then
for implementation, the nearest discrete control value can be
applied. Any error will get propagated to a following control
step, and where it will get corrected. The control is piecewise
4
constant, changing only at the sampling times. Let T
p
be the
prediction horizon, T
c
be the control horizon, T
s
be the control
sampling interval, and N =
T
c
T
s
be the total number of control
steps. The procedure to determine the control strategy at time
t
k
based on MPC is as follows:
Step 1: At time t
k
(i.e. the (k + 1)
th
sampling in-
stant), an estimate of the current state x(t
k
) is
obtained. The nominal power system evolves ac-
cording to Equations (6) and (7). Here, u =
{B
0
m
+
P
k−1
i=0
∆B
i
m1
}
m=M
m=1
is the control variable
(i.e. amounts of shunt capacitors currently in use).
B
0
m
is the amounts of shunt capacitors that exist at
time 0.
P
k−1
i=0
∆B
i
m1
is the amounts of shunt capac-
itors that were added over time [0, t
k
− T
s
]. Time
domain simulation is used to obtain the trajectory
of the nominal system (6) and (7), starting from the
state x(t
k
) at time t
k
to the end of prediction horizon
t
k
+T
p
. At the same time, the trajectory sensitivity of
bus voltages with respect to the shunt capacitors to
be added at instants t
k
+ (n − 1)T
s
, n = 1 . . . N − k
is obtained and denoted as V
kj
B
mn
(t) (see below for
the explanation of notation).
Step 2: At time t
k
, solve the optimization problem over
the prediction horizon [t
k
, t
k
+ T
p
] and the control
horizon [t
k
, t
k
+ T
c
] as stated in (14)-(18). The
objective function is composed of two parts. The
first term is the trajectory deviation, the second
term is the cost of controls. The combination of
the deviation of voltages from nominal values
and the control cost needs to be minimized. The
number of candidate control locations and their
upper limits are determined through a prior planning
step (see for example [15]). The total number
of control variables in the optimization is the
number of candidate control locations times the
number of control steps. The optimization is solved
in Matlab, and it does converge to a global minimum.
Minimize (with respect to ∆B
k
mn
)
Z
t
k
+T
p
t
k
(
b
V
k
(t)−V
ref
)
0
R(
b
V
k
(t)−V
ref
)dt+
X
mn
W
mn
∆B
k
mn
(14)
Subject to
∆B
min
m
≤ ∆B
k
mn
≤ ∆B
max
m
(15)
B
min
m
≤ B
0
m
+
k−1
X
i=0
∆B
i
m1
+
N
X
n=1
∆B
k
mn
≤ B
max
m
(16)
V
kj
min
(t) ≤ V
kj
(t)+
M
X
m=1
min(N ,l)
X
n=1
V
kj
B
mn
(t)∆B
k
mn
≤ V
kj
max
(t)
(17)
∆B
k
mn
≥ 0 (18)
• R is the weight matrix.
b
V
k
(t) is the predicted voltage vector
at the control sampling time t
k
that contains all the bus
voltages in the system at time t. ∆B
k
is the control matrix
calculated at time t
k
.
• W
mn
is the weighted cost of control m to be added at time
t
k
+ (n − 1)T
s
.
• M is the total number of control variables, i.e. the number
of shunt capacitor locations.
• N is the total number of control steps.
• ∆B
k
mn
is the entry ∆B
k
, which is the amount of control
m to be added at time t
k
+ (n − 1)T
s
.
• ∆B
min
m
∈ < is the minimum amount of control m to be
added at any step.
• ∆B
max
m
∈ < is the maximum amount of control m to be
added at any step.
• ∆B
i
m1
is the amount of control m implemented at the
control sampling point t
i
, i = 0, ...k − 1.
• B
min
m
∈ < is the minimum amount of control m that must
be used, typically 0.
• B
max
m
∈ < is the maximum available amount of control m.
• V
kj
(t) ∈ < is the voltage of bus j at time t(t
k
≤ t ≤
t
k
+ T
p
), of the nominal system of time t
k
.
• V
kj
min
(t) is the minimum voltage at bus j desired at time
t
k
≤ t ≤ t
k
+ T
p
.
• V
kj
max
(t) is the maximum voltage at bus j desired at time
t
k
≤ t ≤ t
k
+ T
p
.
• V
kj
B
mn
(t) is the trajectory sensitivity of voltage at bus j at
time t
k
≤ t ≤ t
k
+ T
p
with respect to control m added at
time t
k
+ (n − 1)T
s
.
Step 3: At time t
k
, a solution of the optimization problem
(14)-(18) computes a sequence of controls ∆B
k
mn
.
Add only the first control ∆B
k
m1
at time t
k
and
observe or estimate the system state x(t
k+1
) at time
t
k+1
= t
k
+ T
s
Step 4: Increase k by k + 1 and repeat steps (1)-(3) until
the k = N − 1.
A. Implementation
The functional structure of a real time SPS is shown in
Figure 3. Line flow, bus voltage information, switch status
as well as phase measurement unit (PMU) measurements are
sent to a control center through communication channels of a
SCADA system. These measurements plus a network model
are used by the state estimator (SE) for filtering out the
noise and making best use of the measured data. The results
from the state estimator are used for power flow analysis.
A power flow solution is then used by an on-line dynamic
security assessment program to initialize the state variables
of the dynamic models. Further, it uses system models and
disturbance information to perform the contingency analysis
to evaluate the security margin of the power system. If
a contingency is identified where the system will become
unstable, the MPC based SPS computation will get triggered
at the time an identified critical contingency occurs. The steps
of the MPC computation in the k
th
iteration include:
• Estimate static variables such as voltage magnitudes and
angles at time t
k
as well as the dynamic variables x(t
k
)
such as generator angles, velocities and real and reactive
load recovery.
• Run time-domain simulation to compute the system tra-
jectory given the current state. This step also requires the
knowledge of a complete system model (including both
dynamic and static components).
• Obtain trajectory sensitivities of voltage with respect
to the control variables as a by-product of the time-
domain simulation performed in the previous step. This
5
PMU
measurements
Power and switch
measurement
State estimator
Power flow
analysis
System model
Disturbance
Recorder
On-line dynamic
security program
System
secure?
MPC based SPS
controller
Disturbance
happens?
End
Yes
No
Yes
No
Control signal to
power system
Fig. 3. Structure of a real time SPS
is required for the prediction of system response given a
certain control strategy.
• Solve the quadratic programming optimization problem
and implement the first step of the control.
• Repeat the above steps at each sampling point until the
end of control horizon.
Remark: While we suggested an on-line computation of
MPC based SPS above, it is also possible to do this computa-
tion off-line based on the predicted (rather estimated) values
of the states and trajectory sensitivities.
IV. APPLICATION TO WECC AND TO NEW ENGLAND
SYSTEMS
The proposed method has been applied to the WECC 9-bus
system as well as to the New England 39-bus system. The
exponential recovery load model is used in both cases. The
parameters of the load model are as following:
T
p
= T
q
= 30, α
s
= 0, α
t
= 1, β
s
= 0, β
t
= 4.5.
The parameters in MPC optimization are determined based
on the following considerations. Any voltage instability fol-
lowing a contingency must be stabilized in a certain time
duration (typically the time in which voltage will decrease by
15%). This is the prediction horizon T
p
. The control should
be exercised on a time horizon T
c
, which is shorter than the
prediction horizon, typically the time in which voltage will
decrease by 10% (if no control is applied). A discrete-time
control must be applied within this duration T
c
at a sample-
rate high enough to adequately react to the changing voltage
trajectory, as well as to allow accurate enough predictions
of the voltage trajectory based on the linearization of the
trajectory-sensitivity. This dictates the sampling duration T
s
.
The number of sampling point N is then determined as the
ratio of T
c
and the sampling duration T
s
.
A. WECC 3-generator 9-bus test system
1) System description: Figure 4 is a representation of the
WECC 3-generator 9-bus system. A fourth-order model is
used for modeling each of the three generators. The state
2
7
8
9
5
6
3
4
1
Gen 2
Gen 3
Gen 1
Fig. 4. WECC 3-generator 9-bus test system
variables include the rotor angle δ, the rotor speed ω, the
q-axis transient voltage e
0
q
, and the d-axis transient voltage
e
0
d
. Automatic Voltage Regulator (AVR) defines the primary
voltage regulation for generator 1. The continuously acting
regulator and exciter model [22] is employed in this study. It is
represented by a four-dimensional state equation. The loads at
buses 5, 6 and 8 are taken to be exponential recovery dynamic
load and each load is described by a two-dimensional state
equation. Therefore, the total dimension number of the state
space is 22. At buses 5, 7 and 8, there exist shunt capacitors
for voltage regulation. These are the control variables. Under
normal conditions, all of the shunt capacitors are disconnected.
2) Fault scenario: We consider a three-phase fault at bus
5 at t = 1.0 second, which is cleared at t = 1.2 seconds by
the tripping of the line between bus 4 and bus 5. Based on
the time domain simulation, the voltages at buses 5, 7 and
8 are shown in Fig. 5 and are not satisfactory. At t = 1.0
second, the voltages begin to drop dramatically due to the three
phase to ground fault. At t = 1.2 seconds, the voltages start
to recover since the fault gets cleared. However, the voltages
begin to oscillate. 15 seconds later, voltages begin to decline
gradually. The dynamic load models result in slightly recovery
load consumption, which deteriorate the voltage condition.
These three voltages fall out of the lower limit 0.95 p.u 1
minutes later. According to the system’s operational criteria,
the load bus voltages must be above 0.95 p.u. Therefore, some
control actions are required to satisfy the criterion that the
voltages outlined above remain above 0.95 p.u.
3) Simulation result: In this example, we have chosen
prediction horizon T
p
to be 40 seconds (the time in which
voltage drops by nearly 15% at bus 5). T
c
has been chosen
to be 35 seconds. We found that a sampling duration of T
s
=
7 seconds works well for this example, and so we have the
number of control steps: N =
T
c
T
s
=
35
7
= 5. Model predictive
control approach determines the amounts of shunt capacitors
to be added at each sampling instant so as to recover the local
voltages. Although the capacitors have a positive effect on low
voltage problems, the maximum capacitor to be added at any
step ∆B
max
m
was set to be 0.1 p.u. This is because if large
amounts of capacitors are added at one time, an over-voltage
may occur, which has a bad effect on the electrical devices of
