1 Introduction↩︎

With the increasing integration of distributed renewable generation and power-electronics-based technologies, modern power systems are becoming more dynamic. Today’s grids exhibit diverse and variable subsystem interactions, making simple analytical models inadequate [@Wang2019]. This complexity is compounded by limited data and model sharing among stakeholders: device manufacturers typically withhold proprietary models, grid operators may only have coarse or steady-state system models, and consumers can exhibit significant but largely unknown dynamics (e.g., data centers). In this context, impedance identification offers a promising alternative [@DeMeerendre2020], providing data-driven models of local grid dynamics.

Accurate knowledge of the grid impedance is valuable for a wide range of applications. Its value at the fundamental frequency, \(\omega_g\), can be used to estimate voltage stability margins, maximum power transfer limits, grid strength metrics, etc. When characterized over a wide frequency band, it enables harmonic penetration studies and filter design, characterization of inter-area and subsynchronous oscillations, or model-based control design. It also plays a key role in the optimized operation of grid-connected converters [@Harnefors2007a; @cespedes2014adaptive; @Wang2014; @Wang2018]. Interactions between converters and the grid can degrade power quality and trigger instabilities [@Mollerstedt2000; @Liserre2006; @li2017unstable]. Impedance-based analysis has therefore emerged as an effective tool for assessing small-signal stability [@belkhayat1997; @sun2009small; @Rygg2016; @Chen2024], modeling the system’s small-signal dynamics via an equivalent dynamic impedance \(Z_g(s)\). In three-phase systems, this impedance forms a multi-input multi-output (MIMO) transfer function relating small-signal terminal voltages and currents at the Point of Common Coupling (PCC), typically represented in a synchronous \(dq\)-frame.

During the last two decades, numerous approaches for grid impedance identification have been proposed, primarily in the power electronics literature. They can be classified into passive, quasi-passive, and active methods. A detailed description of most methods is given in [@Stiegler2015] and [@DeMeerendre2020]. Passive methods rely on the small harmonic distortion that naturally exists at the PCC. Their applicability is limited to impedance estimation at \(\omega_g\) and a few harmonic frequencies, [@Gu2012; @Hoffmann2014]. They may not be suitable for tasks that rely on wideband characterizations, such as stability analysis and advanced control design. Quasi-passive methods combine a triggering mechanism [@Garcia2014; @Cobreces2009] with an active method to avoid continuous grid perturbation. Active methods deliberately introduce disturbances by repeatedly switching resistive or capacitive loads [@Girgis1989; @Jordan2018] or by injecting small current or voltage signals through dedicated hardware [@francis2011algorithm; @huang2009small; @Rygg2016]. Examples include frequency-sweep techniques [@huang2009small], which apply perturbations sequentially at discrete frequencies, and wideband excitation methods [@Rygg2016], which excite multiple frequencies simultaneously. Both use idealized steady-state Fourier analysis of the measured voltages and currents to obtain a non-parametric frequency-domain model. However, they require specialized hardware and long measurement times, which limits their practicality.

A more practical approach is to use existing grid-connected converters to excite the grid when needed by superimposing an excitation signal on the reference of an inner control loop. Impulse excitations have been proposed [@cespedes2012online; @liu2020analysis], offering very short perturbation times but at the cost of large disturbances that may jeopardize power quality, excite nonlinear behavior, or trigger protection relays. Alternatively, smaller-amplitude signals such as Maximum-Length Binary Sequences (MLBS) have been employed [@martin2013wide; @riccobono2017noninvasive; @luhtala2018implementation; @roinila2017mimo]; however, these methods either assume \(dq\)-symmetry and neglect cross-coupling effects (i.e., a diagonal \(Z_g(s)\)), or require sequential, linearly independent perturbations, resulting in longer measurement times.

Despite considerable progress, accurate grid impedance identification remains challenging due to inherent trade-offs. A key difficulty is minimizing the perturbation time and amplitude while preserving the accuracy of the identification. Recent efforts have sought to address this: for instance, [@Berg2022] proposed a non-parametric frequency-domain approach using orthogonal binary signals as an alternative to sequential perturbation. Although potentially more robust, its measurement time is the same as MLBS-based methods [@martin2013wide; @riccobono2017noninvasive; @luhtala2018implementation; @roinila2017mimo] and still requires periodic steady-state measurements to prevent spectral leakage in discrete Fourier transform (DFT) analysis. Alternatively, [@Haberle2023] demonstrated the use of parametric time-domain techniques, particularly discrete-time Auto-Regressive Exogenous (ARX) models [@ljung1998system], which eliminate the need for specialized excitation signals and can handle non-periodic or transient measurements. However, their performance is rather sensitive to the chosen parameterization and model order.

To address these challenges, we propose an active non-parametric frequency-domain identification method for \(dq\)-asymmetric grids that does not require sequential perturbations or steady-state measurements. This eliminates the measurement-time limitations of existing methods, leaving only constraints imposed by the required frequency resolution and signal-to-noise ratio (SNR). We leverage complex TFs [@Harnefors2007] to parametrize the grid equivalent impedance using single-input single-output (SISO) complex TFs and show how to reconstruct the full real grid impedance TF matrix once an estimate is obtained. This approach simplifies MIMO impedance parameterization, clearly distinguishes between symmetric and asymmetric cases, and provides an algebraically efficient representation. Furthermore, the SISO complex TF non-parametric estimates can be directly used for stability assessment or control design [@Zhang2018; @Harnefors2020]. They may also be converted to parametric models using vector fitting methods [@Ozdemir2017].

The proposed approach relies on frequency-domain local parametric approximations that are well studied in the system identification literature [@Pintelon2010; @McKelvey2012; @Pintelon2021]. We assume that the frequency response of the complex TFs can be accurately approximated over short frequency intervals by low-order continuous-time ARX models. This is a reasonable assumption, in particular, when the number of samples \(N\) ensures that the 3dB-bandwidth of any resonance spans several spectral lines. No finite global order is assumed, making the estimated model truly non-parametric. Local models are fitted by solving \(N\) small and independent linear least-squares problems. As shown in simulations, the identification accuracy is insensitive to the local model order, unlike parametric methods.

In summary, our approach strikes a balance between fully non-parametric methods that require sequential perturbations, and fully parametric methods that attempt capturing the dynamics with a single high-order parametric model, often with insufficient accuracy. It offers ultrafast (\(<\! 1\) s) identification with a favorable accuracy and data efficiency trade-off.

2 Problem formulation↩︎

2.1 Small-signal model↩︎

The objective is to identify the dynamic small-signal Thévenin equivalent impedance of an AC three-wire, three-phase grid. This is achieved using time-domain samples of terminal voltages \(v_{abc}\) and currents \(i_{abc}\) at the PCC of interest; see Figure 1. No assumptions are made about the topology or strength of the grid, which can include generators, loads, and actively controlled power-electronics systems. Under balanced operation, Park’s transformation at the steady-state frequency of the grid \(\omega_g\) maps the three-phase voltages and currents to constant quantities in synchronous \(dq\)-coordinates, providing a steady-state operating point for small-signal linearization.

The small-signal grid impedance model is given by four SISO real transfer operators that relate the \(dq\) small-signal currents and voltages, \[\begin{bmatrix} \Delta v_d(t)\\ \Delta v_q(t) \end{bmatrix} = \overbrace{\begin{bmatrix} Z_{dd}(\mathrm{\mathstrut p}) & Z_{dq}(\mathrm{\mathstrut p})\\ Z_{qd}(\mathrm{\mathstrut p}) & Z_{qq}(\mathrm{\mathstrut p})\\ \end{bmatrix}}^{=:Z_g(\mathrm{\mathstrut p})} \begin{bmatrix} \Delta i_d(t)\\ \Delta i_q(t) \end{bmatrix},\] with \(Z_g(\mathrm{\mathstrut p})\) being real 2-by-2 transfer operators, \(\mathrm{\mathstrut p}\!=\!\frac{\@ifnextchar^{\DIfF}{\mathop{\mathrm{\mathstrut d}}\nolimits^{}\futurelet\diffarg\let\DiffSpace\!\ifx\diffarg(\let\DiffSpace\relax \else \ifx\diffarg[\let\DiffSpace\relax \else \ifx\diffarg\{\let\DiffSpace\relax \fi\fi\fi\DiffSpace}}{\@ifnextchar^{\DIfF}{\mathop{\mathrm{\mathstrut d}}\nolimits^{}\futurelet\diffarg\let\DiffSpace\!\ifx\diffarg(\let\DiffSpace\relax \else \ifx\diffarg[\let\DiffSpace\relax \else \ifx\diffarg\{\let\DiffSpace\relax \fi\fi\fi\DiffSpace}t}\) is the differential operator, and \(\Delta\) denotes deviations from steady-state values.¹ An alternative equivalent representation can be obtained using complex variables. Define \[\boldsymbol{v}(t) := \Delta v_d(t) +j \Delta v_q(t), \qquad \boldsymbol{i}(t) := \Delta i_d(t) + j \Delta i_q(t).\] Straightforward algebraic manipulations [@Martin2004] then show that \[\label{eq:complex95model} \boldsymbol{v}(t) = \boldsymbol{G}_+(\mathrm{\mathstrut p}) \boldsymbol{i}(t) + \boldsymbol{G}_-(\mathrm{\mathstrut p}) \boldsymbol{i}^\ast(t),\tag{1}\] where \(\boldsymbol{i}^\ast\) is the complex conjugate of \(\boldsymbol{i}\), and \[\label{eq:complex95tfs} \begin{align} \boldsymbol{G}_+(\mathrm{\mathstrut p}) &= 0.5 \big({Z_{dd}(\mathrm{\mathstrut p}) + Z_{qq}(\mathrm{\mathstrut p})} + j ({Z_{qd}(\mathrm{\mathstrut p}) - Z_{dq}(\mathrm{\mathstrut p})}) \big) \\ \boldsymbol{G}_-(\mathrm{\mathstrut p}) &= 0.5 \big({Z_{dd}(\mathrm{\mathstrut p}) - Z_{qq}(\mathrm{\mathstrut p})} + j({Z_{dq}(\mathrm{\mathstrut p}) + Z_{qd}(\mathrm{\mathstrut p})})\big) \end{align}\tag{2}\] are two SISO complex transfer operators. TFs are obtained by applying the Laplace transform to the time-domain transfer operators, under which \(\mathrm{\mathstrut p}\) becomes the Laplace variable \(s\). The frequency response at \(\omega\) is obtained by setting \(s = j\omega\).

If \(Z_g(\mathrm{\mathstrut p})\) is symmetric, i.e. \(Z_{dd}(\mathrm{\mathstrut p}) = Z_{qq}(\mathrm{\mathstrut p}) = G_d(\mathrm{\mathstrut p})\) and \(Z_{qd}(\mathrm{\mathstrut p}) = -Z_{dq}(\mathrm{\mathstrut p}) = G_q(\mathrm{\mathstrut p})\), it holds that \(\boldsymbol{G}_-(\mathrm{\mathstrut p}) = 0\) and the model reduces to \(\boldsymbol{v}(t) = \boldsymbol{G}(\mathrm{\mathstrut p}) \boldsymbol{i}(t)\), with \[\label{eq:single95complex95tf} \boldsymbol{G}(\mathrm{\mathstrut p}) = G_d(\mathrm{\mathstrut p}) + j G_q(\mathrm{\mathstrut p}) = \boldsymbol{G}_+(\mathrm{\mathstrut p}).\tag{3}\] In this special case, \(Z_g(\mathrm{\mathstrut p})\) is represented by one SISO complex TF that can be estimated non-parametrically via one set of periodic measurements. However, in the asymmetric case, \(\boldsymbol{i}^\ast\) is always required, resulting in a double frequency model that cannot be estimated using one set of measurements unless further assumptions are imposed on the frequency response.

2.2 Grid-connected converter↩︎

To demonstrate the approach, we consider a grid-connected voltage-sourced converter with an \(LCL\) filter, as shown in Figure 1. An ideal DC link is assumed. The current control loop is implemented in the \(dq\)-frame, and grid synchronization is performed using a PLL that tracks the voltage on the \(LCL\) filter capacitor. The current controller output serves as the reference for the converter voltage, realized via PWM. Outer power control loops are omitted for the sake of clarity.

The grid is perturbed by adding a wideband excitation signal to the converter voltage references \(u_\mathrm{c,D}^\mathrm{ref}\) and \(u_\mathrm{c,Q}^\mathrm{ref}\). This provides a wider excitation bandwidth than adding it to the current reference. We use a zero-mean random binary signal (RBS) [@ljung1998system], which, unlike MLBS, is non-periodic and can have arbitrary length. Binary signals are popular due to their ideal crest factor, but they are inflexible for spectrum shaping. Excitation spectrum design is beyond the scope of this work.

2.3 Measurement setup↩︎

The measurement setup is shown in Figure 2. Samples of voltage \(v_{abc}(t)\) and current \(i_{abc}(t)\) are recorded, starting at the application instant of the excitation signal, and continue for its entire duration. For accurate identification, \(v_{abc}(t)\) and \(i_{abc}(t)\) are filtered with an anti-alias filter, such as a Chebyshev analog filter, before sampling. In addition, acquisition channels must be synchronized and relatively calibrated to eliminate transducer dynamics (\(S_i\) and \(S_v\)). Sampling is uniform in time (periodic) and should be fast enough for the bandwidth of interest; a possibility is to use the converter’s sampling time \(T_s\), which typically operates at the switching frequency or twice that (tens of kHz). Lastly, Park’s transformation converts the sampled phase measurements to a synchronous \(dq\)-frame.

3 Grid Impedance Identification↩︎

Given a data set \(D_N \!:= \!\{ (\boldsymbol{v}(t_n), \boldsymbol{i}(t_n)), n\! \in\! \{0, \dots, N\!-\!1\} \}\), the goal is to construct a non-parametric estimator of \(\boldsymbol{G}_+(j\omega), \boldsymbol{G}_-(j\omega)\) on a uniform grid of frequencies \(\omega_k\! =\! \frac{2\pi k}{NT_s}\), viz. \(D_N \, \mapsto \, \{ (\, \widehat{\boldsymbol{G}_+(j\omega_k)},\; \widehat{\boldsymbol{G}_-(j\omega_k)}\,), \;\, k \!\in\!\{ 0, \dots, N-1\} \}\). To this end, let \(\boldsymbol{V}\!_k\) and \(\boldsymbol{I}_k\) denote the \(N\)-point DFT of \(\{ \boldsymbol{v}(t_n) \}\) and \(\{ \boldsymbol{i}(t_n) \}\), respectively, defined as \[\label{eq:dft95spectrum} \boldsymbol{V}\!_k = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} \boldsymbol{v}(t_n) \mathrm{e}^{-j\omega_k n T_s}, \quad k \in \{0, \dots, N-1\}, \tag{4}\] and a similar expression for \(\boldsymbol{I}_k\). Recall that the \(N\)-point DFT can be computed very efficiently using fast Fourier transform algorithms, and that the spectra are periodic with period \(N\), namely \(\boldsymbol{I}_{k+N} = \boldsymbol{I}_{k}\). Because we are dealing with complex time-domain signals, the Hermitian symmetry property of the DFT does not hold, but it still satisfies conjugacy and reversal properties, that is, the DFT of \(\{\boldsymbol{i}^\ast(t_k)\}\) is given by \(\{\boldsymbol{I}^\ast_{(N-k)_N} \}\) where \((N-k)_N\) stands for “\(N-k\) modulo \(N\)".

3.1 Relation between voltage and current DFT spectra↩︎

Applying the finite-time (truncated) Fourier transform to 1 over \([0,NT_s]\), gives the model\[\label{eq:dft95relation} \begin{align} \boldsymbol{V}\!_k &= \boldsymbol{G}_+(j\omega_k) \boldsymbol{I}_k + \boldsymbol{G}_-(j\omega_k) \boldsymbol{I}_{(N-k)_N}^\ast + \boldsymbol{T}(j\omega_k),\\ \end{align}\tag{5}\] where \(\boldsymbol{T}(j\omega_k)\) is a transient term decaying at a rate \(\mathcal{O}({N^{\frac{1}{2}}})\). It accounts for sampling, spectral leakage which arises from the mismatch in initial and final conditions, and aliasing effects due to the truncations to the finite time interval \([0, NT_s]\). We emphasize that 5 is an exact relation between the DFT spectra even for arbitrary non-periodic measurements. For more details, the interested reader is referred to [@ljung1998system], [@Pintelon1997] and the Appendix.

This shows that identifying the grid impedance from a single short measurement cycle using DFT spectra presents two main challenges. The first stems from the error introduced by the transient term \(\boldsymbol{T}\). The second is that there are \(2N\) complex unknowns, but only \(N\) complex data equations are available. The first challenge is solved by estimating \(\boldsymbol{T}\) together with the complex TFs. The second challenge is addressed by using local parametric modeling, which provides additional data equations by using neighboring spectral lines, as outlined below.

3.2 Local parametric modeling↩︎

The idea of local parametric modeling (see [@Pintelon2010; @McKelvey2012; @Pintelon2021]) is as follows. To obtain a non-parametric estimate at a frequency \(\omega_k\), we approximate the complex TFs \(\boldsymbol{G}_+(j\omega_k), \boldsymbol{G}_-(j\omega_k)\), and \(\boldsymbol{T}(j\omega_k)\) over a short frequency range around \(\omega_k\) using a parametric model of low order. Estimates at different frequencies are correlated only via raw data, and therefore the method remains truly non-parametric in nature.

By their definition in 2 , \(\boldsymbol{G}_+\) and \(\boldsymbol{G}_-\) have the same poles, which are also poles of \(\boldsymbol{T}\). Therefore, for a local frequency interval \([\omega_{k-\ell}, \omega_{k+\ell}]\) centered on \(\omega_k\), we approximate \[\label{eq:LPM} \begin{align} \boldsymbol{G}_+(j\omega_{k+r}) &\approx \frac{\boldsymbol{B}_k^+(j\omega_{r})}{\boldsymbol{A}_k(j\omega_{r})}, \quad \boldsymbol{G}_-(j\omega_{k+r}) \approx \frac{\boldsymbol{B}_k^-(j\omega_{r})}{\boldsymbol{A}_k(j\omega_{r})}\\ \boldsymbol{T}(j\omega_{k+r}) &\approx \frac{\boldsymbol{C}_k(j\omega_{r})}{\boldsymbol{A}_k(j\omega_{r})}, \end{align}\tag{6}\] where \(r \in \{-\ell, \dots, \ell\}\), \(\boldsymbol{A}_k, \boldsymbol{B}^+_k, \boldsymbol{B}^-_k\), and \(\boldsymbol{C}_k\) are defined as complex polynomials in \(j\omega_r\) of degree \(R\), parameterized as \(\boldsymbol{A}_k(j\omega_r) := \sum_{m = 0}^R a_m(k) r^m\), and similarly for the other polynomials. Although these polynomials can have different degrees, using the same fixed degree for all and keeping it constant over \(k\) provides a simple and effective choice. We normalize \(\boldsymbol{A}_k\) so that \(a_0(k) = 1\) for every \(k\). With this parameterization, the vector of parameters to be estimated is \[\begin{align} &\theta_k :=\left[\begin{matrix} a_1(k) \cdots a_R(k) & b_0^+(k) \cdots\end{matrix}\right. \\ &\left.\begin{matrix} \cdots b_R^+(k) & b_0^-(k)\cdots b_R^-(k) & c_0(k) \cdots c_R(k) \end{matrix}\right]. \end{align}\]

It contains \(4R+3\) unknown complex parameters. From 5 and 6 , the local model linking the DFT spectra becomes \[\begin{align} \boldsymbol{A}_k(j\omega_r) \boldsymbol{V}\!_{k+r} = \,& \boldsymbol{B}_k^+(j\omega_r) \boldsymbol{I}_{k+r} + \boldsymbol{B}_k^-(j\omega_r) \boldsymbol{I}^\ast_{(N-k-r)_N}\\& + \boldsymbol{C}_k(j\omega_r) + \boldsymbol{E}_{k+r} \end{align}\] where \(\boldsymbol{E}_{k+r}\) accounts for the local parametric interpolation error. The parameters for each \(k\) can then be estimated by minimizing the errors in the least squares sense.²

Define the complex column vectors and the data matrix \[\begin{align} \boldsymbol{Y}\!_k &:= \begin{bmatrix} \boldsymbol{V}\!_{k-\ell}\!\! &\dots &\!\!\boldsymbol{V}\!_{k+\ell} \end{bmatrix}^\top\!\!, \quad \boldsymbol{U}^+_k := \begin{bmatrix} \boldsymbol{I}_{k-\ell} \!\!&\dots\!\! &\boldsymbol{I}_{k+\ell} \end{bmatrix}^\top\!\!,\\ \boldsymbol{U}^-_k &:= \begin{bmatrix} \boldsymbol{I}^\ast_{(N-k-\ell)_N} &\dots &\boldsymbol{I}^\ast_{(N-k+\ell))_N} \end{bmatrix}^\top,\\ \Phi_k &:= \begin{bmatrix} \boldsymbol{Y}\!_k \otimes \tilde{\phi}(r) & \boldsymbol{U}^+_k \otimes \phi(r) & \boldsymbol{U}^-_k \otimes \phi(r) & \boldsymbol{1} \otimes \phi(r) \end{bmatrix}, \end{align}\] with \(\tilde{\phi}(r) = \begin{bmatrix} r & r^2& \dots& r^R\end{bmatrix}^\top\) and \(\phi(r) := \begin{bmatrix}1 &\tilde{\phi}(r)^\top\end{bmatrix}^\top\), where \(\otimes\) denotes the Kronecker product. The parameter estimate is then obtained by solving the linear least-squares problem \(\min_\theta\| \boldsymbol{Y}\!_k - \Phi_k \theta \|\), with the closed-form solution \[\label{eq:LS95sol} \hat{\theta}_k = \Phi_k^\dagger \boldsymbol{Y}\!_k.\tag{7}\] The matrix \(\Phi_k^\dagger\) denotes the pseudo-inverse of \(\Phi_k\), which may be computed using a singular value decomposition with an appropriate scaling to improve numerical conditioning. The frequency response estimates at \(\omega_k\) are computed by setting \(r=0\) in 6 and recalling that \(a_0(k) = 1\). This gives \(\widehat{\boldsymbol{G}_+(j\omega_{k})} = \widehat{b_0^+(k)}, \quad \widehat{\boldsymbol{G}_-(j\omega_{k})} = \widehat{b_0^-(k)}\). Notice that the least-squares problems are independent over \(k\), so the computations may be optimized by parallelization.

To guarantee the uniqueness of 7 , the local frequency interval should be wide enough, and the measurement should be sufficiently exciting. A necessary condition is \(2\ell\! +\!1 \geq 4R\! +\! 3\) where \(\ell\) is the radius of the local frequency interval. The measurement is sufficiently exciting at \(\omega_k\) if \(\Phi_k\) is full rank. This imposes a condition on the local spectra \(\boldsymbol{Y}_k\) and \(\boldsymbol{U}_k\), and is satisfied for an RBS excitation signal. Notably, the proposed approach can estimate the equivalent impedance at \(\omega_g\) because it employs local models. This is despite the fact that the RBS does not excite the grid at \(\omega_g\).

3.3 Special cases↩︎

3.3.1 The periodic excitation case↩︎

If the excitation signal is repeated periodically and \(D_N\) contains an integer number of steady-state periods, then \(\boldsymbol{T}(j\omega) = 0\). In this case, it is not necessary to estimate the polynomials \(\boldsymbol{C}_k\), and the corresponding columns in \(\Phi_k\) can be safely removed. However, local parametric modeling is still needed due to the lack of a second measurement.

3.3.2 The \(dq\)-symmetric grid case↩︎

If it is known a priori that the equivalent impedance is \(dq\)-symmetric, the model reduces to \(\boldsymbol{V}\!_k = \boldsymbol{G}(j\omega_k) \boldsymbol{I}_k + \boldsymbol{T}(j\omega_k)\). In this case, it is not necessary to estimate the polynomials \(\boldsymbol{B}_k^-(j\omega)\), and the corresponding columns in \(\Phi_k\) can be safely removed. Although here we only have one SISO complex TF to identify, the effect of leakage errors remains. Therefore, local parametric modeling is still needed to eliminate the spectral leakage.

3.3.3 The periodic excitation and \(dq\)-symmetric grid case↩︎

Without spectral leakage or \(dq\)-asymmetry, an estimate of \(\boldsymbol{G}\) at frequency \(\omega_k\) can be obtained by a single division \(\widehat{\boldsymbol{G}(j\omega_k)} = {\boldsymbol{V}\!_k}/{\boldsymbol{I}_k}\). However, the use of this simple estimator (a.k.a empirical transfer function estimate, ETFE [@ljung1998system]) is justified only if steady-state measurements are possible; otherwise, significant errors would occur [@ljung1998system].

3.4 Extracting the four real TFs↩︎

As pointed out earlier, non-parametric estimates of complex TFs can be used directly for control and stability analysis. However, if needed, they can be mapped numerically to non-parametric estimates of \(Z_g\). For \(\boldsymbol{G}(s)\) in 3 , by noticing that \(\boldsymbol{G}^\ast(s) = [\boldsymbol{G}(s^\ast)]^\ast\), we find that \(G_d(s) = 0.5({\boldsymbol{G}(s) + \boldsymbol{G}^\ast(s)})\), \(G_q(s)\! =\! -0.5j({\boldsymbol{G}(s) \!-\! \boldsymbol{G}^\ast(s)})\). Hence, in the symmetric case \[\begin{align} \widehat{Z_{dd}(j\omega_k)} &= \frac{1}{2} (\widehat{\boldsymbol{G}(j\omega_k)} + [\widehat{\boldsymbol{G}(j\bar{\omega}_k)}]^\ast),\\ \widehat{Z_{qd}(j\omega_k)} &= \frac{1}{2j} (\widehat{\boldsymbol{G}(j\omega_k)} - [\widehat{\boldsymbol{G}(j\bar{\omega}_k)}]^\ast),\end{align}\] where \(\widehat{\boldsymbol{G}(j\bar{\omega}_k)}\) is the estimate of \(\boldsymbol{G}(-j\omega_{k})\), \(\bar{\omega}_k = \omega_{(N-k)_N}\), and \(k \in\{0,\dots, \frac{N}{2}\! -\!1\}\) assuming an even \(N\). Using a similar reasoning for the asymmetric case, 2 gives \[\begin{align} \widehat{Z_{dd}(j\omega_k)} &= \frac{1}{2} \Big( \widehat{\boldsymbol{G}_+(j\omega_k)} + [\widehat{\boldsymbol{G}_+(j\bar{\omega}_k)}]^\ast + \widehat{\boldsymbol{G}_-(j\omega_k)} + [\widehat{\boldsymbol{G}_-(j\bar{\omega}_k)}]^\ast \Big),\\ \widehat{Z_{qq}(j\omega_k)} &= \frac{1}{2} \Big( \widehat{\boldsymbol{G}_+(j\omega_k)} + [\widehat{\boldsymbol{G}_+(j\bar{\omega}_k)}]^\ast - \widehat{\boldsymbol{G}_-(j\omega_k)} - [\widehat{\boldsymbol{G}_-(j\bar{\omega}_k)}]^\ast \Big), \\ \widehat{Z_{dq}(j\omega_k)} &= \!\frac{-1}{2j} \Big( \widehat{\boldsymbol{G}_+(j\omega_k)}\! - [\widehat{\boldsymbol{G}_+(j\bar{\omega}_k)}]^\ast\! - \!\widehat{\boldsymbol{G}_-(j\omega_k)} + [\widehat{\boldsymbol{G}_-(j\bar{\omega}_k)}]^\ast \Big),\\ \widehat{Z_{qd}(j\omega_k)} &= \frac{1}{2j} \Big( \widehat{\boldsymbol{G}_+(j\omega_k)}\! - [\widehat{\boldsymbol{G}_+(j\bar{\omega}_k)}]^\ast \!+ \!\widehat{\boldsymbol{G}_-(j\omega_k)} - [\widehat{\boldsymbol{G}_-(j\bar{\omega}_k)}]^\ast \Big). \end{align}\]

4 Numerical simulation study↩︎

4.1 Experiment Setup↩︎

We demonstrate the performance of the proposed approach using detailed electromagnetic transient simulations in Matlab/Simulink using the Simscape Electrical toolbox.³ The converter is modeled via an IGBT bridge and a discrete space-vector PWM. The ode5 solver is used with a step size \(10^{-7} \si{\second}\).

Figure 3 shows the one-line diagram of the grid. Its analytically computed equivalent impedance exhibits rich dynamics with a few resonances, as shown in Figure 4. The converter, interfaced with the grid using an \(LCL\) filter, is controlled as shown in Figure 1, in addition to an outer PQ loop. The controllers are discretized parallel-form proportional–integral (PI) regulators. All relevant parameters are listed in Table 1.

For clarity, we assume stationary grid operating conditions and do not consider any grid ambient harmonics (no grid disturbances). However, we note that the approach can deal with ambient harmonics efficiently by removing the corresponding spectral lines from the identification process. Also note that although the grid \(dq\)-symmetric, the methods do not assume this symmetry. Instead, they identify a 2-by-2 MIMO model.

To ensure that the validation scenario remains close to a realistic one, the experiment parameters are first fixed to reflect practical limitations. In particular, the total excitation time is fixed at 1 s (and only one measurement cycle), and the sampling time \(T_s =10\,\si{\milli\second}\) is equal to the converter’s switching/control period. The amplitude of the excitation signal is \(\leq 0.05\) p.u. resulting in an acceptable \(\text{THD}_v\) (\(\approx 2.15\%\))⁴ that does not severely compromise power quality. We consider measurement devices of accuracy class 0.5% and model current and voltage noise as discrete-time Gaussian random variables. Note that in this scenario, \(N = 10^4\) and the DFT uniform frequency grid has a resolution of 1 Hz, which allows resonance peaks with a 3dB-bandwidth of 5 Hz to be captured by 5 spectral lines.

4.2 Comparison methods↩︎

Most existing grid impedance identification methods would struggle to provide an accurate estimate in this scenario, mainly due to the short measurement time. Approaches such as frequency sweep and impulse injection methods are not even applicable, as they require multiple measurement cycles and large excitation amplitudes, respectively. While other wideband DFT-based methods (e.g.[@martin2013wide; @riccobono2017noninvasive; @luhtala2018implementation; @roinila2017mimo; @Berg2022]) can be applied, their accuracy is severely compromised by the short measurement time. Discrete-time parametric methods (e.g., [@Haberle2023]), in contrast, are capable of handling short non-periodic measurements, but their accuracy is highly sensitive to the model order.

The proposed approach provides a robust solution to these inherent challenges. To demonstrate this, we compared the proposed approach with a parametric time-domain method using ARX models [@Haberle2023] and a non-parametric sequential perturbation method similar to [@martin2013wide; @riccobono2017noninvasive; @luhtala2018implementation; @roinila2017mimo; @Berg2022]. For the latter, i) the data set \(D_N\) is divided into two equal parts and treated as two different measurements, and ii) the standard Hamming window [@ljung1998system] is applied before computing the DFT spectra to reduce spectral leakage. The proposed approach is tested using local model orders \(2, 4, \dots, 10\), and the radius of the local frequency interval \(\ell = 4R+2\). The parametric discrete-time ARX method is tested using model orders (output lags) \(2, 4, \dots, 10\) and \(20\).

4.3 Evaluation↩︎

To isolate errors caused by measurement noise from those coming from the methods themselves, we simulated two cases: with and without measurement noise. The identification accuracy is reported using the following metrics. For any of the real TFs, let \(Z_{0:k} := \begin{bmatrix} Z(j\omega_0) & Z(j\omega_1) & \dots & Z(j\omega_k) \end{bmatrix}^\top\) be a vector of the true TF evaluated at frequencies \(\omega_0\) to \(\omega_k\), and denote its estimate by \(\widehat{Z}_{0:k}\); the fit metric is then defined as \[\label{eq:fit} \text{Fit\%} := 1 - \frac{\|\widehat{Z}_{0:k} - Z_{0:k} \|_2^2}{\| Z_{0:k} - \text{mean}(Z_{0:k}) \|_2^2} \times 100,\tag{8}\] with \(\text{mean}(Z_{0:k}) = \frac{1}{k+1} \sum_{n = 0}^k Z(j\omega_n)\). Notice that can assume negative values; larger values indicate better estimates. A perfect estimate has a fit of 100%. We also consider the following metric defined for the estimates \(\{\!\widehat{Z_g(\omega_\ell)}\!\}\): \[\label{eq:hinferror} \text{relative } H_\infty \text{ error } := \frac{\max\limits_{\ell\in \{0, \dots, k\}} \bar{\sigma}\left(\widehat{Z_g(j\omega_\ell)} - Z_g(j\omega_\ell) \right)}{ \max\limits_{\ell\in\{0, \dots, k\}}\bar{\sigma}\left(Z_g(j\omega_\ell)\right)},\tag{9}\] \(\{\!Z_g(j\omega_\ell)\!\}\) are the true values, \(\bar{\sigma}\) is the largest singular value.

4.4 Comparison results without measurement noise↩︎

The results for the case without measurement noise are summarized in Table [tab:comparison95without95noise], where the two accuracy metrics 8 and 9 are evaluated in the frequency band \([0,4]\) kHz (i.e. up to 80% of the Nyquist frequency). In this case, the error originates from the method itself; thus, a reliable identification method should exhibit high accuracy.

As the results show, the proposed approach gives an almost perfect estimate for all local orders \(R = 2, 4, \dots, 10\). In contrast, the discrete-time ARX method has large errors for small orders; these are model misspecification errors. In this example, they become negligible if high model orders are used. Recall that here the order of \(Z_g(s)\) is 10; however, in practice a true finite order may not exist. We also remark that the models obtained by the discrete-time parametric ARX method are not always stable. Lastly, the sequential perturbation method incurs significant errors, as expected, due to the spectral leakage despite the use of a Hamming window.

4.5 Comparison results with measurement noise↩︎

The results for the case with measurement noise are summarized in Table [tab:comparison95with95noise], where the two accuracy metrics 8 and 9 are evaluated in the frequency band \([0,2]\, \si{\kHz}\). Because the sequential perturbation method proved inadequate, we only compare the other two methods. The errors in this case originate from two sources: measurement noise and systematic errors of the methods.

The errors of the discrete-time parametric ARX method are relatively large for low orders. In fact, they are comparable to the errors obtained with noise-free data. This indicates that for low orders, the misspecification errors dominate the noise errors. Only with high order models the fit becomes acceptable. Yet, the relative \(H_\infty\) error of the model with the best Fit% (order 20) remains large (about 91%). This is due to the large errors concentrated around the resonance frequencies. In sharp contrast, the proposed approach is able to provide estimates with good accuracy regardless of the chosen local model order. This robustness is one of the main reasons why the proposed local modeling approach is superior to alternative parametric methods.

These observations become evident by inspecting the error plots in Figure 5. The errors of the discrete-time parametric ARX method are concentrated at low frequencies around the resonances (true responses are overlaid in dashed gray). In contrast, the errors of the proposed approach become noticeable only at frequencies higher than 2 kHz. Below this frequency, the estimates are accurate; see Figure 6 where the estimates of \(\boldsymbol{G}_+\) are shown together with the true values. The observed worsening in the accuracy above 2 kHz is due to the limitations imposed by the \(LCL\) filter that reduces the SNR at higher frequencies.

5 Conclusions↩︎

We proposed a non-parametric method for identifying dynamic small-signal \(dq\)-asymmetric grid impedances using grid-connected converters. Our approach avoids assumptions about the grid’s topology or structure and provides estimates over a wide frequency band using short, low-amplitude, non-periodic excitation. The key innovation lies in combining complex transfer function representations of asymmetric systems with local frequency-domain modeling techniques. This strikes a crucial balance between fully non-parametric methods, which require extended excitation times, and fully parametric methods, which are sensitive to the model order. Numerical simulations demonstrated the superior performance and robustness of our approach, with a measurement time of 1 s and a full frequency resolution of 1 Hz. Our future work will cover the experimental validation of this method in more realistic settings, identifying the converter admittance, error analysis, optimizing the excitation signal, and considering unbalanced and multi-converter scenarios.

Define the continuous-time Fourier and convolution integrals \[\mathcal{F}_a^b\{\boldsymbol{v}\}(\omega) \!:=\!\!\int_a^b\! \boldsymbol{v}(t) \mathrm{e}^{-j\omega t} \@ifnextchar^{\DIfF}{\mathop{\mathrm{\mathstrut d}}\nolimits^{}\futurelet\diffarg\let\DiffSpace\!\ifx\diffarg(\let\DiffSpace\relax \else \ifx\diffarg[\let\DiffSpace\relax \else \ifx\diffarg\{\let\DiffSpace\relax \fi\fi\fi\DiffSpace}t, \quad (\boldsymbol{g}\ast \boldsymbol{i})_a^b(t) \!:=\!\! \int_a^b\! \boldsymbol{g}(t-\tau) \boldsymbol{i}(\tau) \@ifnextchar^{\DIfF}{\mathop{\mathrm{\mathstrut d}}\nolimits^{}\futurelet\diffarg\let\DiffSpace\!\ifx\diffarg(\let\DiffSpace\relax \else \ifx\diffarg[\let\DiffSpace\relax \else \ifx\diffarg\{\let\DiffSpace\relax \fi\fi\fi\DiffSpace}\tau,\] and, for clarity, let us drop their arguments, \(\omega\) and \(t\), from the notation. We start by deriving the expressions for the \(dq\)-symmetric case. The general case is established similarly, considering the responses of \(\boldsymbol{G}_+\) and \(\boldsymbol{G}_-\) separately.

The symmetric case↩︎

First, note that \(\boldsymbol{v}(t) = (\boldsymbol{g}\ast \boldsymbol{i})_0^t + (\boldsymbol{g}\ast \boldsymbol{i})_{-\infty}^0\) where \(\boldsymbol{g}\) is the impulse response of the stable and causal TF \(\boldsymbol{G}\) and \(t\in[0,T]\) with \(T=NT_s\) where \(N\) is the number of samples, and \(T_s\) is the sampling time. Then \[\begin{align} \boldsymbol{V}(\omega) &:= \mathcal{F}_{-\infty}^\infty\{ \setbox 0=\boldsymbol{v} \stackon[0pt]{\stackon[1pt]{\boldsymbol{v}}{ \rule{.4pt}{1.1pt}\kern\wd 0\kern-.4pt\kern-.4pt\rule{.4pt}{1.1pt}}} {\rule{\wd 0}{.4pt}}\} = \mathcal{F}_0^T\{\boldsymbol{v}\} \\ &= \mathcal{F}_0^T\{(\boldsymbol{g}\ast \boldsymbol{i})_0^t\} + \mathcal{F}_0^T\{(\boldsymbol{g}\ast \boldsymbol{i})_{-\infty}^0\} \\ & = \mathcal{F}_0^T\{(\boldsymbol{g}\ast \boldsymbol{i})_0^\infty\} -\underbrace{\mathcal{F}_0^T\{(\boldsymbol{g}\ast \boldsymbol{i})_t^\infty\}}_{= \, 0 \text{ by causality of \boldsymbol{G}}} + \mathcal{F}_0^T\{(\boldsymbol{g}\ast \boldsymbol{i})_{-\infty}^0\} \\[0.5em] &= \mathcal{F}_0^\infty\{(\boldsymbol{g}\ast \boldsymbol{i})_0^\infty\} -\mathcal{F}_T^\infty\{(\boldsymbol{g}\ast \boldsymbol{i})_0^\infty\} +\mathcal{F}_0^T\{(\boldsymbol{g}\ast \boldsymbol{i})_{-\infty}^0\}\\ &= \boldsymbol{G}(j\omega) \boldsymbol{I}(\omega) - \boldsymbol{V}_{\!\text{fin}}(\omega) + \boldsymbol{V}_{\!\text{init}}(\omega). \end{align}\] Here, \(\boldsymbol{V}(\omega)\) and \(\boldsymbol{I}(\omega)\) are the spectra of the complex continuous-time signals \(\setbox 0=\boldsymbol{v} \stackon[0pt]{\stackon[1pt]{\boldsymbol{v}}{ \rule{.4pt}{1.1pt}\kern\wd 0\kern-.4pt\kern-.4pt\rule{.4pt}{1.1pt}}} {\rule{\wd 0}{.4pt}}\) and \(\setbox 0=\boldsymbol{i} \stackon[0pt]{\stackon[1pt]{\boldsymbol{i}}{ \rule{.4pt}{1.1pt}\kern\wd 0\kern-.4pt\kern-.4pt\rule{.4pt}{1.1pt}}} {\rule{\wd 0}{.4pt}}\) which are defined as being equal to \(\boldsymbol{v}\) and \(\boldsymbol{i}\) over the finite time interval \([0, T]\) and equal to zero elsewhere. In the third row, the second term is identically zero for a casual \(\boldsymbol{G}\) (because then \(\boldsymbol{g}(t) = 0 \; \forall t <0)\). The last equality follows from the Fourier transform convolution theorem, which implies \(\mathcal{F}_0^\infty\{(\boldsymbol{g}\ast \boldsymbol{i})_0^\infty\} = \boldsymbol{G}(j\omega) \mathcal{F}_0^\infty\{\boldsymbol{i}(t)\}\), and assuming \(\boldsymbol{i}(t) =0 \; \forall t>T\) (this assumption is benign because \(\boldsymbol{G}\) is causal); hence \(\mathcal{F}_0^\infty\{(\boldsymbol{g}\ast \boldsymbol{i})_0^\infty\} = \boldsymbol{G}(j\omega) \boldsymbol{I}(\omega)\). Lastly, we defined the terms \(\boldsymbol{V}_{\!\text{fin}}(\omega):= \mathcal{F}_T^\infty\{(\boldsymbol{g}\ast \boldsymbol{i})_0^\infty\}\), \(\boldsymbol{V}_{\!\text{init}}(\omega) :=\mathcal{F}_0^T\{(\boldsymbol{g}\ast \boldsymbol{i})_{-\infty}^0\}\) as the spectral leakage terms.

The next step is to relate the continuous spectra \(\boldsymbol{V}(\omega)\) and \(\boldsymbol{I}(\omega)\) to the discrete ones obtained by the DFT in 4 . The Dirac delta functions \(\delta(t)\) constitute the classical tool. Recall that these generalized functions are defined using continuous test functions \(\varphi\) via the identity \(\int_{-\infty}^{\infty} \delta(\tau) \varphi(\tau) \@ifnextchar^{\DIfF}{\mathop{\mathrm{\mathstrut d}}\nolimits^{}\futurelet\diffarg\let\DiffSpace\!\ifx\diffarg(\let\DiffSpace\relax \else \ifx\diffarg[\let\DiffSpace\relax \else \ifx\diffarg\{\let\DiffSpace\relax \fi\fi\fi\DiffSpace}\tau = \varphi(0)\), from which all properties of \(\delta\) can be deduced. For example, the impulse train \(\sum_{n=-\infty}^\infty \delta(t+nT)\), which plays an important role in sampling, is a periodic function with a period \(T\), and can be expanded using a Fourier series as \[\sum_{n=-\infty}^\infty \delta(t+nT) = \frac{1}{T} \sum_{n=-\infty}^\infty \mathrm{e}^{jn\omega_1t}, \qquad \omega_1 := \frac{2\pi}{T}.\] By convolving both sides of the last equation with \(\setbox 0=\boldsymbol{v} \stackon[0pt]{\stackon[1pt]{\boldsymbol{v}}{ \rule{.4pt}{1.1pt}\kern\wd 0\kern-.4pt\kern-.4pt\rule{.4pt}{1.1pt}}} {\rule{\wd 0}{.4pt}}(t)\) we get \[\sum_{n = -\infty}^\infty \setbox 0=\boldsymbol{v} \stackon[0pt]{\stackon[1pt]{\boldsymbol{v}}{ \rule{.4pt}{1.1pt}\kern\wd 0\kern-.4pt\kern-.4pt\rule{.4pt}{1.1pt}}} {\rule{\wd 0}{.4pt}}(t+nT) = \frac{1}{T} \sum_{n=-\infty}^\infty \boldsymbol{V}(n\omega_1) \mathrm{e}^{jn\omega_1t}.\] In light of this, the Fourier transform symmetry theorem can be invoked to deduce the following relation between the samples \(\{\boldsymbol{v}(nT_s)\}\) and the continuous spectrum \(\boldsymbol{V}(\omega)\) \[\sum_{n = -\infty}^\infty \boldsymbol{V}(\omega+n\omega_s) = \frac{2\pi}{\omega_s} \sum_{n=0}^{N} {c(n)}{\boldsymbol{v}(nT_s)} \mathrm{e}^{-jnT_s\omega}.\] This is the celebrated Poisson’s sum formula [@papoulis1977signal], where the finite limits of the sum on the right-hand side are because \(\boldsymbol{V}\) is defined by a finite-time Fourier integral over \([0,T]\), and the factor \[c(n) := \begin{cases} \frac{1}{2},\qquad n \in \{0,N\}\\ 1, \qquad \text{otherwise}. \end{cases}\] is needed because of the discontinuity of \(\setbox 0=\boldsymbol{v} \stackon[0pt]{\stackon[1pt]{\boldsymbol{v}}{ \rule{.4pt}{1.1pt}\kern\wd 0\kern-.4pt\kern-.4pt\rule{.4pt}{1.1pt}}} {\rule{\wd 0}{.4pt}}\) at the end points of the observation time window.

Sampling the spectra uniformly over the unit circle, with a frequency resolution \(\frac{2\pi}{T}\), leads to \[\frac{2\pi}{\omega_s} \sum_{n=0}^{N} {c(n)}{\boldsymbol{v}(t_n)} \mathrm{e}^{-jnT_s\omega_k} = \sum_{n = -\infty}^\infty \boldsymbol{V}(\omega_k+n\omega_s)\] where \(t_n := nT_s\), and \(\omega_k :=\frac{2\pi k}{T}\), \(k \in \{0, \dots, N-1\}\). Rearranging and using the definition in 4 , we find that \[\boldsymbol{V}\!_k = \frac{1}{\sqrt{N} T_s}\sum_{n = -\infty}^\infty \boldsymbol{V}\!\left(\omega_n - n \omega_s\right) + \frac{\boldsymbol{v}(0) - \boldsymbol{v}(T)}{2\sqrt{N}}.\] The same relation holds for \(\boldsymbol{I}_k, \boldsymbol{V}_{\!\text{init},k}\), and \(\boldsymbol{V}_{\!\text{fin},k}\). From this we get the model \[\boldsymbol{V}\!_k = \boldsymbol{G}(j\omega_k) \boldsymbol{I}_k + \boldsymbol{T}(j\omega_k)\] where \(\boldsymbol{T}(j\omega_k) = \boldsymbol{V}_{\!\text{init},k} - \boldsymbol{V}_{\!\text{fin},k} + \boldsymbol{\alpha}_k\) is the transient term with \(\boldsymbol{\alpha}_k\) representing the aliasing effects.

The asymmetric case↩︎

The starting point here is 1 , which may be re-written as \[\boldsymbol{v}(t) = \boldsymbol{v}_+(t) + \boldsymbol{v}_-(t),\] with \[\begin{align} \boldsymbol{v}_+(t) &= (\boldsymbol{g}_+\ast \boldsymbol{i})_0^t \;\,+ (\boldsymbol{g}_+\ast \boldsymbol{i})_{-\infty}^0,\\ \boldsymbol{v}_-(t) &= (\boldsymbol{g}_-\ast \boldsymbol{i}^\ast)_0^t + (\boldsymbol{g}_-\ast \boldsymbol{i}^\ast)_{-\infty}^0, \end{align}\] in which \(\boldsymbol{g}_+\) and \(\boldsymbol{g}_-\) are the impulse responses of causal \(\boldsymbol{G}_+\) and \(\boldsymbol{G}_-\), respectively. The above results, from the symmetric case, can be applied directly to \(\boldsymbol{v}_+(t)\) and \(\boldsymbol{v}_-(t)\) to show that the DFTs of \(\{\boldsymbol{v}_+(t_n)\}\) and \(\{\boldsymbol{v}_-(t_n)\}\) are \[\begin{align} \boldsymbol{V}\!_{+k} &= \boldsymbol{G}_+(j\omega_k) \boldsymbol{I}_k + \boldsymbol{T}_+(j\omega_k),\\ \boldsymbol{V}\!_{-k} &= \boldsymbol{G}_-(j\omega_k) \boldsymbol{I}_{(N-k)_N}^\ast + \boldsymbol{T}_-(j\omega_k),\\ \end{align}\] The reversed and conjugated DFT spectrum of \(\boldsymbol{i}\) in the last equation arises because \(\mathcal{F}_0^\infty\{(\boldsymbol{g}\ast \boldsymbol{i}^\ast)_0^\infty\} = \boldsymbol{G}(j\omega) [\boldsymbol{I}(-\omega)]^\ast\), after noticing that \(\boldsymbol{I}^\ast(\omega) := \mathcal{F}_0^\infty\{\boldsymbol{i}^\ast(t)\} = [\boldsymbol{I}(-\omega)]^\ast\). From this we directly get the model \[\boldsymbol{V}\!_k = \boldsymbol{G}(j\omega_k) \boldsymbol{I}_k + \boldsymbol{G}_-(j\omega_k) \boldsymbol{I}_{(N-k)_N}^\ast + {\underbrace{\boldsymbol{T}(j\omega_k)}_{=\boldsymbol{T}_+(j\omega_k)+ \boldsymbol{T}_-(j\omega_k)}}\] which is identical to 5 .

Observe that the transient term \(\boldsymbol{T}(j\omega_k)\) carries two effects:

1. the spectral leakage captured by \(\boldsymbol{V}_{\!\text{init},k}\) and \(\boldsymbol{V}_{\!\text{fin},k}\); these two terms come from unforced decaying responses due to the initial and final conditions, respectively, and therefore they have the same poles as \(\boldsymbol{G}_+\) and \(\boldsymbol{G}_-\),

2. aliasing effect captured by \(\alpha_k\); this term has infinite repetition of poles due to the folding of the spectrum.

---

1. In practice we get these by removing the mean of the time series↩︎

2. This corresponds to fitting a local continuous-time ARX model to the spectra only over \([\omega_{k-\ell}, \omega_{k+\ell}]\).↩︎

3. Code is available at https://doi.org/10.3929/ethz-c-000784827.↩︎

4. Based on 10 cycles, all intraharmonics & up to the 50^th harmonic.↩︎