1

Coupling Load-Following Control with OPF

arXiv:1707.03794v3 [math.OC] 17 Jan 2018

Mohammadhafez Bazrafshan, Student Member, IEEE,
Nikolaos Gatsis, Member, IEEE, Ahmad F. Taha, Member, IEEE, Joshua A. Taylor, Member, IEEE

Abstract—In this paper, the optimal power flow (OPF) problem
is augmented to account for the costs associated with the loadfollowing control of a power network. Load-following control
costs are expressed through the linear quadratic regulator (LQR).
The power network is described by a set of nonlinear differential
algebraic equations (DAEs). By linearizing the DAEs around a
known equilibrium, a linearized OPF that accounts for steadystate operational constraints is formulated first. This linearized
OPF is then augmented by a set of linear matrix inequalities
that are algebraically equivalent to the implementation of an
LQR controller. The resulting formulation, termed LQR-OPF,
is a semidefinite program which furnishes optimal steady-state
setpoints and an optimal feedback law to steer the system to
the new steady state with minimum load-following control costs.
Numerical tests demonstrate that the setpoints computed by
LQR-OPF result in lower overall costs and frequency deviations
compared to the setpoints of a scheme where OPF and loadfollowing control are considered separately.
Index Terms—Optimal power flow, load-following control,
linear quadratic regulator, semidefinite programming.

I. I NTRODUCTION
APACITY expansion and generation planning, economic
dispatch, frequency regulation, and automatic generation control (AGC) are decision making problems in power
networks that are solved over different time horizons. These
problems range from decades in planning to several seconds
in transient control. Although decisions made in shorter time
periods may negatively affect performance over longer time
periods, these problems have traditionally been treated separately. For example, optimization is used for optimal power
flow (OPF) and economic dispatch while feedback control
theory is used for frequency regulation.
This paper aims to integrate two crucial power network
problems with different time-scales. The first problem is the
steady-state OPF, whose decisions are updated every few
minutes (e.g., 5 minutes for real-time market balancing). The
second is the problem of load-following control that spans the
time-scale of several seconds to one minute. Load-following
control, also known as secondary frequency control [1], is
responsible for maintaining system frequencies at nominal
values during normal load fluctuations [2].
For a forecasted load level, optimal steady-state setpoints
that also minimize costs of load-following control are sought.
Steady-state costs account for generator power outputs. Control costs account for the action required to drive the deviation
of frequency and voltage signals from their optimal OPF
setpoints to zero via the linear quadratic regulator (LQR). The
proposed formulation also provides, as an output, a feedback
law to guide the system dynamics to optimal OPF setpoints.

C

This material is based upon work supported by the National Science
Foundation Grants ECCS-1462404 and CCF-1421583.

A. Literature Review
Recent research efforts, organized below in two categories,
have showcased the economic and technical merits of jointly
tackling steady-state and control problems in power networks.
The first category focuses on frequency regulation with a
view towards the economic dispatch [3]–[10]. These works
design controllers based on feedback from frequency measurements and guarantee the stability of system dynamics—
characterized by the swing equation—while ensuring that
system states converge to steady-state values optimal for some
form of the economic dispatch problem.
The designs of these feedback control laws may come from
averaging-based controllers as in [3]–[5]. By leveraging a
continuous-time version of a dual algorithm, generalized versions of the aforementioned controllers are developed in [6],
including the classical AGC as a special case. Similar control
designs come from interpretation of network dynamics as
iterations of a primal-dual algorithm that solves an OPF. For
example, [7] introduces a primary frequency control that minimizes a reverse-engineered load disutility function, and [8]
introduces a modified AGC that solves an economic dispatch
with limited operational constraints.
A frequency control law that further solves economic dispatch with nonlinear power flows in tree networks is devised in
[9] by borrowing virtual dynamics from the KKT conditions. A
novel method to tackle a load-side control problem including
linear equality or inequality constraints for arbitrary network
topologies is presented in [10]. The chief attractive feature
of [3]–[10] is that they afford decentralized or distributed
implementations. Moreover, works such as [6] and [9] develop
controllers accounting for nonlinear power flow equations.
The second category focuses on OPF variations with enhanced stability measures. The goal is to obtain optimal
steady-state setpoints less vulnerable to disturbances, rather
than seeking a stabilizing control law. This goal is achieved
by incorporating additional constraints in OPF that account for
the stability of the steady-state optimal point.
In the context of transient stability, a classical reference
is [11] where system differential equations are converted
to algebraic ones and added to the OPF. More recently,
trajectory sensitivity analysis has been extensively used to
accommodate transient stability specifications in the form of
linear constraints within the OPF. Typically, these constraints
are based on stability margins obtained from the (extended)
equal-area criterion [12]–[15]. In particular, [12] augments
the OPF with automated computationally-efficient rotor angle
constraints around a base stable trajectory. Transient-stability
constrained OPF formulations that are robust to uncertainty
in load dynamics and wind power generation have also been

2

developed [13], [14]. Load-shedding minimization has been
pursued in [15], where in addition to angle stability margins,
trajectory sensitivities have been used to approximate constraints on voltage and frequency security margins.
In the realm of small-signal stability, distance from rotor
instability is guaranteed in [16] by providing a set of stressed
load conditions as a supplement to OPF. The spectral abscissa,
that is, the largest real part of system state matrix eigenvalues,
is upper bounded by a negative number in [17], yielding a
non-smooth OPF. A simpler optimization problem is pursued
in [18] using the pseudo-spectral abscissa as stability measure.
Based on Lyapunov’s stability theorem and the system state
matrix, [19] incorporates small-signal stability constraints into
the OPF. Since the state matrix is a function of the steadystate variables, the overall formulation becomes a nonlinear
and nonconvex semidefinite program (SDP).
B. Paper contributions and organization
The paper contributions are as follows:
• An OPF framework is proposed that solves for optimal
steady-state setpoints and also provides an optimal loadfollowing control law to drive the system to those setpoints.
Optimality of the control law is appraised by an integral cost
on time-varying deviation of system states and controls from
their optimal setpoints. An LQR controller is then applied
which minimizes this cost by providing a feedback law that
is a linear combination of system state deviations. LQR has
previously been used in the context of megawatt-frequency
control [20] and control of oscillatory dynamics [21]. The
proposed framework, in distinction, accounts for LQR costs
from within the steady-state time-scale.
• The setpoints computed by the proposed formulation,
termed LQR-OPF, result in lower overall costs and frequency fluctuations compared to the setpoints of a scheme
where OPF and load-following control are decoupled. The
proposed framework also allows time-varying control costs
to be dependent on steady-state variables. This dependence
enables subsequent regulation pricing schemes similar to
[22] where for example, the cost of frequency regulation
is made dependent on steady-state power generation.
• In comparison to [3]–[10], the formulation in this paper includes voltage dynamics, reactive powers, AC power flows,
and a more realistic dynamical model of the synchronous
generator that distinguishes between generator internal and
external quantities. Costs incurred due to deviations of
system states and controls from their optimal setpoints are
also accounted for. In relation to [11]–[19], the proposed approach incorporates load-following control constraints into
the OPF; but as an extra output, the required control law to
steer the system to stability is also provided.
• The effectiveness of the derived control law is demonstrated
via numerical simulations on the power system described
by nonlinear differential-algebraic equations (DAEs). We
further exhibit that the steady-state setpoints provided by
the proposed LQR-OPF yield smaller overall system costs
even when the traditional AGC is used for load following.
In this paper, in order to obtain a stabilizing feedback law,
nonlinear DAEs are linearized around a known equilibrium

point with respect to the system states and algebraic variables.
It should be emphasized that this linearization is different
than that of trajectory sensitivity analysis. In the latter, the
system dynamics are linearized around a known trajectory
with respect to system initial conditions and parameters [23],
[24]. Formulations in [12]–[14] that use trajectory sensitivity
analysis are suitable for transient stability where it is required
to limit generators’ angle separation subject to contingencies
and large disturbances, e.g., line trips due to three-phase
faults. On the other hand, the load-following or secondary
control, which is considered in this paper, aims at maintaining
frequencies at their nominal value during normal load changes.
The paper is organized as follows. The power system model
is laid out in Section II, followed by a description of a generalized OPF. System linearization is pursued in Section III.
The proposed formulation coupling OPF and load-following
control is detailed in Section IV. Specific generator models,
power flow equations, and connections to the standard OPF
are provided in Section VI. Section VII numerically verifies
the merits of the proposed method. Section VIII provides
pointers for integrating more power system applications into
our proposed framework as future work.
II. P OWER S YSTEM M ODEL
Consider a power network with N buses where N :=
{1, . . . , N } is the set of nodes. Define the partition N = G ∪L
where G = {1, . . . , G} collects G buses that contain generators (and possibly also loads) and L = {G + 1, . . . , G + L}
collects the remaining L load-only buses. Notice that N =
G + L. For a generator i ∈ G with ns states and nc control
inputs, denote by xi (t) ∈ Rns the time-varying vector of state
variables, and denote by ui (t) ∈ Rnc the time-varying control
inputs. For example, adopting a fourth-order model yields
ns = 4 and nc = 2 with xi (t) = {δi (t), ωi (t), ei (t), mi (t)}
and ui (t) = {ri (t), fi (t)} where δi (t), ωi (t), ei (t), mi (t),
ri (t), and fi (t) respectively denote the generator internal
phase angle, rotor electrical velocity, internal EMF, mechanical
power input, reference power setting, and internal field voltage.
Further details are given in Section VI-A.
Denote by ai (t) the vector of algebraic variables. For load
nodes i ∈ L, ai (t) = {vi (t), θi (t)}, where vi (t) and θi (t)
denote the terminal load voltage and phase angle. For generator nodes i ∈ G, ai (t) = {pgi (t), qgi (t), vi (t), θi (t)}, where
pgi (t), qgi (t), vi (t), and θi (t) respectively denote generator
real and reactive power, terminal voltage and phase angle.
For brevity, the dependency of variables xi , ai , and ui on
t is dropped, and the notations x := {xi }i∈G ∈ Rns G ,
u := {ui }i∈G ∈ Rnc G , and a := {ai }i∈N ∈ R2N +2G are
introduced. Finally, let z = {x, a, u} ∈ R(ns +nc +2)G+2N .
The dynamics of a power system can be captured by a set of
nonlinear DAEs
ẋ = g(x, a, u),

(1a)

d = h(x, a),

(1b)

where g : R(ns +nc +2)G+2N → Rns G is given by adopting an appropriate dynamical model of the generator, and
h : R(ns +2)G+2N → R2G+2N includes generator algebraic

3

equations as well as the network power flow equations. Vector
d collects all the network loads as well as leading zero entries
coming from two generator algebraic equations per generator.
A particular example of the mapping g is provided in Section
VI-A; for the corresponding form of the mapping h and the
vector d see Section VI-B.
Given steady-state load conditions deq , the system steadystate operating point is represented by an equilibrium of the
DAEs (1). By setting ẋ = 0 and allowing x, a, and u to reach
steady states xeq , aeq , and ueq , a system of (ns + 2)G + 2N
algebraic equations in (ns +nc +2)G+2N variables is derived:
0 = g(xeq , aeq , ueq ),
d = h(xeq , aeq ).
eq

(2a)
(2b)

Let F (deq ) denote the set of solutions to (2), where the dependency on the load conditions is made explicit. Suppose now
that xeq , aeq , ueq are to be jointly optimized so that a certain
objective function c(xeq , aeq , ueq ) is minimized. This leads to
a generalized OPF [for clarity, the notation (xs , as , us ) is used
to denote optimization variables, and (xeq , aeq , ueq ) is used
to generically denote a DAE equilibrium].
c (xs , as , us )

min

xs ,as ,us

subj. to

s

s

s

s

(3a)
s

0 = g(x , a , u ),

(3b)

s

(3c)
(3d)

d = h(x , a ),
as ∈ A,

where (3d) are the algebraic variable constraints on voltage
magnitudes, line flow limits, and line current capacities. Note
that the parameter vector ds (which includes the constantpower loads) is the input to (3). The term generalized refers
to the fact that (3) considers models of generators within the
OPF, see e.g., [25] for a recent OPF example with machine
models. The connection between (3) and the standard OPF is
explained in Section VI-C. The OPF problem (3) guarantees
optimal steady-state operating costs, but does not provide
minimal control costs. Prior to introducing a formulation that
bridges stability with OPF, linear approximation of the system
dynamics is required which is presented in the next section.
III. L INEAR A PPROXIMATION

OF

S YSTEM DYNAMICS

To obtain an approximate dynamic, (1) is linearized around
a known operating point z0 := (x0 , a0 , u0 ) ∈ F (d0 ). For
example, the point z0 can be a solution of the load-flow
corresponding to an operating point known to the system
operator. The motivation behind this selection is to obtain
tractable constraints to augment (3) as it allows the stability
constraints to take the form of properly formulated linear matrix inequalities (LMIs). The derived control law is eventually
applied to the nonlinear DAEs (1), rather than the linearization
derived in this section.
Consider a generic equilibrium point z0 ∈ F (d0 ), where
0
d is a known load vector. That is, the following holds:
0 = g(x0 , a0 , u0 ),

(4a)

0

(4b)

0

0

d = h(x , a ).

Define ∆ds := ds − d0 as the step-change difference
between ds , the load for which the generalized OPF (3) is
to be solved, and d0 , the generic load that will be used
for linearization. Equations (1a) and (1b) can be linearized
around (x0 , a0 , u0 ) by setting (x, a, u) = (x0 , a0 , u0 ) +
(∆x, ∆a, ∆u):
∆ẋ = gx (z0 )∆x + ga (z0 )∆a + gu (z0 )∆u,
s

0

0

0

0

∆d = hx (x , a )∆x + ha (x , a )∆a,

(5a)
(5b)

where the notation gx defines the Jacobian with respect to
x, and ga , gu , hx , and ha are similarly defined. It is also
understood that ∆x, ∆a, ∆u are functions of time.
Equation (5) represents a set of linear DAEs. Next, (5) is
leveraged in order to develop a proper linear dynamical system
that represents the network without algebraic constraints. In
particular, assuming invertibility of ha (x0 , a0 ), (5b) can be
solved for the algebraic variables as

0
0
s
0
0
∆a = −h−1
(6)
a (x , a ) −∆d + hx (x , a )∆x .

The assumption on invertibility of ha (x0 , a0 ) is very mild in
the sense that it holds for practical networks and for various
operating points; see also [18] and references therein for sufficient conditions in a similar construction. Then, substituting
(6) in (5a) yields
0
0
s
∆ẋ = A(z0 )∆x + B(z0 )∆u + ga (z0 )h−1
a (x , a )∆d , (7)
0
0
0
0
where A(z0 ) = gx (z0 ) − ga (z0 )h−1
a (x , a )hx (x , a ), and
0
0
B(z ) = gu (z ).

IV. C OUPLING L OAD -F OLLOWING C ONTROL WITH OPF
In this section, an optimal control problem using the LQR
is presented to stabilize the linear dynamical system (7) to
an equilibrium point, while ensuring minimal steady-state
operating costs for the equilibrium of the linearized dynamics.
The equilibrium of the linearized dynamics is denoted by
zs = (xs , as , us ) to distinguish it from the equilibrium of
the true nonlinear system in e.g., (2). Section IV-A leverages the linearized dynamical system (7) to obtain a linear
approximation of the constraints (3b) and (3c)—this yields a
linearized OPF. The objective of Section IV-B is to incorporate
stability measures with respect to the dynamical system (7) in
the linearized OPF. The resulting formulation outputs optimal
steady-state values zs = (xs , as , us ) together with the optimal
control law that drives the dynamic variables to stability.
A. Linear approximation of the generalized OPF
The steady state of (5) obtained by setting ∆ẋ = 0 yields
the following system:


0 = gx (z0 ), ga (z0 ), gu (z0 ) (zs − z0 ),
(8a)



s
0
0
0
0
s
s
0
0
∆d = hx (x , a ), ha (x , a ) (x , a ) − (x , a ) (8b)

The equations in (8) are therefore a linear approximation of
constraints (3b) and (3c), which leads to a linearized version
of OPF, formulated as follows:
Linearized OPF :
min s c (xs , as , us ) subj. to (8) and as ∈ A.
s s
x ,a ,u

(9)

4

Notice that (9) has the linearized version of power flow
equations as part of its constraint set.

ds

LQR-OPF (14)

LQR

x̂

Observer

0

A(z )
B(z0 )

B. LQR-OPF formulation
The previous section derived (9), which is the linear approximation of the generalized OPF (3). Likewise, the linear system
in (7) is the linear approximation of the system dynamics (1a)
and (1b) around z0 . This section augments the linearized OPF
(9) with optimal control of the dynamical system in (7).
Writing (7) at its equilibrium zs , that is, setting ∆ẋ = 0,
∆x = ∆xs , and ∆u = ∆us , yields

xeq , ueq

Forecast

z0

Linearized
System (5)

K(x̂ − xeq )
z0

Nonlinear
System (1)

y

Fig. 1. Overall design of the proposed method. The LQR-OPF uses the
forecast information to determine the optimal setpoints for the system. The
feedback control law generated by the LQR drives the system to the optimal
steady-state with minimized control costs.

Matrices Q0 , Q1 , R0 , and R1 are selected so that Q(zs )
and
R(zs ) are positive definite over the bounded range of as
0 = A(z )∆x + B(z )∆u + ga (z )ha (x , a ) ∆d . (10)
prescribed by the set A.
Subtracting (10) from (7) yields the following system with the
The corresponding infinite horizon LQR augmenting the
new state variable ∆x′ := ∆x−∆xs and new control variable linearized OPF (12) can be written as an SDP, as follows:
∆u′ = ∆u − ∆us :
LQR-OPF:
∆ẋ′ = A(z0 )∆x′ + B(z0 )∆u′ .
(11)
Tlqr
min
c(zs ) +
γ
(14a)
s
s
s
s
2
S, z ={x ,a ,u }
The initial conditions of (11) are ∆x′ (0) = ∆x(0) − ∆xs =
subj. to (8a), (8b), as ∈ A
(14b)
x(0) − xs . When ∆x in (5) or (7) converges to ∆xs , ∆x′


s
0 ⊤
−γ
(x
−
x
)
in (11) converges to 0. The objective is to penalize the control
O
(14c)
xs − x0
−S
effort to drive ∆x′ to 0.

 ⊤
To this end, weights on the control actions and state
A S + SA
S
Y⊤ 
deviations are considered. In particular, an optimal LQR con+ BY + Y⊤ B⊤
  O (14d)

troller to drive the system states ∆x′ through a state-feedback


−1 s
S
−Q
(z
)
O
control ∆u′ to their zero steady-state values is considered. By
−1 s
Y
O
−R (z )
augmenting the linearized OPF in (9), this LQR-based OPF is
S  O.
(14e)
written as follows:
The optimal control law generated from the optimal conZtf
Tlqr
(12) is ∆u′ = K∆x′ where
K =
s
′⊤
′
′⊤
′
c(z )+
min
∆x Q∆x + ∆u R∆u dt (12a) trol −1problem
−1
∗
0
s
s s
s
2
z ={x ,a ,u }
−R (z )B(z )⊤ S∗ is the optimal state feedback control
0
∆x′ ,∆u′
gain and (S∗ , z∗ ) is an optimizer of (14). The reader is referred
′
0
′
0
′
subj. to ∆ẋ = A(z )∆x + B(z )∆u
(12b) to [26] for the derivation of this control law.
The LQR-OPF problem (14) is an SDP which solves jointly
∆x′ (0) = x(0) − xs
(12c)
s
for
the new steady-state variables zs and the matrix S while
(8a), (8b), a ∈ A,
(12d)
guaranteeing the stability of the linearized system in (7). The
where Q and R are positive definite matrices penalizing state two problems of OPF and stability control are coupled in (14)
and control actions deviations; Tlqr is a scaling factor to through the LQR weight matrices—with Q−1 and R−1 being
compensate for the time-scale of the stability control problem; affine functions—as well as through the term xs − x0 in
tf is the optimization time-scale for the OPF problem which is LMI (14c). It is worth noting that γ in (14) is equal to the
typically in minutes. Since tf is in minutes, the finite horizon integral in the objective of (12) with tf = ∞.
The optimal steady-state solution zs of the LQR-OPF
LQR problem can be replaced with an infinite horizon LQR
formulation (i.e., tf = ∞) as the solution to the Riccati problem (14) satisfies the linearized steady-state equations (8).
equation reaches steady state [26]. In (12c), if the system In order to obtain an equilibrium for the true nonlinear DAEs
is assumed to operate at x0 , then ∆x(0) = 0, yielding (2), generator setpoints from zs are extracted as follows. For
s
generator i ∈ G, set vieq = vis . If i is nonslack, set peq
∆x′ (0) = x0 − xs .
gi = pgi .
eq
s
Regulating state and control actions can be dependent on the If i is slack, set θi = θi . Then, the system of equations (2)
operating points to encourage smaller steady-state variations. is solved to obtain an equilibrium zeq for the nonlinear DAEs
For example, regulating a generator’s frequency becomes in (1). This task essentially amounts to solving two separate
more costly as the real power generation increases [22]. This sets of nonlinear equalities: first a standard load-flow is solved
coupling is captured by considering Q−1 and R−1 to be affine to obtain the remainder of algebraic variables qgeqi for i ∈ G,
eq
eq
functions of steady-state variable zs . In particular, assume vi for i ∈ L, peq
gi if i is slack, and θi for i nonslack. Second,
−1
−1
by incorporating the algebraic variables in generator equations,
Q and R to be diagonal matrices as follows
the equilibrium states and controls of the generators, that is xeq
−1
s
Q(zs )−1 = Q−1
+
Q
diag(a
),
(13a)
0
1
and ueq , are obtained. This process is detailed in Section VII
−1
−1
s −1
s
R(z ) = R0 + R1 diag(a ).
(13b) that includes simulations.
0

s

0

s

0

0

0 −1

s

5

Figure 1 demonstrates the overall design and the integration
of the proposed LQR-OPF into the power system operation.
Initially, the system is operating at a known equilibrium
point z0 . Using this known operating point, the state-space
matrices A(z0 ) and B(z0 ) are computed according to (7). This
information is then provided to the LQR-OPF block together
with a forecast ds of the demand at the next steady state. The
LQR-OPF calculates the optimal generator control setpoints
included in ueq as well as the optimal state setpoints for
the next steady state xeq . Upon providing the LQR block
with optimal setpoints xeq and ueq , the LQR control law
K(x−xeq )+ueq is then applied to the nonlinear power system
and drives the system to the new desired optimal point zeq .
The purpose of LQR-OPF (14) is to generate steady-state
operating points with desirable stability properties and minimal secondary control effort, so that the optimal zeq in (14) is
cognizant of dynamic stability constraints. The proposed LQROPF facilitates the development of a more general framework
where various dynamical system applications can be integrated
with static operational constraints, leading to the computation
of stability-aware operating points—more concrete pointers to
this end are provided in Section VIII.
Remark 1. In the previous sections, availability of the full
state-vector x(t) (i.e., angles, frequencies, and EMFs for
generators) is assumed for the optimal control law u(t). This
assumption entails that sensors such as phasor measurement
units (PMU) generate x(t) for all generator buses in the
network. However, a PMU is not needed on all generator
buses, as dynamic state estimation tools for power networks
can be leveraged. See e.g., [27], [28] for estimation through
Kalman filters and [29] for deterministic observers. For example, it has been previously demonstrated in [29] that 12 PMUs
are sufficient to estimate the 10th order dynamics of the 16
generators in the New England 68-bus power network.
Remark 2. Problem (14) can be straightforwardly extended
to include convex relaxations of the nonlinear power flow
equations (such as SDP relaxations [30]) in place of the
linearized power flow constraints (8a) and (8b).
V. A PPROXIMATE S OLVER FOR THE LQR-OPF
The LQR-OPF problem presented in (14) has two LMIs
and may be computationally expensive to solve for very
large networks via interior-point methods. In this section, an
algorithm to approximately solve (14) is developed.
To this end, consider the following nonconvex optimization:

s

min
s s

s

z ={x ,a ,u }
P

c(as ) +

Tlqr s
(x − x0 )⊤ P(xs − x0 )
2

subj. to (8a), (8b), as ∈ A
⊤

A P + PA − PBR

−1

s

(15a)
(15b)

s

(z )BP + Q(z ) = O, (15c)

where (15c) is the well-known continuous algebraic Riccati
equation (CARE). Problem (15) acts as a surrogate for (14)
because the optimal values of problems (14) and (15) are
equal. To see this, note first that the objective c(as ) and
constraints (8a), (8b), and as ∈ A are the same for both

Algorithm 1 Approximate solver for LQR-OPF (ALQR-OPF)
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:

Initialize obest ← +∞
Set P(0) to the solution of CARE (15c) for zs = z0
for k = 1 : k max do
Set z(k) = {x(k) , a(k) , u(k) } to the optimal solution
of (15a)–(15b) for P = P(k−1)
Set P(k) to the solution of CARE (15c) for zs = z(k)
T
(k)
if c(a(k) )+ lqr
− x0 )⊤ P(k) (x(k) − x0 )) < obest
2 (x
then
T
(k)
−x0 )⊤ P(k) (x(k) −x0 )
obest ← c(a(k) )+ lqr
2 (x
s
s s
s
(k)
z = {x , a , u } ← z = {x(k) , a(k) , u(k) }
P∗ = P(k)
end if
end for

problems. Then, using the theory in [31, Sec. V], it can be
shown that for any solution P of the CARE in (15c), the
quadratic form (xs − x0 )⊤ P(xs − x0 ) is equal to γ in (14a).
Although (15) is nonconvex, an efficient algorithm that alternates between variables zs and P can be used to approximately
solve (15). Specifically, with zs fixed, the CARE constraint
(15c) can be solved for to obtain P. Then, plugging P into (15)
and removing (15c) yields a quadratic program (QP), which
can be efficiently solved. The process can be repeated as long
as the objective in (15a) is improved, and is summarized in
Algorithm 1. The algorithm is referred to as ALQR-OPF.
The technology for solving the CARE and QPs with very
large number of variables has significantly matured and enables the efficient implementation of Algorithm 1. The numerical tests of Section VII indicate that Algorithm 1 does not
practically compromise optimality and is scalable to networks
with thousands of buses.
VI. G ENERATOR M ODEL , P OWER F LOWS , AND OPF
The developments in this paper and the proposed LQR
formulation are applicable to any generator dynamical model
with control inputs. However, as an example, a specific form
of mappings g and h is given in this section, which will be
used for the numerical tests of Section VII.
A. Generator model
The fourth order model of the synchronous generator internal dynamics for node i ∈ G can be written as
δ̇i = ωi − ωs
(16a)
1
ω̇i =
(16b)
[mi − Di (ωi − ωs ) − pgi ]
Mi 

xdi
1
xdi − x′di
− ′ ei +
ėi =
vi cos(δi − θi ) + fi (16c)
τdi
xdi
x′di


1
1
ṁi =
ri (t) −
(ωi − ωs ) − m(t)
(16d)
τci
Ri
where Mi is the rotor’s inertia constant (pu × sec2 ), Di is the
damping coefficient (pu × sec), τdi is the direct-axis opencircuit time constant (sec), xdi and xqi are respectively the
direct- and quadrature-axis synchronous reactances, and x′di
is the direct-axis transient reactance (pu). Equation (16d) is a

6

simplified model of a prime-mover generator with τci as the
charging time (sec) and a speed-governing mechanism with Ri
as the regulation constant ( Rad×Hz
). The mapping g defined
pu
in (1a) is given by concatenating (16) for i ∈ G.
The following algebraic equations relate the generator real
and reactive power output with generator voltage, internal
EMF, and internal angle and must hold at any time instant
for generator nodes i ∈ G:
ei vi
0 = −pgi + ′ sin(δi − θi )
xdi
x′di − xqi 2
+
v sin[2(δi − θi )]
(17a)
2xqi x′di i
ei vi
x′ + xqi
0 = −qgi + ′ cos(δi − θi ) − di ′ vi2
xdi
2xqi xdi
′
x − xqi
+ di ′ vi2 cos[2(δi − θi )].
(17b)
2xqi xdi

B. Power flow equations
Let Y = G + jB denote the network bus admittance matrix
based on the π-model of transmission lines. Notice that Y
may include transformers, tap-changing voltage regulators, and
phase shifters [32]. The power flow equations are
X
[Gi,j vi vj cos θij
−pli = −pgi + Gi,i vi2 +
j∈Ni

+ Bi,j vi vj sin θij ] , i ∈ G, (18a)
X
2
−qli = −qgi − Bi,i vi +
[Gi,j vi vj sin θij
j∈Ni

− Bi,j vi vj cos θij ] , i ∈ G, (18b)

−pli = Gi,i vi2 +
−qli = −

Bi,i vi2

X

j∈Ni

[Gi,j vi vj cos θij

+ Bi,j vi vj sin θij ] , i ∈ L, (18c)
X
+
[Gi,j vi vj sin θij
j∈Ni

− Bi,j vi vj cos θij ] , i ∈ L, (18d)

where θij := θi − θj ; pli := pli (t) and qli := qli (t) are
respectively the real and reactive power demands at node i
modeled as a time-varying constant-power load, i.e., pli and
qli are not functions of vi . The mapping h in (1b) is given by
concatenating (17) for i ∈ G and (18) for i ∈ N . By defining
plG = {pli }i∈G , qlG = {qli }i∈G , plL = {pli }i∈L , qlL =
{qli }i∈L , the vector d = {02G , −plG , −qlG , −plL , −qlL } is
obtained.
C. Optimal power flow
The standard OPF problem typically only considers the
algebraic variable as and is given as
min
c(as )
s
a

subj. to

(18) and as ∈ A.

(19)

If the cost function c(zs ) of the generalized OPF in (3) is
only a function of algebraic variables as , (3) can be solved by
solving the standard OPF problem (19) to obtain the optimal
as . The variables xs and us can then be found by plugging
in as in equations (17) and (16) while setting ẋ = 0.

VII. N UMERICAL S IMULATIONS
This section provides a numerical assessment of the advantages of the LQR-OPF in comparison to a method where OPF
and load-following control problems are treated separately.
Prior to analyzing the case studies, the general simulation
workflow of Fig. 2 is described.
The decoupled approach, one where OPF and stability are
solved separately, is considered on the left-hand side of Fig. 2.
Initially, the system is in steady-state and operates at z0 . For
a forecasted load demand, ds = d0 + ∆ds , the OPF (19) is
solved yielding optimal steady-state setpoints as , including
generator real and reactive powers (psgi , qgsi ) and generator
voltage magnitudes and phases (vis , θis ). For OPF, it holds
that aeq = as . The algebraic variables are then utilized
to solve (16) upon setting ẋ = 0 together with (17) to
obtain steady-state setpoints of generator states xeq and control
inputs ueq . The next steady-state equilibrium is then simply
represented as zeq = (xeq , aeq , ueq ). The DAEs (1a)–(1b),
upon being subjected to the new load ds , depart from the
initial equilibrium z0 . To steer the DAEs (1a)–(1b) to the next
desired equilibrium point zeq , LQR is performed. Dynamic
performance as well as costs of steady-state and load-following
control are evaluated. The standard OPF (19) is solved using
MATPOWER’s runopf.m.
The proposed methodology is considered on the righthand side of Fig. 2 where the OPF block is replaced by
LQR-OPF (14) followed by a load-flow. LQR-OPF obtains
optimal generator voltage setpoints {vis }i∈G , real power setpoints {psgi }i∈G\{islack } , and θisslack , while accounting for loadfollowing costs that drive the DAEs (1a)–(1b) to the next
desired equilibrium. These obtained setpoints are then input to a standard load-flow (performed by MATPOWER’s
runpf.m) by setting the following for i ∈ G: vieq = vis ,
eq
s
s
peq
gi = pgi for nonslack i, and θislack = θislack . This process
yields the remaining algebraic variables in aeq . Similar to the
previous approach, the DAEs (16) and (17) are solved after
setting ẋ = 0 yielding optimal steady-state setpoints of states
xeq and controls ueq . The next equilibrium is then given by
zeq = (xeq , aeq , ueq ). Once the DAE system departs from
its initial equilibrium point z0 due to ds , LQR is applied to
drive the system to the desired equilibrium zeq . LQR-OPF
is solved using the CVX optimization toolbox [33]. Finally,
a third approach where LQR-OPF in Fig. 2 is replaced by
ALQR-OPF is considered. In Algorithm 1, the CARE is solved
by MATLAB’s care.m and the QP using CVX. Comparisons
of the three approaches (OPF, LQR-OPF, and ALQR-OPF) are
provided in Table I.
The workflow of Fig. 2 is applied to various networks as
described in the next section. It is worth emphasizing that even
though the derived control law required linearized system dynamics, all simulations are performed on the actual nonlinear
dynamics and power flow equations, per the DAEs in (1a)–
(1b). The dynamical system was modeled using MATLAB’s
ode suite. All computations use a personal computer with
32.0 GB RAM and 3.60 GHz CPU processor. MATLAB
scripts that simulate the ensuing case studies are provided
online [34].

7

OPF and LQR decoupled

Start
from z0

OPF (19)
aeq

Proposed approach

LQR-OPF (14)
followed by load-flow
aeq

Solve eq. (16), (17)
upon setting ẋ = 0

Solve eq. (16), (17)
upon setting ẋ = 0

Setpoints
zeq
LQR

Setpoints
zeq
LQR

Feedback gain
Simulation
of nonlinear
DAEs (16)–(18)
Trajectories

Feedback gain
Simulation
of nonlinear
DAEs (16)–(18)
Trajectories

Evaluate performance

Evaluate performance

Total Cost

Total Cost

Fig. 2. Diagram showing the simulation steps.

A. Description of test networks
The workflow of Fig. 2 is applied to 9-, 14-, 39-, and 57bus test networks as well as the 200-bus Illinois network.
Steady-state data are obtained from the corresponding case
files available in MATPOWER [35]. These data include bus
admittance matrix Y, steady-state real and reactive power
demands, that is, p0l := {p0lG , p0lL } and q0l := {q0lG , q0lL }, cost
P
of real power generation c(as ) = i c2 (psgi )2 + c1 psgi + c0 ,
as well as the limits of power generation and voltages in
constraint (3d). Machine constants of (16a)–(16c) and (17)
are taken from the PST toolbox, from case files d3m9bm.m,
d2asbegp.m, and datane.m, respectively for the 9-, 14-,
and 39-bus networks. For the 57- and the 200-bus networks,
as well as the governor model of (16d) for all networks, based
on ranges of values provided in PST [36], typical parameter
values of Mi = 0.2 pu × sec2 , Di = 0 pu × sec, τdi = 5 sec,
xdi = 0.7 pu, xqi = 0.5 pu, x′di = 0.07 pu, τci = 0.2 sec,
and Ri = 0.02 Hz
pu have been selected.
B. Regulation cost matrices Q and R
In accordance with (13), Q−1 and R−1 are selected to be
diagonal with affine entries as follows


psgi
−1
−1
(20a)
Qωi ,δi ,mi = Rri = 1 − α max ,
pgi


qgsi
−1
,
(20b)
Q−1
ei = Rfi = 1 − α max
qgi
where Qωi ,δi ,mi refers to the diagonal entries of Q corresponding to ωi , δi , and mi . Matrices Qei , Rri , and Rfi are
similarly defined. Parameter α is in the interval [0, 1) that
determines the amount of coupling between steady-state quantities and control costs through matrices Q and R. Quantities
are respectively the maximum real and reactive
and qgmax
pmax
gi
i
power limits of generator i ∈ G.

The rationale behind choosing (20a) is that angle and
frequency instability are usually remedied by generating real
power. In this case, increase in steady-state real power generation psgi leads to higher cost of frequency regulation. Similarly,
the rationale behind choosing (20b) is that voltage stability
is typically correlated with reactive power injection. This
choice means that an increase in steady-state reactive power
generation qgsi incurs a higher cost of voltage regulation.
C. Dynamical simulation with LQR
A step load increase of 10% in real power with power factor
0.9 is applied at t = 0. This implies that ∆psli = 0.1p0li
and ∆qlsi = 0.0484ql0i , totaling to a load increase of 31.50 +
j5.56 MVA for the 9-bus system, 25.90 + j3.56 MVA for the
14-bus system, 625.42 + j67.10 MVA for the 39-bus system,
125.08 + j16.28 MVA for the 57-bus, 222.87 + j30.74 MVA
for the 200-bus Illinois system. This load increase drives
the nonlinear dynamics (1) out of the initial equilibrium. By
applying LQR control according to the workflow in Fig. 2, the
dynamics in (1) are steered to arrive at the desired equilibria
obtained from the OPF and LQR-OPF.
With selections of α = 0.6 and Tlqr = 1000, Table I lists the
breakdown of steady-state, control, and total costs, as well as
maximum frequency and voltage deviations from the optimal
equilibrium of the LQR-OPF, ALQR-OPF with two iterations,
and OPF. Column 3 gives the optimal objective of the LQROPF in (14) or the best objective found by the ALQR-OPF
in Algorithm 1. Observe that the optimal objectives of LQROPF and ALQR-OPF are almost identical, which implies
that the approximation in Algorithm 1 does not compromise
optimality.
Steady-state costs in column 4 correspond to c(aeq ) found
by each approach. The LQR step in the diagram of Fig. 2
is solved using MATLAB’s care.m by inputting Q(zeq ) and
R(zeq ). This process yields P and the corresponding feedback
gain K = −R−1 B⊤ P. The estimates of the control costs
T
eq
0 ⊤
eq
0 ⊤
are then calculated as lqr
2 (x − x ) P(x − x ) and are
given in column 5 of Table I. The total estimated costs are
then the summation of steady-state and estimated control costs
and are given in column 6. The computation times of LQROPF, ALQR-OPF, and OPF are listed in column 7. Notice
that for the large 200-bus Illinois network, the LQR-OPF
takes approximately 12 hours, while the ALQR-OPF solves
the problem in less than three seconds and without significant
loss in optimality.
Control
in column 8 of Table I are computed
R tfcosts ′reported
T
⊤
⊤
′
(∆x
Q∆x
+ ∆u′ R∆u′ )dt, that is, through
as lqr
2
0
numerical integration of the trajectories resulting from the
simulation of the nonlinear DAEs. The total cost, given in
column 9, is simply the summation of control and steadystate cost. Between the coupled and decoupled approaches,
OPF exhibits lower steady-state cost but higher control cost.
In terms of total cost, the LQR-OPF and ALQR-OPF show
improved performance. The maximum frequency deviation is
also much lower for LQR-OPF and ALQR-OPF than the OPF.
In Figs. 3 and 4 the dynamical performance of the 39-bus
system under LQR-OPF, ALQR-OPF, and OPF in conjunction

8

TABLE I
C OSTS C OMPARISON BETWEEN LQR-OPF, ALQR-OPF, AND OPF METHODS AT α = 0.6 UNDER LQR

Network

9-bus

14-bus

39-bus

57-bus

200-bus

Method

Obj.

Steady-state
cost
($)

Control
est. cost
($)

Total
est. cost
($)

Comp.
time
(seconds)

Control
cost
($)

LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF

6168.78
6168.78
—
9209.27
9209.28
—
55560.94
55561.20
—
50169.04
50177.06
—
52731.74
54575.99
—

6144.11
6144.25
6113.60
9178.82
9178.57
9127.35
52872.28
52871.80
51386.02
48322.27
48368.33
47199.75
48347.04
50349.65
48271.72

26.30
26.16
329.02
36.20
36.48
248.29
2471.68
2472.94
12504.70
2468.43
2410.75
5829.34
4090.34
2424.85
8724.08

6170.41
6170.41
6442.62
9215.01
9215.04
9375.64
55343.95
55344.74
63890.72
50790.70
50779.08
53029.09
52437.38
52774.50
56995.80

2.57
0.85
0.64
2.32
0.74
0.26
18.95
1.39
0.23
6.01
0.87
0.22
44070.98
2.92
0.72

28.02
27.87
316.01
36.95
37.26
248.13
3120.73
3128.34
13468.42
2306.30
2260.24
5889.98
4089.44
2526.84
9258.04

($)

Max.
freq. dev.
(Hz)

Max.
volt. dev.
(pu)

6172.13
6172.11
6429.61
9215.77
9215.83
9375.48
55993.01
56000.14
64854.44
50628.57
50628.57
53089.73
52436.48
52876.49
57529.76

0.0145
0.0145
0.0202
0.0052
0.0053
0.0120
0.0804
0.0804
0.1333
0.0602
0.0593
0.0944
0.0283
0.0199
0.0468

0.0262
0.0261
0.1240
0.0297
0.0297
0.0297
0.0820
0.0816
0.1157
0.0599
0.0600
0.0637
0.0608
0.0589
0.0745

Total

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 3. Generator frequencies and governor signals using LQR-OPF (left), ALQR-OPF (center), and OPF (right).

with load-following LQR is depicted. Specifically, generator
frequencies and governor reference signals are portrayed in
Fig. 3 where it is observed that quantities, especially frequencies, resulting from OPF undergo higher fluctuations than
those resulting from LQR-OPF and ALQR-OPF. In Fig. 4 the
generator internal EMF, exciter field voltage, and nodal voltages are depicted. Notice that overall, voltages obtained from
LQR-OPF and ALQR-OPF exhibit smaller deviations from the
designated equilibria in comparison with those obtained from
OPF. Deviations of generator angles δ also exhibit similar
behavior. Plots for the remaining quantities (δ and m) and the
corresponding plots for the remaining test networks (including
the 200-bus Illinois network) are available online [34].

D. Dynamical simulation with AGC
The LQR-OPF furnishes a steady-state operating point with
desirable stability properties. After the steady-state operating
point has been computed, one does not have to necessarily use
LQR as a controller, but one could rather implement another
dynamic control law such as AGC or PI-control [37], [38].
The purpose of this section is to examine the control cost to
drive the system to the equilibrium computed by LQR-OPF
or OPF when the control law is AGC.
To this end, the setup of Section VII-C is followed, but
AGC is used to adjust the governor reference signal during the
dynamical simulations instead of LQR. Dynamical equations
describing the AGC for a multi-area power network are
adopted from [37] and [38]. The selection of participation
factors to steer the DAEs to the desired equilibrium follows

9

(a)

(b)

(c)

(d)

(e)

(f)

(g)
Fig. 4.

(h)

(i)

Generator EMF, exciter voltage, and network nodal voltages using LQR-OPF (left), ALQR-OPF (center), and OPF (right).

the suggestions in [39, p. 87] for computer implementations.
The specifics are detailed in the Appendix.
TABLE II
C OSTS C OMPARISON BETWEEN LQR-OPF, ALQR-OPF, AND OPF
METHODS AT α = 0.6 UNDER AGC

Network

9-bus

14-bus

39-bus

57-bus

200-bus

Method

Control
cost
($)

Total
cost
($)

Max.
freq. dev.
(Hz)

Max.
volt. dev.
(pu)

LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF
LQR-OPF
ALQR-OPF
OPF

35.92
35.77
583.90
94.19
94.30
590.84
4819.03
4843.82
17781.72
4507.44
4401.21
12652.76
9214.91
5667.92
14959.59

6180.04
6180.02
6697.50
9273.01
9272.86
9718.19
57691.30
57715.62
69167.74
52829.71
52769.55
59852.51
57561.95
56017.57
63231.30

0.0144
0.0144
0.0161
0.0040
0.0041
0.0096
0.0765
0.0765
0.1321
0.0489
0.0481
0.0760
0.0249
0.0168
0.0358

0.0262
0.0261
0.1240
0.0297
0.0297
0.0297
0.0782
0.0772
0.1062
0.0599
0.0600
0.0637
0.0608
0.0589
0.0741

Table II reports the control and total costs of LQR-OPF
and ALQR-OPF compared with OPF, when AGC is used

to adjust the governor signals ri for i ∈ G. Notice that
the steady-state results are those previously given in Table I
and only quantities pertaining to dynamical simulations have
changed. It is observed that using AGC, setpoints provided
by LQR-OPF and ALQR-OPF result in smaller control costs
compared to OPF.1 The corresponding dynamical performance
for ALQR-OPF and OPF in conjunction with AGC is provided
in Fig. 5–8. The setpoints provided by ALQR-OPF result in
smaller frequency deviation and smaller ACE signal compared
to the setpoints provided by OPF. The performance using the
setpoints of LQR-OPF is similar to that of ALQR-OPF and
has been omitted for brevity. Corresponding plots of other
system quantities for the 39-bus network and the remaining
test networks are also made available online [34].
E. Effect of coupling
Here, the effect of parameter α which couples steadystate variables to load-following control costs through (20) is
studied. When the value of α increases to approach 1, entries
of matrices Q and R increase as the values of psgi and qgsi
1 Table II indicates that total costs of ALQR-OPF are sometimes smaller
than those of LQR-OPF. This happens because the optimality of LQR-OPF
over ALQR-OPF is only guaranteed when the associated controller is LQR.

10

Fig. 5. Generator frequencies using ALQR-OPF in conjunction with AGC.

Fig. 7. Area control signal using ALQR-OPF in conjunction with AGC.

Fig. 6. Generator frequencies using OPF in conjunction with AGC.

Fig. 8. Area control signal using OPF in conjunction with AGC.

approach their respective maximum. It is depicted in Fig. 9 that
as the coupling coefficient α increases, control costs increase
in both OPF and LQR-OPF. However, control costs of LQROPF are significantly lower than the costs incurred by the
scheme where OPF and LQR are solved independently.

the LQR. Numerical tests verify that compared to a scheme
where OPF and load-following control problems are solved
separately, LQR-OPF features significantly improved dynamic
performance and reduced overall system costs.
The proposed framework is general and allows for seamless
incorporation of different power system applications, such as
wind turbine and storage device dynamics—both operating at
different time-scales. Recent work, for instance, demonstrates
the impact of wind power injection on power system oscillations [40]. Future work includes integrating more modern
applications into the proposed framework while investigating
the regulation and cost benefits of this integration.

F. ALQR-OPF on larger networks
This section examines the performance of ALQR-OPF on
larger networks. Table III lists the steady-state costs, estimated
T
eq
0 ⊤
eq
0
control costs (as computed by lqr
2 (x − x ) P(x − x )),
total estimated costs, and computation times. The main observation is that ALQR-OPF yields significantly smaller total
costs than OPF. The computation time of ALQR-OPF is larger
than that of OPF, but it is worth noting that ALQR-OPF has
been solved by a general-purpose solver, while MATPOWER’s
solver is specifically tailored to the OPF problem.
VIII. S UMMARY AND F UTURE W ORK
An OPF framework is presented that in addition to solving
for optimal steady-state setpoints provides an optimal feedback
law to perform load-following control. The costs of loadfollowing control is captured by a classical LQR control that
accounts for deviations of system states and controls from their
optimal steady-state setpoints. A joint formulation of OPF and
load-following control, termed LQR-OPF, is obtained by combining a linearized OPF with an equivalent SDP formulation of

AGC

A PPENDIX
IMPLEMENTATION

Denote by A the set of areas of a power network and by
Aa the set of neighboring areas to area a ∈ A. Further, denote
respectively by paa′ and peq
aa′ , the time-varying and equilibrium
aggregate real power flows on the tie-lines from area a ∈
A into area a′ ∈ A . Notice that paa′ is a function of the
voltage magnitudes and angles at the terminals of all the tielines connecting a to a′ . The area control error (ACE) for area
a ∈ A is then computed as
!
X
1 X
eq
s
ωi − ω
(21)
ACEa =
(paa′ − paa′ ) + ba
|Ga |
′
a ∈Aa

i∈Ga

11

Fig. 9. Effect of coupling between steady-state variables and control cost
matrices Q and R. As coupling coefficient α increases, LQR-OPF results in
more savings than the decoupled OPF and LQR approach.
TABLE III
ALQR-OPF ON LARGER NETWORKS
Network
1354-bus
2383-bus
2869-bus

Method
ALQR-OPF
OPF
ALQR-OPF
OPF
ALQR-OPF
OPF

Steady-state
cost
($)
82612
81688
2235711
2217287
149453
147865

Control
est. cost
($)
2311022
8279013
21174
123999
3255044
14830044

Total
est. cost
($)
2393635
8360702
2256885
2341287
3404498
14977910

Comp.
time
(seconds)
72.65
6.06
325.72
9.86
683.74
52.70

where Ga denotes the set of generators in area a, and ba =
P
1
i∈Ga ( Ri + Di ) is the area bias factor.
AGC uses the ACE signal to provide a command to the governor reference signal. The corresponding dynamical equations
for area a ∈ A are
!
X
eq
(22a)
ẏa = Ka −ya − ACEa +
pgi
i∈Ga

ri = Ki ya ,

i ∈ Ga ,

(22b)

where Ka is an integrator gain, Ki is the participation factors
of generator i, and ri is fed back
The participation
Pinto (16d).
eq
so
that the power
p
/
factors are set to Ki = peq
gi
i∈Ga gi
network will be steered to the desired equilibrium. Notice that
the sum of participation factors per area equals unity.
R EFERENCES
[1] A. Gomez-Exposito, A. J. Conejo, and C. Canizares, Electric Energy
Systems: Analysis and Operation. CRC Press, 2009.
[2] F. Galiana, F. Bouffard, J. Arroyo, and J. Restrepo, “Scheduling and
Pricing of Coupled Energy and Primary, Secondary, and Tertiary Reserves,” Proc. IEEE, vol. 93, no. 11, pp. 1970–1983, Nov. 2005.
[3] M. Andreasson, D. V. Dimarogonas, H. Sandberg, and K. H. Johansson,
“Distributed control of networked dynamical systems: Static feedback,
integral action and consensus,” IEEE Trans. Autom. Control, vol. 59,
no. 7, pp. 1750–1764, July 2014.
[4] Q. Shafiee, J. M. Guerrero, and J. C. Vasquez, “Distributed secondary
control for islanded microgrids—A novel approach,” IEEE Trans. Power
Electron., vol. 29, no. 2, pp. 1018–1031, Feb. 2014.
[5] F. Dörfler, J. W. Simpson-Porco, and F. Bullo, “Breaking the hierarchy:
Distributed control and economic optimality in microgrids,” IEEE Trans.
Control Netw. Syst., vol. 3, no. 3, pp. 241–253, Sept. 2016.
[6] F. Dörfler and S. Grammatico, “Gather-and-broadcast frequency control
in power systems,” Automatica, vol. 79, pp. 296–305, 2017.

[7] C. Zhao, U. Topcu, N. Li, and S. Low, “Design and stability of loadside primary frequency control in power systems,” IEEE Trans. Autom.
Control, vol. 59, no. 5, pp. 1177–1189, May 2014.
[8] N. Li, C. Zhao, and L. Chen, “Connecting automatic generation control
and economic dispatch from an optimization view,” IEEE Trans. Control
Netw. Syst., vol. 3, no. 3, pp. 254–264, Sept. 2016.
[9] X. Zhang and A. Papachristodoulou, “A real-time control framework for
smart power networks: Design methodology and stability,” Automatica,
vol. 58, pp. 43–50, 2015.
[10] E. Mallada, C. Zhao, and S. Low, “Optimal load-side control for
frequency regulation in smart grids,” IEEE Trans. Autom. Control, to
be published.
[11] D. Gan, R. J. Thomas, and R. D. Zimmerman, “Stability-constrained
optimal power flow,” IEEE Trans. Power Syst., vol. 15, no. 2, pp. 535–
540, May 2000.
[12] L. Tang and W. Sun, “An automated transient stability constrained
optimal power flow based on trajectory sensitivity analysis,” IEEE Trans.
Power Syst., vol. 32, no. 1, pp. 590–599, Jan. 2016.
[13] Y. Xu, J. Ma, Z. Y. Dong, and D. J. Hill, “Robust transient stabilityconstrained optimal power flow with uncertain dynamic loads,” IEEE
Trans. Smart Grid, vol. 8, no. 4, pp. 1911–1921, 2017.
[14] Y. Xu, M. Yin, Z. Dong, R. Zhang, D. J. Hill, and Y. Zhang, “Robust Dispatch of High Wind Power-Penetrated Power Systems against
Transient Instability,” IEEE Trans. Power Syst., to be published.
[15] X. Xu, H. Zhang, C. Li, Y. Liu, W. Li, and V. Terzija, “Optimization
of the Event-Driven Emergency Load-Shedding Considering Transient
Security and Stability Constraints,” IEEE Trans. Power Syst., vol. 32,
no. 4, pp. 2581–2592, July 2017.
[16] R. Zarate-Minano, F. Milano, and A. J. Conejo, “An OPF methodology
to ensure small-signal stability,” IEEE Trans. Power Syst., vol. 26, no. 3,
pp. 1050–1061, Aug. 2011.
[17] P. Li, J. Qi, J. Wang, H. Wei, X. Bai, and F. Qiu, “An SQP method
combined with gradient sampling for small-signal stability constrained
opf,” IEEE Trans. Power Syst., vol. 32, no. 3, pp. 2372–2381, May 2017.
[18] E. Mallada and A. Tang, “Dynamics-aware optimal power flow,” in Proc.
52nd IEEE Conf. Decision and Control, Florence, Italy, Dec. 2013, pp.
1646–1652.
[19] P. Li, H. Wei, B. Li, and Y. Yang, “Eigenvalue-optimisation-based optimal power flow with small-signal stability constraints,” IET Generation,
Transmission & Distribution, vol. 7, no. 5, pp. 440–450, May 2013.
[20] C. E. Fosha and O. I. Elgerd, “The megawatt-frequency control problem:
A new approach via optimal control theory,” IEEE Trans. Power App.
Syst., vol. 89, no. 4, pp. 563–577, Apr. 1970.
[21] A. K. Singh and B. C. Pal, “Decentralized control of oscillatory
dynamics in power systems using an extended LQR,” IEEE Trans. Power
Syst., vol. 31, no. 3, pp. 1715–1728, May 2016.
[22] J. A. Taylor, A. Nayyar, D. S. Callaway, and K. Poolla, “Consolidated
dynamic pricing of power system regulation,” IEEE Trans. Power Syst.,
vol. 28, no. 4, pp. 4692–4700, Nov. 2013.
[23] I. Hiskens and M. Pai, “Power system applications of trajectory sensitivities,” in Proc. Power Engineering Society Winter Meeting, vol. 2,
New York, NY, USA, Jan. 2002, pp. 1200–1205.
[24] L. Tang and J. McCalley, “Trajectory sensitivities: Applications in power
systems and estimation accuracy refinement,” in Proc. Power & Energy
Society General Meeting, Vancouver, BC, Canada, July 2013, pp. 1–5.
[25] D. K. Molzahn, “Incorporating squirrel-cage induction machine models
in convex relaxations of OPF problems,” IEEE Trans. Power Syst., to
be published.
[26] S. Boyd, L. E. Ghaoui, E. Feron, and V. Balakrishnan, Linear Matrix
Inequalities in Systems and Control Theory. SIAM, 1994, vol. 15.
[27] Y. Cui and R. Kavasseri, “A particle filter for dynamic state estimation
in multi-machine systems with detailed models,” IEEE Trans. Power
Syst., vol. 30, no. 6, pp. 3377–3385, Nov. 2015.
[28] G. Valverde and V. Terzija, “Unscented Kalman filter for power system
dynamic state estimation,” IET Generation, Transmission & Distribution,
vol. 5, no. 1, pp. 29–37, Jan. 2011.
[29] A. F. Taha, J. Qi, J. Wang, and J. H. Panchal, “Risk mitigation for
dynamic state estimation against cyber attacks and unknown inputs,”
IEEE Trans. Smart Grid, to be published.
[30] R. Madani, S. Sojoudi, and J. Lavaei, “Convex Relaxation for Optimal
Power Flow Problem: Mesh Networks,” IEEE Trans. Power Syst.,
vol. 30, no. 1, pp. 199–211, Jan. 2015.
[31] V. Balakrishnan and L. Vandenberghe, “Semidefinite programming duality and linear time-invariant systems,” IEEE Trans. Autom. Control,
vol. 48, no. 1, pp. 30–41, Jan. 2003.

12

[32] G. B. Giannakis, V. Kekatos, N. Gatsis, S. J. Kim, H. Zhu, and B. F.
Wollenberg, “Monitoring and optimization for power grids: A signal
processing perspective,” IEEE Signal Process. Mag., vol. 30, no. 5, pp.
107–128, Sept. 2013.
[33] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex
programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
[34] M. Bazrafshan, “LQR-OPF,” 2017. [Online]. Available: https://github.
com/hafezbazrafshan/LQR-OPF
[35] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power
systems research and education,” IEEE Trans. Power Syst., vol. 26, no. 1,
pp. 12–19, Feb. 2011.
[36] J. H. Chow and K. W. Cheung, “A toolbox for power system dynamics
and control engineering education and research,” IEEE Trans. Power
Syst., vol. 7, no. 4, pp. 1559–1564, Nov. 1992.
[37] J. Zhang and A. Dominguez-Garcia, “On the Impact of Measurement
Errors on Power System Automatic Generation Control,” IEEE Trans.
Smart Grid, to be published.
[38] J. D. Glover, M. S. Sarma, and T. J. Overbye, Power System Analysis
and Design, 5th ed. Cengage Learning, 2012.
[39] A. J. Wood and B. F. Wollenberg, Power Generation, Operation, and
Control, 3rd ed. John Wiley & Sons, 2012.
[40] S. Chandra, D. F. Gayme, and A. Chakrabortty, “Time-scale modeling
of wind-integrated power systems,” IEEE Trans. Power Syst., vol. 31,
no. 6, pp. 4712–4721, Nov. 2016.

