Handling Initial Conditions in Vector Fitting for Real Time Modeling of Power System Dynamics

This paper develops a predictive modeling algorithm, denoted as Real-Time Vector Fitting (RTVF), which is capable of approximating the real-time linearized dynamics of multi-input multi-output (MIMO) dynamical systems via rational transfer function matrices. Based on a generalization of the well-known Time-Domain Vector Fitting (TDVF) algorithm, RTVF is suitable for online modeling of dynamical systems which experience both initial-state decay contributions in the measured output signals and concurrently active input signals. These adaptations were specifically contrived to meet the needs currently present in the electrical power systems community, where real-time modeling of low frequency power system dynamics is becoming an increasingly coveted tool by power system operators. After introducing and validating the RTVF scheme on synthetic test cases, this paper presents a series of numerical tests on high-order closed-loop generator systems in the IEEE 39-bus test system.


I. INTRODUCTION
D UE to the aggressive deployment of Wide Area Monitoring Systems (WAMS), a deluge of time series data streams are emerging from modernizing smart grids. For instance, as of 2017, Operating Procedure No. 22 specifies that all generation above 100MW in the ISO New England system must provide Phase Measuring Unit (PMU) observability at the point of interconnection [1]. Similarly, the 2018 Nodal Operating Guide specifies that all new generators above 20MVA in the ERCOT (Texas) system must provide PMU observability [2]. In order to both capitalize on these data streams and enforce effective dynamic security assessment (DSA) in the face of a rapidly modernizing energy landscape, the North American Electric Reliability Corporation (NERC) implemented the so-called MOD-033-1 directive in July of 2017 [3]. Not only does "MOD33" mandate the continual development of static and dynamic network planning models, but it also requires that the simulated response of these models be compared to actual time series data collected in the network to validate their accuracy [4]. Associated procedures must also resolve model prediction aberrations. Traditional "staged" testing of generator models can be costly and inconvenient, because the generator must go offline [5]. Thus, performing model tuning and validation online is particularly valuable.
In the electromechanical frequency range, power system dynamics are still dominated by synchronous generators [6]. A wealth of literature exists on the nonlinear physical modeling of these generators and their associated high-order controllers [7]. MOD33-flavored tuning and validation of these models, though, still faces a variety of practical obstacles; these include the exclusion of unmodeled dynamics, measurement noise, uncertain physical parameter values, unknown controller changes, and, in the limiting case of linearized modeling, residual nonlinearities and drifting equilibrium points. In recent years, a variety of "grid-scale" methods have been proposed for the purpose of validating these generator + controller models in the presence of strong network perturbations [3], [5], [8]- [11]. Most of these perform some variation of Kalman filter-based parameter tuning via "play-back" simulation.
The vast majority of the proposed modeling and validation algorithms are characterized by two salient features. First, they are parameterized by a physical prior; that is, they leverage the structure of a given physical generator model can assumably be tuned to explain the full set observed dynamics. Second, these algorithms are typically designed to be deployed in the presence of a sufficiently strong network perturbation (e.g. ground fault). Together, these assumptions allow for a variety of helpful assumptions, such as the neglecting of initial state decay and the availability of a high-fidelity prior. Such algorithmic features, although widespread, can be burdensome if the prior model is incomplete, or if real-time model updates are desired 1 and no recent perturbation events have occurred.
Predictive modeling algorithms which are free from both of these limitations can provide a variety of advantages. If the modeling of generator systems can be performed by a parameterized "black box" approximation, then the associated algorithm will be freed from the constraints of a potentially erroneous, or completely unknown, prior model. Additionally, if this algorithm can be implemented in the face of ambient, rather than severe, operating conditions, then updates can be performed on the fly, and otherwise unused PMU data can be transformed into actionable information.
There is a large collection of literature devoted to the creation of linear dynamical black box models from time series data in power systems, primarily at the "wide-area" level, but also at the individual generator level. Excellent reviews on the topic are provided in [12], [13]. One of the most popular "modal identification" methods, known as Prony analysis, was first applied to power systems in 1990 [14]. It was further shown in [15] that Prony analysis applied to transient data and Wiener-Hopf linear prediction applied to ambient data could both identify very similar system dynamics. Subsequent publications further refined these methods [16], [17], but most of them were focused on identification of system modes (i.e. poles), rather than the full transfer functions. More recent work has focused on building model reduced transfer functions from measured terminal data. Using a simple autoregessive with exogenous input (ARX) model, [18] develops a process for computing a dynamic equivalent transfer function with simulated or measured data. Voltage magnitude and frequency are chosen as inputs, whiles active and reactive power flows across lines are chosen as outputs. MIMO transfer function approaches are developed in [19], where poles, zeros, and model order are iteratively perturbed based on a least-squares numerical procedure inside a standard MATLAB toolbox.
Most recently, vector fitting approaches have been leveraged to construct power system models from measured field data. In [20], authors use the numerical Laplace transform to push the problem into the frequency domain. Next, vector fitting and relative dominant-pole measurements (RDPMs) are iteratively leveraged to build estimations of the system modes. Similarly, [21] combines time domain vector fitting and the popular ring down analysis to develop a novel modal estimation routine. In a more generic application, [22] uses out-of-the-box TDVF to construct a SISO linear model of a turbine governor system from measured PMU; droop gain values are then captured from the derived model.
The black-box identification of dynamic models from time series data goes far beyond power system applications and is regarded as a well established field [23], [24]. As noted, one prominent approach is the Time-Domain Vector Fitting (TDVF) scheme [23], [25], [26]. This approach is applicable for estimating rational models of linear systems, starting from the time domain samples of input and output data. Unfortunately, canonical TDVF has the following restrictions: • The system must be at rest when the data acquisition begins. This guarantees that the input and output data are related by a linear transfer function with no contribution from the zero-input (i.e. initial state decay) response. • For a MIMO system with P ports, only separate modeling of each column of the transfer matrix is possible. This limitation requires that only one of the P input components is acting on the system during the separate modeling periods for each transfer matrix column.
Since online modeling requires the unknown system to be characterized (i) during its dynamic evolution and (ii) in the presence of concurrently acting inputs, the above limitations often make TDVF unsuitable for real-time MIMO applications. This paper introduces a generalized extension of TDVF, known as Real-Time Vector Fitting (RTVF), in order to address these limitations. The remainder of this paper is organized as follows. In Section II, we recall the basic TDVF scheme in order to provide a framework on which the proposed approach is built. We then present our key contributions, namely: 1) In Section III, we remove the restrictions of the basic TDVF scheme by proposing a new problem setting. 2) In Section IV, we develop a first generalization of the basic TDVF scheme by removing the requirement of pure zero-state conditions, hence allowing for the presence of non-vanishing initial conditions. 3) In Section V, we develop a second generalization by allowing all input components to act concurrently in the time series data sets used to train the RTVF model.
In Section VI, we provide implementation details. Numerical results, based on several examples, are given in Section VII.

A. Notation
We denote a generic scalar as x, a generic vector as x or X 0 , and a generic matrix as X. The identity matrix is denoted as I, with size inferred from the context. We use the symbol s for the complex frequency (Laplace) variable; R and C represent the sets of real and complex numbers, respectively.
We consider a possibly nonlinear dynamic system S with input and output signals denoted as u(t) ∈ R P and y(t) ∈ R P , respectively. For later use, we denote with x(t) ∈ RN some unknown system state vector, although we assume no information on the internal system representation. We assume instead that a measurement tool is available that returns K samples of time-domain input-output vectors acquired at sampling rate F s . Without loss of generality, we set t 1 = 0 throughout this paper. All derivations will hold true for t 1 = 0, provided the time variable is redefined as t ← t − t 1 .

B. The Standard Setting for Data-Driven Modeling
When the underlying system S is Linear and Time-Invariant (LTI), Laplace-domain input and output signals are related by An estimate H(s) of the true transfer functionH(s) can be determined from the samples (1) through one of the several available data-driven model order reduction methods.
In particular, the TDVF scheme [23], [25], [26] considered in this work assumes that the system is initially at rest, and the initial conditions vanish identically as This setting guarantees that only the zero-state response contribution is present in the output samples.

C. Time-Domain Vector Fitting
The basic TDVF scheme assumes availability of • time series of each output y ij (t) at port i excited by a single input u j (t) placed at port j and acting alone, with u k =j = 0; this requirement imposes a restriction on the training sequences that can be used for model extraction; • some initial estimate of the dominant system poles {q n , n = 1, . . . , N }. Usually, such poles are initialized as random real or complex conjugate pairs with Re {q n } < 0 and |q n | < Ω, where Ω is the modeling bandwidth of interest [23], [27]. Based on the training data, TDVF constructs the approximation ij and d n are unknown coefficients to be determined via a linear least squares solution. In (4), superscript (n) denotes the result obtained by the singlepole filter (i.e. convolution) via on any arbitrary signal z(t) with z(0) = 0. The regression problem (4) corresponds to the frequency-domain relation which provides an element-wise rational approximation H ij (s) ≈H ij (s) written in barycentric form. The initial poles q n in (6) cancel out, and the actual poles of H ij (s) correspond to the zeros z n of the denominator Problem (4) is solved iteratively, by using these zeros as starting poles for the next iteration via q n ← z n . Iterations stop when the poles and/or the least-squares residual (fitting error) stabilize [23], [28], or alternatively when a maximum number of iterations ν max is reached.

III. A NEW PROBLEM SETTING FOR TDVF
Let us now consider system S operating in real time, with input and output measured signals (1) collected during system operation. Since data recording may start at an arbitrary time instant when system is not at rest, the output samples y(t) may include contributions from both the zero-input and the zerostate response, implying that all initial conditions cannot be considered as vanishing as in (3). Moreover, all input channels are expected to be active concurrently. Therefore, the basic TDVF assumptions do not hold and need a generalization.
From now on, we assume a mildly nonlinear systems S, whose dynamics can be approximated [29] by splitting inputs, outputs and states as where are regarded as non-necessarily vanishing initial conditions. The evolution of the small-signal components u(t), y(t), x(t) can be accurately related by an LTI operator with transfer functionH(s). Our main objective is to devise a numerical scheme that, based on the samples (1), returns an estimate H(s) ≈H(s) of the small-signal transfer function. This objective is attained by exploiting two generalizations of the basic TDVF scheme: 1) We remove the requirements of pure zero-state conditions by allowing for the presence of nonvanishing initial conditions (9). This problem is addressed in Sec. IV. 2) We allow for all input components u j acting concurrently in the training time series, as in common system operation conditions. This problem is addressed in Sec. V.

IV. HANDLING INITIAL CONDITIONS
In this section, we formulate a generalized identification problem that provides an estimate H(s) of the small-signal transfer function, assuming non-vanishing initial conditions. Although initial conditions U 0 , Y 0 are known from the training data, no information on the initial state X 0 is available. Thus, a direct decomposition of y(t) into its zero-state and zero-input contributions is generally not possible. This section provides theoretical justification for the proposed RTVF formulation of Section V and may be skipped at first reading.
In order to characterize the role of the unknown initial state, we first consider a generic LTI system in state-space form: with the only restriction that A should be nonsingular. Inserting the signal decomposition (8) into (10) leads tȯ The output (12) can be equivalently rewritten by splitting the constant and the time-varying small-signal components as Two scenarios are possible: 1) The system is at constant steady-state 2 for t = 0 (equivalently, ∀t ≤ 0). Under this assumption, all small-signal components vanish for t ≤ 0. Therefore, (11) reduces to and provides the initial state condition X 0 . Since A is nonsingular, the system has no poles at the origin and supports constant steady-state operation. Combining (11) with (15), for t > 0, the small-signal components fulfill the standard dynamic equatioṅ Combining (14) and (16) provides the small-signal transfer functionH(s) in terms of the state-space matrices: withD(s) = |sI − A|. Identification of a rational model forH(s) can be performed by subtracting the initial conditions U 0 , Y 0 from the input and output signals and then applying a zero-state identification scheme, such as basic TDVF, to small-signal components u(t), y(t).
2) The second scenario is relevant for our application, and it corresponds to the case where the system is not operating under constant steady-state conditions for t < 0. In this setting, (15) does not hold and Therefore, even if the initial conditions U 0 , Y 0 are removed from the input and output signals, the corresponding small-signal output y(t) still includes a contribution from the initial state. This contribution is analyzed next.

A. Characterization of residual zero-input contributions
Assuming that U 0 and X 0 are known, system evolution in terms of small-signal state components is obtained by integrating the dynamic equation (11) for t > 0 as Simple algebraic manipulations allow us to rewrite (19) as where the contribution of the initial state condition is explicit. Note that, in case of steady-state operation for t < 0, the second term in (20) vanishes since T 0 = 0, and the corresponding solution reduces to the solution of the smallsignal system (16). The term T 0 can thus be considered as the difference between the actual initial state X 0 and the constant state that would be obtained if the system were operating under steady-state conditions excited by the constant input U 0 . Taking the Laplace transform of (20) yields Inserting (21) into the output equation (14) leads to whereH(s) is given by (17). Additionally, whereG(s) is an unknown polynomial vector. Relation (22) is therefore equivalent to The two terms in (24) share the same denominatorD(s) up to a pole at s = 0, which represents the constant contribution of the non-vanishing initial conditions. This observation is the key enabling factor for building a self-consistent vector fitting scheme to estimate model H(s) ≈H(s), as it properly takes into account the presence of the additional term Γ 0 (s) in (22).

V. THE REAL-TIME VECTOR FITTING SCHEME
In this section, we propose our Real-Time Vector Fitting (RTVF) scheme. In particular, our key idea and the proposed generalizations follow: 1) The key idea of our method consists of adding to the standard rational transfer function expression an extra term that shares the same denominator D(s) plus a pole at s = 0, representing the constant contribution of the non-vanishing initial conditions.
2) Since polynomials N(s) and D(s) are typically expanded in a standard barycentric form as in (6) using an initial pole set q n , we propose the use of a similar expansion for the components of the unknown vector G(s): 3) We further propose to account for multiple inputs potentially exciting the system simultaneously by applying linear superposition and expressing each output component Y i (s) in terms of all input components U j (s). The above considerations lead us to parameterize and formulate the fitting condition (25) as Multiplying both sides by the common denominator and taking the inverse Laplace transform leads to the following time-domain fitting condition for t ≥ 0: i Θ (n) (t), ∀i = 1, . . . , P, (28) where Θ(t) is the Heaviside unit step function.

VI. IMPLEMENTATION
In this section, we provide a compact and efficient formulation of the least squares formulation at the heart of the RTVF routine. Writing (28) for t = t k , with k = 1, . . . , K, leads to the following RTVF condition in matrix form: where the vectors collecting the unknown coefficients are and where the regressor matrices collecting the filtered signal time samples (according to (5)) are defined as Further, by defining and collecting all components, (29) reveals the borderedblock-diagonal structure of the RTVF least squares system:  Standard techniques can be employed to avoid the all-zero trivial solution, as explained in [23], [30]. As in standard TDVF, once the set of unknown coefficients is found by solving (35), the zeros z n of the denominator D(s) are computed and used as initial poles for the next iteration. The process is repeated until convergence. Pseudocode for RTVF is provided in Algorithm 1. The final steps (lines 8-10) estimate the residues of a rational approximation based on the fixed poles obtained from the pole relocation process (lines 1-7). In line 10, the vector y i collects all time samples of the i-th small-signal output component.

A. Consistency
We tested the RTVF consistency by running a systematic experimental campaign over a set of synthetic randomly generated LTI reference systems, with the objective of checking Compute filtered signalsỹ Build and solve the least squares problem (35) 5: Compute the zeros zn of denominator D(s) in (7) 6: Set qn ← zn 7: end for 8: Set D(s) = 1 9: Compute filtered signalsũ whether RTVF could provide accurate estimates of all system poles. All modeled systems shared the same dynamic order of N = 10, but had different sizes, with P ranging from 2 to 30. The set of input-output data were generated as colored noise which showed a flat power spectrum up to angular frequency ω max , where the fastest pole of the reference system appeared. The sampling frequency was fixed to F s = 10 ω max /2π , and the total number of collected samples was K = 5000 in all cases. The modeling window started at sample k = 250.
Three metrics were used to assess performance of RTVF: • the consistency of the pole estimates, as measured by the Hausdorff distance 3 d H (P, Q) between the set of exact poles P = {p 1 , . . . , p N } of the true system and the set of numerically computed poles Q = {q 1 , . . . , q N }; • the worst-case time domain output error, computed as • and the RMS-normalized maximum error, computed as The experiments showed that RTVF recovered the system poles almost exactly, with a set distance d H (P, Q) ≤ 10 −10 for all 29 test cases. Similar results were obtained from the 3 The Hausdorff distance between two sets P and Q is defined as output errors: both E ∞ and E RMS ∞ were less than 10 −11 . As Fig. 1 shows, there is no practical difference between the model and the output data samples. In this idealized setting, we conclude that the performance of RTVF is excellent across all investigated metrics.
To further test the consistency of RTVF, we simulated the presence of measurement noise on the input and output signals (for the case where P = 2, N = 10). Signal corruption was performed by adding a vector of zero-mean Gaussian random variables x n to any input or output small-signal vector x as with a prescribed signal to noise ratio In our experiments, we considered increasing levels of SNR, ranging from 10 to 100, with resolution steps of 2. For each level of SNR, we modeled R = 50 different synthetic systems, and we computed the average Signal to Error Ratio (SER) both in time and frequency domain, which is defined as follows. Let z be a vector collecting the samples of either a reference timedomain output signal or a target frequency-domain transfer matrix element, and z M the corresponding response of one of the R models. Then for this signal the SER is defined as For any fixed SNR level, we computed the time-domain TD-SER by averaging the performance induced by (41) over the R models and the two output signals. The frequencydomain FD-SER was computed in the same way, by averaging over the transfer matrix elements. These two metrics are shown in the top panel of Fig. 2.
The meaning of the SER being above or below the black line threshold is that RTVF is either rejecting or amplifying the presence of the noise on the data, respectively. Since the TD-SER is always above the plane bisector (solid black line), RTVF is able to partially reject the presence of measurements noise in the training data. This noise rejection property is expected, since the basis functions involved in the estimation procedure effectively filter the noisy input and output signals via (5). On the other hand, FD-SER follows the bisector almost exactly, confirming also a good frequency-domain accuracy.
The noise-corrupted training input signals for a representative test case for SNR = 16 are depicted in the middle panel of Fig. 2, whereas the corresponding extracted model is validated against the reference time-domain output in the bottom panel. Even with this significant amount of noise, the time-domain prediction capabilities of the model are excellent.

B. IEEE 39-Bus System Results
In order to test the performance of RTVF in a simulated power system setting, we collected data from time domain simulations performed on the IEEE 39-Bus New England system. This system includes 10 generators and 19 ZIP loads; the associated load, network, and generator models and numerical parameters were taken directly from [31]. Accordingly, each generator was modeled as a 6 th order synchronous machine with 3 rd order automatic voltage regulators and power system stabilizers; additionally, each generator was outfitted with a 3 rd order turbine governor (Type I) [32]. Each generator system, with its three controllers, had a total of 15 dynamical states. . All time domain simulations were performed using MATLAB's DAE solver ode23t by setting relative and absolute error tolerances to 10 −7 and 10 −8 , respectively.
Load Perturbations: In order to mimic ambient load fluctuations, we applied an Ornstein-Uhlenbeck (OU) process [33] to the active and reactive power demands at each load. The dynamics of these processes are given by τu p = −u p + η p and τu q = −u q + η q , where η p , η q are zero-mean Gaussian variables, and the "load reversal" time constant τ was set to 50s. We further applied a low-pass filter (LPF), with a cutoff frequency of ∼ 9Hz, to the OU variables such that u = LPF{u}. This filtering operation was applied because dynamics above this frequency range become inconsistent with the quasi-stationary phasor approximation used in modeling the network's dynamics. Therefore, high frequency load behaviour is effectively neglected. Finally, these filtered OU variables were parameterized with time variable t using a cubic spline interpolation and applied to the individual ZIP loads via where V ≡ V/V 0 . As a final step, scalar variable β in (42)-(43) was experimentally tuned to a numerical value of β = 50; this tuning generated load perturbations whose corresponding network voltage perturbations approximately matched those of real PMU data (in terms of signal strength). We then applied RTVF to data measured at the machinenetwork interface in order to model the generator's closedloop dynamics. For modeling purposes, we treated voltage magnitude V(t) and voltage phase θ(t) signals as inputs, and we treated current magnitude I(t) and current phase φ(t) signals as outputs. The RTVF algorithm sought to generate a MIMO model with P = 2 and various reduced orders N .
In order to validate the quality of the model generated in the absence of measurement noise, we refer to the time and frequency domain references provided by the exact machine equations. The results are provided in Fig. 3. The results show the model accuracy is excellent in both the time and frequency domains, even though the reduced model order (in this case N = 9) is less than the machine's true model order (N = 15).
Measurement Noise: In order to further gauge the practical effectiveness of RTVF, we applied measurement noise to the voltage and current signals measured by the PMU. To apply
Top and middle panels of Fig. 4 report the frequency-and time-domain fitting performance of a RTVF model of order N = 7 obtained for SNR = 32dB. The corresponding noisecorrupted voltage magnitude signals are depicted in the bottom panel. Compared to the performance in the noise-free setting (Fig. 3), these results show that the frequency-domain model accuracy is still quite acceptable, and that the accuracy in the time domain seems to be not affected by the presence of noise. Therefore, we conclude that the time prediction capabilities of RTVF models extracted from noisy signals are potentially adequate for power system applications.

VIII. CONCLUSIONS
This work proposes a new procedure, termed Real-Time Vector Fitting (RTVF), which performs real-time predictive modeling of linearized dynamics. While the RTVF algorithm can be applied across a broad range of dynamical engineering systems, it was specifically developed to perform real-time predictive modeling of generator dynamics in the presence of ambient perturbations. Accordingly, RTVF explicitly accounts for the presence of initial state decay in the generator's observed outputs. Also, RTVF can be applied in the face of concurrently active input signals, since this is a practical constraint for any three-phase power system component modeled using a quasi-stationary phasor approximation at its terminals.
Given the widescale deployment of PMUs and the newly mandated MOD-33 modeling/validation ordinances, RTVF provides power system operators with a valuable tool. This tool, for example, can build real-time models of generator systems whose physical prior models are completely unknown. It can also validate or enhance prior models which are supposedly well known. Additionally, since controller parameters (e.g. droop gain) are often changed at the local level without system operator awareness, RTVF can be used to infer the true value of these parameters in real-time. Such capability allows system operators to maintain small-signal stability and ensure that generators are properly adhering to market regulations. To enhance practical effectiveness, future research will couple RTVF with advanced measurement noise filtering mechanisms. Also, physically-aware regularization techniques will be developed to allow operators to more effectively "track" a system's shifting dynamics as operating equilibriums evolve over time.