Electrical Power and Energy Systems 153 (2023) 109279

Contents lists available at ScienceDirect

International Journal of Electrical Power and Energy Systems

journal homepage: www.elsevier.com/locate/ijepes

Modeling microgrid power flow dynamics with covariance based realization

Yangsheng Hu ∗, Raymond A. de Callafon
Department of Mechanical and Aerospace Engineering, University of California San Diego, La Jolla, 92094-0411, CA, USA

ARTICLE INFO ABSTRACT

Keywords: Reducing transient effects of power quality in the distribution network of a microgrid requires a better
Covariance function understanding of the dynamics of power flow and distributed energy resources (DERs) present within the
Microgrid microgrid. Instead of labor intensive and detailed physics-based distribution network and DER modeling, this
Nonlinearity
paper shows how power flow dynamics can also be modeled in a dynamic equivalent with a data-driven
Power flow
approach that exploits the time dependent covariance of power flow at different locations within a microgrid.
System identification
The proposed approach results in multi-input multi-output (MIMO) linear time invariant (LTI) discrete-time
filter concatenated by a static nonlinearity to model both dynamics and power losses. Estimation of the MIMO
LTI dynamics is done by the Covariance Based Realization Algorithm (CoBRA) that can handle noisy perturbed
measurements and provide a low order MIMO LTI filter with a normalized steady-state gain. The additional
static nonlinearity is found by the estimation of the admittance matrix of the distribution network within
the microgrid. The main contribution of this paper is proposing a physics-informed data-driven approach to
modeling microgrid power flow dynamics, where the CoBRA method is implemented to learn parameters from
data while observing a priori physical constraints. The approach is illustrated on experimental data from a
microgrid with multiple DERs and it is shown that an excellent agreement in power flow dynamics is obtained
over a large operating range, which demonstrates its advantage over the N4SID method.

1. Introduction control algorithms with guaranteed stability or robustness margins.

Although physics driven component- and network-based modeling may
Microgrids typically consist of a set of power producing Distributed provide a dynamic model of the power flow in a microgrid, uncertainty
Energy Resources (DERs) interconnected with power consuming as- in component specifications, model complexity and the labor intensive
sets such as lighting, HVAC system, Electric Vehicles (EV) and other modeling work [10] may render the resulting model still unsuitable for
consumer loads. The DERs may include fossil fuel generation, sus- actual feedback control design.
tainable energy generation (wind, solar, etc.), and energy storage sys- Instead of a labor intensive and physics-based distribution network
tem (battery-powered inverters). The microgrid can be either grid- and DER modeling, power flow dynamics can also be modeled with
connected or islanding where the microgrid is disconnected from the a data-driven approach. Early approaches have shown that dynamic
main grid [1]. This flexibility allows autonomous operation of the equivalents [11] are a viable alternative for modeling power flow
microgrid during main grid faults possibly due to events such as ex- dynamics leading to simplified models [12–14]. Although the concept
treme weather, blackouts, and attacks on the power infrastructure to
of formulating a simplified dynamic equivalent of a microgrid is not
provide resilience of the power system [2–4]. Furthermore, microgrids
new, the dynamic modeling is still often limited to modeling the
with modern fast DERs/inverters can help reduce fluctuations in power
effect of active (real) power flow from a single bulk power source
flow and power quality [5–7] and improve operating efficiency of the
or individual DERs to electric values at the PCC [15] via System
(micro)grid [8].
Identification tools [16,17] leading to a scalar linear model. Also,
The ability for fast control of power flow and power quality highly
advanced system identification-based modeling techniques have drawn
depends on real-time feedback control for both active and reactive bi-
much attention in modeling power flow dynamics due to increasing
directional power flow at the point of common coupling (PCC) of a
microgrid. For example, the size of the power flow at the PCC is impor- complexity of power systems and unavailability of physics-based mod-
tant for seamless transition from grid-connected to island operation of els of most inverter-based resources [18]. Among them, the Subspace
the microgrid [9]. For such feedback power control, a dynamic model State Space System Identification (N4SID) method [16,19–23] is popu-
describing the power flow from DERs to the PCC is desired to design lar for identifying a Multi-Input, Multi-Output (MIMO) linear dynamic

∗ Corresponding author.
E-mail address: yah071@ucsd.edu (Y. Hu).

https://doi.org/10.1016/j.ijepes.2023.109279
Received 30 January 2023; Received in revised form 29 April 2023; Accepted 23 May 2023
Available online 8 June 2023
0142-0615/© 2023 Elsevier Ltd. All rights reserved.
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

equivalent. In [19], the N4SID method is used to characterize power

system dynamics by injecting low-level probing signals. A linear dy-
namic equivalent is identified to capture the dominant dynamics for the
electro-mechanical subsystem of microturbines in [20]. In [21], a data-
driven low-order MIMO linear dynamic equivalent of the microgrid
is developed for dynamical analysis and model predictive control for
secondary voltage regulation. The identified models in these work are
usually valid for only a small local operating range due to the nonlin-
earity in nature for most practical systems and the local equilibrium
assumption of the N4SID method.
This paper builds upon the concept of developing a simplified
dynamic equivalent [9], but with the ability to handle both active
and reactive power flow signals from multiple DERs at the same time.
This allows the formulation of a low order MIMO dynamic model
as modeling dynamics for each DER independently does not create a
low order MIMO model. Moreover, the approach in this paper also
models the nonlinear effect of power losses, often ignored in a linear
dynamic equivalent. Access to phasor Measurement Units (PMUs) [24]
is exploited in the proposed data-driven approach for obtaining syn-
chronized phasor data (voltage, current) and frequency and power Fig. 1. The microgrid framework with a multitude of DERs.
flow measurements at different locations in a microgrid. PMUs are
considered to be essential, as traditional SCADA does not provide the
synchronization and measurement update rate needed to model fast
2.2. Identification problem
transient phenomena considered here.
The objective of this paper is to present a physics-informed data-
Referring back to Fig. 1, synchronized power flow measurements are
driven approach that supports modeling of microgrid dynamics with
provided by PMUs that measure the pair of complex-valued (positive
a multitude of DERs operating simultaneously and incorporates power
loss through the microgrid. The end result is a MIMO microgrid dy- sequence) current 𝐈 and voltage 𝐕 phasors and generate the corre-
namic model that consists of a linear dynamic part and a static nonlin- sponding pair of active power 𝑃 and reactive power 𝑄 via 𝑃 = Re{𝐒},
ear mapping that is suitable for (distributed) feedback control design. 𝑄 = Im{𝐒}, where 𝐒 = 𝐕𝐈∗ . For notation convenience we define
For the estimation of dynamic parameters in the model, a preference is ⎡ Pin (𝑡) ⎤ out
⎡ P𝑏 (𝑡) ⎤
given to non-iterative subspace techniques [23,25] and the Covariance ⎢ 𝑏in ⎥ ⎢ out ⎥
Based Realization Algorithm (CoBRA) has been demonstrated to be ⎢Q𝑏 (𝑡)⎥ ⎢Q𝑏 (𝑡)⎥
⎢ in ⎥ ⎢ Pout (𝑡) ⎥
quite efficient in estimating low-order deterministic MIMO models on ⎢ P (𝑡) ⎥
𝐮0 (𝑡) = ⎢ 𝑑in ⎥ , 𝐲(𝑡) = ⎢ 𝑑out ⎥ (1)
the basis of large amount of input/output data heavily contaminated
⎢Q𝑑 (𝑡)⎥ ⎢Q (𝑡)⎥
by noise with unknown and possibly high order spectral character- ⎢ in ⎥ ⎢ 𝑑 ⎥
P
⎢ 𝑠 ⎥ (𝑡) ⎢ Pout
𝑠 (𝑡)
⎥
istics [26–29]. An optimal implementation for CoBRA is proposed ⎢ out ⎥
in [30,31] and shows the best performance among different subspace ⎢Qin (𝑡)⎥ ⎣Q𝑠 (𝑡)⎦
⎣ 𝑠 ⎦
methods under closed-loop settings. It will be shown that the proposed
where 𝐮0 (𝑡) is a vector of DERs active/reactive power set points, while
microgrid dynamic model will capture power flow dynamics in the
𝐲(𝑡) is a vector of active/reactive power of the DERs. Finally, active
microgrid for a large range of operating conditions. The simulation
results will illustrate that the CoBRA method shows better prediction and reactive power flow at the PCC is denoted by P𝑃 𝐶𝐶 (𝑡) and Q𝑃 𝐶𝐶 (𝑡)
accuracy than the N4SID method when the nominal power level moves respectively. Note that the 𝑃 , 𝑄 set points for the load depend on the
far away from the identification range. usage of the consumers and may not be controllable. As such, the load
is regarded as a disturbance or constant during a relatively short period
2. Problem formulation for system identification and not included in the input 𝐮0 (𝑡) and 𝐲(𝑡) in
(1).
2.1. Framework of microgrid network With the definition of these signals, the identification problem can
then be stated as: given the requested DER active and reactive power
To formalize the data-driven approach for microgrid dynamic mod- set points 𝐮0 (𝑡) and PMU-based measurements of the DER active and
eling, consider a microgrid with a network structure, as shown in reactive power output 𝐲(𝑡) in the microgrid, identify a model of the
Fig. 1. The figure illustrates that the microgrid can be either connected active power P𝑃 𝐶𝐶 (𝑡) and reactive power Q𝑃 𝐶𝐶 (𝑡) flow at the PCC. The
or disconnected (islanded) from the main grid through the PCC and objective of the identification problem is to develop a dynamic equiva-
consists of a hybrid collection of DERs interconnected with loads. The lent that models the dynamics of the power flow from the DERs to the
loads include day-to-day residential or industrial appliances that could PCC, while taking into account any nonlinearity associated with power
be modeled as resistive and inductive (R-L) elements. loss within the microgrid. This problem formulation can be stated
DERs and loads within the microgrid are assumed to be tied to
as MIMO system identification or estimation problem, solved by the
individual buses, which are further connected to the main bus at
CoBRA method. The discrete-time nature of the signals is represented
the PCC through distribution lines with certain values of admittance
by 𝑡 = 𝑘𝛥𝑇 , 𝑘 ∈ N here 𝛥𝑇 is the sampling time.
to model power loss. As mentioned in [9], fast switching between
islanding and grid connection is a strict requirement for the microgrid
to survive under main grid or microgrid faults [32]. In addition, a 2.3. Proposed model structure
model based controller design can be used for a quick and smooth
transition from grid-connected mode to islanded mode without the The microgrid dynamic modeling can be performed via a physics-
need for load shedding. The islanded to grid-connected transition can informed data-driven approach by a careful selection of dynamic and
be smoothly realized by referring to phase angle control [33]. Clearly, static components in the dynamic equivalent microgrid model, as sum-
modeling the dynamics of power flow is instrumental in understanding marized in Fig. 2. The model structure is chosen to provide high
grid-connected to islanded power flow transitions. computational efficiency: the computational complexity of the dynamic

2
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

where 𝐴 ∈ R𝑛×𝑛 , 𝐵 ∈ R𝑛×𝑝 , 𝐶 ∈ R𝑚×𝑛 , 𝐷 ∈ R𝑚×𝑝 , 𝐯(𝑡) is a

stochastic noise signal of zero mean and with unknown spectral prop-
erties and uncorrelated with the input 𝐮(𝑡), and 𝐱(𝑡) is the state vector
whose dimension 𝑛 can be determined by low-rank approximations
via the CoBRA method. The objective is to estimate the model pa-
rameters {𝐴, 𝐵, 𝐶, 𝐷} (up to a similarity transformation) as a low-
order deterministic dynamic model for the sampled input/output data
{𝐮(𝑡), 𝐲(𝑡)}.
Fig. 2. Model structure for dynamic microgrid network model.
Since a mixture of both deterministic and stochastic signals is
considered, we denote the (cross) covariance function of two quasi-
stationary vector signals 𝐬(𝑡) and 𝐰(𝑡) as below [16]:
equivalent is mostly determined by the order of the state space model { }
used to represent the linear coupled DERs dynamics. 𝑅𝑠𝑤 (𝜏) = 𝐸̄ 𝐬 (𝑡 + 𝜏) 𝐰(𝑡)T
To explain the proposed model structure, consider the following
1 ∑ {
𝑁−1
}
observations. It is worth recognizing that DERs may have power ramp = lim 𝐸 𝐬 (𝑡 + 𝜏) 𝐰(𝑡)T (3)
𝑁→∞ 𝑁 𝑡=0
constraints, which is modeled by a multi-input rate limiter fed by the
signal 𝐮0 (𝑡) of (1) in the model structure of Fig. 2. The signal 𝐮(𝑡) whose sample counterpart is
denotes the (intermediate signal) after the multi-input rate limiter and
1 ∑
𝑁−1
serves as input for the MIMO linear coupled DER dynamics. 𝑅̂ 𝑠𝑤 (𝜏) = 𝐬 (𝑡 + 𝜏) 𝐰(𝑡)T (4)
As different DERs may have different transient dynamics, the ac- 𝑁 𝑡=0
tive/reactive power signals contained in 𝐲(𝑡) may not reach the desired where 𝑁 is the finite data length and zero padding is assumed. Note
set point immediately. In addition, distribution lines with complex that 𝐸̄ {⋅} is called the extended expectation associated with the normal
admittance within a microgrid may impose coupling between active expectation operator 𝐸{⋅}. The noise effects of 𝐯(𝑡) can be eliminated
and reactive power. The different transient dynamics observed in the by correlating the system signals with an instrumental variable (IV)
ramp-rated DER power inputs 𝐮(𝑡) and the DERs power outputs 𝐲(𝑡) are 𝜉(𝑡) ∈ R𝑝×1 such that
modeled by the MIMO linear coupled DER dynamics.
Finally, without losses and power network dynamics, the PCC active • 𝜉(𝑡) is correlated with the input 𝐮(𝑡)
power flow would simply be the sum of the active power values in • 𝜉(𝑡) is uncorrelated with the noise 𝐯(𝑡)
𝐲(𝑡). In case of power loss, power flow at the PCC can be modeled
by the static network nonlinearity included in the model structure of For open-loop experiments where the input 𝐮(𝑡) excitation is uncorre-
Fig. 2. The effect of the power loss is inherently nonlinear, since it is lated with the noise 𝐯(𝑡), the IV 𝜉(𝑡) = 𝐮(𝑡).
proportional to the current squared. Based on a covariance analysis, we can formulate the following
Following the proposed model structure summarized in Fig. 2, the covariance data equation
estimation methods to find the different parts of the microgrid dynamic ̂ 𝑦𝜉 = O𝑟 𝐑
𝐑 ̂ 𝑥𝜉 + 𝛹 𝐑
̂ 𝑢𝜉 + 𝐑
̂ 𝑣𝜉 (5)
model are presented next. In Section 3, the implementation of CoBRA
is presented to identify the MIMO linear model. Section 4 proposes where
a power loss estimation to account for the nonlinearity caused by [ ( 𝑟−1 )T ]T
the distribution lines. It will be shown that the combined model will
O𝑟 = 𝐶 T (𝐶𝐴)T ⋯ 𝐶𝐴
( )
capture power flow dynamics in the microgrid over a large range of ̂ 𝑥𝜉 = Hankel 𝑅̂ 𝑥𝜉 (𝜏) , 1, 𝜏1 , 𝜏2 − 𝑟 + 1
𝐑
load and generation conditions. ( )
̂ 𝑓 𝜉 = Hankel 𝑅̂ 𝑓 𝜉 (𝜏) , 𝑟, 𝜏1 , 𝜏2 , 𝑓 ∈ {𝑦, 𝑢, 𝑣}
𝐑
( )
3. Linear dynamic model estimation Hankel ℎ (𝜏) , 𝜏0 , 𝜏1 , 𝜏2
( ) ( ) ( )
⎡ ℎ 𝜏1 ℎ 𝜏1 + 1 ⋯ ℎ 𝜏2 − 𝜏0 + 1 ⎤
3.1. Formulation of covariance data equation ⎢ ( ) ( ) ( )⎥
⎢ ℎ 𝜏1 + 1 ℎ 𝜏1 + 2 ⋯ ℎ 𝜏2 − 𝜏0 + 2 ⎥
=⎢ ⎥
In order to obtain 𝐮0 (𝑡), 𝐮(𝑡), 𝑦(𝑡) and P𝑃 𝐶𝐶 (𝑡), Q𝑃 𝐶𝐶 (𝑡) suitable ⎢ ⋮ ⋮ ⋮ ⎥
⎢ ( ) ( ) ( ) ⎥
for modeling the dynamic effects of the DERs, the DERs should be ⎣ℎ 𝜏1 + 𝜏0 − 1 ℎ 𝜏1 + 𝜏0 ⋯ ℎ 𝜏2 ⎦
perturbed around their working equilibrium by small power excitation
and the Toeplitz matrix 𝛹 includes the unknown Markov parameters to
signals without affecting the stability of the microgrid. To guarantee
be identified. Note that the extended observability matrix O𝑟 consists
identifiability of the linear MIMO dynamic model part, the DER power
of the information of the system matrices {𝐴, 𝐶}.
perturbations from different inputs should ideally be uncorrelated. In
practice, uncorrelated input excitation signals can be generated by
3.2. Estimation of matrices 𝐴 and 𝐶 via covariance data
pseudo-random binary signals (PRBS) added to the nominal power level
for each DER set point. Since typically only small power perturbations
̂ 𝑣𝜉 will be asymptotically
As mentioned in [31], the noise effect 𝐑
are allowed, the signal-to-noise ratio (SNR) of the measurement data
may be relatively low due to existence of measurement noise and/or eliminated. In addition, one can easily obtain a shifted version of
possible random load changes. the covariance data equation in (5). The two modified equations are
The CoBRA method serves as a powerful tool to extract a low-order rewritten as follows:
deterministic model from noise contaminated data by first averaging ̂ 𝑦𝜉 = O𝑟 𝐑
𝐑 ̂ 𝑥𝜉 + 𝛹 𝐑
̂ 𝑢𝜉 (6)
the time-dependent covariance between the data and represent corre-
lations as a time-dependent function. To explain the CoBRA method,
̂ → = O𝑟 𝐴𝐑
𝐑 ̂ 𝑥𝜉 + 𝛹
⃗𝐑̂ 𝑢𝜉 (7)
let the dynamic relation between the input/output signals 𝐮(𝑡) and 𝐲(𝑡) 𝑦𝜉
in Fig. 2 be described by a MIMO linear time-invariant minimal system ( )
where 𝐑̂ = Hankel 𝑅̂ 𝑦𝜉 (𝜏) , 𝑟, 𝜏1 + 1, 𝜏2 + 1 is the one-step time-shifted
→
𝑦𝜉
𝐺(𝑞 −1 ) = 𝐶(𝑞𝐼 − 𝐴)−1 𝐵 + 𝐷 with the state space model
{ version of 𝐑̂ 𝑦𝜉 and 𝛹⃗ is the corresponding shifted Toeplitz matrix. As
𝐱 (𝑡 + 1) = 𝐴𝐱 (𝑡) + 𝐵𝐮 (𝑡) will be illustrated in the following derivation, there is no need to pay
(2) ⃗.
𝐲 (𝑡) = 𝐶𝐱 (𝑡) + 𝐷𝐮 (𝑡) + 𝐯 (𝑡) much attention to the detailed structures of both 𝛹 and 𝛹

3
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

With 𝐈 defined as an identity matrix, an orthogonal projection a linear regression optimization based on either the covariance data
[ ]−1 {𝑅̂ 𝑢𝜉 (𝜏), 𝑅̂ 𝑦𝜉 (𝜏)} or the time-domain data {𝑢(𝑡), 𝑦(𝑡)}. In addition, the
̂T 𝐑
𝛱 ⊥̂T = 𝐈 − 𝐑 ̂T
̂ 𝑢𝜉 𝐑 ̂ 𝑢𝜉
𝐑 (8)
𝐑𝑢𝜉 𝑢𝜉 𝑢𝜉 identified model is expected to be applicable for a wide range of
working power levels although the data for identification is obtained
⃗.
is introduced to remove the part containing the information 𝛹 and 𝛹 by exciting the system with small perturbations around a given equi-
Thus, (6) and (7) can be further reduced to librium. Thus, a priori information of the system should be taken into
̂ 𝑥𝜉 𝛱 ⊥
̂ 𝑦𝜉 𝛱 ⊥ = O𝑟 𝐑 account to guarantee a reliable model structure.
𝐑 ̂T ̂T
(9)
𝐑𝑢𝜉 𝐑𝑢𝜉 From the physical perspective, the DERs have unit static gains since
the output active and reactive power will always be adjusted to be
̂ → 𝛱 ⊥ = O𝑟 𝐴𝐑
𝐑 ̂ 𝑥𝜉 𝛱 ⊥ (10)
𝑦𝜉 ̂ T𝐑𝑢𝜉 ̂T 𝐑𝑢𝜉 equal to the set points. For this reason, the unit static gain will be
enforced during the estimation of 𝐵 and 𝐷. On the other hand, the
The matrix 𝐴 can be approximated by solving the optimization problem
mean values will usually be removed beforehand from the covariance
data {𝑅̂ 𝑢𝜉 (𝜏), 𝑅̂ 𝑦𝜉 (𝜏)} to focus on the dynamics. Thus, the raw data
‖ ‖
‖ ̂ →𝛱⊥ ‖
̂ 𝑥𝜉 𝛱 ⊥ − 𝐑 {𝑢(𝑡), 𝑦(𝑡)} will be a better choice for the optimization with a gain
min ‖O𝑟 𝐴𝐑 ̂T ‖
(11)
𝐴 ‖ 𝐑̂T 𝑦𝜉 𝐑 ‖
‖ 𝑢𝜉 𝑢𝜉 ‖𝐹 constraint. The constrained optimization to be solved is then given as
̂ 𝑥𝜉 𝛱 ⊥ are required. Directly solving follows
where an estimation of O𝑟 and 𝐑 ̂T ( )
𝐑𝑢𝜉 ‖ ⎡vec 𝑥0 ⎤ ‖
(11) is typically impossible. However, as mentioned in Section 3.1, ‖ ‖
‖ ‖
min ‖𝐘 − 𝐗 ⎢ vec (𝐵) ⎥ ‖
the estimated model parameters {𝐴, ̂ 𝐵,
̂ 𝐶,
̂ 𝐷}
̂ are just required to be 𝑥0 ,𝐵,𝐷 ‖ ⎢ ⎥ ‖
‖ ⎣ vec (𝐷) ⎦ ‖ (15)
equivalent to the true ones in the sense of being within a similarity ‖ ‖2
transformation. The reason is that both of them are equivalent from s.𝑡. ̂ 𝑛 − 𝐴)
[𝐶(𝐈 ̂ −1 𝐵 + 𝐷]
(𝑘𝑟 ,𝑘𝑐 ) = Gain(𝑘𝑟 ,𝑘𝑐 )
the perspective of the input–output description. Their mathematical
relation is where (𝑘𝑟 , 𝑘𝑐 ) represents the index of a matrix element in row 𝑘𝑟 and
column 𝑘𝑐 , Gain(𝑘𝑟 ,𝑘𝑐 ) is the corresponding predetermined static gain
̂ −1 , 𝐵 = 𝑀 𝐵,
𝐴 = 𝑀 𝐴𝑀 ̂ −1 , 𝐷 = 𝐷
̂ 𝐶 = 𝐶𝑀 from the 𝑘𝑐 th input to the 𝑘𝑟 th output, 𝐴̂ and 𝐶̂ are estimated from
Section 3.2, and
where 𝑀 is a certain nonsingular matrix. Thus, it is sufficient to solve
the following relaxed version of (11) instead ⎡ 𝑦 (0) ⎤ ⎡ 𝜑 (0) ⎤
⎢ ⎥ ⎢ ⎥
‖ ‖ ⎢ 𝑦 (1) ⎥ 𝜑 (1) ⎥
‖̂ ̂̂
min ‖O ⊥ ̂ →𝛱⊥ ‖ 𝐘= ,𝐗 = ⎢
̂T ‖
𝑟 𝐴𝐑𝑥𝜉̂ 𝛱𝐑 −𝐑 (12) ⎢ ⋮ ⎥ ⎢ ⋮ ⎥
𝐴̂ ‖ ‖
̂T 𝑦𝜉 𝐑
‖ 𝑢𝜉 𝑢𝜉 ‖𝐹 ⎢𝑦 (𝑁 − 1)⎥ ⎢𝜑 (𝑁 − 1)⎥
⎣ ⎦ ⎣ ⎦
where Ô = O𝑟 𝑀 and 𝐱̂ (𝑡) = 𝑀 −1 𝐱 (𝑡). One appropriate selection of the [ ∑ ]
{ 𝑟 } 𝜑 (𝑘) = 𝐶̂ 𝐴̂ 𝑘 𝑘−1 T
𝑢(𝑖) ⊗ 𝐶̂ 𝐴̂ 𝑘−𝑖−1 𝑢(𝑘)T ⊗ 𝐈𝑚
pair Ô 𝑟, 𝐑
̂ 𝑥𝜉
̂ 𝛱
⊥ is obtained through the singular value decomposi- 𝑖=0
̂T
𝐑𝑢𝜉
tion of (9) Note that there is no need to specify the static gains for all the
[ ] [ T] subsystems. In our case, only the six diagonal elements of the transfer
[ ] 0 𝑉1
̂ 𝑦𝜉 𝛱 ⊥ = 𝑈1 , 𝑈2 𝛴1
𝐑 function matrix is required to be constrained, i.e.,
𝐑̂T 0 𝛴2 𝑉 T
𝑢𝜉
2
1∕2 1∕2 (13) Gain(𝑗,𝑗) = 1, 𝑗 ∈ {1, 2, … , 6} (16)
≈ 𝑈1 𝛴1 𝛴1 𝑉1T
⏟⏟⏟ ⏟⏟⏟
̂𝑟
O ̂ 𝑥𝜉
𝐑 ⊥ 4. Modeling of the static nonlinearity
̂ 𝛱 T̂
𝐑
𝑢𝜉

where 𝛴1 ∈ R𝑛×𝑛 contains the 𝑛 dominant singular values and thus can 4.1. Admittance matrix of the microgrid network
be used to determine the system order 𝑛. The matrix 𝐴̂ can be computed
by solving (12) and the matrix 𝐶̂ is directly given as the first 𝑚 rows Following the conventional definition of the (positive sequence)
𝑣 𝑖
of O ̂ 𝑟. voltage phasor 𝐕 = 𝑉 𝑟𝑚𝑠 𝑒𝑗𝜙 and current phasor 𝐈 = 𝐼 𝑟𝑚𝑠 𝑒𝑗𝜙 , along
Even if the actual dynamic system is stable, the estimated matrix 𝐴 with the (positive sequence) apparent power
may be unstable due to poor SNR conditions of the experimental data. 𝑣 −𝜙𝑖
𝐒 = 𝐕𝐈∗ = 𝑉 𝑟𝑚𝑠 𝐼 𝑟𝑚𝑠 𝑒𝑗 (𝜙 ) = 𝑃 + 𝑗𝑄 (17)
A stable approximation of matrix 𝐴 can be guaranteed by solving a
modified version of (12) with the pole location constraints where 𝑟𝑚𝑠 denotes the root mean square of a signal and is omitted in
‖ ( )† ‖ the following discussions for notational simplicity. The static coupling
‖̂ ̂ →𝛱⊥ 𝐑 ̂ 𝑥𝜉 ‖
min𝑄,𝑃 ‖O 𝑄−𝐑 𝛱 ⊥̂T 𝑃‖ between phasors and power loss at the PCC can be characterized by
‖ 𝑟 𝑦𝜉 ̂ T
𝐑𝑢𝜉 ̂ 𝐑𝑢𝜉 ‖
‖ ‖𝐹 the admittance matrix of the microgrid network. In particular for the
s.𝑡. 𝑎 ⊗ 𝑃 + 𝑏 ⊗ 𝑄 + 𝑏𝑇 ⊗ 𝑄 T ≥ 0 (14) microgrid framework of Fig. 2, the relation between the phasors is
𝑃 = 𝑃T > 0 given by
T𝑟 (𝑃 ) = 𝑛 ⎡ 𝐕𝑏 ⎤ ⎡ 𝐈𝑏 ⎤
[ ] ⎢ ⎥ ⎢ ⎥
𝐕 𝐈
where 𝑎 = 1 − 𝛿, 𝛿 ∈ [0, 1], 𝑏 =
0 1
, and (⋅)† represents the right 𝐌⎢ 𝑑 ⎥ = ⎢ 𝑑 ⎥ (18)
0 0 ⎢ 𝐕𝑠 ⎥ ⎢ 𝐈𝑠 ⎥
⎢𝐕 ⎥ ⎢ ⎥
pseudo-inverse. Denote the optimum of (14) as {𝑄∗ , 𝑃 ∗ }. Then, matrix ⎣ 𝑃 𝐶𝐶 ⎦ ⎣𝐈𝑃 𝐶𝐶 ⎦
𝐴̂ is estimated by 𝐴 = 𝑄∗ 𝑃 ∗−1 and the eigenvalues of 𝐴̂ are enforced
where 𝐈𝑏 , 𝐈𝑑 , and 𝐈𝑠 are the currents going through the BESS, diesel, and
to be within the disk with radius (1 − 𝛿) in the complex plain. One can
solar/PV respectively, 𝐈𝑃 𝐶𝐶 = −𝐈𝑏 − 𝐈𝑑 − 𝐈𝑠 . The admittance matrix 𝐌
find other different pole location constraints in [34].
in (18) is defined by

3.3. Estimation of matrices 𝐵 and 𝐷 ⎡ 𝑌𝑏 0 0 −𝑌𝑏 ⎤

⎢ ⎥
0 𝑌𝑑 0 −𝑌𝑑
𝐌=⎢ ⎥ (19)
Once the matrices 𝐴 and 𝐶 have been estimated through the ex- ⎢ 0 0 𝑌𝑠 −𝑌𝑠 ⎥
tended observability matrix, the estimation of 𝐵 and 𝐷 is about solving ⎢−𝑌 −𝑌𝑑 −𝑌𝑠 𝑌𝑏 + 𝑌𝑑 + 𝑌𝑠 ⎥⎦
⎣ 𝑏

4
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

Table 1
Nominal values for five scenarios.
Input Initial 1 (0–60 s) 2 (60–70 s) 3 (70–80 s) 4 (80–90 s) 5 (90–100 s)
Pin
b
0 kW 0 kW –20 kW –20 kW –20 kW –20 kW
Qin
b
0 kVAR 0 kVAR 0 kVAR –5 kVAR –5 kVAR –5 kVAR
Pin
d
0 kW 0 kW 0 kW 0 kW 40 kW 40 kW
Qin
d
0 kVAR 0 kVAR 0 kVAR 0 kVAR 0 kVAR 10 kVAR
Pin
s 0 kW 30 kW 30 kW 30 kW 30 kW 30 kW
Qin
s
0 kVAR 0 kVAR 0 kVAR 0 kVAR 0 kVAR 0 kVAR
Output Initial
PPCC 18.3 kW
QPCC 3.45 kVAR

and it should be noted that only the power flows between the main grid the resistance. Thus, 𝑌𝑘 ≈ 𝑅1 by neglecting the imaginary part and (22)
𝑘
and DERs are considered since we are only interested in the dynamic can be further reduced to
modeling from the DERs to the PCC. ⎧( 𝑟𝑒 )2 ( )2
The admittance matrix 𝐌 of the network is dependent on the ⎪ 𝑉𝑘 − 𝑉𝑘𝑟𝑒 𝑉𝑃𝑟𝑒𝐶𝐶 + 𝑉𝑘𝑖𝑚 − 𝑉𝑘𝑖𝑚 𝑉𝑃𝑖𝑚𝐶𝐶 = 𝑅𝑘 𝑃𝑘
⎨ 𝑟𝑒 𝑖𝑚 (24)
structure of the local distribution lines and can be approximated by ⎪𝑉𝑘 𝑉𝑃 𝐶𝐶 − 𝑉𝑘𝑖𝑚 𝑉𝑃𝑟𝑒𝐶𝐶 = 𝑅𝑘 𝑄𝑘
using (18). Only the admittance between PCC and each DER is taken ⎩
into account here. A more general case also considers the admittance The corresponding currents 𝐈𝑘 = 𝐼𝑘𝑟𝑒 + 𝑗𝐼𝑘𝑖𝑚 can be computed from (18)
between the DERs, which results in a similar admittance matrix. How- once (𝑉𝑘𝑟𝑒 , 𝑉𝑘𝑖𝑚 ) is obtained. The active power loss going through each
ever, for the structure (19) used in this paper, each line admittance power channel is then 𝑃𝑘𝑙𝑜𝑠𝑠 = [(𝐼𝑘𝑟𝑒 )2 + (𝐼𝑘𝑖𝑚 )2 ]𝑅𝑘 . The reactive power
can be simply estimated separately instead of computing the whole loss is 𝑄𝑙𝑜𝑠𝑠
𝑘
= [(𝐼𝑘𝑟𝑒 )2 + (𝐼𝑘𝑖𝑚 )2 ]𝑋𝑘 . Note that 𝑄𝑙𝑜𝑠𝑠
𝑘
can be neglected for
admittance matrix as the case of small inductance in the distribution lines and then 𝑄𝑙𝑜𝑠𝑠 = 0.
∑ 𝑘
I𝑘 The total active power loss reflected at the PCC is Ploss = 𝑘 𝑃𝑘𝑙𝑜𝑠𝑠 while
𝑌𝑘 = , 𝑘 = 𝑏, 𝑑, 𝑠 (20) ∑ PCC
V𝑘 − V𝑃 𝐶𝐶 the total reactive power loss Qloss PCC
= 𝑘 𝑄𝑙𝑜𝑠𝑠
𝑘
(= 0).
Denote the initial active and reactive power levels at the PCC as
4.2. Estimation of the power loss Pinit
PCC
and Qinit
PCC
when the three channels are neither producing or
consuming power. Then, the active and reactive power at the PCC
shown in Fig. 2 considering the power loss are estimated by
The distribution line with admittance 𝑌𝑘 (𝑘 = 𝑏, 𝑑, 𝑠) is modeled as
∑
a serial connection of a resistor 𝑅𝑘 and an inductor 𝐿𝑘 , where 𝑌𝑘 = PPCC = Pinit
PCC
− 𝑃𝑘 + Ploss
PCC
(25)
1
𝑅𝑘 +𝑗2𝜋𝑓 𝐿𝑘
. The reactance 𝑋𝑘 = 2𝜋𝑓 𝐿𝑘 and 𝑓 is the fundamental AC 𝑘
power frequency (e.g. 50 or 60 Hz). The total power loss is computed ∑
QPCC = Qinit
PCC
− 𝑄𝑘 + Qloss
PCC
(26)
as 𝐈2𝑘 (𝑅𝑘 + 𝑗𝑋𝑘 ), which is proportional to the through current squared. 𝑘
Since the main grid is working as the voltage source while in grid-
connected mode, one can denote the nominal voltage at the PCC as 5. Validation of modeling algorithm

For validation of the proposed modeling approach, a detailed

𝑣
𝐕𝑃 𝐶𝐶 = 𝑉𝑃 𝐶𝐶 𝑒𝑗𝜙𝑃 𝐶𝐶 = 𝑉𝑃𝑟𝑒𝐶𝐶 + 𝑗𝑉𝑃𝑖𝑚𝐶𝐶 (21) component-based microgrid model is constructed using the toolbox
Simscape Power SystemsTM in MATLAB/Simulink. The resulting Sim-
Similar notations can be defined for 𝑌𝑘 , 𝑉𝑘 , 𝑘 = 𝑏, 𝑑, 𝑠. Then, the
scape model can be used to feed (noisy) data to the proposed estimation
following nonlinear equations can be obtained by combining (17), (18), techniques and compare (noise free) Simscape simulations with the
(20), and (21) simulations created by the estimated model.
⎧ 𝑟𝑒 ( 𝑟𝑒 )2 ( 𝑟𝑒 𝑟𝑒 ) The Simscape model created for this validation study consists of a
⎪𝑌𝑘 𝑉𝑘 − 𝑌𝑘 𝑉𝑃 𝐶𝐶 − 𝑌𝑘𝑖𝑚 𝑉𝑃𝑖𝑚𝐶𝐶 𝑉𝑘𝑟𝑒
hybrid collection of DERs interconnected with distribution lines and
⎪ 𝑟𝑒
( 𝑖𝑚 )2 ( 𝑟𝑒 𝑖𝑚 )
⎪ +𝑌𝑘 𝑉𝑘 − 𝑌𝑘 𝑉𝑃 𝐶𝐶 + 𝑌𝑘𝑖𝑚 𝑉𝑃𝑟𝑒𝐶𝐶 𝑉𝑘𝑖𝑚 = 𝑃𝑘 loads as shown in Fig. 1. The distribution lines are all modeled as
⎨ 𝑖𝑚 ( 𝑟𝑒 )2 ( 𝑟𝑒 𝑖𝑚 ) (22) a serial concatenation of a 1 Ω resistor and 10𝑒−5 H inductor. The
⎪−𝑌𝑘 𝑉𝑘 + 𝑌𝑘 𝑉𝑃 𝐶𝐶 + 𝑌𝑘𝑖𝑚 𝑉𝑃𝑟𝑒𝐶𝐶 𝑉𝑘𝑟𝑒
⎪ corresponding impedance can be approximated as 1 Ω by neglecting the
( )2 ( )
⎪ −𝑌𝑘𝑖𝑚 𝑉𝑘𝑖𝑚 − 𝑌𝑘𝑟𝑒 𝑉𝑃𝑟𝑒𝐶𝐶 − 𝑌𝑘𝑖𝑚 𝑉𝑃𝑖𝑚𝐶𝐶 𝑉𝑘𝑖𝑚 = 𝑄𝑘 inductance. Thus, the admittance 𝑌𝑘 ≈ 1 S for 𝑘 = 𝑏, 𝑑, 𝑠, 𝑙𝑜𝑎𝑑 and the
⎩
reduced version (24) can be used. Note that Qloss
PCC
= 0 in the simplified
where 𝑘 = 𝑏, 𝑑, 𝑠 to represent the battery, diesel, and solar channel case. The nominal voltage at the PCC with the phasor form in (21) is
respectively. Thus, there are actually 3 pairs of nonlinear equations in set as 𝑉𝑃 𝐶𝐶 = 480 V and 𝜙𝑣𝑃 𝐶𝐶 = 10 degree with nominal frequency
total to solve 𝐕𝑘 = 𝑉𝑘𝑟𝑒 + 𝑗𝑉𝑘𝑖𝑚 , 𝑘 = 𝑏, 𝑑, 𝑠. The practical solution can be 60 Hz. Note that, in practice, the frequency within the microgrid will
obtained directly by solving (22) via nonlinear optimization with the not be constant but oscillating around the nominal value regulated by
initial value as (𝑉𝑃𝑟𝑒𝐶𝐶 , 𝑉𝑃𝑖𝑚𝐶𝐶 ). The values of 𝑃𝑘 , 𝑄𝑘 can be obtained from the voltage source. The magnitude and angle of the PCC voltage also
the outputs of the linear DERs system identified using the techniques show a similar behavior. We will still use the nominal PCC values to
in Section 3. In terms of Fig. 2, solve (24) and the estimation result will show that the small oscillation
is not a problem. The load is modeled as the parallel RLC load with 20
⎧ out out
⎪𝑃𝑏 = P𝑏 , 𝑄𝑏 = Q𝑏 kW nominal active power, 4 kVAR nominal inductive reactive power,
⎪ out out and −100 VAR nominal capacitive reactive power. In addition, the rate
⎨𝑃𝑑 = P𝑑 , 𝑄𝑑 = Q𝑑 (23)
⎪ limiter is not considered in the simulation example for easy discussion.
⎪𝑃𝑠 = Pout out
𝑠 , 𝑄𝑠 = Q𝑠
However, the rate limiter will be taken into account in the next section
⎩
where the practical data is used.
It is worthwhile to mention that the problem can be further simpli- Five working scenarios, with 100 s total simulation time and cor-
fied by ignoring the wire inductance if it is much smaller compared to responding to different power levels, are considered and described

5
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

Fig. 3. Validation by comparison of Simscape (true) and linear model equivalent Fig. 6. Validation by comparison of Simscape (true) and linear model equivalent via
via CoBRA without static network nonlinearity to model line losses (uncompensated). N4SID with static network nonlinearity to model line losses (compensated). The first
The first 50 s of data belongs to scenario 1 and is used for identification only. The 50 s of data belongs to scenario 1 and is used for identification only. The vertical
vertical black dashed line separates the time window into the five scenarios of Table 1 black dashed line separates the time window into the five scenarios of Table 1. Obvious
and divergence between true and simulated outputs can be observed during different divergence between true and simulated outputs can be observed only when the nominal
scenarios. power level moves far away from where the linear model equivalent was estimated.
The estimated power loss is also shown.

Fig. 4. Validation by comparison of Simscape (true) and linear model equivalent via
CoBRA with static network nonlinearity to model line losses (compensated). The first
Fig. 7. Comparison of power loss estimation using the CoBRA and N4SID estimated
50 s of data belongs to scenario 1 and is used for identification only. The vertical
linear model equivalents. Obvious discrepancies can be observed when the nominal
black dashed line separates the time window into the five scenarios of Table 1 and
power level moves far away from where the linear model equivalent was estimated.
no divergence between true and simulated outputs can be observed during different
scenarios. The estimated power loss is also shown.

Fig. 8. Schematic overview of Local Area Distribution Network (LADN) for

experimental validation of the physics-informed data-driven modeling approach.

are assumed to be constant during the whole simulation. There will

be only one DER power level step change between two scenarios. For
Fig. 5. Validation by comparison of Simscape (true) and linear model equivalent
via N4SID without static network nonlinearity to model line losses (uncompensated). example, the nominal active power level of the battery inverter system
The first 50 s of data belongs to scenario 1 and is used for identification only. The will drop down to −20 kW from 0 kW at 𝑡 = 60 s while the other set
vertical black dashed line separates the time window into the five scenarios of Table 1 points are not changed.
and divergence between true and simulated outputs can be observed during different
scenarios.
Six uncorrelated Pseudo Random Binary Sequence (PRBS) dispatch
signals with magnitude 3 kW/kVAR are added around the nominal set
points to excite the system. In practice, a PRBS can be generated by
in Table 1. The initial power levels at PCC are Pinit = 18.3 kW and randomly switching the active and reactive dispatch values up/down.
PCC
Qinit = 3.45 kVAR. The nominal power set points of the solar inverter As long as the switching at each DER occurs independently, the dis-
PCC
patch signals will be uncorrelated. The data of the first 50 s under the

6
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

Fig. 9. Determine the rate limiter from the input/output data by calculating the slope between the two selected points denoted by ‘+’.

Fig. 10. The comparison between the model prediction and the true power output at the POI. Compensating the power loss in the model leads to a more accurate prediction.

scenario 1 is used for identification by using the techniques proposed similar to (2). The identified model by using the proposed method
in Sections 3 and 4. The remaining data including all five scenarios shows high accuracy under not only scenario 1 but also the other
is used for validation. The measurements are sampled at 60 Hz and scenarios. This is desirable in practice since one would like the model
contaminated by Gaussian white noise with SNR = 25 dB. identified around one working power level to be applicable for other
To illustrate the need to model distribution line loss, the CoBRA levels.
method is first used to directly identify a linear model from the DER To further illustrate the efficacy of the proposed CoBRA method, the
set points to the output power PPCC and QPCC without a static network N4SID method [16,20,22] is performed on the same simulation data to
nonlinearity. The corresponding results are shown in Fig. 3 and show identify a 14th order state space model for comparison. Similarly, the
that the linear model output is modeling the data of scenario 1 very results without and with considering a static network nonlinearity are
well, but diverges when the nominal power levels are different from shown in Figs. 5 and 6 individually. The power outputs predicted by
scenario 1, as marked by the black dashed lines in the figure at the linear model equivalent via N4SID also diverge when the nominal
specific times. Clearly, the nonlinearity caused by the power loss due power levels move away from the scenario 1 without considering the
to the distribution line impedance should be taken into account. The line losses. Taking the line losses into account may help obtain a better
improvements when also modeling line losses with a static network power prediction when the nominal power levels are not far from
non-linearity are summarized in Fig. 4. Here it can be seen that the scenario 1. However, divergence is observed for scenario 4 and 5. There
dynamic equivalent captures the dynamics and nonlinearity of the are two sources of errors. On one hand, the N4SID method does not
power flow over all scenarios. enforce physical constraints on static gains as (16) in the proposed
The modeling results of Fig. 4 are obtained by choosing the param- CoBRA method. The gain errors during identifying the linear model
eters in the covariance data Eq. (5) as 𝑟 = 30, 𝜏1 = 0, 𝜏2 = 150 and equivalent lead to large prediction errors of the output power for DERs
identifying the linear DERs system as a 14th order state space model when moving to a much different nominal power level. On the other

7
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

Fig. 11. The identified model by CoBRA is able to capture the coupling between Pout
b
and Qout
b
.

hand, the prediction errors in the identified linear DERs system will the apparent power 𝐒 in (17) are used, only cos(⋅) and sin(⋅) of the angle
further propagate to (23) for solving (22) and the power loss will not difference 𝜙𝑣𝑃 𝑂𝐼 − 𝜙𝑖𝑃 𝑂𝐼 are relevant. All PMU data is sampled at 30 H𝑧
be approximated accurately. The comparison of the active power loss and data is collected over a time window of 1 h long. For estimation
estimation is illustrated in Fig. 7. and modeling purposes, only the first 11 min are used as training data
and the rest for validation.
6. Experimental results It is obvious that the network under test is similar to the one in
Fig. 1 and the proposed method can be implemented. Similar to (1),
For further validation of the applicability of the proposed data- the input/output signals for system identification is modified as
driven modeling approach, the estimation algorithms are also applied ⎡ Pin
𝑏 (𝑡) ⎤ ⎡ Pout
𝑏 (𝑡)
⎤
to experimental data from a local area distribution network (LADN). ⎢ in ⎥ ⎢ out ⎥
⎢Q𝑏 (𝑡)⎥ ⎢Q𝑏 (𝑡)⎥
The LADN operates as a grid-connected microgrid where DERs are 𝐮0 (𝑡) = ⎢ in ⎥ , 𝐲(𝑡) = ⎢ out ⎥ (27)
placed behind local distribution transformer. A schematic overview of P
⎢ 𝑠 ⎥ (𝑡) ⎢ P𝑠 (𝑡) ⎥
⎢ in ⎥ ⎢ out ⎥
the LADN with its Point of Interconnect (POI) or PCC is shown in Fig. 8, ⎣Q𝑠 (𝑡)⎦ ⎣Q𝑠 (𝑡)⎦
illustrating the presence of a BESS DER and a PV DER each behind The initial active and reactive power at the POI are computed as
their own distribution transformers. In addition, a non-controllable Pinit = −Pout and Qinit = −Qout . The admittance 𝑌𝑏 , 𝑌𝑠 , and 𝑌𝑙𝑜𝑎𝑑 can be
POI load POI load
but measurable load is present in the LADN behind a distribution determined by (18). As shown in Fig. 9, the rate limiter for BESS and
transformer. The data of this LADN is chosen to illustrate that presence PV DERs can be determined from the slope of two selected points in the
of local distribution transformer can be taken into account in the power output when there is a high rate input change. The intermediate
microgrid dynamic model by estimating the dynamics and power loss signal 𝐮(𝑡) illustrated in Fig. 2 can be then computed. An 11th order
as part of the dynamic model. linear state space model from 𝐮(𝑡) to 𝐲(𝑡) is identified by CoBRA. This
The BESS DER is rated at 250 kW and the PV DER at 500 kW. model can be used to predict 𝐲(𝑡) given any 𝐮(𝑡) without using PMUs.
The BESS DER is excited by random step perturbations on the active The active and reactive power losses, denoted as Ploss and Qloss , are
∑ 𝑙𝑜𝑠𝑠 loss ∑ 𝑙𝑜𝑠𝑠 POI POI
and reactive power, while the PV DER is excited with slowly changing calculated as Ploss
POI
= 𝑃
𝑘 𝑘 , Q POI
= 𝑄
𝑘 𝑘 , where 𝑘 = 𝑏, 𝑠, 𝑙𝑜𝑎𝑑 for
random power variations, curtailing solar power production. In terms this test and (22) will be solved in the process. Similar to (25) and (26),
of notation, these excitation signals correspond to Pin b
, Qin
b
, Pin in
s , Qs as the active and reactive power at the POI can be estimated as
discussed earlier, while there is no diesel generator in the system under
PPOI = Pinit
POI
− (Pout
b
+ Pout loss
s ) + PPOI (28)
test. The network is connected to a main load that can be up to 1.5
MW. The BESS, PV, and load power outputs are measured on bus 2,3,4
respectively and can be denoted as Pout , Qout , Pout out out out QPOI = Qinit − (Qout + Qout loss
s ) + QPOI (29)
b b s , Qs , Pload , and Qload POI b
accordingly. The admittance of the distribution line between bus 1 and The comparison between the model prediction from simulation and
bus 2,3,4 is then denoted as 𝑌𝑏 , 𝑌𝑠 , and 𝑌𝑙𝑜𝑎𝑑 . the true power output at the POI is illustrated in Fig. 10. The model
𝑣
The nominal 𝑉𝑃 𝑂𝐼 = 12 kV in the voltage phasor 𝐕𝑃 𝑂𝐼 = 𝑉𝑃 𝑂𝐼 𝑒𝑗𝜙𝑃 𝑂𝐼 prediction considering the power loss in the distribution line using
at the POI similar to (21). It is worthwhile to note that the positive the experimental data is accurate, as also demonstrated in the earlier
sequence voltage angle 𝜙𝑣𝑃 𝑂𝐼 at the POI is not constant due to the time simulation example. It is worth noting that the CoBRA method can
variations in the AC grid frequency to which the LADN is connected via accurately capture the coupling between different channels as shown
the POI. However, since the active and reactive power components of in Fig. 11.

8
Y. Hu and R.A. de Callafon International Journal of Electrical Power and Energy Systems 153 (2023) 109279

7. Conclusion [10] Popov V, Prykhno V, Prykhno D. Development of the method of determining the
power and electricity losses in distribution network of shop electrical supply. In:
IEEE 6th international conference on energy smart systems. ESS, 2019, p. 104–7.
Synchronized power flow measurements obtained via PMUs with
http://dx.doi.org/10.1109/ESS.2019.8764231.
small simultaneous power perturbations of DERs are used to create [11] Resende FO, Lopes JAP. Development of dynamic equivalents for MicroGrids
a dynamic equivalent of the active and reactive power flow at the using system identification theory. In: 2007 IEEE lausanne power tech. 2007, p.
Point of Common Coupling (PCC) of a microgrid. A model structure 1033–8. http://dx.doi.org/10.1109/PCT.2007.4538457.
combining a linear dynamic system and static nonlinearity is proposed [12] Hu Y, Konakalla SAR, de Callafon RA. Covariance based estimation for reduced
order models of microgrid power flow dynamics. In: The 18th IFAC symposium
to model the multi-input, multi-output (MIMO) active and reactive on system identification–SYSID2018, Vol. 51. IFAC; 2018, p. 903–8.
power flow dynamics of a microgrid with multiple DERs. The linear dy- [13] Shi Y, Hu S, Xua D, Su J, Liu N, Yang X, Dua Y. A simplified microgrid
namic system is identified via a covariance based realization algorithm voltage and frequency response characteristic modelling method based on system
(CoBRA), while the nonlinear part is construed by the estimation of the identification. Int J Electr Power Energy Syst 2020;121:106063.
[14] Conte F, D’Agostino F, Massucco S, Silvestro F, Bossi C, Cabiati M. Experimental
admittance matrix of the power network of the microgrid. Both simula-
validation of a dynamic equivalent model for microgrids. IEEE Trans Ind Appl
tion and experimental results demonstrate an excellent match between 2021;57:2202–11.
measured and simulated power flow, making the resulting dynamic [15] Escobar ED, Manrique T, Isaac IA. Campus microgrid data-driven model
equivalent suitable for not only the current working power levels but identification and secondary voltage control. Energies 2022;15:7846.
also the other operating conditions. With the dynamic model tested and [16] Ljung L. System identification: theory for the user. 2nd ed.. Upper Saddle River,
NJ: Prentice-Hall; 1999.
validated on data, a microgrid power control can dynamically adjust set [17] Verhaegen M, Verdult V. Filtering and system identification: a least squares
point values of DERs to obtain controlled power flow at the PCC while approach. Cambridge University Press; 2007.
taking into account the power loss within the distribution system of [18] Chakraborty R, Jain H, Seo G-S. A review of active probing-based system
the microgrid. Future studies include but are not limited to taking into identification techniques with applications in power systems. Int J Electr Power
Energy Syst 2022;140:108008.
account more physical constraints and investigating the efficacy of the
[19] Zhou N, Pierre JW, Hauer JF. Initial results in power system identification
proposed physics-informed data-driven approach in microgrid power from injected probing signals using a subspace method. IEEE Trans Power Syst
systems controller design. 2006;21(3):1296–302.
[20] Xu X, Li K, Jia H, Yu X, Deng J, Mu Y. Data-driven dynamic modeling of
Declaration of competing interest coupled thermal and electric outputs of microturbines. IEEE Trans Smart Grid
2016;9(2):1387–96.
[21] Escobar ED, Manrique T, Isaac IA. Campus microgrid data-driven model
We have no conflicts of interest to disclose. identification and secondary voltage control. Energies 2022;15(21):7846.
[22] Van Overschee P, De Moor B. N4SID: Subspace algorithms for the identification
Data availability of combined deterministic-stochastic systems. Automatica 1994;30(1):75–93.
[23] Katayama T. Subspace methods for system identification. Springer Science &
Business Media; 2005.
Data will be made available on request. [24] Schweitzer EO, Whitehead DE. Real-time power system control using synchropha-
sors. In: 2008 61st annual conference for protective relay engineers. IEEE; 2008,
References p. 78–88.
[25] Van Overschee P, De Moor B. Subspace identification for linear systems:
[1] Lasseter RH. Microgrids and distributed generation. J Energy Eng theory—-implementation—-applications. Kluwer Academic Publishers; 1996.
2007;133(3):144–9. [26] Miller DN, de Callafon RA. Efficient identification of input dynam-
[2] Weber LG, Nasiri A, Akbari H. Dynamic modeling and control of a syn- ics for correlation function-based subspace identification. IFAC Proc Vol
chronous generator in an AC microgrid environment. IEEE Trans Ind Appl 2011;44(1):6511–6.
2018;54(5):4833–41. [27] Miller DN, De Callafon RA, Brenner MJ. Covariance-based realization algo-
[3] Khodadadi A, Mohammadian M, Hosseinian SH. Robust power sharing control rithm for the identification of aeroelastic dynamics. J Guid Control Dyn
of an autonomous microgrid featuring secondary controller. In: 2016 smart grids 2012;35(4):1169–77.
conference. SGC, IEEE; 2016, p. 1–7. [28] de Callafon RA, Miller DN. CoBRA-software for a covariance based realiza-
[4] Sheikh A, Youssef T, Mohammed O. AC microgrid control using adaptive tion algorithm with convex constraints on pole locations. IFAC-PapersOnLine
synchronous reference frame PLL. In: 2017 ninth annual IEEE green technologies 2015;48(28):763–8.
conference (GreenTech). IEEE; 2017, p. 46–51. [29] Hu Y, de Callafon RA. Covariance based realization algorithm for identification
[5] Jiang T, Costa LM, Siebert N, Tordjman P. Automated microgrid control systems. of linear time periodic systems. In: 2018 annual american control conference.
CIRED-Open Access Proc J 2017;2017(1):961–4. ACC, IEEE; 2018, p. 1102–7.
[6] Bintoudi AD, Zyglakis L, Apostolos T, Ioannidis D, Al-Agtash S, Martinez- [30] Hu Y, de Callafon RA. Optimal weighting for covariance based realization
Ramos JL, Onen A, Azzopardi B, Hadjidemetriou L, Martensen N, et al. Novel algorithm. In: 2017 IEEE 56th annual conference on decision and control. CDC,
hybrid design for microgrid control. In: 2017 IEEE PES asia-pacific power and IEEE; 2017, p. 5274–9.
energy engineering conference. APPEEC, IEEE; 2017, p. 1–6. [31] Hu Y, Jiang Y, de Callafon RA. Variance reduction in covariance based realization
[7] Gwynn BT, de Callafon RA. Robust real-time inverter-based reactive power algorithm with application to closed-loop data. Automatica 2020;113:108683.
compensation. In: 2019 IEEE 58th conference on decision and control. CDC, [32] Abusara MA, Guerrero JM, Sharkh SM. Line-interactive UPS for microgrids. IEEE
IEEE; 2019, p. 6578–83. Trans Ind Electron 2013;61(3):1292–300.
[8] Giannakis GB, Kekatos V, Gatsis N, Kim S-J, Zhu H, Wollenberg BF. Monitoring [33] D’Arco S, Suul JA. A synchronization controller for grid reconnection of islanded
and optimization for power grids: A signal processing perspective. IEEE Signal virtual synchronous machines. In: 2015 IEEE 6th international symposium on
Process Mag 2013;30(5):107–28. power electronics for distributed generation systems. PEDG, IEEE; 2015, p. 1–8.
[9] Konakalla SAR, Valibeygi A, de Callafon RA. Microgrid dynamic model- [34] Miller DN, De Callafon RA. Subspace identification with eigenvalue constraints.
ing and islanding control with synchrophasor data. IEEE Trans Smart Grid Automatica 2013;49(8):2468–73.
2019;11(1):905–15.