IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 39, NO.

1, JANUARY 2024 947

Inertia Estimation Through Covariance Matrix

Federico Bizzarri , Senior Member, IEEE, Davide del Giudice , Member, IEEE,
Samuele Grillo , Senior Member, IEEE, Daniele Linaro , Member, IEEE, Angelo Brambilla , Member, IEEE,
and Federico Milano , Fellow, IEEE

Abstract—This work presents a technique to estimate on-line B. Literature Review

the inertia of a power system based on ambient measurements.
The proposed technique utilizes the covariance matrix of these There are mainly two methods for the on-line estimation
measurements and solves an optimization problem that fits such of the inertia available in a power system: (i) methods based
measurements to the synchronous machine classical model. We on the measurement of frequency and power variations after a
show that the proposed technique is adequate to accurately es- disturbance or by injecting a probing signal and (ii) methods
timate the actual inertia of synchronous machines and also the
virtual inertia provided by the controllers of converter-interfaced
that utilize ambient measurements [3], [4]. Both methods have
generators that emulate the behavior of synchronous machines. advantages and shortcomings.
We also show that the proposed approach is able to estimate the The methods that are based on disturbances can be very accu-
equivalent damping of the classical synchronous machine model. rate [5] provided that the on-set of the disturbance is correctly
This feature is exploited to estimate the droop of grid-following con- detected [6], [7], [8], [9]. This can be challenging in low-inertia
verters, which has a similar effect of the swing equation equivalent
damping. The technique is comprehensively tested on a modified
power systems, which show a “rich” dynamical behavior [10].
version of the IEEE 39-bus system as well as on a dynamic 1479-bus Moreover, the estimation cannot be performed continuously, but
model of the all-island Irish transmission system. only following disturbances [11], [12]. Methods based on signal
probing must inject active power of adequate magnitude and
Index Terms—Inertia estimation, stochastic differential
equations, covariance matrix, ambient noise measurements, require the installation of specific equipment or the modification
synthetic inertia, online estimation. of pre-existing controllers [13].
The methods that utilize ambient measurements are concep-
tually more challenging than those based on disturbances. The
main advantage of these methods is that they employ mea-
I. INTRODUCTION surements obtained in normal operating conditions and hence
the estimation can be performed in a continuous fashion. The
A. Motivation
starting point is to build the dynamical model of the system,
HE inertia constant has become a volatile parameter in
T power systems with high penetration of renewable and
non-synchronous resources [1]. Therefore, the estimation of
which can be either complete or reduced to a small number
of coherent areas [14], [15]. Then, the values of the param-
eters and of the observable state variables are fit to available
the available inertia is a valuable information that can help measurements. It is worth mentioning that the number of areas
system operators to ensure the security of the grid [2]. This and, consequently, the number of equivalent generators, whose
paper addresses this topic and aims at developing an on-line inertia is to be estimated, is relevant. In [16], the authors raise
inertia estimation method based on ambient measurements, their the question whether the center of inertia (COI) of a single-
stochastic behavior and a fair knowledge of the grid model. machine-equivalent power system model could be identified
through ambient measurements. The answer given is that the
combined effects of a well-damped COI dominant mode and of a
Manuscript received 18 June 2022; revised 10 November 2022; accepted 30
low-pass-filter action of the Ornstein-Uhlenbeck (OU) processes
December 2022. Date of publication 12 January 2023; date of current version modelling load behavior make it practically impossible to per-
26 December 2023. The work of Davide del Giudice and Samuele Grillo was form this task. However, this conclusion does not hold when a
supported by Italian MIUR Project under Grant PRIN 2017K4JZEE_006. Paper
no. TPWRS-00896-2022. (Corresponding author: Federico Bizzarri.)
multi-machine system is considered.
Federico Bizzarri is with the Politecnico di Milano, 20133 Milano, Italy, and In the second group of inertia estimation methods, we cite,
also with the Advanced Research Center on Electronic Systems for Informa- for example, [17], where the authors obtain the fitting based on
tion and Communication Technologies E. De Castro (ARCES), University of
Bologna, 41026 Bologna, Italy (e-mail: federico.bizzarri@polimi.it).
a robust Kalman filter. In [13] such a task is performed using
Davide del Giudice, Samuele Grillo, Daniele Linaro, and Angelo Bram- a closed-loop identification method. In [18], [19], [20] inertia
billa are with the Politecnico di Milano, 20133 Milano, Italy (e-mail: da- estimation is performed by observing the step response of the
vide.delgiudice@polimi.it; samuele.grillo@polimi.it; daniele.linaro@polimi.it;
angelo.brambilla@polimi.it).
reduced-order model obtained from the full state-space model
Federico Milano is with the School of Electrical & Electronic Eng., University describing the grid, which is, in turn, obtained through parameter
College Dublin, D04V1W8 Dublin, Ireland (e-mail: federico.milano@ieee.org). identification using ambient measurements. In [21] the power
Color versions of one or more figures in this article are available at
https://doi.org/10.1109/TPWRS.2023.3236059.
grid is reduced to a two-machine system and then the system in-
Digital Object Identifier 10.1109/TPWRS.2023.3236059 ertia is estimated using inter-area oscillation modes. In [22], the

© 2023 The Authors. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see
https://creativecommons.org/licenses/by/4.0/
948 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 39, NO. 1, JANUARY 2024

swing equation of the area single-machine-equivalent generator D. Organization

is transformed in Fourier domain. Then, the electromechanical The remainder of the paper is structured as follows. Section II
oscillation, obtained by applying a fast Fourier transform to the gives the theoretical mathematical background of power sys-
generator estimated electrical power and rotor speed, are used to tems’ modeling and in Section III the mathematical model is
estimate the equivalent inertia. This method requires capturing enriched with the introduction of noise injections to the power
the whole electromechanical oscillation trajectory. The approach system, being this stage instrumental to model load variability
proposed in [23] estimates the inertia constants of each generator and perform the estimation. Section IV presents the proposed
by applying the regression theorem to a set of OU processes. estimation procedure. Then the proposed method is validated
The method utilizes the covariance and correlation matrices of through numerical simulations, which are discussed in Sec-
the state variables to estimate the state matrix of the system. tion V. This section also discusses the challenges posed by the
Then, the inertia constants are obtained solving a least-square proposed approach and presents a solution based on a proper
problem. The main limitation of [23] is the assumption that the filtering of the measurements. Finally, conclusions are drawn in
relationship between the active power variations of synchronous Section VI.
machines and the voltage angle variations, as well as the active
and reactive power variations of converter-interfaced generators II. POWER SYSTEM MODEL
are linear. In a similar way, in [24], the regression theorem of a
multivariate OU process is applied to power systems including The proposed estimation procedure assumes that the gener-
converter-interfaced generators (CIGs), modeled in detail with ators connected to the grid are synchronous machines and fits
their phase locked loops (PLLs) and supposed to work at unity the measurements to this model. In this section, we formulate
power factor. The estimation is carried out separately for CIGs the equations of the system according to this assumption. It
and synchronous generators by identifying the relevant subma- is important to note, however, that every simulation is carried
trices of the system state matrix. out utilizing a full-fledged model representing in detail the
dynamical behavior of all generators, either synchronous or
non-synchronous, and their controllers.
To describe the dynamics of the power system model (PSM),
C. Contribution we consider the following set of differential algebraic equations
The approach proposed in this work falls in the group of (DAEs)
methods that use ambient measurements. We propose to exploit δ̇ = Ω (ω − ω0 ) ,
the colored noise due to random fluctuations of the power
consumption of the loads and their impact on bus voltages, M ω̇ = P m (ω, x, y) − P e (δ, x, v, θ, y) − D(ω − ω0 )
line currents, and power flows across different areas of a power  (δ, ω, x, v, θ, y) ,
T ẋ = F
system. Then, we define analytically the statistical properties of
the colored noise of these electrical quantities as functions of 0 = G (δ, ω, x, v, θ, y) , (1)
the system parameters and of the inertia present in the system.
where, assuming that M is the number of machines, the meaning
Finally, we utilize the variance of the electrical measurements to
of symbols in (1) is as follows:
minimize a non-linear least-squares cost function with respect
to the equivalent inertia constants as well as the equivalent –Ω: the base synchronous frequency in rad/s;
damping coefficients of the generators, either synchronous (see, –ω(t) ∈ RM : the per-unit rotor speeds of the machines;
e.g., equation (3.203) in [25]) or power-electronic based, that –ω0 ∈ R: the per-unit reference synchronous frequency;
are connected to the grid. The solution is the sought value of the –δ(t) ∈ RM : the rotor angles of the machines;
inertia constant or of the damping coefficient that best fits the –M ∈ RM ×M : a diagonal matrix whose entries model the in-
measurements and that satisfies the grid constraints. ertia constants H of the machines. In particular, Mjj = 2Hjj
The main benefit of the proposed approach is that it does (for j = 1, . . . , M );
not require a contingency or a large disturbance to occur to –x(t) ∈ RN : the N state variables of the PSM (ω and δ excluded)
be able to estimate the inertia (or the damping). Only ambient that can influence the dynamics of the machines, and T ∈
measurements are needed, together with a fair knowledge of RN ×N is a mass matrix;
the model of the grid. The price for this is the need to take –y(t) ∈ RS : all the algebraic variables of the PSM but v and θ;
measurements in a given time period. Nevertheless, as shown in –D ∈ RM ×M : a diagonal matrix whose entries djj (for
the case studies, the proposed approach is accurate and provides j = 1, . . . , M ) model the damping factor of the machines;
estimations of the inertia in the same time scale of short-term –v(t) ∈ RP and θ(t) ∈ RP : bus voltages and phases, respec-
dispatch and adjustments markets. The proposed method is tively, where P is the number of buses;
original, as it shows for the first time how to perform inertia –P m (ω, x, y) ∈ RM : the mechanical power regulated by con-
estimation by resorting to the covariance matrix of the power trollers depending on ω, x, and y;
system without needing to estimate the state of the network or to –P e (δ, x, v, θ, y) ∈ RM : the electrical power exchanged by
derive (or measure) the rotor speed or any other internal variable machines;
of both synchronous machines and other devices performing  (·) accounts for regulators and other dynamics included in
–F
frequency control. the system; and
BIZZARRI et al.: INERTIA ESTIMATION THROUGH COVARIANCE MATRIX 949

–G(·) models algebraic constraints such as lumped models of the solution of (6) can be written as ξ = ξ o + ξ η and ζ = ζ o +
transmission lines, transformers, and static loads. ζ η , that is, the effects of perturbations are assumed as additive. It
is thus possible to linearise (6) around ξ o and obtain the random
Equation (1) can be conveniently rewritten as
ordinary differential equation [32]
 (ξ, ζ) ,
Λξ̇ = F  ζ G−1 Ξ η,
Λξ̇ η = F ξ ξ η − F ζ (8)
0 = G(ξ, ζ), (2)
 ζ , and Gζ matrices are computed at (ξ o , ζ o ).
where the F ξ , F
T T
where ξ(t) ≡ [ω, δ, x] , ζ(t) ≡ [v, θ, y] , and Λ is a non- The small-signal variations of the algebraic variables are given
singular diagonal matrix including M . G and F  are assumed by
to be continuously differentiable in their definition domain and
 ζ,
 ξ, F ξη
matrices of their partial derivatives are referred to as F ζ η = −G−1
ζ [ Gξ | Ξ ] . (9)
   η
Gξ , and Gζ . As an example, F  ξ = ∂F  j /∂ξk , for j, k =
jk E
1, . . ., 2M + N .
The implicit function theorem guarantees that if G(ξ ∗ , ζ ∗ ) = Augmenting (8) with (7) allows obtaining the linear SDEs (in
0, provided that Gζ (ξ ∗ , ζ ∗ ) is not singular, a unique and smooth narrow sense [26]) governing the overall dynamics of the lin-
function Γ : R2M +N → R2P +S exists so that ζ ∗ = Γ(ξ ∗ ). If earised noisy PSM. Adopting the typical formalism of the SDEs,
the conditions of the implicit function theorem are satisfied, (2) the complete set of linearised SDEs reads
can be rewritten as ⎡ J ⎤
ξ η ⎣ −1  
 ζ G−1 Ξ⎦ ξ η dt + 0 dW t ,
 (ξ, Γ(ξ)) ≡ F (ξ),
Λξ̇ = F (3) d = Λ F ξ −Λ−1 F (10)
η ζ η Σ
   0 −Υ 
the equilibrium points of which, say ξo , satisfy the condition Xt    B
A
F (ξ o ) = 0, (4)
The solution of (10) with a normally distributed (or constant)
with the Jacobian matrix initial condition is a (2M + N + Z)−dimensional Gaussian

J (ξ o ) = Λ−1F ξ = Λ−1 (F  ζ G−1Gξ )
ξ− F . (5)
stochastic process.
ζ
ξ=ξo ,ζ=Γ(ξo )
IV. PROPOSED TECHNIQUE TO ESTIMATE THE INERTIA
III. INCLUSION OF NOISE
The mean and the covariance matrix of process (10) can
We assume stochastic noise is injected into the PSM as be derived by solving two sets of linear ordinary differential
 (ξ, ζ) ,
Λξ̇ = F equations (ODEs) [26]. Since (10) is stable under the hypothesis
that (3) is stable at ξ o , these ODEs reveal that, at steady state,
0 = G(ξ, ζ) + Ξ η, (6) the mean of X t is zero, and its K Xt covariance matrix derives
from the solution of the following Lyapunov equation
where η is a vector of Z independent Ornstein-Uhlenbeck (OU)
processes [26], and Ξ ∈ R(2P +S)×Z is a constant matrix. AK Xt + K Xt AT + BB T = 0. (11)
The OU processes are defined through a set of stochastic
differential equations (SDEs), as follows The diagonal elements of K Xt are the (steady-state) variances
of the components of the X t process. In particular, the last Z
dη = −Υη dt + Σ dW t , (7) diagonal elements, corresponding to the η sub-vector of X t , are
where the drift Υ ∈ R Z×Z
and diffusion Σ ∈ R Z×Z
are di- the (steady-state) variances of the Z independent OU processes
agonal matrices with positive entries, W t ∈ RZ is a vector introduced in (7). Hence, these terms can be written as σz2 /(2υz ),
of Z Wiener processes, and the differentials rather than time for z = 1, . . ., Z, where σz and υz are the diagonal elements of
derivatives are utilized to account for the idiosyncrasies of the Σ and Υ, respectively [26]. The remaining 2M + N diagonal
Wiener processes. The OU processes are characterized by a elements of K Xt refer to the ξ η sub-vector of X t and are
mean-reversion property and show bound standard deviation. influenced by η through the sub-matrix −Λ−1 F  ζ G−1 Ξ in (10).
ζ
Moreover, these processes show a spectrum that is an accurate According to (9), ζ η can be written as a linear combination
model of the stochastic variability of power loads [27], [28], of the entries of X t through E, thus being a multidimensional
[29], [30], [31]. In (6), we assume that noise injection models Gaussian process, too. Hence, it is possible to write the K ζη
the effect of the consumption randomness of Z loads. covariance matrix of the small-signal algebraic variables as the
Note that, in (6), the vector of stochastic processes η appears quadratic form [33], [34]
only in the algebraic equations. This is assumed for simplicity K ζη = EK Xt E T . (12)
but without lack of generality. The interested reader can find for
instance in [29] a stochastic PSM model where the noise perturbs The proposed approach exploits (12) and the fact that K ζη
both differential and algebraic equations. (and, thus, the variance of the algebraic variables) depends
We also assume that the action of η can be safely modeled as through K Xt (and, thus, A) on the elements of Λ, a subset
a small-signal perturbation around an equilibrium point. Hence, of which are the inertia constants of the synchronous machines.
950 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 39, NO. 1, JANUARY 2024

Fig. 1. The flow-diagram of the proposed inertia estimation procedure.

Let us assume one wants to estimate a subset M of these

inertia constants and let Z be a set of PMU measurements of bus Fig. 2. The single-line diagram of the IEEE 39-bus power system.
voltages and line currents flowing through a given number of
transmission lines and/or transformers.
The above derivations show that, if the model and the param- the devices, which can be attained by minimizing a cost function
eters of the grid are known, one can compute the σZ2 variances for a subset of damping coefficients D, in a similar way as for
of the measured quantities, which are the diagonal elements of M in (13),
the K ζη matrix. However, it is also possible to solve an inverse  (
σZ2 − σZ2 (D))2
problem: knowing through measurements the variance of the C(D)|F ξ ,F ζ ,Gξ ,Gζ = . (14)
σZ4 (D∗ )
elements of Z, the values of the elements of M can be deter- Z
mined by minimizing the following non-linear least-squares cost
function (see also Fig. 1) V. CASE STUDIES
 (
σZ2 − σZ2 (M))2 In this Section, we consider two grids, the well-known IEEE
C(M)|F ξ ,F ζ ,Gξ ,Gζ = , (13) 39-bus system and a 1479-bus dynamic model of the all-island
σZ4 (M∗ )
Z Irish transmission system AIITS. The IEEE 39-bus system allows
where σ Z2 are the variances of the measured variables, σZ2 (M) us illustrating the various features and challenges of the proposed
are the variances computed in run-time during the minimization method. With this aim, we first consider the conventional system
procedure w.r.t. the elements of M, while σZ2 (M∗ ) are the with synchronous machines and their regulators (Section V-A1).
corresponding variances of the same variables obtained with Then, we evaluate the behavior of the proposed method when
the model for a reference set of M values, which are used to the system includes grid-forming (Section V-A2) and grid-
normalize the cost function. During the optimization process following converters (Section V-A3). The all-island Irish system
the estimation of the elements of M makes necessary to update, (Section V-B), on the other hand, serves to demonstrate the
iteratively, Λ and, consequently, A.1 robustness of the proposed method when applied to a real-world
We note that the procedure above is agnostic w.r.t. to the complex network.
devices that are connected to the grid. That is, the estimation of
the inertia is obtained by fitting the synchronous machine model A. IEEE 39-Bus System
to the measurements. If there exist non-synchronous devices, We use as a benchmark the IEEE 39-bus power system shown
such as converter-interfaced generation, their “equivalent” iner- in Fig. 2 [35]. This is a simplified model of the transmission
tia constants can be obtained using (13) regardless of the fact system in the New England region in the Northeastern United
that the converter is set up to be grid forming or grid following. States and is composed of 10 generators, 34 lines, 19 loads,
The procedure is also agnostic in terms of the parameters to be and 12 transformers. The network models and parameters are
estimated, provided that such parameters do have an effect on those adopted and directly extracted from its DIgSILENT Pow-
the frequency variations of the system. In particular, we use this erFactory implementation [36] and derived from [37]. The G1
property to estimate also the equivalent damping coefficients of generator, which represents the connection of the New England
System to the rest of the US and Canadian grid, is modeled with a
1 Actually if, for instance, some specific implementation of virtual syn- constant excitation [37], viz. there is no automatic voltage regu-
chronous generators are considered, it may happen that the equivalent inertia lator (AVR) connected to it. On the contrary, the other generators
constant of those devices directly affects some of the entries of the F ξ Jacobian are equipped with AVRs (specifically, rotating excitation systems
matrix instead of Λ. This situation, despite being more involved, can also
be handled by updating A during the optimization process. For the sake of of IEEE Type 1 according to [37]). Governors are considered as
simplicity, we will not dwell further on this point here. IEEE Type G1 (steam turbine) for G2 –G9 , and IEEE Type G3
BIZZARRI et al.: INERTIA ESTIMATION THROUGH COVARIANCE MATRIX 951

TABLE I
SYNCHRONOUS GENERATORS H AND PRAT

(hydro turbine) for G10 . Finally, in Table I we report the value

of inertia constants and power rating of the 10 generators.
The active power consumption of each of the Z = 19 loads
is altered by one of the ηz (t) OU processes (z = 1, . . . , Z) in
(10), which has zero mean, and variance σz2 /(2υz ). All loads Fig. 3. Violin plots of the inertia estimations of each generator of the IEEE
are assumed to incorporate stochastic power fluctuations [29] 39-bus system. The black solid circle markers correspond to the median value
 γ of the optimization results, whereas the empty square markers represent the
|Vz (t)| nominal values of the inertia constants. The magenta bars represent the IQR,
Lz (t) = ηz (t) PL0z , (15) viz. the spread difference between the 75th and 25th percentiles of the data. The
V 0z green solid circle markers represent the upper adjacent value (i.e., the largest
where (for z = 1, . . . , Z) PL0z is the nominal active power of observation that is less than or equal to the third quartile plus 1.5 × IQR) and
the lower adjacent value (i.e., the smallest observation that is greater than or
the load, V0z is the load voltage rating, Vz (t) is the bus voltage at equal to the first quartile minus 1.5 × IQR). Time window: 15 min. Optimization
which the load is connected and γ governs the dependence of the performed on the unfiltered currents.
load on bus voltage. In the simulation below, we assume υz =
0.5 and σz is defined in such a way that the standard deviation of
ηz (t) is 5% of PL0 . The zero mean implies that stochastic load The H  inertia constant estimates were obtained with the MAT-
power fluctuations do not perturb, on average, the operating point LAB Global Optimization Toolbox and the fgoalattain
of the system. Furthermore, we assumed γ = 0. function. The search interval of the optimization procedure
To carry out the simulations discussed below, the numerical was lower-bounded to zero since inertia constants cannot be
integration of the multi-dimensional OU process in (6) was negative. We performed 20 independent trials each one being
based on the numerical scheme proposed by Gillespie in [38]. a 24 h simulation of the grid and an optimization problem is
Furthermore, the second-order trapezoidal implicit weak scheme solved every 15 min of simulated time, thus collecting 1920
for stochastic differential equations with colored noise, available estimates per synchronous generator. The violin plots3 in Fig.
in the simulator PAN [39], [40], [41], was adopted [42]. 3 (one for each generator) summarize the obtained results. The
1) Conventional Power System: The objective of the first median of the estimated value of the inertia constants (black solid
scenario considered in this case study is to estimate the inertia dots in Fig. 3) are in good agreement with the corresponding
constants of the 10 generators of the IEEE 39-bus system across nominal values (empty squared markers in Fig. 3) listed in
a time period spanning one day. The target is to obtain an Table I. Despite the amplitude of the interquartile range (IQR)
estimation of the inertia constants every 15 min (estimation intervals (magenta vertical bars in Fig. 3) are reasonably narrow,
window) of the working day. To this end, we compute the the upper and lower adjacent values (green solid dots in Fig.
variance of the currents flowing through the generators, using 3) are quite far from the median values, and the same holds
a 15 min-long moving window. The values of variance are the for the tails of the distribution of the estimated values. The
elements of the M set introduced in Section IV. The sampling maximum value of the relative percentage error, computed as
period is 25 ms, i.e., a sampling rate of 40 Hz.2 H 
|Median(H)−H|
We also assume that the matrices F̃ ξ , F̃ ζ , Gξ , Gζ are % = 100 · H , is obtained for G5 and amounts to
3.87%. The 0.44% minimum value is achieved for G3 .
constant, i.e., that the parameters of the IEEE 39-bus system
The quality of the results can be improved by filtering the
do not vary during the simulation. This assumption simplifies
considered currents before computing their moving variance.
the analysis but is not a binding requirement of the proposed
This is done by resorting to a steep bandpass filter acting from
procedure. If the operating point and/or the topology of the grid
0.1 Hz to 1.5 Hz. This bandwidth was chosen to remove the
change, one just needs to update those matrices.
low frequency contribution of the OU process. This is done since
these components account for the very slow fluctuations of the
2 The steady-state variance of a given time-domain stochastic signal corre-
currents that heavily affect the trend of their moving variance
sponds to the integral of its power spectral density (PSD) in the frequency domain. w.r.t. the steady-state values of the variance itself. The effect of
When computing the variance of a given set of samples of such a signal, one must
choose a sampling frequency and a sampling window such that the corresponding this filter is shown in Fig. 4 for the i8 current (i.e., the current
PSD is accurate enough w.r.t. the continuous-time original signal. By doing so,
the steady-state variance computed in time domain from a finite set of samples
will be accurate enough. According to several numerical tests, we derived that 3 A violin plot is a combination of a box plot and a kernel density plot.
a sampling period of 25 ms in a sampling window of 15 min is a good trade-off Specifically, it starts with a box plot and then adds a rotated kernel density
between accuracy and computational burden. plot to each side of the box plot.
952 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 39, NO. 1, JANUARY 2024

Fig. 4. The E% percentage relative error of the moving variance of the i8

current, calculated over a sliding window of 15 min, with respect to the steady-
Fig. 6. Upper panel: time evolution of the λ coefficient used to overload the
state variance of i8 . The black and the red curves refer to the unfiltered and
IEEE 39-bus system. Lower panel: the overall active power supplied by the solar
filtered case, respectively.
plants that were added to the IEEE 39-bus system.

synchronous machines and are thus expected to provide a virtual

inertia to the system. The controllers of grid-forming converters
are based on [43]. The objective of this section is twofold: (i)
showing that the proposed inertia estimation approach works
correctly under variable generation/load conditions and (ii) il-
lustrating the ability of the proposed estimation technique to
capture the effect of grid-forming converters.
For what regards point (i), we first modified the IEEE 39-bus
system by mimicking a typical daily variation of the loads. To
do so, we multiplied the power absorbed by the loads and the
(active) power generated by the synchronous machines by a
coefficient λ. To alter the nominal value of the power loads, λ
Fig. 5. Violin plots of the inertia estimations of each generator of the IEEE was varied continuously in time according to the profile shown in
39-bus system. The optimization was performed on the filtered currents, on a the upper panel of Fig. 6. This profile was sampled every 15 min:
time window of 15 min. Note the different scale of the vertical axes w.r.t. Fig. 3.
this same time basis was used to update the set points of the gen-
erators. The effect of energy production by three (aggregated)
solar plants connected to bus8 , bus24 , and bus27 was emulated
flowing through G8 ) in terms of the absolute value of the relative by reducing the absorbed active power of the loads connected at
percentage error between the moving variance and its steady- those buses. Load power was decreased by the same active power
state value over 24 h in one of the 20 trials. In particular, the supplied by those (aggregated) solar plants (each one supplying
black and red curves refer to the unfiltered and the filtered case, one third of the electrical power reported in the lower panel
respectively, and the improvement over the unfiltered case is of Fig. 6). To balance generation, the power supplied by these
evident. plants was subtracted from the power of synchronous generators.
The state equations governing the dynamics of the filters were In particular, the active power supplied by G3 , G7 , and G8 was
added to the set of ODEs reported in (3) and this allowed to derive varied in proportion to their power ratings. Concerning point
the steady-state variance of the filters output solving (11). The (ii), the power rating of each one of the virtual synchronous
inertia constants of the IEEE 39-bus generators was estimated generators (i.e., the solar plants connected to bus8 , bus24 , and
from the filtered currents and the results are reported in Fig. 5. bus27 ) is 100 MW. We regulated these generators so that each
The effect of this signal processing is a significant reduction of of them provided, through the controllers of their inverters, a
the amplitude of the interval spanned by the upper and lower total virtual inertia equivalent to H = 10 s only from 8:15 am
adjacent values, thus guaranteeing a more reliable estimate of to 5:30 pm.
the inertia constants. The H % relative percentage error spans For what regards point (ii), results are shown in Fig. 7, using
from 0.14% (for G2 ) to 1.78% (for G9 ). These results illustrate again violin plots. The left panel refers to the time interval in
well that the proposed approach accurately estimates the inertia which the solar plants are not connected and thus the inertia
constants of the ten synchronous generators of the system. constants of their grid-forming converters is zero. The right panel
2) Inclusion of Grid-Forming Converters: This section il- corresponds to the 8:15 am – 5:30 pm time interval. In the high-
lustrates the behavior of the proposed estimation technique inertia case, the H % percentage error for GFM1 , GFM2 , and
for a scenario where the IEEE 39-bus system is modified to GFM3 is equal to 0.29%, 1.23%, and 0.60%, respectively. In the
include a variable generation/load profile and grid-forming con- low-inertia case, H % is almost the same for the three grid-forming
verters. These devices are assumed to mimic the behavior of converters and amounts to 0.02%.
BIZZARRI et al.: INERTIA ESTIMATION THROUGH COVARIANCE MATRIX 953

Fig. 7. Violin plots of the inertia constants estimates of the grid-forming

converters GFM1 , GFM2 , and GFM3 (all of which were set to 10 s from
8:15 am to 5:30 pm). Left and right panels refer to the low- and high-inertia
cases, respectively. The time window is 15 min.

The estimation is accurate even though we used only the

Fig. 8. From top to bottom: inertia constant estimates over time of the grid-
measurements of the current flowing through the lines connect- forming converters GFM1 , GFM2 , and GFM3 (all of which were set to 10 s
ing bus1 − bus39 , bus9 − bus39 , bus3 − bus4 , bus16 − bus17 , from 8:15 am to 5:30 pm). The results of the 50 estimation procedures are
and bus14 − bus15 , (i.e., the branches indicated with dashed depicted in black. In red we report the median value computed every 15 min
over the 50 trials.
lines in Fig. 2), which are far away from the buses connecting
the solar power plants. As reported in the literature, these lines
split the IEEE 39-bus system in areas. TABLE II
The measurements used in this case highlight that we do not ESTIMATION OF THE EQUIVALENT INERTIA AND OF THE DAMPING/DROOP FOR
THE GRID-FOLLOWING CONVERTERS IN THE MODIFIED IEEE 39-BUS SYSTEM
necessarily need to put PMUs in front of synchronous generators
or grid-forming/following converters. However, we do not have
a systematic method that clearly indicates the optimal locations
where PMUs should be connected to provide measurements used
during the estimation process. This is an open issue that we will
cope with in future work.
Moreover, we observe the birth of additional “swing equation”
modes due to the presence of this kind of the virtual synchronous
generators [44]. These, in fact, induce inter-area oscillations with inclusion of the same solar power plants considered in the
that are clearly identified by peaks in the frequency responses previous section. However, we assume that these plants are now
by frequency scan. Modifying the inertia constant of these equipped with conventional frequency droop controllers.
elements changes significantly the variance of the measured Table II shows that the estimation of the damping/droop and
electrical quantities and, thus, helps increasing the accuracy of of the inertia are strongly correlated. An increase in the droop
the proposed technique. coefficients of the frequency controllers leads to a lower inertia
Fig. 8 shows the estimation of the inertia of the system across constant estimation. This result had to be expected as a lower
the full day of measurements. This illustrates the effectiveness droop makes the response of the controller slower and thus its
of the proposed technique to continuously estimate the inertia time scale tends to overlap with that of the inertial response. The
by stochastic fluctuations. In Fig. 8, the black traces are single estimation of the inertia constant and of the damping coefficient
trials while the red one is the median value over 50 trials. All are performed independently, since two separate least-squares
realisations give accurate results with only a few time instants problems are solved: one fits the inertia constant and a second
being affected by a relatively higher error. one fits the damping coefficient.
3) Inclusion of Grid-Following (GFL) Converters: In the The main conclusion that can be drawn from the results shown
third and last scenario, we discuss the effect of the frequency in Table II is that the value of the inertia constant alone is
droop control of grid-following converters on the estimation not sufficient to define the ability of a device to regulate the
of the inertia with the proposed technique. As is well-known, frequency. If only the inertia constant is estimated, in fact, one
the frequency droop control of the converter is equivalent to a would conclude that the case with higher droop coefficient is
second order model, and is thus similar to that of synchronous the one with lower frequency containment, which is not. On
machines, except for the fact that it shows a large damping the one hand, if both parameters, namely inertia and damping,
and approximately zero inertia [45]. In this scenario, thus, we are estimated, then one can conclude, correctly, that the system
repeat the simulations of the modified IEEE 39-bus system with higher droop provides “more” fast frequency control than
954 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 39, NO. 1, JANUARY 2024

TABLE III
SYNCHRONOUS GENERATORS H AND PRAT OF THE AIITS

inertial response. On the other hand, lower droop values increase

the ability of the system to have an inertial response but lead,
as expected, to a lower capability of providing a service as fast
Fig. 9. Violin plots of the inertia constants estimates for some generators of
frequency controllers. the AIITS, normalized w.r.t. their effective value.

B. All-Island Irish Transmission System (AIITS)

proposed method when applied to a real-world complex system
In this section, the proposed method is applied to a real-world with several hundreds of buses and dynamic devices.
complex system. The model of the AIITS considered here consists
of 1479 buses, 1851 transmission lines or transformers, 245
VI. CONCLUSION
loads, 22 conventional synchronous power plants with AVRs and
turbine governors, 6 power system stabilizers, and 169 wind This work proposes a technique to estimate the inertia based
power plants. Note that the secondary frequency control of the on the variance of measurements. Assuming a given model of the
AIITS is implemented manually and, thus, no AGC is considered in power system and the dynamics of the generators, the technique
the model. Both load fluctuations and wind speeds are modeled fits this model to the measurements by solving a least-squares
as OU processes. The resulting set of DAEs for the AIITS includes problem. Simulation results show that the proposed technique,
2118 state variables (663 of which are stochastic processes) and given an adequate filtering of the measurements, is accurate
6175 algebraic variables. Both the active and the reactive power even in a real-world complex system and can take into account
of the Z = 245 loads is altered by an OU process characterised both conventional and virtual synchronous machines, as well as
by vz = 0.1, while σz is defined in such a way that the standard the effect of frequency droop controllers. However, the correct
deviation of those processes is 0.5% of the nominal active or interpretation of the effect of droop controllers is consistent only
reactive power of each load. Finally, the standard deviation and if also the damping/droop coefficients are estimated, which the
drift term of the OU processes modeling the wind speeds are set proposed technique duly allows to obtain.
to 0.002 and 0.02, respectively. The time domain simulations of In particular, we believe that the proposed technique can
the AIITS are obtained with the software tool Dome [46]. be useful to reward the frequency support provided by non-
To test the applicability of the proposed method, we estimated synchronous devices. In turn, we recommend that the reward
the H inertia constants of 11 conventional synchronous power should be evaluated based on the effect that the controllers of
plants chosen to be representative of non-homogeneous H con- these devices have on the system, not on the actual implemen-
stants and power ratings (see Table III). An optimization problem tation of the control itself.
was solved every 15 min of simulated time by collecting 100 Future work will focus on further investigating the theoretical
estimates per synchronous generator. Making reference to the and practical aspects of the proposed technique. For example,
diagram in Fig. 1, the subset Z contains only the magnitude relevant questions are how to improve the accuracy of the
of the voltage at the buses where the generators are connected. estimated inertia constants through (i) an optimal number and
Hence we did not resort to the phase of such voltages, to any location of the set of measurements and (ii) of the filtering of
current of the network, or to any rotor speed of the synchronous the measurements.
generators.
The H % relative percentage error spans from 0.22% (for G7 ) REFERENCES
to 5.42% (for G10 ). The quality of the results, measured through
H [1] A. Ulbig, T. S. Borsche, and G. Andersson, “Impact of low rotational
% , is similar to that obtained for the IEEE 39-bus system. inertia on power system stability and operation,” IFAC Proc. Volumes,
For the sake of readability, the violin plots in Fig. 9 report the vol. 47, no. 3, pp. 7290–7297, 2014. [Online]. Available: https://www.
 estimated inertia values normalized w.r.t. the
statistics of the H sciencedirect.com/science/article/pii/S1474667016427618
[2] F. Milano, F. Dörfler, G. Hug, D. J. Hill, and G. Verbič, “Foundations and
H nominal inertia of each generator. So doing, for each violin challenges of low-inertia systems (invited paper),” in Proc. Power Syst.
plot in Fig. 9 the reference value is 1. The distance of the upper Comput. Conf., 2018, pp. 1–25.
and lower adjacent values (green solid dots in Fig. 9) from the [3] B. Tan, J. Zhao, M. Netto, V. Krishnan, V. Terzija, and Y. Zhang, “Power
expected normalized value of H  is lower than 0.22 for all the 11 system inertia estimation: Review of methods and the impacts of converter-
interfaced generations,” Int. J. Elect. Power Energy Syst., vol. 134, 2022,
synchronous generators, which confirms the robusteness of the Art. no. 107362.
BIZZARRI et al.: INERTIA ESTIMATION THROUGH COVARIANCE MATRIX 955

[4] S. C. Dimoulias, E. O. Kontis, and G. K. Papagiannis, “Inertia estimation [28] C. Nwankpa and S. Shahidehpour, “Colored noise modelling in the relia-
of synchronous devices: Review of available techniques and comparative bility evaluation of electric power systems,” Appl. Math. Model., vol. 14,
assessment of conventional measurement-based approaches,” Energies, no. 7, pp. 338–351, 1990.
vol. 15, no. 20, 2022, Art. no. 7767. [29] F. Milano and R. Zárate-Miñano, “A systematic method to model power
[5] U. Tamrakar, N. Guruwacharya, N. Bhujel, F. Wilches-Bernal, T. M. systems as stochastic differential algebraic equations,” IEEE Trans. Power
Hansen, and R. Tonkoski, “Inertia estimation in power systems using Syst., vol. 28, no. 4, pp. 4537–4544, Nov. 2013.
energy storage and system identification techniques,” in Proc. Int. Symp. [30] C. Roberts, E. M. Stewart, and F. Milano, “Validation of the Ornstein-
Power Electron., Elect. Drives, Automat. Motion, 2020, pp. 577–582. Uhlenbeck process for load modeling based on µPMU measurements,” in
[6] P. Wall and V. Terzija, “Simultaneous estimation of the time of disturbance Proc. Power Syst. Comput. Conf., 2016, pp. 1–7.
and inertia in power systems,” IEEE Trans. Power Del., vol. 29, no. 4, [31] H. Hua, Y. Qin, C. Hao, and J. Cao, “Stochastic optimal control for energy
pp. 2018–2031, Aug. 2014. internet: A bottom-up energy management approach,” IEEE Trans. Ind.
[7] D. Zografos, M. Ghandhari, and R. Eriksson, “Power system inertia esti- Informat., vol. 15, no. 3, pp. 1788–1797, Mar. 2019.
mation: Utilization of frequency and voltage response after a disturbance,” [32] H. Xiaoying and P. Kloeden, Random Ordinary Differential Equations and
Elect. Power Syst. Res., vol. 161, pp. 52–60, 2018. Their Numerical Solution (Probability Theory and Stochastic Modelling
[8] P. M. Ashton, C. S. Saunders, G. A. Taylor, A. M. Carter, and M. E. Series). Berlin, Germany: Springer, 2017.
Bradley, “Inertia estimation of the GB power system using synchrophasor [33] S. Provost and A. Mathai, Quadratic Forms in Random Variables: Theory
measurements,” IEEE Trans. Power Syst., vol. 30, no. 2, pp. 701–709, and Applications (Statistics: Textbooks and Monographs Series). New
Mar. 2015. York City, NY, USA: Marcel Dekker, 1992.
[9] D. del Giudice and S. Grillo, “Analysis of the sensitivity of extended [34] M. Adeen et al., “On the calculation of the variance of algebraic variables
Kalman filter-based inertia estimation method to the assumed time of in power system dynamic models with stochastic processes,” IEEE Trans.
disturbance,” Energies, vol. 12, 2019, Art. no. 483. Power Syst., early access, doi: 10.1109/TPWRS.2022.3226076.
[10] E. Heylen, F. Teng, and G. Strbac, “Challenges and opportunities of inertia [35] T. Athay, R. Podmore, and S. Virmani, “A practical method for the direct
estimation and forecasting in low-inertia power systems,” Renewable analysis of transient stability,” IEEE Trans. Power App. Syst., vol. PAS- 98,
Sustain. Energy Rev., vol. 147, 2021, Art. no. 111176. no. 2, pp. 573–584, Mar./Apr. 1979.
[11] R. K. Panda, A. Mohapatra, and S.C. Srivastava, “Online estimation of [36] PowerFactory, “39 bus New England system,” DIgSILENT GmbH, Tech.
system inertia in a power network utilizing synchrophasor measurements,” Rep. r1338, 2014. [Online]. Available: https://qdoc.tips/39-bus-new-
IEEE Trans. Power Syst., vol. 35, no. 4, pp. 3122–3132, Jul. 2020. england-system-pdf-free.html
[12] G. Cai, B. Wang, D. Yang, Z. Sun, and L. Wang, “Inertia estimation [37] M. A. Pai, Energy Function Analysis for Power System Stability (Kluwer
based on observed electromechanical oscillation response for power International Series in Engineering and Computer Science, Power Elec-
systems,” IEEE Trans. Power Syst., vol. 34, no. 6, pp. 4291–4299, tronics & Power Systems). Boston, MA, USA: Kluwer, 1989.
Nov. 2019. [38] D. T. Gillespie, “Exact numerical simulation of the Ornstein-Uhlenbeck
[13] J. Zhang and H. Xu, “Online identification of power system equivalent in- process and its integral,” Phys. Rev. E, vol. 54, pp. 2084–2091, Aug. 1996.
ertia constant,” IEEE Trans. Ind. Electron., vol. 64, no. 10, pp. 8098–8107, [39] F. Bizzarri and A. Brambilla, “PAN and MPanSuite: Simulation vehicles
Oct. 2017. towards the analysis and design of heterogeneous mixed electrical sys-
[14] G. R. Moraes et al., “Measurement-based inertia estimation method con- tems,” in Proc. New Gener. CAS, Genova, Italy, Sep. 2017, pp. 1–4.
sidering system reduction strategies and dynamic equivalents,” in Proc. [40] F. Bizzarri, A. Brambilla, G. Storti Gajani, and S. Banerjee, “Simulation
IEEE Milan PowerTech, 2019, pp. 1–6. of real world circuits: Extending conventional analysis methods to circuits
[15] K. Tuttelberg, J. Kilter, D. Wilson, and K. Uhlen, “Estimation of power described by heterogeneous languages,” IEEE Circuits Syst. Mag., vol. 14,
system inertia from ambient wide area measurements,” IEEE Trans. Power no. 4, pp. 51–70, Oct.–Dec. 2014.
Syst., vol. 33, no. 6, pp. 7249–7257, Nov. 2018. [41] D. Linaro, D. del Giudice, F. Bizzarri, and A. Brambilla, “Pansuite: A
[16] A. Gorbunov, J. C. -H. Peng, J. W. Bialek, and P. Vorobev, “Can center- free simulation environment for the analysis of hybrid electrical power
of-inertia model be identified from ambient frequency measurements?,” systems,” Electric Power Syst. Res., vol. 212, 2022, Art. no. 108354.
IEEE Trans. Power Syst., vol. 37, no. 3, pp. 2459–2462, May 2022. [42] G. N. Milshtein and M. V. Tret’yakov, “Numerical solution of differential
[17] J. Zhao, Y. Tang, and V. Terzija, “Robust online estimation of power equations with colored noise,” J. Stat. Phys., vol. 77, no. 3, pp. 691–715,
system center of inertia frequency,” IEEE Trans. Power Syst., vol. 34, 1994.
no. 1, pp. 821–825, Jan. 2019. [43] B. Barać, M. Krpan, T. Capuder, and I. Kuzle, “Modeling and initialization
[18] F. Zeng, J. Zhang, G. Chen, Z. Wu, S. Huang, and Y. Liang, “Online estima- of a virtual synchronous machine for power system fundamental frequency
tion of power system inertia constant under normal operating conditions,” simulations,” IEEE Access, vol. 9, pp. 160116–160134, 2021.
IEEE Access, vol. 8, pp. 101426–101436, 2020. [44] S. Hadavi, S. P. Me, B. Bahrani, M. Fard, and A. Zadeh, “Virtual syn-
[19] V. Baruzzi, M. Lodi, A. Oliveri, and M. Storace, “Analysis and improve- chronous generator versus synchronous condensers: An electromagnetic
ment of an algorithm for the online inertia estimation in power grids with transient simulation-based comparison,” CIGRE Sci. Eng., vol. 2022,
RES,” in Proc. IEEE Int. Symp. Circuits Syst. 2021, pp. 1–5. no. 24, pp. 1–19, 2022.
[20] V. Baruzzi, M. Lodi, A. Oliveri, and M. Storace, “Estimation of inertia [45] S. D’Arco and J. A. Suul, “Equivalence of virtual synchronous machines
in power grids with turbine governors,” in Proc. IEEE Int. Symp. Circuits and frequency-droops for converter-based MicroGrids,” IEEE Trans.
Syst., 2022, pp. 1–5. Smart Grid, vol. 5, no. 1, pp. 394–395, Jan. 2014.
[21] D. Yang et al., “Ambient-data-driven modal-identification-based approach [46] F. Milano, “A Python-based software tool for power system analysis,” in
to estimate the inertia of an interconnected power system,” IEEE Access, Proc. IEEE Power Energy Soc. Gen. Meeting, 2013, pp. 1–5.
vol. 8, pp. 118799–118807, 2020.
[22] B. Wang, D. Yang, G. Cai, Z. Chen, and J. Ma, “An improved electrome-
chanical oscillation-based inertia estimation method,” IEEE Trans. Power
Syst., vol. 37, no. 3, pp. 2479–2482, May 2022. Federico Bizzarri (Senior Member, IEEE) was born
[23] J. Guo, X. Wang, and B.-T. Ooi, “Online purely data-driven estimation in Genoa, Italy, in 1974. He received the Laurea
of inertia and center-of-inertia frequency for power systems with vsc- (M.Sc.) five-year degree (summa cum laude) in elec-
interfaced energy sources,” Int. J. Elect. Power Energy Syst., vol. 137, tronic engineering and the Ph.D. degree in electrical
2022, Art. no. 107643. engineering from the University of Genoa, Genoa,
[24] J. Guo, X. Wang, and B.-T. Ooi, “Estimation of inertia for synchronous Italy, in 1998 and 2001, respectively. Since October
and non-synchronous generators based on ambient measurements,” IEEE 2018, he has been an Associate Professor with the
Trans. Power Syst., vol. 37, no. 5, pp. 3747–3757, Sep. 2022. Electronic and Information Department, Politecnico
[25] P. Kundur, Power System Stability and Control. New York, NY, USA: di Milano, Milan, Italy. He is currently a Research
McGraw-Hill, 1994. Fellow with the Advanced Research Center on Elec-
[26] L. Arnold, Stochastic Differential Equations. Hoboken, NJ, USA: Wiley, tronic Systems for Information and Communication
1974. Technologies E. De Castro (ARCES), University of Bologna, Italy. From 2012 to
[27] R. H. Hirpara and S. N. Sharma, “An Ornstein-Uhlenbeck process- 2015, he was an Associate Editor for the IEEE TRANSACTIONS ON CIRCUITS AND
driven power system dynamics,” IFAC-PapersOnLine, vol. 48, no. 30, SYSTEMS — PART I. He was recipient as one of the 2012–2013 Best Associate
pp. 409–414, 2015. Editors for IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS — PART I.
956 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 39, NO. 1, JANUARY 2024

Davide del Giudice (Member, IEEE) was born in Angelo Brambilla (Member, IEEE) received the
Milan, Italy, in 1993. He received the M.S. degree Dr. Ing. degree in electronics engineering from the
and the Ph.D. degree in electrical engineering from University of Pavia, Pavia, Italy, in 1986. He is
the Polytechnic of Milan, Milan, Italy, in 2017 and currently a Full Professor with the Dipartimento di
2022, respectively. He is currently a Researcher with Elettronica, Informazione e Bioingegneria, Politec-
the Department of Electronics, Information and Bio- nico di Milano, Milano, Italy. His research interests
engineering, Polytechnic of Milan. His main research include the areas of circuit analysis, simulation, and
include simulation techniques for electric power sys- modeling.
tems with a high penetration of converter-interfaced
elements, such as high voltage direct current sys-
tems, electric vehicles, and generation fuelled by
renewables.

Samuele Grillo (Senior Member, IEEE) received

the Laurea degree in electronic engineering, and the
Federico Milano (Fellow, IEEE) received the M.E.
Ph.D. degree in power systems from the University of
and Ph.D. degrees in electrical engineering from the
Genova, Genoa, Italy, in 2004 and 2008, respectively.
University of Genoa, Genoa, Italy, in 1999 and 2003,
He is currently an Associate Professor with the Dipar-
respectively. From 2001 to 2002, he was with the
timento di Elettronica, Informazione e Bioingegneria,
University of Waterloo, Waterloo, ON, Canada, as
Politecnico di Milano, Milan, Italy. Since 2018, he
a Visiting Scholar. From 2003 to 2013, he was with
has been a contributor to CIGRÉ Working Group
B5.65 Enhancing Protection System Performance by the University of Castilla-La Mancha, Ciudad Real,
Spain. In 2013, he joined the University College
Optimising the Response of Inverter-Based Sources.
Dublin, Dublin, Ireland, where he is a full profes-
His research interests include smart grids, integration
sor. In 2022, he was appointed Co-Editor-in-Chief
of distributed renewable sources, and energy storage devices in power networks,
optimization, and control techniques applied to power systems. of the IET Generation Transmission & Distribution.
In 2023, he was appointed Senior Editor of the IEEE Transactions on Power
Systems. His research interests include power system modeling, control, and
stability analysis.

Daniele Linaro (Member, IEEE) received the M.Sc.

degree in electronic engineering and the Ph.D. de-
gree in electrical engineering from the University
of Genoa, Genoa, Italy, in 2007 and 2011, respec-
tively. Since 2018, he has been an Assistant Pro-
fessor with the Department of Electronics, Informa-
tion Technology and Bioengineering, Polytechnic of
Milan, Milan, Italy. His main research interests in-
clude the area of circuit theory, nonlinear dynamical
systems with applications to electronic oscillators
and power systems, computational neuroscience, and
biophysically-realistic single-cell models of neuronal cells.

Open Access provided by ‘Politecnico di Milano’ within the CRUI CARE Agreement