1

Complex Frequency
Federico Milano, IEEE Fellow

Abstract— The paper introduces the concept of complex fre- the modeling approaches that go beyond the classical phasor
quency. The imaginary part of the complex frequency is the representation or that focus on analysis of non-sinusoidal
variation with respect of a synchronous reference of the local signals (see, for example, [3] for a state-of-the-art survey on
bus frequency as commonly defined in power system studies.
The real part is defined based on the variation of the voltage this topic and the several references therein). On the contrary,
magnitude. The latter term is crucial for the correct interpre- the proposed concept of “complex frequency” is compatible
tation and analysis of the variation of the frequency at each with the approaches that have been proposed in the literature as
bus of the network. The paper also develops a set of differential
arXiv:2105.07769v2 [eess.SY] 16 Jul 2021

it allows interpreting angle and magnitude variations as com-

equations that describe the link between complex powers and plementary components of the same phenomenon, provided
complex frequencies at network buses in transient conditions. No
simplifications are assumed except for the usual approximations that one accepts to extend the domain of frequency to the
of the models utilized for the transient stability analysis of power complex numbers.
systems. A variety of analytical and numerical examples show the This paper focuses on electro-mechanical transients in high-
applications and potentials of the proposed concept. voltage transmission systems. Thus, the starting point is sim-
Index Terms— Power system dynamics, converter-interfaced ilar to that of [4]–[7], that is, the transient conditions during
generation, frequency control, low-inertia systems. which the magnitude and the phase angle of bus voltage
phasors change according to the inertial response of syn-
chronous machines and the frequency control of synchronous
I. I NTRODUCTION
and non-synchronous devices. On the other hand, harmonics,
A well-known and accepted definition of the frequency unbalanced conditions and electro-magnetic transients are not
of a signal x(t) = Xm (t) cos(ϑ(t)) is given in the IEEE taken into consideration.
Std. IEC/IEEE 60255-118-1 [1], as follows: The resulting formulation is exact, in the measure that power
1 1 system models based on the dqo transform for voltage and
f (t) = ϑ̇(t) = θ̇(t) + fo , (1) angle stability analysis are exact; general, as it provides a
2π 2π
framework to study the dynamic effect of any device on the
where θ is the phase difference, in radians, between the angular
local frequency variations at network buses; and systematic,
position ϑ, also in radians, of the signal x(t) and the phase
because it provides with the tools to determine analytically
due to the reference nominal frequency fo , expressed in Hz.
the impact of each device on bus frequencies.
If the magnitude Xm of the signal is constant, this definition
The remainder of the paper is organized as follows. Section
is adequate. However, if Xm changes with time, the definition
II provides the background for the proposed theoretical frame-
of the frequency in (1) does not provide a meaningful way
work. Section III provides the formal definition of complex
to separate the effects of the variations of ϑ and Xm . Thus,
frequency and its link with complex power injections, voltages
the definition of frequency in the most general conditions is a
and currents and network topology. The special cases of
highly controversial concept that has been discussed at length
constant power and constant current injections as well as con-
in the literature (see the interesting discussion in [2] and the
stant impedances are discussed in Section III-B. Section III-C
references therein).
discusses a variety of relevant approximated expressions that
This paper provides a novel interpretation of “frequency” as
link the complex frequency to bus power injections. Section
complex quantity, i.e., composed of a real and an imaginary
IV illustrates some applications of the analytical expressions
part. This complex frequency takes into account the time
derived in Section III to simulation, state estimation and
dependency of both ϑ and Xm . The focus is on the theory,
control. Finally, Section V draws conclusions and outlines
modeling, simulation and some application aspects of the
future work.
proposed definition. The proposed complex frequency allows
a neat and compact representation as well as a consistent
interpretation of frequency variations in ac power systems. II. BACKGROUND
The complex frequency is also capable of explaining the The starting point is the set of equations that describe the
interactions among active and reactive power injections at complex power injections, in per unit, at the n network buses
buses and flows in network branches. It is important to note of the system, say s̄ ∈ Cn , as follows:
that the proposed approach does not attempt to substitute
s̄(t) = p(t) + q(t) = v̄(t) ◦ ı̄∗ (t) , (2)
F. Milano is with the School of Electrical & Electronic Engineering,
University College Dublin, Belfield, Ireland. E-mail: [email protected] where p ∈ Rn×1 and q ∈ Rn×1 are the active and reactive
This work was supported by Science Foundation Ireland, by funding power injections at network buses, respectively; v̄ ∈ Cn×1 and
F. Milano under project AMPSAS, Grant No. SFI/15/IA/3074; and by the
European Commission by funding F. Milano under project EdgeFLEX, Grant ı̄ ∈ Cn×1 are the voltages and current injections at network
No. 883710. buses; ∗ indicates the conjugate of a complex quantity; and ◦
 2

is the Hadamard product, i.e. the element-by-element product Rphase angles referred to a stationary reference and θo =
of two vectors.1 ω dt is the angle of the rotating dq-axis reference frame and
t o
In steady-state, balanced conditions, (2) expresses the well- ωo is the angular frequency in rad/s of the dq-axis reference
known power flow equations. However, it is important to note frame.
that, in (2), all quantities are assumed to be time dependent. From (1), the time derivative of θ gives:
The elements of the voltage and current vectors v̄ and ı̄ that
appears in (2), in fact, are not to be interpreted as conventional ω(t) = θ̇(t) = ϑ̇(t) − ωo (t) , (8)
phasors, but as dynamic quantities, which in some references where ω is the vector of frequency deviations with respect
are called Park’s vectors [8], [9]. A Park’s vector is a complex to the reference frequency at the network buses. In [1], it is
quantity obtained from the dq-axis components of the well- assumed that ωo = 2πfo is constant and equal to the nominal
known dqo transform. For example, for the voltage, one has: angular frequency of the grid, e.g., ωo = 2π 60 rad/s in North
v̄(t) = v d (t) + v q (t) . (3) American transmission grids. Note that ωo being constant is
not a requirement of the transform in (5). However, to carry out
where the components vd,k and vq,k of the k-th element of the derivations presented in Section III, ωo is assumed to be
the vector v̄ are calculated as follows: constant when it is utilized to calculate the values of reactances
   
vd,k (t) va,k (t) and susceptances. As for the reactive power, this is again a
vq,k (t) = P(t) vb,k (t) , (4) widely-accepted approximation utilized in RMS models for
vo,k (t) vc,k (t) angle and voltage stability analysis and consists in assuming
that the link between current injections and voltages is given
where
by
cos(θo0 (t)) cos(θo00 (t))
 
q cos(θo (t)) ı̄(t) ≈ Ȳ v̄(t) , (9)
P(t) = 2  sin(θ (t))
3 o sin(θo0 (t)) 00
sin(θo (t))  , (5)
√1 √1 √1 where Ȳ = G + B ∈ Cn×n is the well-known admit-
2 2 2
tance matrix of the network. It is important not to confuse
and θo is the angle between the phase a and the q-axis, with (9) with the conventional relationship between current and
θ˙o = ωo , and θo0 = θo − 2π 00 2π
3 and θo = θo + 3 . The same voltage phasors (in which case (9) is an exact equality). ı̄
transformation (4) is applied to the abc currents. Since no and v̄ are Park’s vectors, i.e., complex quantities with time-
assumption is made on the abc quantities, the d- and q- varying real and imaginary parts and, hence, (9) represents
axis components of the Park’s vectors v̄ and ı̄ and, hence, an approximation of the dynamics of the grid. In turn, to
(2) can be assumed to be valid in transient conditions. It is obtain (9), it is assumed that, for network inductances and
important to note that the reactive power is not well defined for capacitances the relationships between voltages and currents
non-sinusoidal signals.2 However, it is a common assumption, can be approximated with:
which effectively underpins the vast majority of studies on the
transient stability of power systems [11], to approximate the d
v̄ = Lı̄˙ = L( + ωo )ı̄ ≈ ωo Lı̄ = Xı̄ ,
reactive power as in (2). dt (10)
The vo,k is the o-axis or zero component and is null for d
ı̄ = C v̄˙ = C( + ωo )ı̄ ≈ ωo C v̄ = Bv̄ ,
balanced systems. If the system is not balanced and the o- dt
axis components are not null, then the vector p in (2) does where dt d
is the time derivative relative to the Park rotating
not represent the total active power injections at network frame; ωo is the term due to the rotation of the Park reference;
buses as it does not include the term v o ◦ ıo . The hypothesis and L, C, X, B are the inductance, capacitance, reactance and
of balanced system is not necessary for the developments susceptance, respectively. The quantities in (10) are assumed in
presented below. However, since the focus is on high-voltage absolute values. In turn, the approximation above assumes that
transmission systems, in the remainder of this paper, balanced, electro-magnetic transients in the elements of the transmission
positive sequence operating conditions are assumed. lines and transformers are fast and can be assumed to be in
For the purposes of the developments given below, it is Quasi-Steady-State (QSS). The approximation (10) is applied
convenient to rewrite (3) in polar form: also to the equations of the circuits of the devices connected to
v̄(t) = v(t) ◦ ∠θ(t) , (6) the grid, e.g., the equations of the synchronous machine (see
(71) in Section IV-B). The focus of this paper is, in fact, on the
where v = |v̄|, ∠θ = cos(θ) +  sin(θ) and time scales of electro-mechanical and primary frequency and
θ(t) = ϑ(t) − θo (t) , (7) voltage control transients, which are a few orders of magnitude
slower than electro-magnetic dynamics.
namely, θ is the vector of bus voltage phase angles referred to Merging (2) and (9) becomes:
the rotating dq-axis reference frame, ϑ are the bus voltage
s̄(t) = v̄(t) ◦ [Ȳ v̄(t)]∗ . (11)
1 The Hadarmard product of two column vectors x and z can be also written
as x ◦ z = diag(x) z, where diag(x) is a diagonal matrix whose element These equations resemble the well-known power flow equa-
(i, i) is the i-th element of the vector x. tions except for the fact that the voltages are Park’s vectors,
2 Several attempts have been done to try to overcome this issue. Among
these works, we mention the monograph on the “instantaneous power theory” not phasors, and, hence, the power injections at buses are, in
[10]. general, time-varying quantities.
 3

A. Time Derivative of Algebraic Equations III. D ERIVATION

An important aspect of the developments discussed in the For the sake of the derivation, it is convenient to drop
next section is whether (11) can be differentiated with respect the dependency on time and rewrite (11) in an element-wise
to an independent variable and, in particular, with respect to notation. For a network with n buses, one has:
time. With this aim, observe that (9) leads to the well-known Pn
ph = vh k=1 vk [Ghk cos θhk + Bhk sin θhk ] ,
QSS model for power system angle and voltage transient Pn (16)
stability analysis, as follows [11], [12]: qh = vh k=1 vk [Ghk sin θhk − Bhk cos θhk ] ,
ẋ = f (x, y) , where Ghk and Bhk are the real and imaginary parts of the
(12) element (h, k) of the network admittance matrix, i.e. Ȳhk =
0 = g(x, y) ,
Ghk + Bhk ; vh and vk denote the voltage magnitudes at
where f ∈ Rnx +ny 7→ Rnx are the differential equations; g ∈ buses h and k, respectively; and θhk = θh − θk , where
Rnx +ny 7→ Rny are the algebraic equations, x ∈ X ⊂ Rnx θh and θk are the voltage phase angles at buses h and k,
are the state variables; and y ∈ Y ⊂ Rny are the algebraic respectively. Equations (16) and all equations with subindex h
variables. Equations (12) are smooth in Rnx +ny except for a in the remainder of this section are valid for h = 1, 2, . . . , n.
finite number of points at which discrete events occur, such Equations (16) can be equivalently written as:
as line outages and faults. In practice, these events can be Pn Pn
modeled as if-then rules that modify the structure of f and g. ph = k=1 phk , and qh = k=1 qhk , (17)
There exist, of course, more sophisticated approaches to model
where
events. These often involve the definition of additional vari-
ables and equations such as the hybrid automaton described phk = vh vk [Ghk cos θhk + Bhk sin θhk ] ,
(18)
in [13], or the approach based on Filippov theory discussed in qhk = vh vk [Ghk sin θhk − Bhk cos θhk ] .
[14]. For simplicity, however, but without lack of generality,
model (12) is considered in the remainder of this work. Differentiating (16) and writing the active power injections
The implicit function theorem indicates that, if the Jacobian as the sum of two components:
matrix ∂g/∂y is not singular, there exists a function φ such Pn ∂ph Pn ∂ph
that: dph = k=1 dθhk + k=1 dvk ≡ dp0h +dp00h , (19)
∂θhk ∂vk
y = φ(x) . (13)
In (19), dph is the total variation of power at bus h; dp0h is the
Equation (13) is often utilized to reduce the set of Differential- quota of the active power that depends on bus voltage phase
Algebraic Equations (DAEs) in (12) into a set of Ordinary angle variations; and dp00h is the quota of active power that
Differential Equations (ODEs) that depends only on x and z. depends on bus voltage magnitude variations. From (11) and
In this work, however, (13) is utilized the other way round, using the same procedure that leads to (19), one can define
i.e., to guarantee that it is possible to define the time derivative two components also for the reactive power, as follows:
of y except for the finite number of points at which the events Pn ∂qh Pn ∂qh
occur. This condition leads to: dqh = k=1 dθhk + k=1 dvk ≡ dqh0 + dqh00 . (20)
 −1 ∂θhk ∂vk
∂φ ∂g ∂g
ẏ = ẋ = ẋ From (17), it is relevant to observe that:
∂x ∂y ∂x
 −1 (14) ∂ph ∂qh
∂g ∂g = −qhk , and = phk , (21)
= f (x, φ(x)) . ∂θhk ∂θhk
∂y ∂x
The condition (14) implies that the set of DAEs in (12) is which leads to rewrite dp0h and dqh0 as:
assumed to be index 1 [15], which is the form of DAEs Pn
dp0h = − k=1 qhk dθhk ,
that describes most physical systems, including power systems n (22)
dqh0 = k=1 phk dθhk ,
P
[12].
The voltage magnitudes v and phase angles θ that appears Recalling that θhk = θh − θk , one has:
in (11), and hence also the real and imaginary parts of the Pn
dp0h = −qh dθh + k=1 qhk dθk ,
complex power s̄, are algebraic variables in the conventional Pn (23)
formulation of QSS models. Thus, the assumption of index- dqh0 = ph dθh − k=1 phk dθk ,
1 DAEs allows rewriting the current injections at bus and, where the identities (17) have been used.
hence, the complex power s̄ as functions of state and discrete In the same vein, from (17), (19) and (20), dp00h and dqh00
variables, as well as of the bus voltages v̄, namely s̄(v̄, x), can be rewritten as:
which are smooth, except at the points at which the events ph Pn phk
occur. Then, the time derivatives of s̄ can be computed with dp00h = dvh + k=1 dvk ,
vh vk
the chain rule as: qh Pn qhk (24)
∂s̄ ∂s̄ dqh00 = dvh + k=1 dvk .
s̄˙ = v̄˙ + ẋ . (15) vh vk
∂ v̄ ∂x Let us define the quantity:
The next section of this paper elaborates on (15) and deduces
an expression that involves the concept of complex frequency. uh ≡ ln (vh ) , (25)
 4

where both the logarithm and its argument are dimensionless. the elements of s̄ are the inputs or boundary conditions at
With this aim, vh is in per unit, namely it is the ratio of two network buses and depend on the devices connected to grid,
voltages (thus without units). The differential of (25) gives: whereas S̄ depends only on network quantities.
dvh Another way to write (32) is by splitting s̄˙ into its compo-
duh = . (26) nents s̄˙ 0 and s̄˙ 00 . According the definitions of s̄0 and s̄00 , s̄˙ 0
vh
does not depend on %, whereas s̄˙ 00 does not depend on ω, as
Then, (24) can be rewritten as: follows:
dp00h = ph duh +
Pn
phk duk , s̄˙ 0 = s̄ ◦ ω − S̄ ω ,
Pnk=1 (27) (35)
dqh00 = qh duh + k=1 qhk duk . s̄˙ 00 = s̄ ◦ % + S̄ % .

Equations (23) and (27) can be expressed in terms of complex It is important to note that, in general, the expressions of s̄˙ 0
powers variations: and s̄˙ 00 are not known a priori. These components, however,
Pn can be determined using (32), (35) and:
ds̄0h = dp0h + dqh0 = s̄h dθh −  k=1 s̄hk dθk ,
n (28) s̄˙ = s̄˙ 0 + s̄˙ 00 . (36)
ds̄00h = dp00h + dqh00 = s̄h duh + k=1 s̄hk duk ,
P

and, finally, defining the complex quantity: With this regard, Section IV-B explains through an example
how to calculate ω and % based on (32). In the following, (32),
ζ̄h ≡ uh +  θh , (29) (35) and (36) are utilized to discuss relevant special cases.
the total complex power variation is given by:
ds̄h = ds̄0h + ds̄00h A. Alternative Expressions
= dp0h + dp00h + (dqh0 + dqh00 ) (30) We now derive (32) in an alternative and more compact
Pn formulation as a function of the currents. First, observe that:
= s̄h dζ̄h + k=1 s̄hk dζ̄k∗ .
Equation (30) can be written in a compact matrix form, as v̄˙ = v̄ ◦ η̄ . (37)
follows: The proof of (37) is given in Appendix I. Then, recalling (2),
∗
ds̄ = s̄ ◦ dζ̄ + S̄ dζ̄ , (31) the time derivative of s̄ with respect to the dq-axis reference
where S̄ ∈ Cn×n is a matrix whose (h, k)-th element is s̄hk . frame can be written as:
The expression (31) has been obtained in general, i.e., as- d
s̄˙ = (v̄ ◦ ı̄∗ )
suming a differentiation with respect to a generic independent dt
parameter. If this parameter is the time t, (31) can be rewritten = v̄˙ ◦ ı̄∗ + v̄ ◦ ı̄˙ ∗ (38)
as a set of differential equations: = v̄ ◦ η̄ ◦ ı̄∗ + v̄ ◦ ı̄˙ ∗
s̄˙ − s̄ ◦ η̄ = S̄ η̄ ∗ (32) = s̄ ◦ η̄ + v̄ ◦ ı̄˙ ∗ ,

where where the commutative property of the Hadamard product has

η̄ = ζ̄˙ = u̇ + θ̇ , (33) been used. Substituting (38) into (32), one obtains:

and recalling the derivative of θ given in (8) and defining v̄ ◦ ı̄˙ ∗ = S̄ η̄ ∗ (39)
% ≡ u̇, one obtains:
Yet another way to write (32) can be obtained by differenti-
η̄ ≡ % +  ω (34) ating with respect of time (9), or, which is the same, dividing
each row h of S̄ by the corresponding voltage v̄h in (39).
We define η̄ as the vector of complex frequencies of the buses Implementing this division and using (37) lead to (see also
of an AC grid. Note that both real and imaginary part of (34) footnote 1):
have, in fact, the dimension of s−1 , as u is dimensionless
and ω is expressed in rad/s. In (34), the imaginary part is ı̄˙ = Ȳ [v̄ ◦ η̄] = Ȳ diag(v̄) η̄ , (40)
the usual angular frequency (relative to the reference ωo ). On
and, defining Ī = Ȳ diag(v̄), one obtains:
the other hand, it is more involved to determine the physical
meaning of the real part of (34), %. From the definition of uh , ı̄˙ = Ī η̄ (41)
one has vh = exp(uh ), that is, the magnitude of the voltage
is expressed as a function whose derivative is equal to the As per (32), the right-hand sides of (39) and (41) depend
function itself. This concept is key in the theory of Lie groups exclusively on network quantities, whereas the left-hand side
and algebra, which defines the space of linear transformations is device dependent. While equivalent, the relevant feature of
of generalized “rates of change” [16]. (39) and (41) with respect to (32) is that the complex frequency
Equation (32) is the sought expression of the relationship vector only appears once.
between frequency variations and power flows in an ac grid. It is worth noticing that one can also proceed the other
It contains the information on how power injections of the way round, namely obtain (32) from (41). In fact, taking the
devices connected to the grid impact on the frequency at their conjugate of (41) and multiplying both sides by (v̄ ◦) and
point of connection as well as on the rest of the grid. In (32), using the identity (38), one re-obtains (32), as expected.
 5

Finally, we note that (41) requires less calculations than and, defining Ȳh,tot = Ȳho + Ȳhh , one obtains:
(32). Thus, in a software tool where currents are modelled ∗ Pn
−Ȳh,tot vh2 η̄h∗ = h6=k s̄hk η̄k∗ , (53)
as state variables and, hence, their first time derivatives are
available as a byproduct of the integration of (12), (41) can which indicates that the two components of the complex
be an efficient alternative to (32) for the calculation of η̄. frequency, %h and ωh at the bus h of a constant admittance load
are linear combinations of the %k and ωk at the neighboring
B. Special Cases buses. This, in turn and as expected, means that constant
Three cases are considered in this section, namely, con- admittance loads are passive devices and cannot modify the
stant power injection, constant current injection, and constant frequency at their point of connection but rather “take” the
admittance load. These cases illustrate the utilization of the frequency that is imposed by the rest of the grid. This
formulas deduced above, namely (32), (39) and (41). conclusion generalizes the results of the appendix of [7] that
1) Constant Power Injection: We illustrate first an appli- considers the simplified case of an admittance load connected
cation of (32) for a constant power injection, say s̄h = s̄ho . to the rest of the system through a lossless line.
From (36), the boundary condition at the h-th bus is: 3) Constant Current Injection: This last example shows an
application of (41). For a constant current injection, we have
s̄˙ h = 0 ⇒ s̄˙ 0h = −s̄˙ 00h , (42) two cases. If the magnitude and phase angle of the current are
and, from (35): constant, then ı̄˙h = 0 and, from (41), one obtains immediately:
Pn
Pn Pn
s̄ho ωh −  k=1 s̄hk ωk = −s̄ho %h − k=1 s̄hk %k , (43) 0 = k=1 ı̄hk η̄k , (54)

and, from (32): where, ı̄hk is the (h, k) element of Ī. Equivalently, from (39),
Pn the condition ı̄˙h = 0 leads to:
−s̄ho η̄h = k=1 s̄hk η̄k∗ . (44) Pn
0 = k=1 s̄hk η̄k∗ . (55)
2) Constant Admittance Load: This case illustrates an
application of (39). For a constant admittance load, the current On the other hand, it is unlikely that a device is able to impose
consumption at the h-th bus is: the phase angle of its current injection independently from the
phase angle of its bus voltage. More likely, a device imposes
ı̄h = −Ȳho v̄h , (45) a constant magnitude and power factor. In this case, the phase
where the negative sign indicates that the current is drawn angle of the current will depend on the phase angle θh of the
from bus h. From (39) and (45) and recalling that v̄˙ h = v̄h η̄h voltage at bus h (see Appendix II). This means that, according
(see Appendix I), one obtains: to (28), the term ds̄00h = 0 and, from (36), s̄˙ 00h = 0 and s̄˙ h = s̄˙ 0h .
Pn This result generalizes the one obtained in [17]. From (41),
∗ 2 ∗
−Ȳho vh η̄h = k=1 s̄hk η̄k∗ . (46) a current injection with constant magnitude and power factor
leads to:
The same result can be obtained also from (32). The power Pn
ı̄˙h =  k=1 ı̄hk ωk ,
consumption at bus h can be written as: Pn (56)
0 = k=1 ı̄hk %k ,
s̄h = v̄h ı̄∗h = −Ȳho
∗ 2
vh , (47)
and, since the angular frequency of the current of a constant
which indicates that the power consumption s̄h in (47) does power factor load is the same of the angular frequency ωh of
not depend on θh . This means that, according to (28), the term the voltage at the load bus (see Appendix II and (85)), the first
ds̄0h = 0 and, from (36), s̄˙ h = s̄˙ 00h and s̄˙ 0h = 0. From the first equation of (56) can be rewritten as:
equation of (35), one has: ı̇h

Pn ı̄hk

∗ 2 Pn Re{ξ¯h } = = − ωh − k=1 ωk . (57)
−Ȳho vh ωh = s̄h ωh = k=1 s̄hk ωk . (48) ıh ı̄h
Then, from the time derivative of (47) and the second equation
C. Approximated Expressions
of (35), one has:
∗ ∗ 2 Pn The mathematical developments carried out so far have
−2Ȳho vh v̇h = −Ȳho vh %h + k=1 s̄hk %k . (49) assumed no simplifications except for neglecting the electro-
Observing that, from (26), %h = v̇h /vh , then: magnetic dynamics of network branches. All formulas that
have been deduced are thus accurate in the measure that
vh v̇h = vh2 %h , (50) the effect of electro-magnetic transients are negligible. It is,
and, hence, (49) can be rewritten as: however, relevant to explore whether the expressions (32), (39)
Pn and (41) can be approximated while retaining the information
∗ 2
−Ȳho vh %h = s̄h %h = k=1 s̄hk %k , (51) on the relationship between power injections and frequency
variations at network buses.
The expressions (48) and (51) have the same structure and
Except during faults and some post-fault transients, it is not
can be merged into (46). Moreover, in the summations on the
uncommon the case for which one can assume that vh ≈ 1
right-hand-sides of (46), the term s̄hh is given by:
pu and that bus voltage phase angle differences are small,
∗ 2
s̄hh = Ȳhh vh , (52) hence sin(θh − θk ) ≈ θh − θk and cos(θh − θk ) ≈ 1. These
 6

assumptions, which, in turn, are the well-known approximation A. Ratio between ω and %
utilized in the fast decoupled power flow method [18], lead to:
This first example illustrates the transient behavior of the
∗
s̄hk ≈ Ȳhk . (58) two components of the complex frequency. Figure 1 shows
Equation (58) allows simplifying (32) as: the transient behavior of the components of the complex
∗
frequency at a generator and a load bus, bus 2 and bus 8,
s̄˙ − s̄ ◦ η̄ ≈ Ȳ η̄ ∗ , (59) respectively, for the well-known WSCC 9-bus system [11].
and (39) and (41) as: The estimation of ωh and %h at the buses of the grid is obtained
through a synchronous-reference frame Phase-Locked Loop
ı̄˙ ≈ Ȳ η̄ . (60) (PLL) model [20]. Simulation results show that |ωh |  |%h |.
One can further simplify the expressions above for high- This inequality holds for all networks and scenarios that we
voltage transmission systems, for which Ȳ ≈ B: have tested for the preparation of this work.
It is important to note that the inequality |ωh |  |%h | does
s̄˙ − s̄ ◦ η̄ ≈ −B η̄ ∗ , (61) not imply s̄˙ 0h  s̄˙ 00h . As a matter of fact, taking as an example
and (39) and (41) as: constant impedance loads, s̄˙ 0h = 0 and s̄˙ 00h = s̄˙ h . The rationale
behind this observation can be explained by rewriting (35)
ı̄˙ ≈ B η̄ . (62) using (17) and an element-by-element notation:
∗
Then, approximating the term s̄ ◦ η̄ ≈ Ȳdiag η̄,
where Ȳdiag Pn
is a matrix obtained using the diagonal elements of Ȳ, and s̄˙ 0h =  k=1 s̄hk (ωh − ωk ) ,
Pn (67)
splitting the real and imaginary part of η̄, (61) leads to: s̄˙ 00h = k=1 s̄hk (%h + %k ) ,
ṗ0 ≈ B0 ω , (63) which indicates that s̄˙ 0h is proportional to the difference of the
elements of ω, whereas s̄˙ 00h is proportional to the sum of the
0 0
Pn
where Bhk = −Bhk and Bhh = h6=k Bhk are the elements
of B0 , and elements of %.
q̇ 00 ≈ B00 % , (64)
00 00
where Bhk = −Bhk and Bhh = −2Bhh are the elements
00 B. Implementation of Equation (32)
of B . Equation (63) is the expression deduced in [7] and
that, with due simplifications, leads to the frequency divider The implementation of (32) in a software tool for the
formulas presented in [5]. simulation of power systems can be useful to determine the
Considering the resistive parts of the network branches, the “exact” frequency variations at network buses in a RMS model
following dual expressions hold: for transient stability analysis. This topic has been discussed
and solved under various hypotheses in [4]–[6]. In particular,
ṗ00 ≈ G00 % , q̇ 0 ≈ G0 ω , (65)
the Frequency Divider Formula (FDF) proposed in [5] is based
0
where the elements of GP and G00 are defined as G0hk = on (63), which is an approximation of (32).
00 0 n 00
Ghk = −Ghk , Ghh = h6=k Ghk , and Ghh = −2Ghh . The expression (32) states the link between the complex
Interestingly, from (65), it descends that, in lossy networks, frequency and the rest of the variables of the system. Thus,
the reactive power can be utilized to regulate the frequency. in (32), the unknowns are η̄ or, equivalently, ρ and ω. One
Combining together (63), (64) and (65) leads to the follow- can, of course, calculate these quantities by differentiating with
ing approximated expressions: respect to time u and θ, respectively. However, since u and
0∗
s̄˙ 0 ≈ Ȳ ω , s̄˙ 00
00
≈ Ȳ % . (66) θ are algebraic variables, (12) does not provide directly the
values of their time derivatives. So either one has to use some
IV. E XAMPLES sort of numerical differentiation (with the issues that arise at
discontinuities); use some sort of estimation (e.g., the PLL
The examples presented below apply the theory developed
utilized in the previous example); or solve (32). The latter is
in the previous section and discuss applications to power sys-
the approach utilized in this example.
tem modeling, state estimation and control. In particular, three
devices are discussed, namely, the synchronous machine; the The procedure implemented in the simulations presented in
voltage dependent load; and a converter-interfaced generator this section is as follows.
with frequency and voltage control capability. The objective • The system variables x and y are initialized using (12).
is to illustrate the methodological approach discussed above • The DAEs defined by (12) are integrated using a con-
and showcase some problems that the proposed definition of ventional numerical scheme (in this case, the implicit
complex frequency and the expression (32) make possible to trapezoidal method with fixed time step).
solve. In the following, all simulation results are obtained • Discrete events are also handled with a conventional
using the software tool Dome [19]. With this regard, note that, approach. Specifically, an approach based on if-then rules
in this software tool, ωo is set to be equal to the frequency of that switch the equations on the right-hand-side of (12)
the Center of Inertia (CoI), nevertheless, as usually assumed has been utilized.
in transient stability analysis, the values of reactances and • Finally, (32) – or, alternatively, (39) or (41) – and the
susceptances are kept constant. values x, y and ẋ obtained at each step of the time
 7

0.4
respect to η̄, so it remains to determine the dependency of the
0.3
elements of s̄˙ − s̄ ◦ η̄ ∗ – or, alternatively, of v̄ ◦ ı̄˙ ∗ from (39)
0.2
or ı̄˙ from (41) – on η̄ and, eventually, on the time derivatives
0.1
ωh [Hz]

of the state variables ẋ.

0 We illustrate the procedure using a conventional 4th order
−0.1 model of the synchronous machine and (32). The stator voltage
−0.2 Bus 2 of the machine with respect to its dq-axis reference frame is
−0.3 Bus 8 linked to the grid voltage vh with the following equations [11]:
−0.4
0 2.5 5 7.5 10 12.5 15 17.5 20 v̄s = vs,d +  vs,q = v̄h ∠( π2 − δr ) , (68)
Time [s]
0.4 where δr is th rotor angle of the machine. The time derivative
0.3 of (68) gives (see Appendix I for the derivative of v̄h ):
0.2
0.1 v̄˙ s = v̄s (η̄h −  ωr ) , (69)
%h [Hz]

0 where ωr is the deviation in rad/s of the angular speed of the

−0.1 machine with respect to ωo . The real and imaginary terms of
−0.2 Bus 2 (69) can be written as:
−0.3 Bus 8
v̇s,d = vs,d %h − vs,q ωh + vs,q ωr ,
−0.4 (70)
0 2.5 5 7.5 10 12.5 15 17.5 20 v̇s,q = vs,q %h + vs,d ωh − vs,d ωr .
Time [s]
Then, the stator electrical and magnetic equations are:
(a) Outage of the load at bus 5.
Xd0 ıs,d + Ra ıs,q = −vs,q + e0r,q
2 (71)
Ra ıs,d − Xq0 ıs,q = −vs,d + e0r,d
1.5
1 where ı̄s = ıs,d +  ıs,q is the stator current. From (71), one
0.5 can thus obtain the expressions of ıs,d and ıs,q as a function
ωh [Hz]

0 of vs,d , vs,q and the state variables e0r,d and e0r,q . In this case,
−0.5 these relationships are linear, so, one can obtain the current
−1 Bus 2 explicitly as:
−1.5 Bus 8
−1 
−vs,q + e0r,q
   0 
ıs,d Xd Ra
−2 = , (72)
0 2.5 5 7.5 10 12.5 15 17.5 20
ıs,q Ra −Xq0 −vs,d + e0r,d
Time [s] or, equivalently:
2
1.5
ıs,d = kv,dd vs,d + kv,dq vs,q + ke,dd e0r,d + ke,dq e0r,q ,
(73)
1 ıs,q = kv,qd vs,d + kv,qq vs,q + ke,qd e0r,d + ke,qq e0r,q ,
0.5 where the k-coefficients are constant and depend on Ra , Xd0
%h [Hz]

0 and Xq0 . Differentiating (73) with respect to time, substituting

−0.5 the expressions of v̇s,d and v̇s,d from (70) and applying the
−1 Bus 2 chain rule, the time derivatives of the components of the
−1.5 Bus 8 current can be written in the form:
−2 ∂ıs,d ∂ıs,d ∂ıs,d ∂ıs,d 0 ∂ıs,d 0
0 2.5 5 7.5 10 12.5 15 17.5 20 ı̇s,d = ∂%h %h + ∂ωh ωh + ∂ωr ωr + ∂e0r,d ėr,d + ∂e0r,q ėr,q ,
Time [s]
∂ıs,q ∂ıs,q ∂ıs,q ∂ıs,q 0 ∂ıs,q 0
(b) Fault at bus 7 cleared by tripping line 5-7. ı̇s,q = ∂%h %h + ∂ωh ωh + ∂ωr ωr + ∂e0r,d ėr,d
, + ∂e0r,q ėr,q
(74)
Fig. 1: Transient behavior of the components of the complex fre-
quency for the WSCC 9-bus system. where, in the right-hand side, %h and ωh are unknown and all
other variables and terms are calculated based on the current
solution of (12). Then, from (39), one obtains:
domain integration are utilized to calculate η̄.3 Re{v̄h ı̄˙∗h } = Re{v̄s ı̄˙∗s } = vs,d ı̇s,d + vs,q ı̇s,q ,
With regard to the last step, i.e., the determination of η̄, note (75)
Im{v̄h ı̄˙∗h } = Im{v̄s ı̄˙∗s } = vs,q ı̇s,d − vs,d ı̇s,q .
that (32) can be rewritten as a set of 2n real equations, with
unknowns (%, ω). The right-hand side of (32) is linear with Finally, substituting (74) into (75), one obtains the expressions
of the left-hand side of (39) at the buses of the synchronous
3 The procedure described above is an open-loop, i.e., (32) are solved “off- machines.
line” based on the solution obtained by integrating (12). In a closed loop, The very same procedure can be applied to any device of
i.e., if the elements of the complex frequency are utilized as inputs to some
controllers or one wants to take into account the dependence on the frequency the grid and, in the vast majority of the cases, the resulting
in the model of a device, then (12) has to be extended to include (32). expressions are linear with respect to (%, ω). Hence, at each
 8

0.015
time step of a time domain simulation one need to solve a
CF
problem of the type Aχ = b, where χ = [%T , ω T ]T , and: 0.01
PLL
0.005

ω bus 1 [Hz]
 
[H + diag(p) − P% ] [K − diag(q) − Pω ] FDF
A= (76)
[K + diag(q) − Q% ] −[H − diag(p) + Qω ] 0

and   −0.005
Px ẋ −0.01
b= , (77)
Qx ẋ
−0.015
where H = Re{S̄} and K = Im{S̄}; P% , Pω , Q% and Qω 0 0.5 1 1.5 2 2.5 3
Time [s]
are the Jacobian n × n matrices of p and q with respect to %
(a) 4th order machine models
and ω, respectively; and Px and Qx are the Jacobian n × nx
matrices of p and q with respect to the state variables x, 0.015
respectively. CF
0.01
Figure 2 shows the imaginary part of the complex frequency PLL
0.005

ω bus 3 [Hz]
as obtained for the WSCC 9-bus system using a 4th and a
2nd order model of the synchronous machines, as well as 0
for a scenario where the machine at bus 3 is substituted −0.005
for a Converter-Interfaced Generator (CIG). The model of
the CIG is shown in Fig. 3. The results obtained using −0.01
PLLs matches well the “exact” results obtained with proposed −0.015
method (indicated with CF in the legends) except for the 0 0.5 1 1.5 2 2.5 3
Time [s]
numerical spikes that following the load disconnection and
(b) 4th order machine models & noise
a small delay that is due to the control loop of the PLL. The
proposed formula is also robust with respect to noise as shown 0.015
in Fig. 2b. CF
0.01
The events create spikes in the estimation of ω obtained with PLL
0.005
ω bus 1 [Hz]

the PLL. The behavior of ω as estimated with the complex FDF

frequency, on the other hand, depends on the dynamics of the 0
device connected to the bus where the frequency is estimated. −0.005
It is smooth for the 2nd order synchronous machine and show a
−0.01
small jump for the 4th order synchronous machine. The latter
behavior is due to the dependence of the frequency on the −0.015
0 0.5 1 1.5 2 2.5 3
voltage magnitude, which is modeled as an algebraic variable
Time [s]
in (12). On the other hand, the behavior of the ω estimated
(c) 2nd order machine models
with the complex frequency at the CIG bus shows a relatively
fast transient after the load disconnection (see Fig. 2d). This 0.06
is due to the dynamics of the currents of the converter. CF
0.04
Figure 2 also shows the results obtained using the FDF PLL
0.02
ω bus 3 [Hz]

proposed in [5]. The FDF matches closely the results of FDF

the proposed method when considering the simplified 2nd 0
order model of the machines but introduces some errors if −0.02
the machine model includes rotor flux dynamics or when
−0.04
considering the CIG.
−0.06
0 0.5 1 1.5 2 2.5 3
C. Voltage Dependent Load (VDL) Time [s]
(d) CIG connected at bus 3
The power consumption of a VDL is:
Fig. 2: Comparison of bus frequencies using PLL, FDF and the
γ γ
s̄h = ph + jqh = −po vhp −  qo vhq , (78) proposed approach based on the complex frequency. The plots refers
to the WSCC 9-bus system following the disconnection at t = 1 s
then: of the load at bus 5.
γ −1 γ −1
s̄˙ h = −po γp v̇h vhp −  qo γq v̇h vhq
(79)
= (γp ph +  γq qh ) %h ,
As an application, we utilize (79) to estimate the parameters
where it has been assumed that the exponents γp and γq are γp and γq of a VDL using a similar technique as the one
constant and that po and qo vary “slowly” with respect to vh . proposed in [17]. From (32) and (79), one obtains:
If γp = γq = γ, then s̄˙ h = γ s̄h %h , which generalizes the Pn
results obtained in Section III-B.2. (γp ph +  γq qh ) %h = k=1 [s̄hk (η̄h + η̄k∗ )] , (80)
 9

2.6 2.6
s̄ref
h Current ı̄ref
h 1 ı̄h 2.4 2.4 Paper [8]
set point 2.2 2.2 Equation (79)
controller 1 + sT 2 2

γ̂p

γ̂q
1.8 1.8
Converter 1.6 1.6
v̄h Paper [8]
1.4 1.4
1.2 Equation (79) 1.2
ω ref + Kp pref
h
vhref + Kq qhref 1 1
2 6 10 14 2 6 10 14
− 1 + sTp − 1 + sTq Time [s] Time [s]
ωh vh (a) Ratio of transmission lines: R/X  1.

Fig. 3: Simplified scheme of a CIG and its controllers. 2.6 2.6

2.4 2.4 Paper [8]
2.2 2.2 Equation (79)
2 2
and, splitting real and imaginary parts:

γ̂p

γ̂q
1.8 1.8
Pn 1.6 1.6
γp = (ph %h )−1 Re { k=1 s̄hk (η̄h + η̄k∗ )} , 1.4
Paper [8]
1.4
Pn (81) Equation (79)
γq = (qh %h )−1 Im { k=1 s̄hk (η̄h + η̄k∗ )} , 1.2 1.2
1 1
2 6 10 14 2 6 10 14
where the right-hand sides can be determined based on mea- Time [s] Time [s]
surements. The fact that %h → 0 in steady-state can create (b) Ratio of transmission lines: R/X ≈ 1.
numerical issues, which can be solved, as discussed in [17],
using finite differences over a period of time ∆t, namely Fig. 4: Estimation of the exponents of the VDL connected at bus 8
following the disconnection of 15% of the load at bus 5 at t = 1 s
η̄h ≈ ∆ζ̄h /∆t, η̄k∗ ≈ ∆ζ̄k∗ /∆t, and %h ≈ ∆uh /∆t, as follows: for the WSCC 9-bus system.
Pn
γ̂p ≈ (ph ∆uh )−1 Re ∗
k=1 s̄hk (∆ζ̄h + ∆ζ̄k ) ,
n (82)
γ̂q ≈ (qh ∆uh )−1 Im ∗
 P
k=1 s̄hk (∆ζ̄h + ∆ζ̄k ) .
The simplified converter-interfaced generator model shown
Equations (81) and (82) generalize the empirical formulas to in Fig. 5, with differential equations:
estimate γp and γq proposed in [17]. The latter, in fact, can
be obtained from (82) by approximating s̄hk ≈ −jBhk . Tp ṗh = Kp (ω ref − ωh ) − ph ,
(83)
Figure 4 shows the results obtained for the WSCC 9-bus Tq q̇h = Kq (v ref − vh ) − qh .
system where the bus connected at bus 8 is a VDL with γp = 2
and γq = 1.5. The results show that (82) is, as expected,
more precise than the approximated estimations based on the ω ref + Kp pref
h
ω ref + Kq qhref
expression proposed in [17]. Equation (82) is, in fact, an − 1 + sTp − 1 + sTq
exact expression and its accuracy depends exclusively on the ωh ωh
accuracy of the measurements of ζ̄h and ζ̄k , which can be
Fig. 5: Alternative frequency control through active and reactive
obtained, for example, with PMUs. If the R/X ratio of the powers for the CIG of Fig. 3.
transmission lines of the system is changed to resemble that
of a distribution system, (82) appears also numerically more
robust than its approximated counterparts (see Fig. 4b).
V. C ONCLUSIONS
The paper introduces a new physical quantity, namely, the
D. Converter-Interfaced Generation
complex frequency. This quantity allows writing the differen-
This last example illustrates an application of the approx- tial of power flow equations in terms of complex powers and
imated expression (65). We focus in particular on the link voltages at network buses. The most significant property of the
between ω and q̇ through the resistances of network branches. newly defined complex frequency is its ability to give a more
Using again the WSCC 9-bus system and substituting two robust and clean indication of frequency than what generally
synchronous machines with CIGs, we compare the dynamic accepted in the literature, especially to describe the behavior
response of the system following a load variation using the of the frequency at buses close to a disturbance. For example,
control scheme of Fig. 3 (Control 1), and the control scheme it is well known that many other methods result in meaningless
shown in Fig. 5 (Control 2). The latter regulates the frequency spikes in frequency at the inception and clearing of faults and
by both the active and reactive powers of the CIGs. other sudden disturbances. The proposed definition provides a
Figure 6 shows that Control 2 is more effective than Control solution to this issue.
1 to reduce the variations of the frequency. This result indicates The proposed definition of complex frequency is in accor-
that the relationship between q and ω is not weak and can dance with the commonly accepted definition of frequency and
be exploited to improve the frequency response of low-inertia generalizes it. In fact, if the magnitude of the voltage is con-
systems. This result is predicted by (65) and hence by (32). stant, then ρ = 0 and the magnitude of the complex frequency
 10

60.1
A PPENDIX I
60 T IME -D ERIVATIVE OF THE VOLTAGE
59.9 This appendix provides the proof of (37). With this aim, let
ωCoI [Hz]

59.8
us consider the h-th element of (37), namely v̄h = vh ∠θh =
vh (cos θh +  sin θh ). The time derivative of v̄h with respect
59.7 to the dq-axis reference frame gives:
Control 1
59.6 Control 2 v̄˙ h = v̇h ∠θh +  vh ωh ∠θh
59.5 = vh (%h ∠θh +  ωh ∠θh )
0 2 4 6 8 10
Time [s] = vh ∠θh (%h +  ωh )
1.025 = v̄h η̄h ,
1.02
where the following identities hold:
vbus 5 [pu(kV)]

1.015
d
∠θh = ωh (− sin θh +  cos θh )
1.01 dt
1.005 =  ωh (cos θh +  sin θh )
Control 1 =  ωh ∠θh .
1 Control 2
0.995 A PPENDIX II
0 2 4 6 8 10
Time [s] C OMPLEX F REQUENCY OF C URRENT I NJECTIONS
Fig. 6: Frequency of the CoI and voltage at bus 5 following the Let us assume that ı̄h = ıh ∠βh , then, one has:
connection of p = 0.25 pu at bus 5 at t = 1 s for the WSCC 9-bus  
ı̇
system with high penetration of CIG. ˙ı̄h = ı̄h h + β̇h = ı̄h ξ¯h , (84)
ıh
where ξ¯h is the complex frequency of the current injection at
bus h which is determined using the same procedure that leads
coincides with (1). The common definition of frequency is to (37) and that is described in the Appendix.
thus a special case of the complex frequency. On the other A relevant special case is that of constant power factor
hand, if the magnitude of the voltage is not constant, then, devices (either generators or loads), for which βh = θh + ϕho ,
the complex frequency extends the concept of frequency by where ϕho is the constant power factor angle. Then:
including an additional term, ρ.
β̇h = θ̇h + ϕ̇ho = ωh . (85)
Then, the paper shows how the complex frequencies are On the other hand, Re{ξ¯h } = %h only if ıh = kvh . Since a
related to the time derivatives of the voltages and to the constant impedance has constant power factor and its current is
power/current injections at the network buses. In this vein, proportional to the voltage, it satisfies the condition ξ¯h = η̄h ,
noteworthy expressions are (32), (39) and (41). It is also shown ∀t.
that these expressions are a generalization of the FDF proposed
in [5].
ACKNOWLEDGMENTS
The proposed framework can be used, for example, for eval- The author wishes to thank Dr. Ioannis Dassios, Univer-
uating how good is the estimation of the frequency obtained sity College Dublin, for carefully revising the mathematical
with a PLL and also appears as an useful tool to design derivations given in the paper.
novel and efficient controllers. The several analytical and
numerical examples discussed in Sections III and IV show the R EFERENCES
prospective applications of the proposed theoretical approach.
[1] “IEEE/IEC International Standard - Measuring relays and protection
The newly introduced term ρ can be also expected to have equipment - Part 118-1: Synchrophasor for power systems - Measure-
practical applications. An example is given in Section IV-C. ments,” pp. 1–78, 2018.
[2] H. Kirkham, W. Dickerson, and A. Phadke, “Defining power system
Future work will focus on further elaborating the expression frequency,” in IEEE PES General Meeting, 2018, pp. 1–5.
(32) and exploiting its features for the control and state [3] M. Paolone, T. Gaunt, X. Guillaud, M. Liserre, S. Meliopoulos,
A. Monti, T. Van Cutsem, V. Vittal, and C. Vournas, “Fundamentals
estimation of power systems. A relevant question worth of of power systems modelling in the presence of converter-interfaced
further study is how to take into account fast dynamics, generation,” Electric Power Systems Research, vol. 189, p. 106811, 2020.
such as transmission line dynamics, in the evaluation of the [4] J. Nutaro and V. Protopopescu, “Calculating frequency at loads in
simulations of electro-mechanical transients,” IEEE Trans. on Smart
complex frequency. It appears relevant also to combine the Grid, vol. 3, no. 1, pp. 233–240, 2012.
proposed complex frequency approach to some of the tech- [5] F. Milano and Á. Ortega, “Frequency divider,” IEEE Trans. on Power
niques described in [3]. Particularly interesting, for example, is Systems, vol. 32, no. 2, pp. 1493–1501, 2017.
[6] H. Golpı̂ra and A. R. Messina, “A center-of-gravity-based approach to
the “beyond-phasor” approach based on the Hilbert transform estimate slow power and frequency variations,” IEEE Trans. on Power
proposed in [21]. Systems, vol. 33, no. 1, pp. 1026–1035, 2018.
 11

[7] F. Milano and Á. Ortega, “A method for evaluating frequency regulation
in an electrical grid Part I: Theory,” IEEE Trans. on Power Systems,
2020, in press.
[8] F. Saccomanno, Electric Power Systems - Analysis and Control. New
York: John Wiley & Sons, 2003.
[9] F. Milano and Á. Ortega, Frequency Variations in Power Systems:
Modeling, State Estimation, and Control. Hoboken, NJ: Wiley, 2020.
[10] H. Akagi, E. H. Watanabe, and M. Aredes, Instantaneous Power Theory
and Applications to Power Conditioning, 2nd ed. New York: Wiley
IEEE Press, 2017.
[11] P. W. Sauer and M. A. Pai, Power System Dynamics and Stability. Upper
Saddle River, NJ: Prentice Hall, 1998.
[12] F. Milano, Power System Modelling and Scripting. London, UK:
Springer, 2010.
[13] I. Hiskens, “Power system modeling for inverse problems,” IEEE
Transactions on Circuits and Systems I: Regular Papers, vol. 51, no. 3,
pp. 539–551, 2004.
[14] M. A. A. Murad, M. Liu, and F. Milano, “Modeling and simulation of
variable limits on conditional anti-windup pi controllers for vsc-based
devices,” IEEE Transactions on Circuits and Systems I: Regular Papers,
vol. 68, no. 7, pp. 3079–3088, 2021.
[15] U. M. Ascher and L. R. Petzold, Computer Methods for Ordinary
Differential Equations and Differential-Algebraic Equations, 1st ed.
USA: Society for Industrial and Applied Mathematics, 1998.
[16] H. Stephani, Differential Equations: Their Solution Using Symmetries,
M. MacCallum, Ed. Cambridge University Press, 1990.
[17] Á. Ortega and F. Milano, “Estimation of voltage dependent load models
through power and frequency measurements,” IEEE Trans. on Power
Systems, vol. 35, no. 4, pp. 3308–3311, 2020.
[18] B. Stott and O. Alsaç, “Fast decoupled load flow,” IEEE Trans. on Power
Systems, vol. PAS-93, no. 3, pp. 859–869, 1974.
[19] F. Milano, “A Python-based software tool for power system analysis,”
in IEEE PES General Meeting, 2013, pp. 1–5.
[20] Á. Ortega and F. Milano, “Comparison of different PLL implementations
for frequency estimation and control,” in ICHQP, 2018, pp. 1–6.
[21] A. Derviškadić, G. Frigo, and M. Paolone, “Beyond phasors: Modeling
of power system signals using the hilbert transform,” IEEE Transactions
on Power Systems, vol. 35, no. 4, pp. 2971–2980, 2020.

Federico Milano (F’16) received from the Univer-

sity of Genoa, Italy, the ME and Ph.D. in Electrical
Engineering in 1999 and 2003, respectively. From
2001 to 2002, he was with the Univ. of Water-
loo, Canada. From 2003 to 2013, he was with the
Univ. of Castilla-La Mancha, Spain. In 2013, he
joined the Univ. College Dublin, Ireland, where he
is currently Professor of Power Systems Control and
Protections and Head of Electrical Engineering. His
research interests include power systems modeling,
control and stability analysis.