Description
Modeling, Analysis, and Control of Renewable Energy Systems
Problem 1. Consider an AC electric power network with synchronous generators indexed in set G
and transmission lines indexed in set E. For each generator g ∈ G, let δg, ωg, P
m
g
, P
e
g
, and P
r
g denote
the electrical angular position, angular speed, turbine mechanical power, electrical power output, and
reference power input, respectively. Dynamics of generator g can be described by the swing equations
augmented with a simplified turbine-governor model, as follows:
˙δg = ωg − ωs, (1)
Mgω˙ g = P
m
g − P
e
g
, (2)
τgP˙ m
g = P
r
g − P
m
g −
1
Rgωs
(ωg − ωs), (3)
where Mg, τg, and Rg denote its inertia constant, governor time constant, and droop constant, respectively; and ωs = 2π60 [rad/s] is the synchronous frequency of the system. The electrical power
output, P
e
g
, is a function of bus-voltage magnitudes and phase angles as solved from the network power
balance. The generator g reference power input is given by
P
r
g = P
?
g + αg(ξ −
X
j∈G
P
?
j
), (4)
where P
?
g
is the setpoint from economic dispatch and αg is the AGC participation factor with P
g∈G αg =
1. Moreover, ξ is the AGC state whose evolution is dictated by
˙ξ = −ξ − ACE +X
g∈G
P
e
g
, (5)
where ACE denotes the area control error that accounts for deviations in frequency from the synchronous value. For a single-area system with no tie-line flows, the area control error is formulated as
ACE = b(ω − ωs) where b > 0 is the area bias factor and ω the prevailing frequency of the system.
(Without loss of generality, we assume ω is computed with real-time measurements as the average
electrical frequency of all generators.) Let Cg(Pg) denote the cost function for generator g, and collect
Pg, g ∈ G, in PG. We will assume the economic dispatch problem takes the following form:
min.
Pg, g∈G
X
g∈G
Cg(Pg) (6a)
s.t.
X
g∈G
Pg = Pload + Ploss(PG), (6b)
where Pload is the look-ahead net load and Ploss(PG) is the system loss modeled as a function of PG.
Part (i). Let us denote values taken by variables at a new steady-state operating point by X (corresponding values at the initial operating point are denoted by X). Suppose the system net load Pload
changes by amount ∆Pload, so that the new net load is Pload = Pload + ∆Pload. In steady state after
the load change, show that the dynamical system (1)–(5) converges to the following operating point:
ωg = ωs, ∀ g ∈ G, (7)
P
m
g = P
?
g + αg(∆Pload + ∆Ploss), ∀ g ∈ G, (8)
ξ = Pload + ∆Pload + ∆Ploss, (9)
1
EE 8744: Modeling, Analysis, and Control of Renewable Energy Systems Spring 2021
where ∆Ploss is the change in system loss due to the change in operating point.
Part (ii). Show that the generator g steady-state electrical power output with the new net load is
P
e
g = P
?
g + αg(∆Pload + ∆Ploss), ∀ g ∈ G. (10)
Part (iii). Show that the Karush–Kuhn–Tucker (KKT) conditions for the economic dispatch problem (6) are given by:
C
0
g
(P
?
g
) −
λ
?
Λ?
g
= 0, ∀g ∈ G, (11)
where Λ?
g
is the loss penalty factor for generator g evaluated at P
?
G
; defined as
Λ
?
g
:=
1 −
∂Ploss(P
?
G
)
∂Pg
−1
. (12)
Part (iv). In the lecture modules, we have seen that a good choice for the AGC participation factors,
αg, g ∈ G from an economic point of view is given by
αg =
(C
00
g
(P
?
g
))−1
P
j∈G
(C00
j
(P
?
j
))−1
, ∀ g ∈ G. (13)
We will now explore the impact of this choice choice of participation factors on system dynamic performance. Suppose the system initially operates in steady state with net load, Pload, and generator
electrical power outputs corresponding to a KKT point, P
?
g
, g ∈ G, of the economic dispatch problem (6). Note that a KKT point is a solution that satisfies the KKT conditions but may not be a
global/local minimum of the optimization problem. (Every minimum/maximum is a KKT point but
not vice versa.) Consider AGC action triggered by a change in net load to Pload = Pload + ∆Pload.
Denote by P
?
g
, g ∈ G, a KKT point of the economic dispatch problem solved with the new net load,
Pload. Further, denote the loss penalty factors corresponding to P
?
g and P
?
g
, by Λ?
g and Λ
?
g
, respectively.
Consider the following assumptions to hold true
[A1] The cost function Cg(Pg) is twice differentiable.
[A2] The generator-cost and system-loss functions are such that for the KKT points P
?
g
, P
?
g
:
Λ
?
gC
0
g
(P
?
g
) − Λ
?
gC
0
g
(P
?
g
) = (P
?
g − P
?
g
)C
00
g
(P
?
g
). (14)
Under assumptions [A1]–[A2], prove that if the AGC participation factors are picked per (13), then
the steady-state generator electrical power outputs following the net-load change correspond to KKT
point, P
?
g
, i.e.,
P
e
g = P
?
g
, ∀ g ∈ G. (15)
Comment on the implication/significance of the above result.
2
EE 8744: Modeling, Analysis, and Control of Renewable Energy Systems Spring 2021
Part (v). Can you provide specific cases (types of cost functions, types of networks) under which
assumptions [A1]–[A2] are automatically satisfied?
Problem 2. Notation. The matrix transpose will be denoted by (·)
T, complex conjugate by (·)
∗
,
real and imaginary parts of a complex number by Re{·} and Im{·}, respectively, magnitude of a
complex scalar by | · |, and j := √
−1. A diagonal matrix formed with entries of the vector x is
denoted by diag(x); diag(x/y) forms a diagonal matrix with the `th entry given by x`/y`
, where x`
and y` are the `th entries of vectors x and y, respectively; and diag(1/x) forms a diagonal matrix
with the `th entry given by x
−1
`
. For a vector x = [x1, . . . , xN ]
T, cos(x) := [cos(x1), . . . , cos(xN )]T
and sin(x) := [sin(x1), . . . ,sin(xN )]T. We will routinely decompose the complex-valued vector x ∈ CN
(complex-valued matrix X ∈ CN×N ) into its real and imaginary parts as follows: x = xre + jxim
(X = Xre + jXim, respectively). The spaces of N × 1 real-valued and complex-valued vectors are
denoted by RN and CN , respectively; TN denotes the N-dimensional torus. With 0N and 1N , we
denote N-dimensional column vectors with all entries equal to 0 and 1, respectively. The N × N
identity matrix is denoted by IN×N and the N × N matrix with all 0 entries is denoted by 0N×N .
Power System Model. Consider a power system with N + 1 buses collected in the set N . We model
loads as the parallel interconnection of constant impedance and a constant power component. Without
loss of generality, the slack bus is fixed to be the N + 1 bus, and its voltage is denoted by V◦e
jθ◦
. Let
V = [V1, . . . , VN ]
T ∈ CN , where V` = |V |`
6 θ` ∈ C represents the voltage phasor at bus `. In
subsequent developments, we will find it useful to define the vectors |V | = [|V1|, . . . , |VN |]
T ∈ RN
>0
and θ = [θ1, . . . , θN ]
T ∈ TN . Given our focus on rectangular coordinates, we will also routinely
express V = Vre + jVim, where Vre, Vim ∈ RN denote the real and imaginary components of V . Let
I = [I1, . . . , IN ]
T, where I` ∈ C denotes the current injected into bus `. Kirchhoff’s current law for
the buses in the power system can be compactly represented in matrix-vector form as follows:
I
IN+1
=
”
Y Y
Y
T
y
#
V
V◦e
jθ◦
, (16)
where V◦e
jθ◦
is the slack-bus voltage, IN+1 denotes the current injected into the slack bus, and the
entries of the admittance matrix have the following dimensions: Y ∈ CN×N , Y ∈ CN , and y ∈ C\{0}.
We will decompose the matrix Y and its inverse Y
−1 as follows:
Y = G + jB, Y −1 = R + jX. (17)
Denote the vector of complex-power bus injections by S = [S1, . . . , SN ]
T, where S` = P` + jQ`
. By
convention, P` and Q` are positive for generators and negative for loads.
Part (i). Using (16), show that the complex-power bus injections can be compactly written as
S = diag (V )
Y
∗V
∗ + Y
∗
V◦e
−jθ◦
. (18)
Part (ii). Express V = V
nom + ∆V , where V
nom is some a priori determined nominal voltage vector,
and entries of ∆V capture perturbations around V
nom. With the choice
V
nom = −Y
−1Y V◦e
jθ◦
, (19)
show that ∆V can be solved from the following linear equation (once second-order terms are neglected)
diag ((V
nom)
∗
) Y ∆V = S
∗
. (20)
3
EE 8744: Modeling, Analysis, and Control of Renewable Energy Systems Spring 2021
Part (iii). Define the following matrices
K = diag(V
nom)Y
∗
, J =
Re(K) Im(K)
Im(K) −Re(K)
. (21)
Show that the real and imaginary parts of ∆V , ∆Vre and ∆Vim are given by the solution of the
following linear equations
∆Vre
∆Vim
= J
−1
P
Q
. (22)
Part (iv). Consider a complex scalar η such that |η| << 1. Show that
|1 + η| ≈ 1 + Re(η), (23)
6 (1 + η) ≈ Im(η). (24)
Part (v). Leveraging the approximations in (23)–(24) and the expressions in (21)–(22), show that
|V | ≈ |V
nom| + [diag(cos(θ
nom)) diag(sin(θ
nom))]J
−1
P
Q
. (25)
Part (vi). Leveraging the approximations in (23)–(24) and the expressions in (21)–(22), show that
θ ≈ θ
nom + diag(|V
nom|)
−1
[−diag(sin(θ
nom)) diag(cos(θ
nom))]J
−1
P
Q
. (26)
Part (vii). Consider a two-bus power system composed of a single P Q load connected to bus 1 and
bus 2 representing the slack bus. At buses 1 and 2, there are shunt admittances of 0.02j; and the
impedance of the line connecting bus 1 and 2 is 0.01 + 0.05j. The real power drawn by the load at
bus 1 is 2 p.u., and the reactive power drawn is 0.75 p.u. The voltage of the slack bus, Vo
6 θo = 16 30o
.
Compute the voltage magnitude and phase angle at bus 1 using the Newton-Raphson approach run
out for 5 iterations. (Choose the bus-2 voltage as the initial guess for the voltage magnitude and phase
angle.) Verify the accuracy of the approximations in (25)-(26) with the result you obtain from the
Newton-Raphson method. Include your code in the solution.
4
