Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

1

Hybrid Controller for Wind Turbine Generators to
Ensure Adequate Frequency Response in Power
Networks

arXiv:1802.03429v1 [cs.SY] 9 Feb 2018

Yichen Zhang, Student Member, IEEE, Kevin Tomsovic, Fellow, IEEE, Seddik M. Djouadi, Member, IEEE
and Hector Pulgar-Painemal, Member, IEEE

Abstract—Converter-interfaced power sources (CIPS) are hybrid control systems as they may switch between multiple
operating modes. Due to increasing penetration, the hybrid
behavior of CIPS, such as, wind turbine generators (WTG), may
have significant impact on power system dynamics. In this paper,
the frequency dynamics under inertia emulation and primary
support from WTG is studied. A mode switching for WTG to
ensure adequate frequency response is proposed. The switching
instants are determined by our proposed concept of a region
of safety (ROS), which is the initial set of safe trajectories.
The barrier certificate methodology is employed to derive a
new algorithm to obtain and enlarge the ROS for the given
desired safe limits and the worst-case disturbance scenarios. Then
critical switching instants and a safe recovery procedure are
found. In addition, the emulated inertia and load-damping effect
is derived in the time frame of inertia and primary frequency
response, respectively. The theoretical results under critical cases
are consistent with simulations and can be used as guidance for
practical control design.
Index Terms—System frequency response, deadband, hybrid
system, inertia emulation, primary frequency support, safety
verification, sum of squares decomposition, semidefinite programming, wind turbine generator.

I. I NTRODUCTION
Ybrid behaviors in complex power networks have not
been carefully studied. However, with the increasing
connection of converter-interfaced power sources (CIPS), such
as, wind turbine generators (WTGs), into the power grid,
complex switching behaviors have been introduced as CIPS
can operate in many different modes such as grid-feeding,
grid-forming and grid-supporting [1]. The complex hybrid
behaviors from integrated CIPS will have more significant
impact on the traditional grid due to increasing penetration [2].
Analysis of the existing modes and corresponding switchings
is important to understand system limits and guidelines for
control design.
Wind power is a dominant source among all renewable
sources. Variable speed wind turbine generators (WTG) are
mechanically decoupled from the grid and do not automatically respond to frequency changes. With increasing penetra-

H

This work was supported in part by the National Science Foundation under
grant No CNS-1239366, in part by the National Science Foundation under
grant No ECCS-1509114, and in part by the Engineering Research Center
Program of the National Science Foundation and the Department of Energy
under NSF Award Number EEC-1041877.
The authors are with the Min H. Kao Department of Electrical Engineering
and Computer Science, The University of Tennessee, Knoxville, TN 37996
USA (e-mail: yzhan124@utk.edu).

tion of wind power, the natural frequency support of traditional
synchronous generators is decreasing. Configuring the WTG
controls to participate in frequency control could benefit power
system dynamics by reducing frequency excursions [3]. This
participation can be realized by adding supplementary control
loops to the normal maximum power point tracking (MPPT)
mode of the WTG as shown in Fig. 1 [4].
Wind farm should operate at MPPT mode during most
times for efficient energy extraction, but provide enough active
power to form a synthetic inertia during certain events to
ensure system frequency stay within safe limits to avoid triggering protection [5]. Such performance guaranteed control
concept have proposed as a new objective for highly controllable converters [6] [7]. The physical component corresponds
to this hybrid dynamics is a deadband. It is necessary for
efficient operation by guaranteeing more power extraction
from the wind and less mechanical stress on the gearbox;
however, a large deadband may limit the opportunity for the
WTG to provide sufficient inertia during a disturbance [8].
This is a crucial trade-off between economics and reliable
operations.
The aforementioned issues lead to the following two questions: Under a certain disturbance, can the designed inertia
emulation preserve the desired frequency limits? If so, what
is the largest deadband that preserves these limits? These
questions arise from actual power system operations faced
by transmission system operators (TSO) such as, HydroQuébec (HQ) [9]. However, as pointed out in [7], available
time (equivalent to deadband setting) for CIPS to maintain
bounded frequency is usually unknown. Few methods have
been proposed to answer the above questions beyond extensive
simulations.
In this paper, we propose a systematic theoretical analysis
to find precise answers to the above questions by considering
hybrid dynamic models. Based on selected modal analysis
(SMA) [10], a computationally truncable reduced-order model
is obtained and the above questions become tractable and
solvable based on the barrier certificate approaches for hybrid
system safety verification [11]. In addition, the synthetic
inertia and damping provided by CIPS is derived based on
the reduced-order model.
A. Related works
1) Hybrid system safety verification: A hybrid system consists of continuous dynamic subsystems and discrete events

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

w

Pref

MPPT with Full Load

Frequency Support
with Partial Load

Frequency Support
with Full Load

Saturated Frequency
Support with Partial
Load

Saturated Frequency
Support with Full
Load

Pref

r

w

r

(Optimal Tracking)

Dw

MPPT with Partial
Load

2

d
dt

K ie

Support
Deactivation

K pc
D w , D Pm , D Pv , D w

Emergency
Shunt down
r

Fig. 2. Mode transitions of a WTG.
Fig. 1. Active power control of wind turbine generator with inertia emulation
and primary frequency control.

that capture interactions between them. Safety verification of
hybrid systems combines the automatic verification techniques
for finite state concurrent systems (so-called model checking
techniques) and computation of reachable sets for continuous
dynamic systems [12]. The mode transitions of controller
in Fig. 1 is shown in Fig. 2. Some switchings take place
autonomously due to physical limits (solid lines), while others
(dash lines) are designed for specific purposes, which are in the
scope of this paper. Due to the relatively simple transition map,
safety verification reduces down to reachable sets computation
under different vector fields.
The different approaches in reachability analysis of continuous dynamic system can be categorized into Lagrangian
and Eulerian methods [13]. Lagrangian methods seek efficient
over-approximation of the reachable set by propagating certain
initial sets represented usually by polygons or ellipses under
system vector field. Lagrangian methods are computationally
feasible for high-dimensional systems, and have been successfully applied to large scale power systems [14] [15] [16].
These approximations lack accuracy when the shape of the
reachable set is not a polygon or an ellipse. On the other
hand, the goal of Eulerian methods (also known as level set
method) is to calculate as closely as possible the true reachable
set by computing a numerical solution to the Hamilton-Jacobi
partial differential equation (HJ PDE), where the initial sets
are implicitly represented by zero sublevel sets of an appropriate function. This is known as convergent approximation.
Transient [17] [18] and voltage [19] stability can be precisely
analyzed with the help of this method. To obtain numerical
solutions to the HJ PDEs, one needs to discretize the state
space, which leads to exponentially increasing computational
complexity and limits its application to systems with no more
than four continuous variables [14].
If the system dynamics and safety specifications can be
represented as polynomials, references [20] and [11] propose a
passivity-based approach that formulates safety verification as
a sum of squares (SOS) optimization problem. As long as the
SOS program is feasible, the safety property can be verified

and a polynomial barrier certificate is obtained such that no
trajectory of the system starting from the initial set can cross
this barrier to reach an unsafe region. The formulation in [11]
leads to an arbitrary barrier certificate, which can not represent the maximum safety preserving capability. The safety
supervisors for wind turbine emergency shutdown is designed
in [21] where the conservatism is reduced by maximizing
the volume of an elliptical reference shape inside the barrier
certificate. Still, no trajectory evaluation is attempted and
because of the shape limitation of the ellipsoid, it is difficult
to attain satisfactory results. In addition, the formulations
in [11] and [21] require a specific initial set. In a hybrid
system, the initial conditions after switching depend on the
switching instant. This requires a framework that builds the
condition without specific initial sets. Despite these issues,
the framework provides useful flexibility between accuracy
and computational complexity by choosing an appropriate
polynomial order.
2) Frequency control with participation of WTG considering deadband and safety limits: Accurate modeling of
deadband and other thresholds as hybrid systems can be done
using piecewise linear approximation [22]. GE has proposed
a detailed WTG model with inertia emulation and primary
frequency controllers, where the deadband is included. A
simplified structure is shown in Fig. 1 [8]. Recommended
values are 0.15 Hz for inertia emulation and 0.24 Hz for
primary frequency control, respectively, but the justification for
these values is not clear. Based on this control configuration,
a unified deadband of 0.1 Hz has been studied through
simulation [23]. In [9] Hydro-Québec identified the need for
quantifying inertia emulation with deadband under the safety
requirement. In their approach, a frequency excursion limit of
1.5 Hz is set up to prevent load shedding, then a certain amount
of inertia emulation and deadband are determined based on
simulation studies. A second drop of grid frequency will occur
once the primary frequency control of a partially loaded WTG
is deactivated [4] [8] [24]. A soft recovery procedure, i.e.,
deactivating primary frequency control of different WTGs at
different times, is proposed in [25] [26]. The safety region has
not been considered in this scenario.

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

B. Our Contributions
We propose the very concept region of safety to clarify the
hybrid system safety verification. A convergent approximation
algorithm is developed to estimate the largest region of safety
instead of achieving an arbitrary barrier certificate containing
a specific initial set. From the practical engineering point of
view, the established framework is applied to solve the hybrid
mode synthesis problem. The switching criteria or equivalent
deadbands for bounded response are built and explained.
The framework and analysis in this paper provide a general
guideline for the unclear safe switching problem.

C. Outline
The outline of the paper is as follows. In Section II, the
barrier certificate methodology and our proposed algorithm
is introduced. In Section III-A, the concept of representing
deadband as a hybrid systems is presented. Then selected
modal analysis (SMA) based model reduction is employed.
The controller gain is preserved under certain assumptions and
their equivalent emulated inertia and load-damping coefficient
is calculated. A case study is presented in Section IV and
followed by the conclusions in Section V.

3

Theorem 3: Let the system ẋ = f (x, d), and the sets X ,
XI , XU and D be given, with f continuous. If there exists a
differentiable function B : Rn −→ R such that
B(x) ≤ 0

∀x ∈ XI

(2)

B(x) > 0

∀x ∈ XU

(3)

∂B
f (x, d) < 0
∀(x, d) ∈ X × D s.t. B(x) = 0 (4)
∂x
then the safety of the system in the sense of Definition 1 is
guaranteed [11].
B(x) is called a barrier certificate. The zero level set
of B(x) defines an invariant set containing XI , that is, no
trajectory starting in XI can leave. Thus, XI is a region of
safety (ROS) due to the existence of B(x). Eq. (4) relaxes
the passivity condition from the state space to the zero level
set of B(x) and thus, reduces conservatism. The other source
of conservatism is the initial set XI usually represented by
a ball containing the equilibrium point. However, this set
could change under disturbances. Based on this observation,
we propose to solve the following problem.
Problem 4: Let ẋ = f (x, d), X , XU and D be given. The
region of safety XI is obtained by solving:
max

XI ,B(x)

Volume(XI )

subject to
II. P RINCIPLE

OF

S AFE M ODE S WITCHING S YNTHESIS

Safety denotes the property that all system trajectories stay
within given bounded regions, thus, equipment damage or
relay trigger can be avoided (Note this is similar, but not
identical, to what is called security in power industry but
for purposes of this paper we will assume satisfying safety
conditions ensures secure operation). Consider the dynamics
of a power system governed by a set of ordinary differential
equations (ODEs) as
ẋ(t) = f (x(t), d(t))

(1)

where x(t) ∈ Rn denotes the vector of state variables and
d(t) ∈ Rm denotes certain disturbances, such as, generation
loss or an abrupt load change. Such a disturbance may be
assumed to be piecewise constant and bounded in the set D.
Let X ⊆ Rn be the computational domain of interest, XI ⊆ X
be the initial set, XU ⊆ X be the unsafe set, X(XI , t, d(t))
be the set of trajectories initialized in XI . Then the formal
definition of the safety property is given as follows.
Definition 1 (Safety): Given (1), X , XI , XU and D, the
safety property holds if there exists no time instant T ≥ 0 and
no piecewise constant bounded disturbance d : [0, T ] −→ D
such that X(XI , t, d(t)) ∩ XU 6= ∅ for any t ∈ [0, T ].
Definition 2 (Region of Safety): A set that only initializes
trajectories with the property in Definition 1 is called a region
of safety.
Within the bounded set D, the safety property above is
defined in the worse-case scenario as well as the region of
safety since there is no further limits on disturbance value.
Then safety can be verified by the following theorem.

B(x) ≤ 0 ∀x ∈ XI
B(x) > 0 ∀x ∈ XU
∂B
f (x, d) < 0 ∀(x, d) ∈ X × D
∂x

s.t.

B(x) = 0

Remark 5: The importance of introducing the concept of
ROS and Problem 4 stems from the fact that we have need to
work with the initial sets instead of the invariant sets. Consider
the invariant sets {x ∈ Rn : Bi (x) ≤ 0} and ROS XI,i (green
regions) with i = 1, 2 calculated under different vector field (or
modes) f1 (x) and f2 (x) with the same safety limits and D as
shown in Fig. 3. Consider set up an emergency alert for f1 (x).
Based on Theorem 3 once the trajectory crosses B1 (x) = 0,
the alert is triggered. The dynamics after the alert is not our
concern (black dash lines in Fig. 2). The safety supervisor
of wind turbine shutdown in [21] is based on this principle.
Now consider a transition from f1 (x) to f2 (x), since the state
variables are continuous, they will evolve according to f2 (x)
beginning at the last points before transition. As a result, safety
can be guaranteed only if this initial value belongs to the ROS
under f2 (x). In Fig. 3, only the point b is safe under f1 (x)
and after switching. The point a is safe under f1 (x) but unsafe
after the transition to f2 (x). The point c is a safe switching
point but not safe under f1 (x).
Proposition 6: In a hybrid system with several modes, safe
switching to mode i is guaranteed if the trajectory of the
current mode belongs to the ROS of mode i. Moreover, if ROS
is represented by some sublevel set of a continuous function
in terms of system states, then this function represents a safe
switching guard.
Proposition 6 illustrates the fundamental principle for
switching analysis in a hybrid system. Hence, the purpose

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

of Problem 4 is to lessen conservatism. Before we introduce
the proposed iterative algorithm to approximate the solution to
Problem 4, the computational techniques for Theorem 3 under
polynomial data is introduced. Polynomial data denotes that
all sets are basic closed semi-algebraic sets (hence defined by
finitely many polynomial inequalities and equality constraints)
and vector fields are polynomial [27]. Then the property that
the polynomials are non-negative on the given semi-algebraic
sets can be checked by sums of squares decomposition, which
can be further converted to semidefinite programming (SDP)
[20]. From now on all functions are assumed to be polynomial
unless specified otherwise. Then the conditions in Theorem 3
can be written into a sums of squares programming (SOSP)
problem. First let us denote by Σ 2 [x] the space of SOS
polynomials, and by Σ 2 [x]p the space of SOS polynomials
of degree at most 2p.
x = f1 ( x )

Vector Field Switching

x = f2 ( x)

c

c

I ,1

a

b

Region of safety: c

I ,1

Safety boundary : gU ( x ) = 0
Zero level set of barriers : Bi ( x ) = 0

(0)

−B (0) (x) − λI,i (x)gI,i (x) ∈ Σ 2 [x]
(0)

B (0) (x) − ǫ − λU (x)gU (x) ∈ Σ 2 [x]
−

•

∂B (0)
(0)
(x)f (x, d) − λD (x, d)gD (d)
∂x
(0)
−λX (x, d)gX (x) − rB (0) (x) ∈ Σ 2 [x]

Iteration k
(a) Fix the barrier certificate B (k−1) (x) from k − 1
step, solve the SOS optimization for multiplier
(k )
λB a (x, d):
∂B (k−1)
(k )
(x)f (x, d) − λD a (x, d)gD (d)
∂x
(k )
(k )
−λX a (x, d)gX (x) − λB a (x, d)B (k−1) (x) ∈ Σ 2 [x]

c

c

a

Initialization Let xi0 for i = 1, · · · , N be several initial points with safety verified, and XI,i =
{x ∈ Rn : gI,i (x) ≥ 0} represent a small ball centered at
xi0 . Choose λB (x, d) equal to a sufficiently small positive
real number r and solve the following SOS optimization
for i = 1, · · · , N :

−

B2 ( x ) = 0
B1 ( x ) = 0

•

4

I ,2

(b) Fix the barrier certificate B (k−1) (x) from k − 1 step,
(k )
the multiplier λB a (x, d) from k (a) step, solve the
following SOS optimization for B (k) (x):

b

Region of safety: c

(k)

−B (k) (x) − λI (x)B (k−1) (x) ∈ Σ 2 [x]

I ,2

(k)

B (k) (x) − ǫ − λU (x)gU (x) ∈ Σ 2 [x]

Equilibrium Points
Region of safety: c I ,i

∂B (k)
(k)
(x)f (x, d) − λD (x, d)gD (d)
∂x
(k )
(k)
−λX (x, d)gX (x) − λB a (x, d)B (k) (x) ∈ Σ 2 [x]
−

Fig. 3. Region of safety and corresponding barriers under different vector
fields.

Problem 7: Let X = {x ∈ Rn : gX (x) ≥ 0}, XI =
{x ∈ Rn : gI (x) ≥ 0}, XU = {x ∈ Rn : gU (x) ≥ 0}, and
D = {d ∈ Rm : gD (d) ≥ 0}, which are represented by the
zero superlevel set of the polynomials gX (x), gI (x), gU (x),
and gD (d), respectively, and some small positive number ǫ be
given. Then
−B(x) − λI (x)gI (x) ∈ Σ 2 [x]

(5)

B(x) − ǫ − λU (x)gU (x) ∈ Σ 2 [x]

(6)

∂B
(x)f (x, d) − λD (x, d)gD (d)
∂x
−λX (x, d)gX (x) − λB (x, d)B(x) ∈ Σ 2 [x]

(7)

−

with multipliers λI (x), λU (x), λX (x, d), λD (x, d) and
λB (x, d) SOS polynomials.
Conversion of Problem 7 to SDP has been implemented in
solvers such as SOSTOOLS [28] or the sum of squares module
[29] in YALMIP [30]. Then the powerful SDP solver MOSEK
[31] can be employed. Now let us introduce the algorithm to
approximate the solution of Problem 4.
Algorithm 8: Let X = {x ∈ Rn : gX (x) ≥ 0}, XU =
{x ∈ Rn : gU (x) ≥ 0}, D = {d ∈ Rm : gD (d) ≥ 0}, which
are represented by the zero superlevel of the polynomials
gX (x), gU (x) and gD (d), respectively, some small positive
number ǫ, initial order 2p and maximal order 2pmax for barrier
certificate computation be given.

(c) If step k (b) is feasible, then let k = k + 1. If
infeasible, then increase the polynomial order of
B (k) by two, i.e., 2p = 2p + 2. If p = pmax but
step k (b) is
 still infeasible, then the algorithm stops
and XI∗ = x : B (k−2) (x) ≤ 0 with B (k−1) (x) the
barrier.
The key idea of the proposed algorithm is to use the zero
level set of a feasible barrier certificate as an initial condition
and to search for a larger invariant set. Once feasible, this
initial condition becomes the ROS due to the existence of
corresponding invariant sets. A judicious choice of the initial
points in the initialization step can reduce the number of
iterations, and also helps to have a precise estimate in certain
sub-dimensions, if a full dimensional estimate is hard due to
computational complexity.
III. H YBRID R EDUCED - ORDER M ODEL OF WTG
SMA- BASED M ODEL R EDUCTION

VIA

Consider the active power control in Fig. 1. The WTG is
assumed to operate at partial loaded condition, which means
max
Pgen < Pgen
and ωrmin < ωr < ωrmax . And no pitch angle
control is considered. The control signal sent to the voltagesource converter (VSC) is given by
Pref = Copt ωr3 + Kie (ωgrid )∆ω̇ + Kpc (ωgrid )∆ω

(8)

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

where Copt ωr3 , Kie (ωgrid )∆ω̇, and Kpc (ωgrid )∆ω provide functionalities of maximum power point tracking (MPPT), inertia
emulation (IE) and primary frequency control (PFC), respectively. Kie and Kpc are equal to zero when | ∆ωgrid |<
∆ωdeadband , or equal to certain pre-set values otherwise. The
deadband effect is equivalent to the switching guard Gdb (x) in
the hybrid system illustrated in Fig. 4. Gda1 (x) and Gda2 (x)
represent signals (dash line in Fig. 1) for deactivating the
corresponding support in order to have a faster frequency
restoration after safety is preserved. It is worth mentioning
that if ROSs of corresponding modes are employed for those
guards, then by Proposition 6, safety can be guaranteed.
A. Model Reduction
As shown in [32], there is a trade-off between system
dimension and the order of polynomials for set representation.
Since the active power variations and frequency dynamics
in power system are dominantly governed by mechanical
dynamics and modes, a reduced-order system is desired in
analysis so that higher order polynomials can be used for
better estimation of the ROS. The selective modal analysis
(SMA) based model reduction has proven to be successful in
capturing active power dynamics of WTG [10] and is chosen
for our study.

where
Asys = As − Bs Ds−1 Cs

Csys = Es − Fs Ds−1 Cs

Bsys1 = Ms1 − Bs Ds−1 Ns1

Dsys1 = −Fs Ds−1 Ns1

Bs Ds−1 Ns2

Dsys2 = −Fs Ds−1 Ns2

Bsys2 = Ms2 −

The WTG rotor speed ∆ωr dynamic is closely related to
its active power output, and the mode where ∆ωr has the
highest participation would capture the relevant active power
dynamics. Therefore, ∆ωr is considered as the most relevant
state, and the other states are less relevant and denoted as z(t).
Eq. (12) can be rearranged as

 


∆ω˙r
A11 A12
∆ωr
=
ż
A21 A22
z




(13a)
Br2
Br1
∆ω
∆ω̇ +
+
Bz2
Bz1


∆ωr
∆Pgen = [Cr Cz ]
z
(13b)
+ Dsys1 ∆ω̇ + Dsys2 ∆ω
The less relevant dynamics are:
ż = A22 z + A21 ∆ωr + Bz1 ∆ω̇ + Bz2 ∆ω

MPPT

IEPFC

IE

(15)

In (15), z can be represented by the following expression:
Z t
A22 (t−t0 )
eA22 (t−τ ) A21 ∆ωr (τ )dτ
z(t) = e
z(t0 ) +
t0
{z
}
|

MPPT

K ie
K pc

response without control input

+

Z

t

t0

Fig. 4. Modelling deadband as a hybrid transition system.

Consider a type-3 WTG connected to a reference bus. The
detailed model and the meaning of each variable is given in
Appendix B. Linearize (33), (34), (35) and (36) about the
equilibrium point given in Appendix A. Keep ∆ω̇ and ∆ω
as input variables and omit the variations of vwind . Then the
following model is obtained:



 
 
∆ẋ
∆x
As Bs
Ms1
∆ω̇
=
+
0
∆y
Cs Ds
Ns1


(9a)
Ms2
∆ω
+
Ns2


∆x
∆Pgen = [Es Fs ]
(9b)
∆y
where
x = [EqD , EdD , ωr , x1 , x2 , x3 , x4 ]T
y = [Vqr , Vdr , Iqr , Idr , Pgen , Qgen , Ids , Iqs , VD , θD ]

(14)

Thus, the most relevant dynamic is described by:
∆ω̇r = A11 ∆ωr + A12 z + Br1 ∆ω̇ + Br2 ∆ω

Gda 2 ( x )

Gda1 ( x )

Gdb ( x )

5

(10)
T

(11)

Using Kron reduction [33] on Eq. (9) yields the following
state-space model:
∆ẋ = Asys ∆x + Bsys1 ∆ω̇ + Bsys2 ∆ω

(12a)

∆Pgen = Csys ∆x + Dsys1 ∆ω̇ + Dsys2 ∆ω

(12b)

+

|

eA22 (t−τ ) Bz1 ∆ω̇(τ )dτ
{z
}

(16)

response under inertia emulation
t
eA22 (t−τ ) Bz2 ∆ω(τ )dτ
t0

Z
|

{z

}

response under primary frequency control

Using the most relevant mode, ∆ωr (τ ) can be expressed as
[10]:
∆ωr (τ ) = cr vr eλr τ
(17)
where λr is the relevant eigenvalue, vr is the corresponding
eigenvector and cr is an arbitrary constant. The accuracy of
(17) is guaranteed by the dominant term of ∆ωr , which can be
used in solving the first integral in (16). Since A22 is Hurwitz
and its largest eigenvalue is much smaller than λr , the natural
response will decay faster and can be omitted. The essential
reason is that A22 represents electrical dynamics which are
faster than the electro-mechanical dynamic represented by
λr . Then the response without control input in (16) will
approximately equal to the forced response represented as
follows:
Z t
A22 (t−t0 )
(18)
eA22 (t−τ ) A21 ∆ωr (τ )dτ
e
z(t0 ) +
t0
{z
}
|
response without control input

≈ (λr I − A22 )−1 A21 ∆ωr

(19)

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

The rate of change of frequency (RoCoF) ∆ω̇ and the stabilized frequency deviation ∆ω are assumed to be fixed during
the time window of interest, then the other two integrals are
easily calculated as
Z t
eA22 (t−τ ) Bz1 ∆ω̇(τ )dτ ≈ (−A22 )−1 Bz1 ∆ω̇
(20)
t0
|
{z
}
response under inertia emulation
t
eA22 (t−τ ) Bz2 ∆ω(τ )dτ
t0

Z

{z

|

≈ (−A22 )−1 Bz2 ∆ω

(21)

}

Finally, the reduced-order WTG model with control inputs is
∆ω̇r = Ard ∆ωr + Brd1 ∆ω̇ + Brd2 ∆ω

(22a)

∆Pgen = Crd ∆ωr + Drd1 ∆ω̇ + Drd2 ∆ω

(22b)

where
Ard = A11 + A12 (λr I − A22 )−1 A21
Crd = Cr + Cz (λr I − A22 )−1 A21
Brd1 = Br1 + A12 (−A22 )−1 Bz1
Drd1 = Dsys1 + Cz (−A22 )−1 Bz1
Brd2 = Br2 + A12 (−A22 )−1 Bz2
Drd2 = Dsys2 + Cz (−A22 )−1 Bz2
B. Quantification of Frequency Support from WTG
Consider the swing equation with the active power increment ∆Pgen from WTGs
ωs
D
∆ω)
(∆Pm + ∆Pgen − ∆Pe −
2H
ωs

t0

(24)
During the time window of inertia response Th =
{t : 0 ≤ t ≤ th }, the RoCoF is approximately fixed. Then ∆ω̇
can be pulled outside the integral and integrating (24) with
t0 = 0 yields
(25)

Substituting (25) back into (23) and rearranging the state yields
∆ω̇ =

D
ωs
∆ω)
(∆Pm − ∆Pe −
2H + 2He (t)
ωs

(26)

where
Ard t
He (t) = 0.5ωs (−Drd1 + Crd A−1
)Brd1 )
rd (I − e

t0

(28)
After the frequency is stabilized by the governor, i.e., t ∈ Tp =
{t : tp ≤ t ≤ ts }, the term ∆ω can be pulled outside of the
integral and integrating (29) with t0 = tp yields
Ard (t−tp )
∆Pgen (t) = (Drd2 − Crd A−1
)Brd2 )∆ω (29)
rd (I − e

∆ω̇ =

(D + De (t))
ωs
∆ω)
(∆Pm − ∆Pe −
2H
ωs

(30)

where
Ard (t−tp )
De (t) = ωs (−Drd2 + Crd A−1
)Brd2 )
rd (I − e

(31)

Eq. (30) and (31) have clearly illustrated that PFC is actually
emulating load-damping characteristic.
Remark 9: The controller gain is preserved under the
approximations that the RoCoF ∆ω̇ and the stabilized frequency deviation ∆ω is constant within the time window of
interest. The comparison to the full order model will show
that this approximation is accurate within these windows.
Then the emulated inertia He and load-damping coefficient
De are expressed in term of Brd1 , Drd1 and Brd2 , Drd2 , which
correspond to Kie and Kpc . Note that He (t) and De (t) are
time-varying and their values are considered to be accurate
only within the corresponding time windows, i.e., Th and Tp .
IV. C ASE S TUDY

(23)

The WTG active power output in (23) due to the signal ∆ω̇
and ∆ω will influence the values of H and D independently.
To evaluate the emulated inertia, the terms related to the
PFC, i.e., Brd2 and Drd2 in (22), are set to zero. The explicit
forced output response of (22) due to ∆ω̇ is given by
Z t
∆Pgen (t) = Crd
eArd (t−τ ) Brd1 ∆ω̇(τ )dτ + Drd1 ∆ω̇(t)

Ard t
∆Pgen (t) = (Drd1 − Crd A−1
)Brd1 )∆ω̇
rd (I − e

to zero. The explicit forced output response of (22) due to ∆ω
is given as
Z t
∆Pgen (t) = Crd
eArd (t−τ ) Brd2 ∆ω(τ )dτ + Drd2 ∆ω(t)

Substituting (29) into (23) yields

response under primary frequency control

∆ω̇ =

6

(27)

To evaluate the emulated load-damping effect, the terms
related to inertia emulation, i.e., Brd1 and Drd1 in (22), are set

Consider the four-bus system with a 600 MW thermal plant
made up of four identical units in Fig. 5. The frequency
dynamics of the above system can be represented by the classic
system frequency response (SFR) model [34] as follows:
D
ωs
∆ω)
(∆Pm − ∆Pe −
2H
ωs
1
∆Ṗm =
(∆Pv − ∆Pm )
τch
1
1
∆Ṗv = (−∆Pv − ∆ω)
τg
R

∆ω̇ =

(32a)
(32b)
(32c)

The power flow equation is ∆Pe = ∆Pd − ∆Pgen , where ∆Pd
denotes a large disturbance, such as, generation loss or abrupt
load changes, and ∆Pgen given in (22b) represents the active
power variation due to the frequency control loop. Although
a one-area case is studied, the SFR model has the potential
to describe system frequency response in a complex power
network as shown in many recent studies [35] [36] [37]. In
this one-area system, as shown in Fig. 6, the response of the
SFR model and two-axis nonlinear model are the same under
the same disturbance.
The wind farm is assumed to be an aggregation of 200
individual GE 1.5 MW WTGs with rated speed of 450 rad/s
(or 72 Hz) and rated output of 300 MW. Under the operating
condition given in Appendix A, the reduced-order WTG model
can be obtained with Ard = −0.0723 and Crd = 0.0127. Brd1 ,

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

6.5
6

3.5
Mode 4
Mode 5
Emulated Load−damping

Emulated Inertia

Frequency (Hz)

5.5
59.8

59.7

59.6

5
4.5
4
3.5

0

Mode 2
Mode 3

3

Frequency Deviation (Hz)

3rd-order SFR model
Two-axis model

60

59.9

2.5
2
1.5
1

10

15

20

25

30

Time (s)

2.5
0

−0.2
−0.3
−0.4

0.5

3
59.4
5

Mode 1
Mode 2
Mode 3
Mode 4
Mode 5
Nadir Limit

−0.1

−0.5

59.5

Fig. 5. Case study: four-bus system
[10].

7

0.5
1
Time (s)

0

1.5

0
10

20
Time (s)

1

2

3

4

5

Time (s)

30

Fig. 6. Response comparison of sys- Fig. 7. Time-varying emulated inertia Fig. 8. Frequency response of different modes under the worst-case scetem frequency response model and and load-damping coefficient.
nario: 150 MW generation loss.
two-axis model.

as the true one. The comparison shows that the proposed
algorithm successfully reduces conservatism in the estimate
for the corner effect in the polynomial-based set study.
0.5
0.4
Mechanical Power (P. U.)

Drd1 and Brd2 , Drd2 with the corresponding Kie and Kpc are
listed in Table I. The emulated He and De are time-varying
and shown in Fig. 7.
The worst-case scenario is assumed to be the loss of
one unit (150 MW), which occurs at 1 s. The safety limit
is set to be a 0.5 Hz deviation [38] to avoid triggering
load shedding [5]. The frequency response of all modes
under this scenario is given in Fig. 8. The inertia emulation
effect can be observed as the RoCoF becomes slower from
the response of Modes 1-3. The ROS (safety switching
guard) is calculated under the reduced-order model in Eq.
(32), but the full-order linearized model in Fig. 5 is used
∆Pm , ∆Pv , ∆ωr ] and x =
for verification. Denote′ xrd = [∆ω,
′
∆ω, ∆Pm , ∆Pv , ∆EqD , ∆EdD
, ∆ωr , ∆x1 , ∆x2 , ∆x3 , ∆x4
for theoretical analysis and simulation verification,
respectively.

MPPT
IE
IE
IEPFC
IEPFC

1
2
3
4
5

Kie
0
-0.10
-0.20
-0.10
-0.20

Brd1

Drd1

Kpc

Brd2

Drd2

0
0.6246
1.2492
0.6246
1.2492

0
-0.10
-0.20
-0.10
-0.20

0
0
0
-0.03
-0.06

0
0
0
0.1874
0.3748

0
0
0
-0.03
-0.06

0.1
0
−0.1
−0.2

−0.4
−0.5

−0.4

−0.2
0
0.2
System Frequency (Hz)

0.4

0.6

Fig. 9. Iteration in calculating ROS with different initializations.

0.6

Mechanical Power Deivation (PU)

Number

0.2

−0.3

TABLE I
G AIN OF F REQUENCY S UPPORT M ODE AND C ORRESPONDING M ATRIX
VALUE
Mode

0.3

ROS by massive simulations
ROS by proposed Algorithm 8

0.4
0.2
0
−0.2
−0.4
−0.6

A. Model and Algorithm Validation
To validate the reduced-order model, consider the worstcase scenario above. The four state variables ∆ω, ∆Pm , ∆Pv ,
∆ωr between reduced-order and full-order model of Mode 2-5
in Table I are compared in Fig. 11. The excellent agreement in
mode behaviour ensures the reduced-order model based ROS
should be sufficient to find the switching for the full-order
dynamics.
With the given safety limit, the ROS calculation for Mode
1 under no disturbance can be projected onto the plane ∆ω∆Pm as illustrated in Fig. 9 with two different initializations.
The iteration sequence indicates that if more initial guess
points are used, the fewer iterations needed and a better
estimation can be achieved (as shown in the blue case). The
final result is shown in Fig. 10. The green region is the
ROS obtained by extensive simulations and can be regarded

−0.4

−0.2

0

0.2

0.4

0.6

System Frequency (Hz)

Fig. 10. ROS of Mode 1 under normal condition obtained by proposed
Algorithm 8 and extensive simulations.

B. IE Mode Only
The ROS under the worst-case scenario of Modes 1-3 are
calculated with representation of polynomials in terms of xrd
up to degree 8. Denote these regions as
Worst-case ROS 1: Sd1 = {xrd : Bd1 (xrd ) ≤ 0}
Worst-case ROS 2: Sd2 = {xrd : Bd2 (xrd ) ≤ 0}
Worst-case ROS 3: Sd3 = {xrd : Bd3 (xrd ) ≤ 0}
where Bd1 (xrd ), Bd2 (xrd ) and Bd3 (xrd ) serve as safety switching guards.

4

5

Reduced−order
Full−order

−0.1
−0.2
−0.3
−0.4
−0.5

0

1

2
3
Time (s)

4

5

−0.5

0

1

2
3
Time (s)
Mode 5

4

−0.2
−0.3
−0.4
−0.5

0

1

2
3
Time (s)

4

−2

0

1

2
3
Time (s)
Mode 4

4

Reduced−order
Full−order

0

−2

−4

0

1

2
3
Time (s)

4

(a)

Mode 2

Mode 3

−4

0

1

2
3
Time (s)
Mode 5

4

2

0

5
Time (s)
Mode 4

Reduced−order
Full−order

0

−4

0

1

2
3
Time (s)

4

5

0

5
Time (s)
Mode 5

Reduced−order
Full−order

0.1

0

5
Time (s)

(b)

10

0.1

0

5
Time (s)
Mode 4

Reduced−order
Full−order

0.1

0

5
Time (s)

10

0.1

0

0

5
Time (s)
Mode 5

10

0.3
Reduced−order
Full−order

0.2

0.1

0

Reduced−order
Full−order

0.2

10

0.3

0.2

0

0.2

0

10

0.3

0.2

0

0.1

0

10

0.3

−2

−6

0.1

0

5

0.2

Reduced−order
Full−order

Valve Position

−2

0.2

Reduced−order
Full−order

Valve Position

5

2

Mode 3

Reduced−order
Full−order

Mechnical Power

5

Reduced−order
Full−order

−0.1

5

−1

0

Mechnical Power

5

−0.4

Mode 2

Reduced−order
Full−order

Valve Position

2
3
Time (s)
Mode 4

0

Mode 3

8

2

Valve Position

1

Reduced−order
Full−order

1

Mechnical Power

0

−0.3

2

Mechnical Power

−0.5

−0.2

Mode 2

WTG Rotor Speed (Rad/s)

−0.4

Reduced−order
Full−order

−0.1

WTG Rotor Speed (Rad/s)

−0.3

WTG Rotor Speed (Rad/s)

−0.2

Mode 3

WTG Rotor Speed (Rad/s)

Reduced−order
Full−order

Frequency Deviation (Hz)

Mode 2
−0.1

Frequency Deviation (Hz)

Frequency Deviation (Hz)

Frequency Deviation (Hz)

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

0

5
Time (s)

0.1

0

10

(c)

Reduced−order
Full−order

0.2

0

5
Time (s)

10

(d)

Fig. 11. Dynamics between full-order and reduced-order model: (a) Frequency deviation; (b) WTG rotor speed; (c) Turbine-governor mechanical power; (d)
Turbine-governor valve position.

To determine if the safety limits can be preserved, one
needs to check whether the intersection between Sd1 and the
pre-disturbed operating point x0 is empty. In our case, the
fact that Sd1 ∩ {x0 } = ∅ is graphically shown in Fig. 15
and mathematically verified by Bd1 (x0 ) > 0. According to
Proposition 6, the safety of frequency cannot be preserved
without inertia emulation as shown in Fig 8.
To verify the largest deadband or equivalently the critical
switching instant from Mode 1 to Mode 2 or 3, the values of
Bd2 (xrd ) and Bd3 (xrd ) with respect to the disturbed trajectory
of Mode 1, denoted as Xd1 (dash line in Fig. 12), is calculated.
Note that Xd1 is from the full-order model and only relevant
states X̄d1 = [Xd1 (1), Xd1 (2), Xd1 (3), Xd1 (6)] are substi-

0

1

1.5
Time (s)

2

2.5

3

Certificate Value

0.4

0
−0.2
−0.4

−0.2
0
0.2
System Frequency (Hz)

0.4

ROS of Mode 2 and 3 under given

−50

−0.3

−0.4

−0.5

0.5

1

1.5
Time (s)

2

2.5

3

Fig. 13. Value of Bd2 (x) (upper) and Bd3 (x)
(lower) w.r.t the disturbed trajectory Xd1

←

Critical DB 0.42 Hz

−0.5
0

1

2
3
Time(s)

4

5

5
Time(s)

10

−0.2
Mode 1
Mode 2
Mode 1(Unsafe)

−0.25

zero−crossing at 2.02s

0

1

2

3

4

5

500
zero−crossing at 2.06s

0

1

2

3

4

5

Time (s)

Fig. 16. Value of Bd1 (x) w.r.t. the disturbed
trajectories Xd12 (upper) and Xd13 (lower).

Mode 1
Mode 3
Mode 1(Unsafe)

−0.25
−0.3

−0.3

−500
0

0

Fig. 14. Frequency dynamics of full-order model
under calculated critical deadband.

500

Time (s)
Certificate Value

← Critical DB 0.3 Hz

−0.3

← Recommended
DB 0.15 Hz
−0.2

zero−crossing at 1.44s

−500
0

Fig. 15. ROS of Mode 1 under the given scenario.

−0.2

−0.4

0

0

0.6

Certificate Value

Fig. 12.
scenario.

−0.4

50

Frequency(Hz)

0.2

−0.1
← Recommended
DB 0.15 Hz

Frequency(Hz)

0.5

Mode 1
Mode 3 (Critical)
Mode 3 (Recom.)

0

−0.1

zero−crossing at 1.29s
−20
0

0.6

Mode 1
Mode 2 (Critical)
Mode 2 (Recom.)

0

20

Frequency(Hz)

1
0.8

40

Frequency(Hz)

ROS of Mode 2
ROS of Mode 3
Disturbed Trajectory of Mode 1

1.2
Mechanical Power (P. U.)

Certificate Value

1.4

tuted into the guards. The zero-crossing point from negative
to positive values denotes the critical switching instant tcr ,
or equivalently largest deadband with the value ∆ω(tcr ). As
shown in Fig. 13, the largest deadband (critical switching
instant) is 0.30 Hz (1.29 s) if Mode 2 is used, and 0.42 Hz
(1.44 s) if Mode 3 is used. Simulation of each scenario with
the suggested largest deadband as well as the recommended
value from GE (0.15 Hz) is carried out and shown in Fig. 14.
As seen the system safety is preserved, but the recommended
values are conservative especially when Mode 3 is activated.
On the other hand, the fact that the largest frequency excursion
point is extremely close to the limit indicates that the estimated
ROS is not overly conservative.

−0.35
−0.4

−0.35
−0.4

−0.45

−0.45

−0.5

−0.5

1

2
3
Time(s)

4

1

2
3
Time(s)

4

5

Fig. 17. Frequency dynamics under critical deadband and threshold.

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

The additional PFC loop will artificially create additional
load-frequency sensitivity and the maximum frequency excursion will decrease. The deadband analysis procedure is
similar and will not be repeated. However, when it comes to
support deactivation, due to the additional constant frequency
deviation, a safe switching time window appears. Thus, the
PFC mode needs to be deactivated before a critical time. The
mechanism is illustrated in Fig. 18.
The WTG attempts to draw the total energy that is pulled
out during support mode to restore the rotor speed so the
upper area and lower areas have to be equal. When the PFC is
deactivated, Pgen decreases to satisfy this equal area criterion.
This sudden shortage of active power, if large enough, will
lead to an unsafe trajectory. Let Mode 5 be designed with
a deadband of 0.3 Hz and substitute the disturbed trajectory
X̄d15 into Bd1 (xrd ), then the critical switching window can
be observed as in the upper plot of Fig. 19, where the critical
deactivating instants suggested by the guard is 15.2 s. A
deactivation at 22 s leads to an unsafe trajectory. Frequency
dynamics of both cases are shown in the upper plot of Fig.
20.

WTG Active Power Output

Switching before tcr
Switching after tcr
No Support Control
Maximum Power Outage for Safety

2

3

4

5

6

7

Certificate Value

50
0

zero−crossing at 1.85s
zero−crossing at 15.20s

−50

Certificate Value

5

10

15

20
25
Time (s)

30

35

40

35

40

20
10
zero−crossing at 30.21s

0
−10
−20
−30
0

5

10

15

20
25
Time (s)

30

Fig. 19. Upper: Value of Bd1 (x) and Bd3 (x) w.r.t. the disturbed trajectory
Xd15 .

−0.3

Mode 1
Mode 5
Mode 1(Unsafe)

−0.35
−0.4
−0.45
−0.5
0

5

10

15

20
Time(s)

25

30

35

40

Mode 1
Mode 5
Mode 3
Mode 3(Critical)

−0.3
−0.35
−0.4
−0.45
−0.5
0

5

10

15

20
Time(s)

25

30

35

40

Fig. 20. Frequency dynamics under 0.3 Hz deadband and critical deactivation:
normal sequence Mode 1-5-1 (Upper) and safe recovery sequence Mode 1-53-1 (lower)

V. C ONCLUSION

Switching before tcr

1

100

−100
0

Frequency(Hz)

C. IEPFC Mode and Safety Recovery

instant from Mode 5 to 3 is suggested to be 30.21 s by the
guard in Fig. 19 and verified by simulation in Fig. 20.

Frequency(Hz)

Beyond safety, the earliest support deactivation (earliest
switching instant) is established so that the emulated inertia
can be shed to obtain faster frequency restoration. Let the
trajectories in Fig. 14 be denoted as Xd12 (left) and Xd13
(right) and substituted into Bd1 (xrd ) to find the zero-crossing
point from positive to negative, which occurs approximately
at 2 s for both cases. Negativity of the safety switching guard
Bd1 (X̄d1i ) guarantees safe switching from Mode i to Mode 1.
Early switching when Bd1 (X̄d1i ) > 0 will lead to an unsafe
trajectory. Both cases are shown in Fig. 17.

9

8

Time (s)

Fig. 18. An equal-area criterion in support-deactivation procedure.

To extend the critical deactivating instant, we propose a
safety recovery procedure illustrated by Fig. 4. When the PFC
mode is deactivated, the corresponding IE mode is kept to
manage the sudden shortage of active power. By checking the
value of Bd3 (xrd ) with respect to the trajectory X̄d15 (lower
plot of Fig. 19) this procedure extends the critical deactivating
instant by 15 s. The original unsafe switching from Mode 5
directly to Mode 1 at 22 s is now safely switched to Mode 3
as shown in the lower plot of Fig. 20. The critical switching

This paper proposes a principle of WTG mode switching
synthesis for safe frequency response. Given desired safety
limits of grid frequency and worst-case scenarios, an SOS
optimization based algorithm is developed to enlarge the
estimation of ROS, and then the critical switching instants,
equivalent to largest deadband, for safety preservation are
obtained. Simulation shows that the critical switching instants
are not overly conservative. A switching sequence for safe
recovery of WTG rotor speed is proposed as well. In addition,
the emulated inertia and load-damping effect is derived in
the time frame of inertia and primary frequency response,
respectively.
Future study will focus on analyzing larger power networks.
From a methodology point of view, instead of choosing
monomials as a dense basis to represent the set of continuous
functions on compact sets, other representations such as the
Handelman representation can be employed, which will lead
to linear relaxations of polynomial positivity rather than SOS
relaxations [39]. Coordination of multiple renewable sources
as a redundant actuator set under a severe contingency can be
considered as well [40].

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

A PPENDIX A
S YSTEM PARAMETERS AND O PERATING C ONDITION
Type-3 wind turbine generator parameters
Xm = 3.5092, Xs = 3.5547, Xr = 3.5859, Rs = 0.01015,
Rr = 0.0088, HD = 4(s), p = 4, ρ = 1.225(kg/m3 ), Rt =
38.5(m), SbD = 1(MVA), Copt = 3.2397 × 10−7 (s3 /Hz3 ),
k = 1/45, KP 1 = KP 2 = KP 3 = KP 4 = 1, KI1 = KI2 =
KI3 = KI4 = 5.
System frequency response model parameters
ωs = 60(Hz), D = 1, H = 4(s), τch = 0.3(s), τg = 0.1(s),
R = 0.05.
Network parameters and operating condition
Xt = 0.07, V̄ = 1, θ̄ = 0(rad), θ̄t = 0(rad), v̄wind = 12(m/s),
ω̄r = 72(Hz), P̄gen = 0.3, Q̄set = 0, Sb = 1000(MVA)

A PPENDIX B
W IND T URBINE G ENERATOR M ODEL
The standard wind turbine model is given as follows [10]
2kωr Rt
pvwind
−1

0.035
1
− 3
λi =
λ + 0.08θt θt + 1


116
− 12.5
− 0.4θt − 5 e λi
Cp =0.22
λi
3
1 ρπRt2 ωb Cp vwind
Tm =
2
Sb ω r
λ=

(33a)

Xm
1
′
(EqD
+ (Xs − Xs′ )Ids ) + ωs
Vdr
′
T0
Xr
′
− (ωs − ωr )EdD
1
Xm
′
′
ĖdD
= − ′ (EdD
+ (Xs − Xs′ )Iqs ) + ωs
Vqr
T0
Xr
′
− (ωs − ωr )EqD
ωs
′
′
(Tm − EdD
Ids − EqD
Iqs )
ω̇r =
2HD
ẋ1 =KI1 (Pref − Pgen )
ẋ2 =KI2 (KP 1 (Pref − Pgen ) + x1 − Iqr )
ẋ3 =KI3 (Qref − Qgen )
ẋ4 =KI4 (KP 3 (Qref − Qgen ) + x3 − Idr )

0 =KP 2 (KP 1 (Pref − Pgen ) + x1 − Iqr )
+ x2 − Vqr

(35a)

0 =KP 4 (KP 3 (Qref − Qgen ) + x3 − Idr )
+ x4 − Vdr

(35b)

′
′
2
2
0 = − Pgen + EdD
Ids + EqD
Iqs − Rs (Ids
+ Iqs
)

− (Vqr Iqr + Vdr Idr )

(35c)

′
′
2
2
0 = − Qgen + EqD
Ids + EdD
Iqs − Xs′ (Ids
+ Iqs
) (35d)
′
EqD
Xm
+
Ids
(35e)
0 = − Idr +
Xm
Xr
E′
Xm
0 = − Iqr − dD +
Iqs
(35f)
Xm
Xr
where Vdr , Vqr , Idr , Iqr are rotor d q axis voltage and current,
respectively. Vds , Vqs , Ids , Iqs are stator d q axis voltage and
current. Pgen and Qgen are WTG active and reactive power
output. VD and θD are voltage magnitude and angle of the
bus which WTGs are connected to.
The network algebraic equations are given as follows
′
′
EqD
− jEdD
= (Rs + jXs′ )(Iqs − jIds ) + VD

(36a)

jθ

(36b)

VD e

jθD

= jXt (Iqs − jIds − IGC )e

jθD

+Ve

where
Vdr Iqr + Vdr Idr
VD

ACKNOWLEDGMENT
(33c)
(33d)

The type-3 wind turbine generator differential equations are
given as follows [10]
′
ĖqD
=−

The type-3 wind turbine generator algebraic equations are
given as follows [10]

IGC =
(33b)

10

(34a)

(34b)
(34c)
(34d)
(34e)
(34f)
(34g)

′
′
where EdD
, EqD
and ωr are d q axis voltage and rotor speed of
wind turbine generator, respectively. x1 to x4 are proportionalintegral regulator induced states. And Pref = Copt ωr3 , Qref =
X2
r
and Xs′ = Xs − Xmr .
Qset , T0′ = ωX
s Rr

This work was supported in part by the National Science
Foundation under grant No CNS-1239366, in part by the
National Science Foundation under grant No ECCS-1509114,
and in part by the Engineering Research Center Program of the
National Science Foundation and the Department of Energy
under NSF Award Number EEC-1041877.
R EFERENCES
[1] J. Rocabert, A. Luna, F. Blaabjerg, and P. Rodriguez, “Control of power
converters in ac microgrids,” IEEE Trans. Power Electron., vol. 27,
no. 11, pp. 4734–4749, 2012.
[2] “Global wind report, annual market update,” Global Wind Energy
Council, Tech. Rep., 2014.
[3] W. Lu and D. G. Infield, “Towards an assessment of power system frequency support from wind plant – modeling aggregate inertial response,”
IEEE Trans. Power Syst., vol. 28, no. 3, pp. 2283–2291, 2013.
[4] J. Morren, S. W. De Haan, W. L. Kling, and J. Ferreira, “Wind turbines
emulating inertia and supporting primary frequency control,” IEEE
Trans. Power Syst., vol. 21, no. 1, pp. 433–434, 2006.
[5] “Western electricity coordinating council off-nominal frequency load
shedding plan,” Western Electricity Coordinating Council, Tech. Rep.,
2011.
[6] A. Ulbig, T. S. Borsche, and G. Andersson, “Impact of low rotational
inertia on power system stability and operation,” in Proc. 19th IFAC
World Congr., 2014, pp. 7290–7297.
[7] F. M. Uriarte, C. Smith, S. VanBroekhoven, and R. E. Hebner, “Microgrid ramp rates and the inertial stability margin,” IEEE Trans. Power
Syst., vol. 30, no. 6, pp. 3209–3216, 2015.
[8] K. Clark, N. W. Miller, and J. J. Sanchez-Gasca, “Modeling of GE
wind turbine-generators for grid studies, version 4.5,” General Electric
International, Inc., Tech. Rep., 2010.
[9] J. Brisebois and N. Aubut, “Wind farm inertia emulation to fulfill hydro–
qubec’s specific need,” in Proc. 2011 IEEE Power and Energy Soc.
General Meeting, 2011, pp. 1–7.

Published on IEEE JOURNAL ON EMERGING AND SELECTED TOPICS IN CIRCUITS AND SYSTEMS in September (DOI: 10.1109/JETCAS.2017.2675879)

[10] H. A. Pulgar-Painemal, “Wind farm model for power system stability
analysis,” Ph.D. dissertation, Univ. of Illinois at Urbana-Champaign,
Champaign, IL, 2010.
[11] S. Prajna, A. Jadbabaie, and G. J. Pappas, “A framework for worst-case
and stochastic safety verification using barrier certificates,” IEEE Trans.
Autom. Control, vol. 52, no. 8, pp. 1415–1428, 2007.
[12] C. J. Tomlin, I. Mitchell, A. M. Bayen, and M. Oishi, “Computational
techniques for the verification of hybrid systems,” Proc. IEEE, vol. 91,
no. 7, pp. 986–1001, 2003.
[13] J. N. Maidens, S. Kaynama, I. M. Mitchell, M. M. Oishi, and G. A.
Dumont, “Lagrangian methods for approximating the viability kernel in
high-dimensional systems,” Automatica, vol. 49, no. 7, pp. 2017–2029,
2013.
[14] M. Althoff and B. H. Krogh, “Reachability analysis of nonlinear
differential-algebraic systems,” IEEE Trans. Autom. Control, vol. 59,
no. 2, pp. 371–383, 2014.
[15] M. Althoff, “Formal and compositional analysis of power systems using
reachable sets,” IEEE Trans. Power Syst., vol. 29, no. 5, pp. 2270–2280,
2014.
[16] Y. C. Chen and A. D. Domı́nguez-Garcı́a, “A method to study the effect
of renewable resource variability on power system dynamics,” IEEE
Trans. Power Syst., vol. 27, no. 4, pp. 1978–1989, 2012.
[17] Y. Susuki, T. J. Koo, H. Ebina, T. Yamazaki, T. Ochi, T. Uemura, and
T. Hikihara, “A hybrid system approach to the analysis and design of
power grid dynamic performance,” Proc. IEEE, vol. 100, no. 1, pp.
225–239, 2012.
[18] L. Jin, R. Kumar, and N. Elia, “Reachability analysis based transient
stability design in power systems,” Int. J. Elect. Power Energy Syst.,
vol. 32, no. 7, pp. 782–787, 2010.
[19] Y. Susuki and T. Hikihara, “Predicting voltage instability of power
system via hybrid system reachability analysis,” in Proc. Amer. Control,
2007, pp. 4166–4171.
[20] P. A. Parrilo, “Structured semidefinite programs and semialgebraic
geometry methods in robustness and optimization,” Ph.D. dissertation,
Calif. Inst. Technol., Pasadena, CA, 2000.
[21] M. Laurijsse, R. Wisniewski, and S. Weiland, “Optimizing safety supervisors for wind turbines using barrier certificates,” in Proc. 19th IFAC
World Congr., 2014, pp. 11 450–11 457.
[22] I. A. Hiskens, “Dynamics of type-3 wind turbine generator models,”
IEEE Trans. Power Syst., vol. 27, no. 1, pp. 465–474, 2012.
[23] A.-M. Bucurenciu, “Primary frequency response by mtdc offshore
grids,” Master’s thesis, Delft Univ. of Technology, 2014.
[24] Y. Sun, Z. Zhang, G. Li, and J. Lin, “Review on frequency control of
power systems with wind power penetration,” in Proc. 2010 International Conference on Power System Technology (POWERCON), 2010,
pp. 1–8.
[25] J. F. Conroy and R. Watson, “Frequency response capability of full
converter wind turbine generators in comparison to conventional generation,” IEEE Trans. Power Syst., vol. 23, no. 2, pp. 649–656, 2008.
[26] P.-K. Keung, P. Li, H. Banakar, and B. T. Ooi, “Kinetic energy of wind–
turbine generators for system frequency support,” IEEE Trans. Power
Syst., vol. 1, no. 24, pp. 279–287, 2009.
[27] D. Henrion, J. B. Lasserre, and C. Savorgnan, “Approximate volume
and integration for basic semialgebraic sets,” SIAM Rev., vol. 51, pp.
722–743, 2009.
[28] A. Papachristodoulou, J. Anderson, G. Valmorbida, S. Prajna,
P. Seiler, and P. A. Parrilo, SOSTOOLS: Sum of squares
optimization toolbox for MATLAB, http://arxiv.org/abs/1310.4716 ,
2013,
available
from
http://www.eng.ox.ac.uk/control/sostools,
http://www.cds.caltech.edu/sostools
and
http://www.mit.edu/˜parrilo/sostools.
[29] J. Löfberg, “Pre-and post-processing sum-of-squares programs in practice,” IEEE Trans. Autom. Control, vol. 54, no. 5, pp. 1007–1011, 2009.
[30] ——, “YALMIP: A toolbox for modeling and optimization in
MATLAB,” in Proc. IEEE CCA/ISIC/CACSD Conf., Taipei, Taiwan,
2004. [Online]. Available: http://users.isy.liu.se/johanl/yalmip/
[31] M. ApS, The MOSEK optimization toolbox for MATLAB
manual. Version 7.1 (Revision 28)., 2015. [Online]. Available:
http://docs.mosek.com/7.1/toolbox/index.html
[32] D. Henrion and M. Korda, “Convex computation of the region of
attraction of polynomial control systems,” IEEE Trans. Autom. Control,
vol. 59, no. 2, pp. 297–312, 2014.
[33] P. Kundur, N. J. Balu, and M. G. Lauby, Power system stability and
control. McGraw-hill New York, 1994.
[34] X. Yu and K. Tomsovic, “Application of linear matrix inequalities for
load frequency control with communication delays,” IEEE Trans. Power
Syst., vol. 19, no. 3, pp. 1508–1515, 2004.

11

[35] M. S. N. W. Miller and S. Venkataraman, “California iso (caiso)
frequency response study,” General Electric International, Inc., Tech.
Rep., 2011.
[36] V. Singhvi, P. Pourbeik, N. Bhatt, and D. Brooks, “Impact of wind active
power control strategies on frequency response of an interconnection,”
in Proc. IEEE Power Energy Soc. Gen. Meeting (PES GM), 2013, pp.
1–5.
[37] V. Gevorgian, Y. Zhang, and E. Ela, “Investigating the impacts of wind
generation participation in interconnection frequency response,” IEEE
Trans. Sustain. Energy, vol. 6, no. 3, pp. 1–9, 2014.
[38] “Frequency control concerns in the north american electric power
system,” Oak Ridge National Laboratory, Tech. Rep., 2002.
[39] M. A. B. Sassi, S. Sankaranarayanan, X. Chen, and E. Ábrahám, “Linear
relaxations of polynomial positivity for polynomial lyapunov function
synthesis,” IMA J. Math. Contr. Inform., pp. 1–34, 2015.
[40] M. E. Raoufat, K. Tomsovic, and S. M. Djouadi, “Virtual actuators for
wide-area damping control of power systems,” IEEE Trans. Power Syst.,
no. 99, pp. 1–9, 2016.

