Locating the Source of Oscillation With Two-Tier Dynamic Mode Decomposition Integrating Early-Stage Energy

This article presents an improved method to locate the sustained oscillation source using measurements, e.g., voltage measurements from phasor measurement units (PMUs). Based on dynamic mode decomposition (DMD), the proposed method exploits the spatiotemporal patterns of dominant dynamics from the measurements to detect the oscillation and identify the source. It significantly enhances the overall performance by incorporating the initial stats as input such that the state matrix represents the post-mortem dynamics. We also propose a two-tier structure to further improve the robustness and accuracy: full-time analysis over the entire time frame of interest to identify the problematic modes and candidate sources and peak-time analysis to identify the actual source. These analyses provide insights into the oscillatory energy distribution among the monitored buses at the early stage of the oscillation and normalized contribution factors pointing to the source robustly. Studies on the simulation data in WECC 179 and 240 bus systems and those on a real case in ISO-NE validate the effectiveness of the proposed method.

these oscillations and their sources (or causes) are increasing due to the rapid increase of inverter-based renewable sources along with power electronics-based transmission or distribution controllers widely deployed in the power system [5], [6].
Although there have been various efforts to mitigate the oscillation events through internal or external controls [7], [8], disconnection of the oscillation source from the power system, if allowed, has been a practical and immediate solution [2], [3], [9].Either way, locating the oscillation source is the first and crucial step of the mitigation strategy.Accordingly, various oscillation source location methods have been proposed, which are well categorized in [3].Among various categories, recent studies are mainly focusing on the damping torque-based method [2] and energy-based method [9], [10], [11], [12], [13].The former method estimates the damping coefficient of the generator model, and the latter one calculates the net energy flow and considers the machine with the largest dissipating energy as the source.Especially, the energy-based method is built upon a solid mathematical formulation and provides credible explanations for the source.It demonstrated good performance on both forced and natural oscillation cases in real power systems [10], [11], [13].However, the energy-based method can mislocate the forced oscillation source under diverse operating conditions, as reported in [9].Furthermore, the energy-based method requires the measurements of all the branches to ensure its accuracy.Although both damping torque and energy-based methods are expected to provide accurate location and support the theoretical reasoning, accurate modeling of the target system should be a prerequisite.Additional information other than the measurements from PMUs to construct the system model should thus be required for the implementation and best performance.
Data-driven or measurement-based methods can overcome the limitations of traditional model-based methods.In [14], robust principal component analysis (RPCA) is adopted to decompose the voltage measurement matrix, and the forced oscillation source is located based on the resonance-free component.Meanwhile, [15] extracts the power system dynamical model from the voltage and frequency measurements using the sparse identification of nonlinear dynamics (SINDy).Both methods show high accuracy in locating the forced oscillation cases only with practical PMU measurements, even if the information about the network topology and model parameters is not given.However, both methods are not proven to be effective in natural oscillation cases nor the forced oscillation case triggered by non-generator sources.
This article notices the dynamic mode decomposition (DMD), an extension of the singular value decomposition (SVD), as an effective tool to capture the system dynamics from the PMU measurements.Although SINDy has been regarded as the upgraded version of DMD, it requires higher model complexity and longer computation time [16].Therefore, we have chosen to adopt DMD instead of SINDy, considering that the location algorithm needs to be implemented online.DMD has the advantage of providing the spatiotemporal distribution of oscillation modes.There have been various efforts to analyze the oscillation events by estimating the mode shapes using the DMD [1], [17], [18].In [18], DMD-based analysis presented both the mode shapes and the electro-mechanical coherency of the generators accurately.The method was shown to be more efficient than Prony and Koopman analysis in the computational sense.Randomized DMD [19] demonstrated higher reconstruction accuracy than the other methods.However, there has not been a sufficient effort to derive a DMD-based location method.Although the DMD can show the overall tendency of the signal, the original DMD-based method may fail to locate sources, especially for the forced oscillations [3].It is also noteworthy that less attention has been paid to oscillation cases such as maximum oscillation magnitude at a non-source bus, rectangular forced signal injection, and the oscillation case caused by a non-generator source, etc., even though it is crucial yet challenging to locate them because their patterns differ from the usual cases.
This article thus proposes a fully data-driven oscillation source location method using DMD, specifically addressing the aforementioned challenging cases as well as the natural oscillation cases.Building on the concept of DMD with control (DMDc) [20], we have reformulated the DMD by substituting the control input matrix of DMDc with the initial state matrix of the measurements, which allows the state matrix to focus on the oscillatory dynamics.Furthermore, the proposed method only uses the voltage magnitude and voltage angle as variables to align with the practical measurements of PMUs.Finally, in order to heighten the location accuracy and adaptability to unusual cases, the method employs a two-tier DMD approach consisting of full-time DMD and peak-time DMD.Full-time DMD decomposes data across the entire time span of interest, providing potential oscillation source insights.Peak-time DMD operates after the oscillation energy at the desired frequency reaches its initial high and low peaks.By synthesizing the high and low peak results, the method assesses the impact of each potential source bus on the early stage of ambient oscillations.The contributions of this article can be summarized as follows: r Introduction of an improved DMD formulation that en- hances the clarity of oscillation dynamics in the state matrix.
r Proposal of online oscillation detection and source location algorithms based on spatiotemporal energy distribution derived from the improved DMD; These algorithms prioritize location accuracy and computational efficiency.r Utilization of a two-tier structure in the source location method to significantly improve accuracy across various oscillation scenarios.Even in cases where measurements are unavailable at the source bus, the method provides valuable information for system operators.The case study results validate the accuracy and feasibility of the proposed method for two different simulation cases and one real oscillation data, which can be downloaded from [21].The objective of the source location is limited to below the substation level.The rest of this article is organized as follows.Section II briefly reviews DMD from the viewpoint of SVD, which is well introduced in [22], and the oscillation analysis through DMD.Section III presents the mathematical formulation of the improved DMD and the source location method using the two-tier structure on full-time and peak-time.Case study results and analysis are given in Section IV.The conclusion and future works are drawn in Section V.

A. Dynamic Mode Decomposition
Dynamic mode decomposition (DMD) has been expressed and interpreted into two perspectives: Arnoldi and singular value decomposition (SVD)-based approaches [22], [23].From the viewpoint of SVD, DMD is the extension of SVD to characterize the nonlinear system without direct computation of A, s.t.X = AX.Since the mathematical formulation of DMD is well represented in [22] and [24], this section briefly introduces the mathematical derivation of DMD from SVD.
SVD is the generalized eigen-decomposition of a square matrix into two unitary matrices and one diagonal matrix.Based on the decomposed structure, key features of the correlation matrix can be drawn from large data; thus, SVD can be regarded as the data dimensionality reduction method [25].For a measurement matrix X, the overall structure of various SVDs, including full SVD, economy SVD, and truncated SVD, is shown in Fig. 1.The mathematical formulation of full SVD can be given by where n is the number of measuring units and m is the number of measurements.Under the assumption that n m > r, economy SVD and truncated SVD can be derived by taking m and r rows of U as Economy SVD : X = Û ΣV * , Truncated SVD : X = Ũ Σ Ṽ * . (2) If we assume that X and X are future and present measurement matrices, an objective system can be linearized as X = AX for the best-fit linear operator A. According to Eckart-Young's theorem [26], the best rank-r approximation of X becomes the rank-r truncated SVD of X, given by argmin X X − X ≈ Ũ Σ Ṽ * s.t.rank( X) = r.(3) If we substitute (3) into the linearized system, X can be represented as Supposing Ã is the projection of A onto Ũ, we can calculate Ã without A using (4) as If we let the eigenvalues and eigenvectors of Ã as Λ and W, i.e., ÃW = WΛ, X can be calculated from (5) as where the DMD mode Φ = ŨW denotes the spatio-temporal coherency and Γ = W −1 Σ Ṽ * represent the temporal distribution of the mode amplitudes.It is worth noticing that the eigenvalues of X also become Λ, and the DMD mode is equal to the projected eigenvectors of A. Each column vector of the measurement matrices can be expressed separately using the spectral decomposition as where k = 1, 2, . .., m.

B. Oscillation Analysis via DMD
Analyzing the oscillation using DMD starts from two assumptions: 1) power system states are dominated by a timeindependent linear system, and 2) power system states can be expressed by the linear combination of previous states for sufficiently long data.These assumptions can work even for the forced oscillation cases since the impact of the forcing signal would have been reflected in the states.By combining the assumptions, we can determine the Frobenius matrix of power system states [18], leading to the estimated power system model in (4).The state matrix A can be estimated through (5); thereby, the state dynamics can be represented using DMD modes achieved from the measurements based on (6) and (7) without examining the power system dynamics.
The states to conduct DMD can be rotor speed, voltage, and current of measuring buses.In [18], the rotor speed of generators is adopted, while voltage magnitude is used in [17].Regardless of which state is used, the power system state at any measuring bus k can be expressed with (7).When we assume that λ j s are the diagonal components of Λ (empirical Ritz values) for j = 1, 2, . . ., n, they exist as complex conjugate pairs for the unstable oscillation cases [27].In addition, the oscillation frequency and damping ratio can be acquired based on the analysis of the empirical Ritz values as [28] The spatial and temporal distributions of the oscillation at f j can be obtained from φ j and γ j .Especially, [18] identifies the coherency of the generators by using the spatial distributions.

III. PROPOSED OSCILLATION SOURCE LOCATION METHOD
The proposed oscillation source location method includes the oscillation detection and the source location, as represented in Fig. 2. The detection algorithm continuously calculates the oscillatory energy of the measurements within a relatively short time span.Once the oscillation is detected, the location algorithm gathers longer data and analyzes the spatiotemporal energy distribution, and determines the oscillation source.This article further improves DMD by including the initial state as an additional input of the dynamic equation so that the oscillatory energy can be estimated more accurately.Both detection and location algorithms adopt the improved DMD, although they have different evaluation matrices and different interesting timeframes.Improved DMD, detection, and location methods are explained in the following subsections.

A. Improved Dynamic Mode Decomposition
In the original DMD, the state matrix A should take every dynamics that occurred in the target system.DMDc puts the current control signals as the additional input to divide the impact of the external control from the system dynamics.The dynamic equation with control signal Υ changes from (4) to Motivated by the concept of DMDc, we improve DMD to focus on the oscillation dynamics itself.Therefore, we use the initial state of the measurement matrix as an additional input instead of Υ.The initial state of the measurement matrix is a n × 1 vector, where X ini = X 1 .The decomposition process of the improved DMD is similar to one of DMDc, where the mapping between the control input and future measurements is unknown [20].Improved DMD can be formulated by replacing Υ of ( 9) by X ini as Assuming that f (k) denotes a function describing the oscillation dynamics for every discrete time of the measurements k ∈ (1, 2, . .., m − 1), we can modify (10) as Since we use the RMS values for the voltage and current magnitudes, ( 11) is reasonable as the RMS value will be constant without oscillation.For . If we substitute this into (11), the matrix A represents the linear relationship between the change of the states and change of the oscillation dynamics: Therefore, the state matrix A can focus on the post-mortem modified dynamics, while B takes the basic dynamics driven by the initial states.Unlike the original DMD, the improved DMD helps analyze the oscillation dynamics aside from the initial state, whether the oscillation is forced or natural.In addition, note that B in (10) plays a different role from the input matrix of the original system state space equation.The overall process of the oscillation analysis using the improved DMD is represented in Fig. 3. Considering that the present PMUs do not measure the rotor speed and angle of the generators, we mostly observe the voltage magnitude and angle, and the current magnitude is also adopted for real oscillation cases.Furthermore, in order to focus on the oscillation signal, the deviations of voltage magnitude and angle are adopted instead of the original data.Therefore, the current measurement matrix can be presented as X = [ΔV Δθ] T ; thus, the input matrix Ω ∈ R 4l×(m−1) for l measured buses.Similar to the original DMD, we compute the truncated SVD of the input matrix as Ω = Ũ Σ Ṽ * .If we substitute the SVD results for Ω, the approximates of A and B can be found as To find the singular vector, which can be utilized to represent the subspace of the system state, another SVD needs to be computed on X .Although truncated-SVD is adopted in DMDc, full-SVD is conducted for the improved DMD since the number of measured buses is generally less than the number of measured timesteps.Therefore, the second SVD result can be represented as X = UΣV * and the best approximations based on both SVDs become According to the eigen-decomposition of Ã, the DMD mode and temporal distribution can be derived as Since we perform the full SVD for X , the rank r should be selected as r = 2l for the number of measuring buses l to reconstruct the matrix X .Based on (15), the future measurement matrix can be decomposed into Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Accordingly, the size of the decomposed matrices can be determined as Φ ∈ C 2l×2l and Γ ∈ C 2l×(m−1) .We can calculate the oscillation frequency and damping ratio based on (8).Temporal distribution of energy for each oscillation frequency can be achieved from Γ, while Φ represents how the energy is spatially distributed among the measuring buses.Therefore, we can determine when the oscillation starts by analyzing Γ and which bus contributes the largest energy to the oscillation from Φ.

B. Online Oscillation Detection
The proposed oscillation detection algorithm continuously performs the improved DMD for online measurements.For the number of measurements for oscillation detection m detect , the future measurement matrix has the form of 2l × (m detect − 1).For the truncated SVD, the length of the measurement matrix should be larger than the rank r; thus, the number of measurements should be selected as m detect − 1 >= r = 2l.Note that the size of the screening window is equal to m detect .The window size needs to be selected as short as possible to guarantee online detection, but it should be long enough to detect the low-frequency oscillation.If we assume that l is sufficiently large enough, the number of measurements for the online oscillation detection is chosen as m detect = 2l + 1.
At some time t and fixed timestep Δt, the overall timeline of the proposed method is represented in Fig. 4(a).As a result of the improved DMD computation from t − (2l + 1)Δt to t, we can get Λ| t and Γ| t .Here, | t denotes that the corresponding values are achieved from the measurements until t.Then, a specific oscillation frequency, f t i , existing at that time can be calculated from λ i | t .If there is no significant oscillation during the time, the norm of the non-oscillatory energy would be larger than about a third of the energy of any other frequencies as follows: ∀i ∈ (1, 2, . .., 2l), If ( 17) is satisfied, the proposed method concludes that oscillation does not exist, and DMD is computed for the next measurements.The number of measurements to be gathered is determined considering the computation time of the improved DMD and timestep of the measurements, neglecting the communication delay.In our study, we update the dataset for every measurement because the average computation time is about 0.025 s, which is shorter than the timestep of 30 Hz measurements.When the computation results do not meet (17), our method concludes that there exists oscillation and moves to the source location process at t detect .Although evaluating the presence of the oscillation based on the norm can delay the discovery of the oscillation, it can prevent the false alarm by data outliers as well as enhance the detection accuracy.

C. Two-Tier Oscillation Source Location Method
After the oscillation is detected in step I, the proposed method moves on to step II-1, gathers the measurements longer than the dataset used for the oscillation detection, and loads the previous measurements.Note that the dataset in orange color is wider than the screening window in blue in Fig. 4(a).The oscillation source location method performs the improved DMD for the extended oscillation source location (OSL) dataset from t c to t f , where the critical time t c < t detect and the final measurement time t f > t detect .For the stable and accurate calculation of the improved DMD, t f should be selected as t f − t detect ≥ 2lΔt, and we mostly choose t f = t detect + 4lΔt in our experiments for the accurate and fast location.The critical time t c should be selected differently according to the existence of the disturbance in the measurements.Since we only take account of the momentary disturbances in this article, disturbances in longer term such as cyclic load, harmonic distortion are not considered.The disturbance type of interest includes momentary sag, swell, and surge caused by short circuits, major equipment shutdown, etc.These disturbances appear as impulse signals, distorting the dominant frequency and corresponding oscillatory energy.Fig. 4(a) illustrates the fault occurrence and clearing as a representative example of the momentary disturbance.In the case of the oscillation with the disturbance, the critical time is chosen as the time when the disturbance signal is removed; thus, the critical time for the oscillation with disturbance, t c,d , is equal to the fault clearing time in Fig. 4(a).On the other hand, if there is no disturbance in the measurements, the critical time should be selected based on the detection time.Similar to the final time, the critical time without the disturbance is set as t c,nd = t detect − 4lΔt so that the extended dataset can give sufficient information for the DMD analysis.
Then, the improved DMD is computed for the full-time of the extended OSL dataset.Although the length of the dataset increases, the number of Ritz values remains the same with 2l.Assuming that i is one of the Ritz values, i ∈ (1, 2, . .., 2l), we can examine the total energy of each oscillation frequency for the full time from ||γ i | t f ||.Note that | t f represents that the norm value is achieved from the measurements until t f .According to the magnitude of norm values, the dominant component of the oscillation, d, and the according dominant oscillation frequency, f d , can be determined.Subsequently, the norm value of DMD mode ||φ d || represents how much energy each bus contributes to the dominant oscillation.Therefore, the buses with a large value of ||φ d || would be the source candidates of the oscillation.However, because the contributions on both voltage magnitude and angle should be assessed and the phase of the mode needs to be considered, the oscillation contribution factor (termed as OCF in this article) for the full-time DMD is organized as follows: where the subscript ft means the full-time, and the mode magnitudes are separated into MM d,v and MM d,θ .The former magnitude corresponds to one of the upper half of φ d , while the latter takes the remains.In this article, the dot product • represents the element-wise multiplication of the matrices.Aside from whether the mode phase is leading or lagging, the buses having large OCF ft values are considered as the oscillation source candidates.
Although DMD shows high accuracy in the reconstruction of the measurement data, it has a limitation in accurately determining the oscillation sources.Because the full-time DMD shows the participation over the whole time, it cannot estimate the temporal participation of each bus.Furthermore, due to the inborn characteristics of DMD, the length of the measurement matrix should be larger than the height, which prevents DMD from observing only the early part of the oscillation.In other words, the DMD-based method cannot guarantee location accuracy if there is a non-source bus whose oscillatory magnitude is larger than the source bus.In order to fully exploit the strength of DMD and accurately locate the source, this article proposes a two-tier oscillation source location method, which can be divided into the full-time DMD and the peak-time DMD.As mentioned earlier, full-time DMD handles the whole-time measurements after removing the oscillation trigger, e.g., fault clearing and injection of the forcing signal.On the other hand, peak-time DMD is performed for the measurements after the oscillation energy of the dominant frequency reaches the initial high and low peak values.Fig. 4(b) denotes the illustrative example of the two-tier DMDs and how to select the peak times.By comparing the results from high and low peak time DMDs, the energy and contribution factor of the potential sources at the early stage of the oscillation can be observed.Furthermore, notice that the temporal energy distribution of the dominant frequency is generally shaped as the blue line in Fig. 4(b).The temporal energy ||γ t d | t f || initially increases to a certain level, repeating increase and decrease since it represents the energy flow of the entire observed system, not the energy of the source bus.Therefore, even if the forcing signal is injected continuously into the system, the energy exhibits a fluctuating pattern, rather than a continuous increase, because other power system components, such as loads and generators, primarily respond to the oscillation by reducing its energy.After the energy takes the highest level, it gradually stabilizes to the lower level.The early stage is mostly formed before the energy reaches the highest energy, within 3 s from the critical time in our experiments.
The step II-2 process of the proposed method can be expressed as follows.As the results of the improved full-time DMD discussed in the previous subsection, the potential oscillation source buses with the number of p can be selected as k = [k 1 , k 2 , . . ., k p ] for the dominant frequency, f d .In order to choose the initial peak time, the temporal oscillation energy needs to be considered instead of the averaged energy ||γ d | t f ||; thus, the magnitude of energy at each timestep, ||γ t d | t f ||, is considered.As represented in Fig. 4(b), the peak times are chosen to represent the first swing of the temporal energy for the dominant frequency.The high peak time can be easily selected when Δ||γ t d | t f || becomes negative for the first time or right before the time.Then, the low peak time is chosen as the lowest point of ||γ t d | t f || until the energy magnitude exceeds the high peak level.However, all the energy norm does not have the same form as the general shape presented in Fig. 4(b).Therefore, we provide the guidance for the selection of the peak times as follows: r The energy norm of the non-disturbance oscillation cases rapidly increases from t c , similar to the general shape of the energy norm.The high peak time can be chosen around the time when Δ||γ t d | t f || becomes negative for the first time.The low peak time is selected at the lowest point of energy until the energy surpasses the energy at the high peak time.Notice that Δ||γ t d | t f || changes from positive to negative at the high peak time and vice versa at the low peak time.
r If the time span of the first swing is not sufficiently long, the high and low peak times are determined in the next swing of the energy.In this article, the first swing is neglected when the time, which meets Δ||γ t d | t f || < 0 until the energy exceeds the high peak time energy, is shorter than twice the time step of the PMU.
r Unlike the non-disturbance oscillation, the energy norm of the oscillation with disturbance takes various form case by case.Specifically, the energy norm is highly influenced by the voltage recovery process after the fault.Hence, the energy norm decreases below the initial energy level for the first several times as the impact of the fault decreases.
r If the energy norm records the lowest level in a short time, mostly 1-2 seconds, and re-increases after the impact of the fault is reduced, the peak times are determined after the lowest level.However, if the energy reduction maintains longer, the peak times should be determined before the energy norm reaches the lowest level.The practical choice of peak times is further discussed in the next section with examples.
After the peak times are determined, the high and low peak time DMDs are computed for the measurements from t hp to t f Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and from t lp to t f .For computational efficiency, the peak-time DMD oscillation contribution factors are only calculated for the buses included in k.To equivalently evaluate the participation on both high peak and low peak, the contribution factor of bus k should be normalized as It is worth noticing that the contribution factors of high and low peak-time DMDs are determined by the subtracted results between the full-time and peak-time mode magnitudes; thus, the factors can assess the contribution from t c to t hp or t lp , the ambient stage of the oscillation.This helps the proposed location method to discover the proper source by neglecting the other subsequent effects like the control reaction.By aggregating the high and low peak-time factors, the proposed final oscillation contribution factor can be given by where the subscript tt denotes the two-tier.The reason why both high and low peak contribution factors are used is to check which bus consistently intensifies the oscillation at the early stage, irrelevant to whether the total oscillation energy of f d increases or decreases.As the bus has a larger value of OCF tt , the larger and more consistent oscillation energy would be radiated from the bus.Therefore, the bus with the largest OCF tt is considered the oscillation source of the dominant frequency.
Remark: The mode magnitudes are calculated from the norm of DMD modes for the interested time from t c or peak times to t f .Therefore, the mode magnitudes represent the cumulative energy throughout the interested time, and OCF hp,lp values denote the energy accumulation during the early stage of the oscillation from t c to peak times.We design OCF tt equal to the multiplication of two normalized peak OCF values so that the index can focus on whether each bus consistently contributes on the oscillation, even when the total oscillation energy decreases.

A. Data Description
The proposed source location method is applied for the simulated cases (Sim cases) on WECC 179 bus system, the 2021 NASPI contest cases (NASPI cases) on WECC 240 bus system, where all the cases are conducted by TSAT simulation [21].In order to validate the effectiveness of the proposed method in the presence of noise and outlying data, a real oscillation case in ISO-NE has been further adopted.According to the characteristics, we initially select the representative cases among the simulated and NASPI cases and divide them into five categories as follows:  I .In addition, the estimated frequency and source through the proposed method with its elapsed time are represented.Oscillation location is simulated on MATLAB R2022a environment, implemented on a computer system with Intel Core i7-10700 2.90 GHz CPU and 32 GB RAM.Corresponding results on each category will be provided in the following subsection to validate the proposed method.Then, the validated method is applied to the real oscillation data in the one after the next subsection.

B. Application on the Simulated Oscillation Cases
As the first example, the proposed method is applied to normal condition cases.The first case, Sim 2F, is the forced oscillation case in which the oscillation is caused by the forcing signal without the system disturbance, while the others are the natural oscillation cases, where the short circuit or fault provokes the oscillation.For Sim cases, the truncation ranks r is 59, considering the number of measuring buses.Based on the proposed detection algorithm, the oscillatory energy of a specific frequency, 0.87 Hz, exceeds the stationary energy at 5.8 s.Therefore, the source location algorithm activates at t f =17.5 s for the data from t c = 0 to t f .The overall process of the location method on the Sim 2F case is shown in Fig. 5.As the result of the full-time DMD, the dominant frequency is determined as 0.8567 Hz.Then, the high and low peak times are identified through the temporal oscillation of the dominant frequency.As represented in Fig. 5(b), the high peak time can be chosen when Δ||γ d | t f || becomes negative for the first time and the low peak time as it returns to the positive.In addition, the buses whose OCF ft values are two times larger than the average value are selected as the potential oscillation sources, as k = [30, 35, 70, 77, 79].After the computation of the high and low peak time DMD, it can be observed that the frequencies of the mode with the highest magnitude have similar values with 0.8557 and 0.8555 Hz each.Finally, considering the magnitude of OCF tt , bus 79 is decided as the oscillation source.Since the location algorithm is computed for 0.41 s, we can find out the source at 17.91 s by adding the computation time and t f .
The simulation results on NASPI 7 and Sim 6ND cases are represented in Fig. 6.In NASPI 7 case, the fault is cleared at 27.06 s, which we set as t c .Since the NASPI cases have more measuring buses than Sim cases, r is set as 116, and the corresponding t f is about 67 s.In the Sim 6ND case, the fault clearing time is 0.05 s with t f = 21 s.The dominant frequency can be easily determined as 0.3799 Hz based on the norm of the oscillation energy.Similarly, 1.402 Hz shows the largest norm in the Sim 6ND case.Both cases have higher ||γ d | t=0 t f || than the next few seconds because the signal is affected by  the fault clearing.For these cases, we select the peak times after the energy norm records the minimum value.Based on OCF ft values, the potential oscillation sources are selected as k = [6553, 6404, 2634, 2604, 2438, 6507] for NASPI 7 case and k = [36, 159, 4, 9, 45] for Sim 6ND case.Notice that two buses show comparable OCF tt values for both cases; thus, we need to consider the possibility of a multiple-source oscillation case.In NASPI 7 case, bus 2634 can be determined as a single source bus since both nOCF hp and nOCF lp values of bus 2634 are higher than those of bus 2604.Bus 2604 records a high OCF tt value because it is directly connected to bus 2634.On the other hand, bus 45 has a larger nOCF hp value than one of bus 159, while nOCF lp value is lower.Furthermore, buses 45 and 159 are not adjacent since bus 159 belongs to zone 1-B and bus 49 to zone 2-B [29].Therefore, the Sim 6ND case is the multiple sources oscillation case whose sources are bus 45 and 159.It is worth noticing that network configuration and components should help validate and reinforce the performance and accuracy data-driven method.For example, information about the network topology can help to distinguish the single and multiple source oscillations with certainty.
Then, we applied the method to the second category, where the maximum oscillation magnitude exists on the non-source bus.The high and low peak times of all cases in the second category are shown in Fig. 7(a).According to the basic criteria, the initial high peak and corresponding low peak time can be easily selected in Sim 3F and NASPI 2 cases.In the Sim 1ND case, there is a previous high peak point before the selected high peak time.However, the previous peak time is neglected because Δ||γ t d | t f || becomes a positive value directly.The energy norm in NASPI 8 case continuously decreases until about 32 s.The high peak time is selected after the energy hits the initial low value, and the low peak time is selected before the energy increases larger than the high peak value.Fig. 7(b) shows OCF ft values of the second category cases.The blue line denotes two times the average OCF ft values, and the bus number with the red circle represents the oscillation source.Although the oscillation source is involved in the potential sources, the non-source bus has a higher OCF ft value than the source bus, except for the Sim 1ND case.Therefore, if locating the source is performed based only on the full-time DMD, the adequate source could not be specified.Furthermore, when comparing OCF ft values of NASPI 8 case based on the original DMD to those based on the improved DMD, it can be found that the source bus has a very small value; thus, the source bus would be removed from the source candidates when the original DMD is adopted.Based on the values of OCF tt shown in Fig. 7(c), the proposed method finds the source adequately for all the cases.Especially, notice that in Sim 3F case, OCF tt of bus 45, which is the source, is higher than one of bus 159, while OCF ft values are the opposite.For Sim 1ND and NASPI 8 cases, the source buses, i.e., bus 77 and bus 6333, have the highest OCF tt among the potential source buses.In NASPI 2 case, bus 2604 and 2634 have similar OCF tt values.However, bus 2634 is selected as a single source bus, based on the same reason as NASPI 7 case.
The third category includes the cases where the rectangular forcing signal causes oscillations.Although the rectangular forcing signal has an oscillating frequency similar to the sinusoidal signal, it also contains the harmonics of the basic oscillating frequency.The oscillation source can be easily specified if the basic oscillating frequency is dominant compared to the harmonics like the Sim 6F3 case in Fig. 8.The dominant frequency of the Sim 6F3 case is obtained as 0.4008 Hz, which is 0.2% different from the real value.Then, the oscillation source can be determined to bus 79 based on the two-tier method.However, NASPI 12 case has a quite different norm energy distribution from the Sim 6F3 case.Note that the dominant frequencies of the conjugate pairs include 1.481, 0.7444, and 0.3723 Hz, which are the integer multiple of 0.3723 Hz.Therefore, we determine 0.3723 Hz as the dominant frequency despite the largest norm magnitude observed at 0.7444 Hz.In general, when the rectangular forcing signal is injected, the frequencies with high ||γ i | t f || values are composed of the integer multiples of a specific frequency f s as nf s for natural number n.For these cases, the dominant frequency is determined as f s , even if the other integer multiple frequency has a higher energy norm than f s .Another feature of the NASPI 12 case is that the maximum oscillation magnitude is not at the source; it exists at bus 4131.However, the proposed method finds out the actual source, bus 6335, using the OCF tt values.
The fourth category includes NASPI 3 case, whose voltage data at the source is not monitored.Fig. 9 represents the reduced WECC 240 bus system and the simulation results on NASPI 3 case.In Fig. 9(a), monitored buses are represented as the black box and the non-monitored source bus as the red dotted box.According to the full-time DMD results, six buses are classified as potential sources.Among them, only buses 1431 and 1034 contain generators, but they do not have considerable OCF tt magnitudes.The largest value is observed at bus 1002, which is the non-generator bus.Although the load in bus 1002 can be the source of the oscillation, OCF tt value is too low to be considered as the source.In addition, the monitored buses around bus 1002 have a lower OCF tt value than bus 1002.Therefore, we can conclude that the unmonitored buses, which are close to bus 1002, i.e., bus 1131 and bus 1032, are the oscillation source.However, the proposed method cannot designate the source between the two because of the lack of information.
The oscillation of the NASPI 13 case in the final category is derived from the HVDC controls.Hence, it is hard to accurately pick when the oscillation starts, as shown in Fig. 10(a).Based on the proposed scheme, the start point of the oscillation can be easily designated using the energy level.The oscillation is detected at around 45 s, then t c and t f are selected around 30 s and 70 s each.The dominant frequency is observed as 0.6149 Hz, with its flow depicted in Fig. 10(b).Since the energy decreases below the start level in the beginning, the high peak time is chosen after the lowest level of the energy, and the low peak time can be easily chosen before the energy surpasses the high peak energy.According to the NASPI contest solution, both terminal buses of HVDC 4010 and 2619 are deemed as oscillation sources, while the source is bus 2619.Notice that the proposed method can clearly discover the source bus 2619 in Fig. 10(c).

C. Application on the Real Oscillation Cases
As the representative data of the real oscillation cases, ISO-NE case 3 is used in this article.This case belongs to category 4, where the source is not monitored.Unlike the simulated cases, the real-world measurements have a high possibility that the noise is added to the original signal and outlier data exists.Notice that the original data of the ISO-NE case in Fig. 11(a) contains variable frequency components because of the noise.Therefore, an adequate preprocessing process is necessary for the accurate detection of the oscillation and location of the source with noisy data.Missing data is replaced with the average of the previous and next data.The low pass filter (LPF) is the first step to reducing the noise.However, if we apply LPF to the voltage magnitude data, the noise will re-arise after computing the voltage deviation.Therefore, considering the general range of the low-frequency oscillation, we apply LPF with a 3 Hz cut-off frequency on the voltage deviation data.We apply the smoothing function based on the moving average to the filtered data to further lessen the high-frequency component.The power spectra after the filter and smoothing applications for the real case are shown in Fig. 11(a).Notice that the frequency components above 5 Hz are removed by using LPF, and irrelevant components below 5 Hz are lessened by the smoothing function.Based on the proposed detection algorithm on the noisereduced data, the oscillation is detected at around 69 s.If the measurements did not have noise, data from 58 s to 70 s   would be required to compute the DMD.However, we use the extended data from 58 s to 72 s because of spike generation after the filtering.Based on the full-time DMD results using the noise-reduced data, the dominant frequency of the ISO-NE 3 case can be determined as 1.1507 Hz, which differs 1.83% from the actual frequency of 1.13 Hz.The energy norm values per mode are shown in Fig. 12(a).Contribution factors on the dominant frequency can be found in Fig. 12(b), where the largest factor is observed at both line 2 and line 4 of substation 2 with the same magnitude.Hence, we can infer that the source would be located on the bus connected at both line 2 and line 4, i.e., substation 2 or another unknown bus.However, notice that the OCF tt value of line 3, connected to the west of substation 2, is less than those of lines 2 and 4, connected to the east.Therefore, the source may be located east of substation 2, and the oscillation has been propagated from east to west.

V. CONCLUSION AND FUTURE WORKS
This article proposes an improved DMD-based approach to locate the oscillation source beyond the conventional usage of DMD to identify the oscillation modes.Motivated by the structural limitation of the original DMD in locating the forced oscillations with unusual patterns, we have reformulated the DMD with control (DMDc) by replacing the control input with the initial state of the measurements and successfully demonstrated the performance.Furthermore, the two-tier structure and corresponding normalized contribution factor help to enhance the accuracy and robustness by analyzing the spatial distribution of the energy in the initial stage of the oscillation.Case studies on various oscillation data represent the promising results of the proposed method, even with partial observability.
However, the proposed method has been implemented with the assumption that all data are synchronized and treated in a centralized manner.The distributed implementation of the algorithm and the loss of data synchronization need to be considered in the future.Furthermore, significant efforts need to be exerted to investigate the sub-synchronous oscillations, especially with the fast deployment of inverter-based resources.

Fig. 2 .
Fig. 2. Flow chart of the proposed oscillation source location framework.

Fig. 3 .
Fig. 3. Overall process of the oscillation analysis via the improved DMD.

Fig. 4 .
Fig. 4. Schematic diagram of (a) the timeline for the proposed method and (b) selection of high and low peak times.

r Category 1 :
Normal oscillation cases r Category 2: Maximum oscillation magnitude is observed at non-source bus r Category 3: A rectangular forcing signal is injected into the source r Category 4: The voltage signal at the oscillation source is not measured r Category 5: Other cases The cases within each category and their specifications are shown in Table

Fig. 6 .
Fig. 6.Determination of the dominant frequency and peak times, and the oscillation contribution factors of the proposed method for (a) NASPI 7 case and (b) Sim 6ND case.

Fig. 8 .
Fig. 8. Schematic diagram of the simulation results for the third category cases.The dominant frequency for (a) Sim 6F3 and (b) NASPI 12 case and OCF tt values for (c) Sim 6F3 and (d) NASPI 12 case.

Fig. 10 .
Fig. 10.Schematic diagram of (a) the voltage measurements, (b) selection of the peak times, and (c) OCF tt values for the NASPI 13 case.

Fig. 11 .
Fig. 11.Comparison of (a) power spectra for original, filtered, and filtered and smoothed data and (b) voltage deviation signals of the original and noise-reduced voltage deviation for ISO-NE case 3.