Locating the source of forced oscillations in transmission power grids

Forced oscillation event in power grids refers to a state where malfunctioning or abnormally operating equipment causes persisting periodic disturbances in the system. While power grids are designed to damp most of perturbations during standard operations, some of them can excite normal modes of the system and cause significant energy transfers across the system, creating large oscillations thousands of miles away from the source. Localization of the source of such disturbances remains an outstanding challenge due to a limited knowledge of the system parameters outside of the zone of responsibility of system operators. Here, we propose a new method for locating the source of forced oscillations which addresses this challenge by performing a simultaneous dynamic model identification using a principled maximum likelihood approach. We illustrate the validity of the algorithm on a variety of examples where forcing leads to resonance conditions in the system dynamics. Our results establish that an accurate knowledge of system parameters is not required for a successful inference of the source and frequency of a forced oscillation. We anticipate that our method will find a broader application in general dynamical systems that can be well-described by their linearized dynamics over short periods of time.


INTRODUCTION
The power grid is indubitably one if not the greatest engineering achievement of the past century, as recognized by the National Academy of Engineering [1]. It is an intricate and complex network that has the size of a continent with thousands of constituents that require constant supervision. An important aspect of this surveillance is to ensure that voltage frequencies remain within a narrow band (50 ± 0.05Hz in Europe [2] and 60 ± 0.05Hz in the U.S. [3]), as the violation of this requirement can cause significant damage to vital assets and results in blackouts [4]. Voltage frequency fluctuations are primarily caused by real-time imbalances between production and consumption, such as variations in the charging of electric vehicles, or a sudden gust at a wind farm. The power grid is designed to damp these random electromechanical oscillations appearing during standard operating conditions before it creates a resonance with one of the normal modes of the system. However, malfunctioning equipment or abnormal operating conditions can cause periodic disturbances that would persist over time, creating an undesirable transfer of energy across the system, an effect referred to as forced oscillations.
Whereas most forced oscillations are localized to a particular area, some may be close in frequency to one of the dominant normal modes, resulting in a system-wide response and significant energy transfers [5]. Potential impacts of these wide-area oscillations include equipment failure, inadvertent tripping or control actions, and problems with the automatic generation control. This is why fast and reliable location of the source of forced oscillations is crucial in ensuring the safety and reliability of power grids. However, it remains an outstanding challenge, even when forced oscillation events are detected on the network. For instance, on January 11, 2019, a forced oscillation event happened across the entire Eastern Interconnection on the U.S power grid that was promptly noticed by the reliability coordinators. Nonetheless, existing tools were ineffective at identifying the source location, and a wide-area operator action did not contribute to mitigating the event [6]. The root cause was later fortuitously identified as a faulty input from a steam turbine at a combined-cycle power plant in Florida which forced the system to oscillate for around 18 minutes before local plant personnel removed the unit from service. The forced oscillations created by the faulty turbine had a peak-to-peak amplitude of 200 MW at the generating unit, with power swings of about 50 MW observed as far as the New England area [6], which shows how these disturbances have the power to affect an entire continent. In the case of the November 29, 2005 Western American Oscillation event, a forced oscillation with amplitude 20MW originating from Alberta in Canada created a resonance effect across the entire Western Interconnection. This led to oscillations of amplitude 200 MW registered on the California-Oregon Interface, thousands of miles away from the source [7]. This ability of forced oscillations to cause disturbances at long distances and to be amplified by the grid dynamics seriously complicates the search of their source. Forced oscillations pose a permanent threat to the power grid with more than 20 large-scale events in the past 30 years documented in the U.S., with some of them still lacking a well-identified root cause [8].
The increasing deployment of time-synchronous and distributed frequency sensors in the grid, such as Phasor Measurement Units (PMUs) [9], presents an opportunity for developing advanced data-driven and automated detection and localization techniques of forced oscillations. In the majority of cases, detection of forced oscillations poses a limited challenge as they can be directly ob-served from the Fourier peaks in the signal spectrum [10]. However, forced oscillations need to be differentiated from weakly damped normal modes, or free oscillations [11][12][13]. Weakly damped modes are typically analyzed through Prony analysis [14], and are mitigated via preventative measures by power system operators. On the other hand, the localization of forced oscillations constitutes a much tougher challenge. A complete and perfect knowledge of the system and its dynamics would allow one to locate the source of a forcing [15][16][17]. However, an accurate instantaneous knowledge of the power grid dynamics appears as too strong of an assumption given that the system parameters can fluctuate on the scale of tens of minutes due to local feedback control or temperature variations [18], while the details on the system topology may be unavailable outside of the zone of responsibility for reliability coordinators. Different methods have been proposed to circumvent this lack of information about the system. Some techniques are based on local physical properties such as the monitoring of front arrival times [19], the evaluation of energy flow [20, 21], signal decomposition [22], or verification of the linear relation between voltage and current by multiple PMUs [23,24]. Unfortunately, these methods can be very sensitive to modeling errors such as an inaccurate assessment of the fluctuation propagation speed, or can fail to localize perturbations that are amplified by normal modes. Blackbox machine learning methods [25-27] have been developed with the aim of being fully model-agnostic, but suffer from a prohibitive requirement in training examples of forced oscillations events. These various shortcomings motivated recent calls from system operators and regulators to develop robust tools for performing forced oscillation analysis and localization [6,28], which led to a further exploration [29,30]. Nevertheless, a correct localization of the source in the case where dominant system modes are excited due to the resonance phenomenon remains an outstanding challenge.
In this paper, we propose a new principled method of detection and localization of forced oscillations which is agnostic to the knowledge of the system topology and parameters, fully capable of identifying the source of distant normal modes excitation, and which does not rely on any offline training. Our method operates within a much broader framework that extends rather universally to any dynamical system that can be well-described over short period of time by its linearized dynamics. This makes it a method of choice for the detection and localization of energy transfers created by small disturbances on a large class of complex networks. The key insight in our method consists in leveraging random frequencies fluctuations naturally present in the system to build in real time an effective dynamic model of the network. The dynamic model identification and source localization are performed at the same time using a principled maximum likelihood approach. We illustrate the performance of our approach on a number of examples where forcing excites the natural system modes and create a resonance phenomenon, as well as on a real-world PMU data set.

CHALLENGES AND SOLUTIONS FOR LOCATING THE OSCILLATION SOURCE
Learning of the dynamics of linear stochastic systems lies at the foundation of our approach. In the ambient regime of fluctuations, a complex system such as power grid with a general non-linear dynamics is typically found in a state which is close to a stable equilibrium. In this ambient regime of small perturbations, the dynamics of the system is well described by a linear stochastic equation corresponding to a celebrated model of coupled harmonic oscillators. Without surprise, this family of models includes the popular swing equations describing the ambient dynamics of generators in a transmission power grid, which are derived under the assumption of small deviations from the steady state [32, 33] (see Supplementary Information, section S1 [31] for more details).
We consider the linear stochastic network dynamics with an additional forcing at a single node l, indexed by e l , the canonical basis vector with nonzero lth component. Mathematically, this dynamics is described by the following linear stochastic differential equation, where x t represents the network state variables (for a power grid, they correspond to deviations of phases from the steady state values), p t dt = dx t is the generalized momentum (deviation of frequencies from the steady state values in a power grid), M , D, L are generalized mass, damping, and network coupling parameters, γ, f , and ϕ are the amplitude, frequency and phase of the forcing, respectively, and W t is a Wiener noise process describing random fluctuations. Our goal is to reconstruct the oscillation frequency and the location of the forcing source from the measured time series of x t and p t of length T at a sequence of N discrete time steps t ∈ {t 1 , ..., t N }. We do not assume any knowledge of the system parameters or topology, which represents a realistic scenario in power grids: instantaneous awareness of system parameters are almost never available to system operators [18]. Hence, we do not suppose any knowledge on the inertia M , damping D, or grid Laplacian matrix L. Obviously, we also do not assume any knowledge of the parameters related to the forcing, i.e., γ, l, f , or ϕ. All these parameters of interest need to be recovered from the noisy data {x t } and {p t }.
In Figure 1, we illustrate with a toy example of a threenode network the key mechanism that renders the localization of forced oscillations elusive to standard signal processing analysis. Namely, a forcing with a frequency close to a natural frequency of the system may excite the natural modes of the system peaked at a different node. As a result, neither the correct forcing frequency nor the source of the oscillations are evident from the Fourier spectrum in Figure 1 (d). This effect is reminiscent of forced oscillation events in power grids such as those of November 29, 2005 and January 11, 2019 discussed above, where large perturbations can be seen far from the source and at a very different frequency.
To address this challenge, we develop a principled method for determining the oscillation frequency f and locating the source of forced oscillations based on a maximum likelihood approach (see Methods and Supplementary Information, section S3 [31] for a detailed derivation). A naïve real-space estimator leads to a complicated non-convex optimization problem in frequency, as we explain in the Supplementary Information, section S4 [31].
A key insight which leads to an efficient solution consists in realizing that the finite length of the time series imposes a finite resolution on the frequencies. This leads to a tractable formulation featuring the Fourier transformed quantities, which can be efficiently solved with the stateof-the-art interior-point method solvers (see Methods).
Knowledge of the number of sources leads to a discrete formulation of the problem, where optimization is run for every node (or every pair of nodes if two sources of forcing are present, etc.). We refer to this method as to the System Agnostic Localization of Oscillations (SALO) algorithm. SALO approach is fully parallelizable over all nodes in the network, making it the method of choice for a high performance computing system. An example of an application of the SALO framework is given in Figure 1 (e) for our toy resonance example: both the source and the frequency of forced oscillations are unambiguously and correctly identified. For streaming applications where a faster identification is desired, we consider a computationally advantageous relaxation of the problem. In this SALO-relaxed version of the algorithm, all nodes are formally allowed to be a source with a respective amplitude γ i for each node i (see Methods for more details).
In what follows, we benchmark the method on a number of simulated and real use cases.

TESTS ON SYNTHETIC SYSTEMS AND REAL DATA
An effect where interaction between the forcing frequency and one of the natural modes of the system leads to a peak response away from the forcing source may arise in networks with a much more complex structure compared to the three-node system illustrated in Figure 1. Such a resonance behavior represents an outstanding challenge for detection algorithms whereby the peaks in the Fourier spectrum may not only involve nodes far away from the source, but also point to frequencies related to the natural modes of the system rather than to the frequency of the forcing. We showcase this phenomenon on a synthetically generated data set according to the model (1) on a network inspired by the UK highvoltage grid, see Figure 2. This test case demonstrates some of the features observed in the 2005 Western Interconnection and in the 2019 Eastern Interconnection oscillation events. In particular, the largest amplitudes in the Fourier components can be located very far away from the source, at distances comparable to the diameter of the network. On the other hand, we see that the maximum likelihood based SALO method identifies the correct forcing frequency and the correct source, for which the Fourier signals are otherwise completely hidden among the responses of other nodes in the network, see Figure 2 (b), (c), and (d).
In the Supplementary Information, section S5 [31], we further illustrate the challenges of Fourier-based source localization under the resonance phenomenon on a standard IEEE test case topology with 57 nodes. Yet in this case again, as shown in the Figure S3, SALO method precisely and unequivocally identifies the correct frequency and location of the forcing source, without exploiting any prior knowledge on the system topology and parameters. We also used this test case to demonstrate the performance of the SALO-relaxed version of the algorithm, which accurately points to the correct location of the source of forced oscillations, while benefiting from the computational complexity of a single source verification under the full maximum likelihood approach. This shows that the SALO-relaxed version can be used for a quicker assessment of the forced oscillations once they have been detected in the system, prior to running parallelized computations under the SALO framework. We further show this computational advantage on a series of synthetic instances of increasing size in Figure S4, whereby the ratio of run-times of SALO and SALO-relaxed scales linearly with the system size.
In the tests described above, we assumed that the time series have been produced using the model (1), albeit the system parameters are not known to the reconstruction algorithm. In particular, it is assumed that there is a single and observed source node l, and that the forcing of the type cos(2π(f t + ϕ)) is associated with a single frequency f . In Figure 3 we look at the results produced by SALO in situations where these assumptions are violated. In Figure 3 (a), the source node is outside of the observable system. We see that in this case, the SALO algorithm is still able to correctly identify the forcing frequency, and points to the immediate neighbors of the source node inside the observable system. The Figure 3 (b) demonstrates the case where two sources of oscillations at different nodes and frequencies are simultaneously present in the system. Notably, two peaks corresponding to both forcings appear in the rescaled likelihood score. Finally, in Figure 3 (c), we show the rescaled likelihood scores for the input forcing signals of different types. Even in this case, SALO correctly identifies the source location showing up as several peaks in the likelihood at different harmonics of the input forcing signal. This makes the algorithm remarkably robust to the assumptions behind the model. In the Supplementary Information, section S6 [31], we show an example of an application of the SALO algorithm to real data which display a combination of features observed in Figure 3.
Another possible misspecification is the assumption of linearity of the dynamics. For instance, real power systems do not exactly follow the linear model of the type (1) with constant system parameters. Instead, a linear swing model which falls within the class (1) as discussed in the Supplementary Information, section S1 [31], represents an approximation to a general non-linear dynamics of generators under small deviations from the steady state, valid over finite periods of time [18,34]. However, with the understanding that (1) only serves as an effective model providing an adequate description of the complex system dynamics over short time scales, we show that the SALO method is still applicable to data from real transmission power grids with the purpose of identifying the conditions pertaining to forced oscillations. Here, we benchmark the SALO algorithm on an instance of real PMU data from a U.S. transmission power grid with 200 nodes, where a presence of sustained oscillations has been suggested in previous studies. A Fourier-type analysis of the PMU data from a U.S. Independent System Operator showed a presence of sustained oscillations with a frequency in the 4-6[s −1 ] range, responsible for the emergence of correlations between several node clusters [35]. The analysis of these time series with the SALO method, see Figure 4, confidently points to a single source of oscillations, with a frequency close to the range previously identified in [35]. Incidentally, although the ground truth for this system is unknown, the identified node is consistently pointed to as the most likely source of sustained oscillations even for PMU data separated by a time interval of about a month, see Figure 4 (b) and (c). An extended range of candidate frequencies is likely to be connected to the fundamental limits on the data resolution, as exemplified in the Supplementary Information, section S4 [31], where shorter time series lead to a wider log-likelihood objective function in the frequency domain (see Figure S1).

CONCLUSIONS
We proposed a rigorous maximum-likelihood-based framework enabling a simultaneous system identification and localization of the source of forced oscillations at the same time, without any prior knowledge of the system topology or parameters. In particular, our method is able to perfectly locate the oscillation source, even in the case where the forcing excites one of the natural modes of the system and creates an amplitude peak at a faraway node and at a different frequency. This scenario is reminiscent of some of the real-world historical events, such as the 2005 Western Interconnection event, which turned out to be the most challenging from the oscillation source localization perspective. The ease of parallelization and a relaxed version of the algorithm makes the method scalable to large network instances and multiple sources, while robustness to modeling assumptions makes the algorithm applicable to data produced by real-world dynamical systems.
The SALO algorithm can be naturally adapted to the situation where additional prior information is available. For instance, the matrices M , D, and L do not need to be reconstructed from data if a prior knowledge on the network structure and parameters is available. In the Supplementary Information, section S7 [31], we show that under this scenario, the SALO algorithm is able to identify the correct source and frequency of the forcing using significantly less data compared to the most challenging setting considered in this work, where no prior information on the system is available. Similarly, knowledge of the disturbance type or number of sources can also be directly incorporated into the algorithm. In the context of power grid applications, it could be practical to extend the SALO method to the case of partial sensor coverage, as well as to a source localization under the model which specifically takes into account different types of generators, buses, and PMUs in the grid. For instance, the primary setting considered in our work assumes that all buses have non-zero inertia, which is true for generators in the power grid, but not necessarily for other buses, e.g. those representing loads.
Due to a wider applicability of our approach to general stochastic linear dynamics and coupled harmonic oscillators, we anticipate that our methods will find a broader range of applications beyond power grids, e.g., vehicular platoon subjected to malfunctioning elements or malicious attacks, and more broadly multi-agent systems where automated units must reach an overall consensus [36], as well as forced oscillations in wave propagation dynamics. (1) can be reformulated as a first-order linear dynamic system

Model. The dynamic equation
where we used p t dt = dx t , and denoted X t = (x t , p t ). Discretizing Eq. (2) with an Euler-Maruyama approximation scheme over N points yields the finite difference equation, for j = 0, ..., N − 1, where ∆ tj = (X tj+1 − X tj )/τ for τ = T /N , we assumed t j = jτ , and the frequency relates to the integer 0 < k <  N/2 as k = f T . As discussed in the problem formulation above, the measurements {X jτ } j=0,...,N −1 are assumed to be available, but the parameters of the system A, γ, l, f , and ϕ are unknown, and need to be estimated from the data.
Reconstruction algorithm. Given that on finite time intervals, the contribution of the Wiener process is an i.i.d. Gaussian variable with zero mean, i.e., ξ tj ∼ N (0, τ ), the negative log-likelihood reads, Note that a discretized set of frequencies appearing in the forcing term as a result of finiteness in T and N is crucial, because an optimization over a continuous variable f leads to a hard nonlinear optimization problem, as exemplified in the Supplementary Information, section  S4 [31]. For a fixed frequency k and node l, the joint min-imization over A, γ, and ϕ remains challenging to solve as the negative log-likelihood possesses numerous local minima. Nevertheless, the global minimization of (4) over the phase ϕ, as a function of the remaining parameters, can be performed exactly, leading to the following expression for the partial negative log-likelihood: see Supplementary  However, we have observed that this cost function, unlike its naive counterpart (4), seems to always possess a single minimum that can be efficiently found using state-ofthe-art interior point methods. In this work, we used the optimization software Ipopt [37] within the Julia/JuMP modeling framework for mathematical optimization [38]. Note that the optimization is fully parallelizable over k and l. We refer to the minimization of the partial negative log-likelihood as to the SALO method. We also consider an accelerated method referred to as the SALO-relaxed method, which exploits the spatial relaxation of the problem in the variable γ. Under this variant of the algorithm, all nodes are formally allowed to be a source, and their relative likelihood to be a source is encoded by a respective amplitude γ i for each node i (so that the regression now runs over the vector γ): We benchmark this method in the Supplementary Information, section S5 [31], where we also discuss the computational speed-up of SALO-relaxed compared to the SALO estimator.

DATA AVAILABILITY
All data that support the plots within this paper and other findings of this study are available from the authors on reasonable request.

CODE AVAILABILITY
Code implementing SALO and SALO-relaxed algorithms in Julia can be found at the following Github repository: https://github.com/lanl-ansi/SALO. Michael Chertkov from University of Arizona for fruitful discussions.

AUTHOR CONTRIBUTIONS
All authors designed and performed the research, wrote the manuscript, reviewed and edited the paper.

COMPETING INTERESTS
The authors declare no competing interests.

Supplementary information is available for this paper [31].
Correspondence and requests for materials should be addressed to MV.
[20] L. Chen, Y. Min, and W. Hu, An energy-based method for location of power system oscillation source, IEEE

Supplementary information
Appendix A: Linear second-order dynamics and relation to swing equations In this section, we show that the linearized swing equations describing the voltage phase dynamics in the ambient regime [32,33], falls in the class of linear dynamic models [Eq. (1) in the main text] that we consider in this work. When the system remains in the vicinity of a steady-state fixed point, the linearized equations give a fair approximation of the swing dynamics [32, Sec. 14.3] at node i ∈ {1, ..., n}: where θ i and ω i :=θ i are respectively phase and frequency deviations from the fixed point, m i and d i are respectively the effective inertia and damping, b ij is the susceptance of the line between nodes i and j (zero if they are not connected), and η i accounts for a additive disturbances, including noise (typically considered white and Gaussian) and a potential forcing. Aggregating the phases and frequencies in a vector X t := (θ ⊤ t , ω ⊤ t ) ⊤ , we write Eq. (A1) as the stochastic differential equation, for t ∈ [0, T ], where the matrices M and D are the diagnonal matrices of inertias and dampings respectively, L is the Laplacian matrix of the grid, weighted by the susceptances, γ, f , and ϕ are the amplitude, frequency and phase of the forcing respectively, I is the identity matrix of appropriate size, e l is the canonical basis vector with nonzero lth component, and W t is a Wiener process, accounting for unpredictable disturbances. Therefore, this model is a particular case of the dynamic model [Eq. (1) in the main text] considered in this work.

Appendix B: Parameters of the studied synthetic test cases
Three-node test case in Fig. 1 (main text). The dynamic state matrix A (see Eq. (A2)) used for the example of Fig. 1 (main text) is given by: It corresponds to the dynamics of three nodes, symmetrically coupled, with damping parameters varying over two orders of magnitude. Choosing these parameters and selecting a forcing frequency close to the imaginary part of the eigenvalues of matrix A, Eq. (B1) allows one to have the largest amplitude of oscillation at a node that is not the source of the forcing and at a frequency that is not exactly the one of the forcing, but instead corresponding to the imaginary part of the selected system mode.
UK high-voltage network test case in Fig. 2 (main text). The network topology used to generate the time series in the results of Fig. 2 is a model of the UK high-voltage transmission grids. It contains 120 nodes and 165 lines. Its adjacency list is given in Table I.
The time series are obtained by numerically solving Eq. (1) in the main text, including a forcing at one bus.

Appendix C: Derivation of the SALO and SALO-relaxed algorithms
In this section, we provide the derivation of Eq. (5) from the Methods section of the main text. Starting with the negative likelihood expression (4), we have Define the discrete Fourier transforms of the data Then where [Y ] l is the lth component of a vector Y . Let us consider each of the terms in this expression separately. In (C5), the term 1 N N −1 j=0 ∆ ⊤ tj ∆ tj is independent of the variables being optimized, and hence can be dropped from the optimization. Using the definitions the relevant terms in (C5) can be simply written as Tr(A ⊤ AΣ 0 ) − 2Tr(AΣ 1 ). Regardless of the values of ϕ in (C6), it can be shown that this term is simply equal to 1 2 γ 2 . Indeed, where in the last line we used the fact that k is an integer. Finally, notice that the last term (C7) is the only one that depends on the phase ϕ, and the minimization over the phase ϕ can be performed independently of other variables. The minimum of (C7) over ϕ is given by where M l,· denotes the lth row of the matrix, and we defined Therefore, using all the transformation above and defining a new objective function where minimization of L A, γ, l, k | {X tj } N j=1 over ϕ has been performed, we finally obtain the expression for the SALO objective function given in the Methods section of the main text: In order to obtain the SALO-relaxed version of the algorithm, we notice that the objective (C17) can be rewritten in an equivalent form, where the candidate source node l and the amplitude γ are encoded as the lth component of the one-hot vector γ that spans all nodes in the system: s.t. ∥γ∥ 0 = 1.
A spatially relaxed version of the algorithm over the variable γ is obtained by dropping the one-hot-encoding constraint ∥γ∥ 0 = 1, a version that we refer to as the SALO-relaxed algorithm: In this version of the algorithm, the source location is determined from the indices of the largest components of the vector γ.
Appendix D: Challenges for the optimization of the continuous likelihood function The length and the sampling rate of the available time series impose a natural limit to the frequency resolution, which sets a limit on a precision with which the forcing frequency f can be identified. This idea lies at the heart of the discretization in the frequency domain explained in the Methods section of the main text. Both SALO and SALO-relaxed methods are then evaluated for each discrete value of possible frequencies k. To provide an additional insight into our approach, in this section we illustrate why a direct optimization over the continuous frequency f results in a hard-to-solve non-convex optimization problem. Consider a direct optimization over a continuous variable f of the negative log-likelihood from N measurements {X jτ } j=0,...,N −1 taken at times t j = jτ over the observation window of length T , so that the interval τ = T /N . For simplicity, let us fix variables A, γ, l, and ϕ to their ground-truth values. In this case, we can plot the objective function (D1) for a test problem as a one-dimensional function of f , see Figure 5. We see that with the increasing observation window T , the objective function is oscillating around a certain constant value for most values of f , with a sharp and narrow peak at the correct frequency value. It is possible to understand this behavior using a toy example presented in Figure 6. When parameters A, γ, l, and ϕ are fixed to their ground-truth values, the optimization problem in (D1) becomes essentially equivalent to fitting of a noisy periodic function cos(νt) + ξ t , where ξ t is a white Gaussian noise, with a function cos(f t). For any value of f which is not equal to ν and a sufficiently long time series, a small deviation of f from ν makes the fit eventually diverge to the point where the two curves appear in counterphase. This effect is responsible for the form of the objective function presented in the Figure 5. This type of landscapes with a large number of local minima and a single narrow global minimum is particularly challenging for optimization solvers. This challenge is alleviated in our approach by switching to the discrete frequencies k and turning the minimization into a mixed-integer optimization problem.  The SALO-relaxed algorithm has an improved computational time compared to the SALO algorithm given that it does not need to be run for all possible candidate sources in the network. In this section, we empirically benchmark the SALO-relaxed algorithm against the SALO method. For all the test cases used throughout the main text, we observed that the SALO-relaxed algorithm points to the same source, consistently with SALO. Here, we use additional test cases to illustrate that while the SALO-relaxed algorithm provides an improvement in computational time, it does not sacrifice the accuracy of the source localization.
In Figure 7, we compare the outcome of both methods on the network structure extracted from the standard IEEE 57-bus test case [? ]. Similarly to the test case presented in the Figure 2 in the main text, we modified the parameters of the network to produce the challenging resonance conditions, where the largest response seen in the Fourier spectrum appears on nodes that are far away from the source. Importantly, under these challenging conditions, SALO-relaxed method agrees with the SALO algorithm, and confidently points both to the correct source of forced oscillations, as well as to the ground-truth forcing frequency.
In Figure 8, we run a scaling experiment to compare the computational speed of SALO-relaxed and SALO methods, without using parallelized computation for the latter one. For this purpose, we use random Watts-Strogatz networks of increasing sizes. At every network instance, the source and frequency of forced oscillations is correctly identified by both algorithms, however, SALO-relaxed algorithm enjoys a computational advantage over SALO that scales linearly with the system size.
We provide further results of testing of the SALO-relaxed in the situation where some of the modeling assumptions are violated in the next section F. In the Figure 3 of the main text, we showed the results produced by the SALO method in the case of multiple sources injecting oscillations at different frequencies, potentially with sources lying outside of the observed system. In this section, we apply the SALO method to real data, which shows a combination of features observed in Figure 3 of the main text. We use the data collected under an oscillatory event in the U.S. Eastern Interconnection from FNET/GridEye, a GPS-synchronized wide-area frequency measurement network [? ]. We have analyzed the subset of the time series without missing data. The ground-truth location of the sources of forced oscillations is not known in the considered data set. The network topology or parameters are also not known, which calls for the application of the SALO method. The known locations of a subset of measurement devices is depicted in Figure 9 (a). The analysis of the time series with our SALO method points to the most likely sources of forced oscillation, highlighted with the orange, blue, and purple colors in the rescaled log-likelihood scores in Figure 9 (b). A precise location of the node corresponding to the likelihood score highlighted in orange is not available, however the other likely sources are shown with the blue and purple colors in Figure 9 (a). The features shown in the log-likelihood scores are reminiscent of the ones observed in the synthetic study Figure 3 of the main text. In particular, both the orange and the blue nodes indicate the presence of oscillations at the exact same frequency f = 4.54 [s −1 ]. This feature is similar to the response of the algorithm observed in the case where the source is located outside of the observed system, see Figure 3 (a) in the main text. Although the geographical location of the orange node is not available, the blue node is indeed located at the edge of the observed area, which makes this scenario plausible. The analysis with SALO also points to a presence of the second source of oscillations at a different frequency f = 1.35 [s −1 ], similarly to the algorithm response observed in a synthetic case with two oscillation sources shown in Figure 3 (b) of the main text. The respective source node found by SALO is located inside the bulk of the system, as shown in purple in Figure 9 (a). In Figure 10, we use a synthetic network example to verify that the combination of these features -two sources of oscillations, inside and outside the observed systemindeed leads to a picture corresponding to a combination of features observed in Figures 3 (a) and 3 (b) of the main text. Interestingly, the nature of the leading peaks in the rescaled log-likelihood observed in Figure 10 shows a lot of resemblance with the results produced by the SALO algorithm in Figure 9 (b). FIG. 10. Rescaled log-likelihood score obtained by SALO for the scenario with one observed and one hidden sources of forced oscillations. The algorithm correctly identifies the location and forcing frequency for the visible source, and points to the neighbors of the hidden source inside the visible system for the detected oscillation at the other forcing frequency. As before, the envelope of scores for non-highlighted nodes is shown in gray.
Finally, we test the performance of the SALO-relaxed algorithm under the scenarios of violated modeling assumptions studied in Figures 3 (a) and 3 (b) of the main text, as well as for the combination of hidden and multiple sources examplified in Figure 10. The results of these tests are reported in Figure 11. In all three cases, the outcome produced by the SALO-relaxed algorithm is in complete agreement with the results shown by SALO. . Rescaled log-likelihood and amplitude vector inferred using the SALO-relaxed version of the algorithm for three scenarios of model mismatch, illustrated on a network with 20 nodes. In the top panel, forced oscillation is injected at the unobserved gray node (a). SALO-relaxed identifies the correct forcing frequency (b) and points to the immediate neighbors (c) of the gray node as the most likely sources. The center panel studies the situation where two independent forcings are present at different nodes (e). For two frequency peaks identified in the log-likelihood score (f ), we show two respective amplitude vectors, with the largest components pointing to the correct sources. The bottom panel shows the combination of the previous two scenarios, with two sources, one observed (purple) and one unobserved (gray) (g). Both forcing frequencies (h) and the respective locations from the largest components of the amplitude vector (i) are correctly identified. All results are consistent with the outcome produced by SALO, see Figure 3   Rescaled log-likelihood score obtained by the SALO algorithm from time series data with three different lengths, illustrated on a network with 20 nodes (a) with the forced oscillation injected at node 20 at a frequency about 0.08 [s −1 ]. For shorter time series with length of (b) 10 and (c) 23 [s], the system-informed version of SALO algorithm is capable of correctly identifying the forcing, whereas the system-agnostic version requires more data since it simulataneous needs to estimate the parameter matrix A. The time series with length of (d) 39 [s] are enough for both the system-agnostic and system-informed versions of the method to recover both location and frequency of the forced oscillation. Rescaled log-likelihood score and location amplitude vector obtained by the SALO-relaxed algorithm from time series data with three different lengths, illustrated on a network with 20 nodes (a) with the forced oscillation injected at node 20 at a frequency about 0.08 [s −1 ]. For shorter time series with length of (b) 10 and (c) 23 [s], the system-informed version of SALO-relaxed algorithm is capable of correctly identifying the forcing, whereas the system-agnostic version requires more data since it simulataneous needs to estimate the parameter matrix A. The time series with length of (d) 39 [s] are enough for both the system-agnostic and system-informed versions of the method to recover both location and frequency of the forced oscillation.