Three-Stage Data-Driven Phase Analysis to Reveal Generator-Site Origin Source of Forced Oscillations Under Resonance

This paper proposes a novel three-stage near real-time phase-driven procedure to locate the generator-type source of forced oscillations, inspired by the concept of Transient Energy Function (TEF), independent of the amplitude of any signal. In the commencement of the process, each generation bus is assigned with its active power and angular velocity angles reached from the reduced power system graph. Next, the difference in the synchronous generators’ active power angles and the difference in their angular velocity angles is exploited as the graph branches’ weight coefficients. Then, three exclusive and first-proposed decision functions are applied sequentially. First, the weighting coefficients of branches and nodes are fed to the first decision function. Regardless of the grid size, the result of this step is limited to up to four generators. The next step focuses exclusively on the first step outcomes, which results in two generators as recommended items. The last step, relying on the output of the second one, reveals the target generator accurately. The methodological procedure has been validated on the New England 10-machine 39-bus benchmark power system modelled in the Real-Time Digital Simulator (RTDS/RSCAD) and then scrutinised in the MATLAB environment. In each study scenario, the results are compared with the conventional transient energy method. The simulation results revealed that the presented approach reliably releases all sources, including limit cycle and turbine governing reference signal modulation, under the most intense parametric resonance.


I. INTRODUCTION
Forced oscillations (FOs) mechanism in the power systems can be investigated by monitoring the angular manner of the network's strategic points. Under non-resonant conditions, the study of this issue is almost straightforward. The problem becomes significantly more complex when the frequency of the periodic driving force is in the vicinity or equal to the frequency of the system's electromechanical modes [1]. This perilous situation is called resonance, which its presence, as a precursor, can potentially interfere with the overall power system variables [2]. Each point may experience a unique phase shift and amplitude intensification depending on the network The associate editor coordinating the review of this manuscript and approving it for publication was Xiaodong Liang . topology and essence. We are interested in finding openings by revealing patterns of behaviour that are not affected by resonance. Achieving this goal requires a deep understanding of the phenomenon of FOs and the concept of resonance, and how it impacts [3], [4].
FOs emerge as a major threat to the security and stability of the power networks following the occurrence of an external periodic driving force [5], or a mistuned generator controller [6]. Since the process of modelling a real-world power system is always prone to errors, then FOs may even be imposed due to the boldness of unmodeled dynamics [7]. These persistent destructive oscillations can occur at any frequency. However, the most challenging is when they arise in the system's local and inter-area modes' frequency range (i.e., 0.2 Hz to 2 Hz). Although the history of the study of VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ FOs in the power systems dates back to the 1960s, when the power system response to cyclic loads was first studied [8], it has attracted remarkable attention and efforts over the past decade. Due to the dimensions enlargement of modern power systems, we are witnessing widespread blackouts following this phenomenon. References [9], [10] are two comprehensive survey papers on the FOs, which have categorized the matter and location methods from different aspects. FOs location methods are divided into two major categories; the first group uses the system model and measurement data at the same time, such as hybrid methods [11], and mode estimation methods [12]. The second category is exclusively based on measurement data. Energy-based methods [13]- [17] fall into this category.
Although these undesirable sustained oscillations may have other origins, such as cyclic loads [8], the most common and challenging type goes back to the power plant side. Some of the most challenging generation-site motivations which immerse the power system in the FOs are the unsatisfactory nonlinear performance of an apparatus such as synchronous generator AVR [18], over-excitation limiters (OEL) [19], that create stable limit cycles induced by supercritical Poincare-Andropov-Hopf bifurcation [20], modulated steam turbine extraction-control valve reference [21], synchronous generator mechanical power modulation [22], inadequately designed power system stabilizer (PSS) [23], the electro-hydraulic regulating system of a steam turbine [24], and nuclear power plant accelerators [25].
In this paper, we propose a novel phase-driven methodology to disclose the exact source location of the FOs while concentrating on the sources with generator-site origin. The proposed algorithm establishes three sets of decision functions for the first time and performs quite reliable, even under the most severe resonance conditions. This method is inspired by the transient energy methodology but with a completely new concept. Unlike the conventional transient energy method, this approach has no dependence on the amplitude of any signal; Like a function of forced oscillation's frequency, it depends only on the phase angle of the generator's active power and angular velocity. After detecting the system being caught in the trap of the FOs, the power system is reduced to the internal buses of its generators. Phasor measurement units (PMUs) data are used to estimate the phase angles of the angular velocities and active powers of the synchronous generators. Subsequently, the proposed approach demonstrates its performance in three stages. In each step, a set of decision functions are established. The results of the first and second stages affect the continuance of the route. Regardless of the dimension and complexity of the problem, the first step will be able to reduce the number of candidate generators to a maximum of 4. The second step then targets precisely two of them. The real source location is finally revealed at the end of the third step. The algorithm's strength is evaluated by the most challenging scenarios, including the limit cycle and turbine governing reference signal modulations. A source with a frequency equal to one of the natural modes is injected in each scenario to engage the whole benchmark power system in challenging resonance circumstances. The continuation of this paper is organized as follows: Section II discusses how resonance influences the phase angles of different parts of the network. The novel three-stage phase-driven idea of this work is incorporated in Section III. Then, simulations and validation of the proposed methodology are given in Section IV. Finally, Section V is devoted to a summary and conclusions.

II. PHASE-SHIFT DEPENDENCE ON DRIVING FREQUENCY
In this section, the primary purpose is to show the dependence of the angular shift on the frequency, which prompted us to come up with a novel idea (In Section III ) to overcome this dependency. Assume the swing equation of a second-order system to which an external alternating force is applied, as follows [4], [26], [27] The right-hand side of equation (1) is related to the applied alternating external force with frequency ω f . The general solution is derived below, The first term on the right hand side is the particular solution to intermittent driving force. In equation (2), µ = β/2α, and ω d = ω n 2 − µ 2 , ω n = √ γ /α, and f d = ω d /2π are the angular frequency of the damped oscillation, eigen angular frequency, and free oscillations frequency, respectively. The particular response amplitude is related to the external source frequency as follows, The dependence of the angular phase on the frequency of the alternating source, which is the main interest and focus of this work, is thus revealed, Finally, the phase shift is obtained as follows, It should be noted that θ ω f is the phase lag over the excitation force. Equation (5) actually shows how the angular phase shift depends on forcing frequency, especially in full resonance conditions. In the following, while extracting a codified pattern of behaviour, we try to overcome this dependence in order to reliably locate the source of FOs.

III. PROPOSED THREE-STAGE ANGULAR PHASE DATA-DRIVEN METHODOLOGY
This section discusses the proposed methodology in detail. The work is based on the angular position of generators' active powers and angular velocities. During FOs, all network variables oscillate at a frequency equal to the external periodic driving force frequency. Suppose this frequency is far enough from the system's natural modes, i.e., non-resonant conditions. In that case, we will see smaller oscillation amplitudes and more angular differences as the distance from the source location increases. The main challenge is when resonance occurs. Under such circumstances, the situation is not so straightforward and will be very complicated. This complexity stems from the fact that the resonance overshadows network behaviour significantly. In such situations, the network operator needs a reliable way to deal with the problem in a nearly real-time fashion. In this regard, a novel phasedriven approach is proposed in this paper. It can quickly and accurately detect the location of the source under the most severe FOs behaviours. The proposed algorithm shows its performance in three consecutive stages. The following are the concepts and formulations of this method.

A. CONSTRUCTING THE FIRST GROUP OF DECISION FUNCTIONS
The power system is reduced to the generator buses in the first step. The dimensions of the problem are equal to the number of network generators. Therefore, if we consider it as a graph, its vertices are the generators' internal buses. Here we suppose the number of network generators as m. Assuming that the phase of the active power and the angular velocity of the i th generator are θ P i ω f and θ ω i ω f , respectively, we have where θ P i,j ω f and θ ω i,j ω f represent the active powers angle difference and the angular velocities angle difference of the i th and j th generator, respectively. It is necessary to mention that By extending the set of equations (6) to all the network generators in pairs, and then aggregate them into square matrices, we get We call P ω f and ω ω f the generators' active powers and angular velocities phase difference matrices, respectively, whose i, j and j, i arrays belong to the generators i and j. Next, by combining the set of equations (7) matrices, we establish two new functions as follows We define new vectors using the arrays of Π 1 ω f and Π 2 ω f matrices defined in the set of equations (8) as follows, We are trying to reveal a latent pattern by focusing on the set of equations (9) elements. For this purpose, we define the following four decision functions: It is worh mentioning that Ψ 1 and Ψ 2 are 1 × m vectors. The decision arrays ψ 1 Max ω f in the set of equations (10) must ultimately point to one individual generator. Otherwise (if there are duplicate values for each of them), it must be ignored, and in the continuation of the process, it should not be used. Regardless of the grid size, a maximum of four candidate generators are selected. If only one generator results, the work is done, and there is no need to continue. If two generators are selected, it enters the third step directly. In the cases where there are three or four candidates, the second and third steps are taken one after the other. The source generator is among the candidates, and therefore from this moment on, the other generators are ignored. With these considerations in mind, the results of this stage will be re-evaluated in later steps.

B. CONSTRUCTING THE SECOND GROUP OF DECISION FUNCTIONS
In this phase, two new decision functions are inaugurated as follows (11a) and (11b), as shown at the bottom of the page.
The decision functions ϒ Max ω f and ϒ Min ω f are generally defined and each only mentions a unique generator. In practice, the set of equations (11) only applies to the nominee generators acquired in the prior step. The number of candidates has been reduced from a maximum of four generators to strictly two generators. It is noteworthy that if exactly two generators are output in the first step, there is no necessity to apply the second step, and it can jump directly to the third stage.

C. CONSTRUCTING THE THIRD GROUP OF DECISION FUNCTIONS
The third stage begins with the hypothesis that there are precisely two candidate generators. These two generators are delivered directly from the first step or have passed through the second stage of filtration. Then the third group of decision functions is characterised as follows The whole algorithm is just a closed block that receives the network graph information in the first step and reveals the source generator as the third step's output. Despite the fact that the functions in the set of equations (12) are generally defined, since there is no information about the results of the previous steps in practice, then on this basis two, three, or all of them are implemented. The special case occurs when the third step takes its inputs directly from the first one, and only two candidates are the result of the first round, in a way that, at the same time, one or both of them are assigned two values. Following the application of the third step, only two of these decision functions are implemented intelligently. In the next step, relying on the set of equations (12), we define 1 and 2 as In the last step, the precise location of the periodic external driving force source is released as follows.
where η ω f refers to that generator in which the occurrence of an intermittent disturbance has involved the whole power system in FOs. Fig. 1 illustrates the flowchart of the proposed algorithm, where N represents the number of candidate generators.

IV. SIMULATION RESULTS AND DISCUSSIONS
This section tests the proposed method on the New England 10-machine 39-bus benchmark power system modelled in the RTDS Simulator (Real-time Simulator CAD/RSCAD) at Gerald Mainman power systems dynamics and control laboratory, University of Manitoba, to provide the required data. The system's one-line diagram is represented in Fig. 2; Further details can be achieved in [28]. All generator terminal buses have been equipped with PMUs. In practice, however, the states can be estimated by the optimal placement of the units. The simultaneous rating of the PMUs is 60 samples per second. All the PMUs contain Gaussian noises (SNR=40 dB) to follow real-world conventions. Different scenarios of FOs occurrence, including limit cycle and turbine governing reference signal modulation, have been investigated. In each scenario, a specific generator is considered the external driving force with a frequency equal to one of the benchmark power system's predominant natural frequencies.
The preceding data-processing procedures were performed  Toolbox). Following a small signal disturbance in the network, PMU datasets have been utilized to estimate critical dominant oscillation information. An instantaneous pulse with an amplitude of 0.05 per unit has been modulated with the AVR reference signal of generator 5. Spectrum analysis is performed on the angular velocity signals using the Welch algorithm with a 95% confidence interval for 15 seconds. The power spectral density (PSD) results reveal that three electromechanical modes dominate the benchmark system, including an inter-area mode (with a frequency of 0.667 Hz) and two local modes (with frequencies of 1.068 Hz and 1.333 Hz), as can be witnessed in Fig. 3. In the continuous, the most important and probable scenarios that may occur in practice have been investigated.

A. GENERATOR 5 TURBINE GOVERNING REFERENCE SIGNAL MODULATION
The turbine governing reference signal of generator 5 is modulated as τ m 5 (t) = τ m 5 o (t) + 0.05sin(2π(1.068)t). The 1.068 is one of the dominant natural frequencies of the system. So we expect a full resonance to happen. The main component (i.e.,1.068 Hz) and the second to fifth harmonics (i.e., 2.136, 3.2, 4.273, and 5.34 Hz) of the voltage amplitude signal of generator 5 terminal bus are shown in Fig. 4.
Due to this local mode, generator 5 oscillates against generators 6 and 7 [29]. The ideology of such a choice has been to study circumstances where the source of FOs is   located precisely where the oscillation frequency is equal to the local mode of that generator. The generators' angular velocities and active powers envelopes estimated by PMUs are displayed in Fig. 5.
The result of applying the proposed method to this scenario is tabulated in Table 1. As shown in Fig. 6, the energy of the fifth generator shows a inconsistent oscillation pattern. This confirms the results documented in Table 1.    As can be seen from Table 1, ψ 1 Min ω f = 1.068 = 141.29, ψ 1 Max ω f = 1.068 = 1287.68, and ψ 2 Max ω f = 1.068 = 1174.25 belong to generators 5, 9, and 10, respectively. In the case of the ψ 2 Min ω f = 1.068 , since the value of 724.87 is repeated for generators 4, 6, 7, and 9, it is ignored. Therefore, ψ 2 Min ω f = 1.068 does not exist in this scenario. By the end of first stage, three candidate generators 5, 9, and 10 have been proposed. The source of the FOs will be among these cases. In the second stage we get ϒ Max ω f = 1.068 = 5.077 and ϒ Min ω f = 1.068 =   Table 2 also illustrates the results of the method presented in [30]. Based on what is seen in Table 2, the largest value is 144.985, which belongs to the fifth generator, the one with the most energy increase.

B. GENERATOR 6 TURBINE GOVERNING REFERENCE SIGNAL MODULATION
The turbine governing reference signal of generator 6 is modulated as τ m 6 (t) = τ m 6 o (t) + 0.05sin(2π(1.333)t). In this local mode, generator 6 oscillates against generator 7 [29]. The reason for choosing such a scenario was to cover a situation where the source of FOs is located exactly in a place whose oscillation frequency is equal to the local mode of the generator. The generators' angular velocities and active powers envelope estimated by the PMUs are portrayed in Fig. 7.
The main component (i.e.,1333 Hz) and the second and third harmonics (i.e., 2.67 and 4.0 Hz) of the voltage amplitude signal of generator 6 terminal bus are shown in Fig. 8. Table 3 details this scenario.
The sixth generator is identified as the source location in both the proposed and conventional methods. The results of the energy method are shown in Fig. 9.
Following the application of the first step, it results that    are connected to generators 6, 9, and 10, respectively. ψ 2 Min ω f = 1.333 = 577.36. The value ψ 2 Min ω f = 1.333 = 577.36 is repeated in all cases 2, 3, 5, 7, and 9, so it is avoided. So far, three generators 6, 9, and 10 are potential candidates that need further investigation. In other words, the rest of the network generators are removed from the list. Now, with this assumption, we enter the second step. As can be seen, the values ϒ Min ω f = 1.333 = 0.107, and ϒ Max ω f = 1.333 = 0.881 refer to generators 6 and 9, respectively. As a result, generator 10 was eliminated from  the competition. Now, relying on the results of the previous two steps, we reach the third stage. According to the last row of Table 3, it can be seen that, ξ Min  Table 4 also represents the results of the method presented in [30]. Based on what is seen in Table 4, The infinite value in the table belongs to the sixth generator, the one who has the most energy increase.

C. OCCURRENCE OF LIMIT CYCLE IN GENERATOR 3
Generator's 3 exciter model has been changed to BBSEX1, while the rest are the IEEE type DC1 excitation system. We have reduced the generator three exciter's upper bound to a specific value of 2.2847871875 (Obtained by the Bisection Method) and set its AVR gain to 10. A three-phase symmetric fault has occurred in the meantime with a 0.16-second duration. The generators' angular velocities and active powers envelope estimated by PMUs are shown in Fig. 10. The frequency of FOs observed at the origin of this limit cycle is 1.78 Hz. The process of formation of the limit cycle on generator three angular velocity and active power have been illustrated in Fig. 11.
The main component (i.e.,1.8 Hz) and the second to eighth harmonics (i.e., 3.6, 5.4, 7.2, 9.0, 10.8, 12.6, and 14.4 Hz) of the voltage amplitude signal of generator 3 terminal bus can be seen in Fig. 12. The simulation outcome of this scenario is available in Table 5. Let us now analyze this fascinating scenario. By applying the decision functions of the first step, we have that ψ 1 Min ω f = 1.8 = 308.92, ψ 1 Max ω f = 1.8 = 815.00, and ψ 2 Max ω f = 1.8 = 898.32, which target generators 10, 9 and 3, respectively. The VOLUME 10, 2022

minimum value of ψ 2
Min ω f = 1.8 = 516.18 is repeated in all cases 1, 4, 5, 7, 9, and 10, so it is ignored. So far, candidate generators 3, 9, and 10 will advance to the next stage and the rest will be removed from the list. At this point, the second step is applied. According to the results we have that ϒ Min ω f = 1.8 = 0.872, and ϒ Max ω f = 1.8 = 1.546, which correspond to generators 3 and 9, respectively. So generator 10 is also ignored. In the third stage ξ Max 938. Subsequently we get the value η ω f = 1.8 = 2.938, which pinpoint generator 3 as the source of the FOs. As can be seen in Fig. 13, similar results are seen in the potential energy method, in which generator 3 exhibits different behaviors in terms of energy. Table 6 also shows the results of the method presented in [30]. Based on what is seen in Table 6,  The infinite value in the table belongs to generator 3; The generator with the most energy reduction. *** A critical note: The result of this scenario, as reported in the Table 5, may be misleading as long as the proposed methodology is not deeply understood. It should be noted that both the values of 2.938 and 2.044 are related to ξ Max 2 ω f = 1.8 and ξ Max 1 ω f = 1.8 , respectively, both of which are in the 1 ω f = 1.8 function. In fact, it is better to make it clear that the maximum of these two numbers is the result of 1 ω f = 1.8 . This means that there is no 2 ω f = 1.8 that we want to compare it with

D. GENERATOR 9 TURBINE GOVERNING REFERENCE SIGNAL MODULATION
The turbine governing reference signal of generator 9 is modulated as τ m 9 (t) = τ m 9 o (t) + 0.1sin(2π(0.667)t). This is another dominant network natural frequency, whereby a powerful resonance pervades the benchmark. The third harmonic (i.e., 2.0 Hz) of the voltage amplitude signal of generator 9 terminal bus can be caught in Fig. 14. The results are detailed in Table 7.
In this inter-area mode, coherent generators (G4, G5, G6, and G7) oscillate against coherent groups (G2 and G3) and (G8, G9, and G10) [29]. The reason for concentrating on this study scenario has been to examine the particular circumstances in which one of the generator coherent groups   includes the generator locating exactly where the source of the forced oscillations is while not having the most significant contribution to this natural oscillation mode. The generators' angular velocities and active powers envelope estimated by PMUs are displayed in Fig.15.
Next, we will examine the results of applying decision functions. After applying the first step, we have ψ 1 Min ω f = 0.667 = 173.12, ψ 1 Max ω f = 0.667 = 292.06, and ψ 2 Max ω f = 0.667 = 253.27, which represent generators 8, 10 and 9, respectively. The other minimum  value is ψ 2 Min ω f = 0.667 = 155.39 in all generators 2, 3, 4, 5, 6 and 7, which we pass easily. With these descriptions, generators 8, 9 and 10, as possible candidates for the source location, will advance to the second stage. At this point, the decision functions of the second stage come into play. The result illustrate ϒ Min ω f = 0.667 = 0.069, and ϒ Max ω f = 0.667 = 0.484, which belong to generators 9 and 8, respectively. From now on, generator 10 will be removed from the race scene. Instead, generators 8 and 9 go to the next stage. After performing the third step, we conclude ξ Min 1 ω f = 0.667 = 6.206, and ξ Max 2 ω f = 0.667 = 3.300, and correspondingly 1 ω f = 0.667 = 3.300 and 2 ω f = 0.667 = 6.206 which results η ω f = 0.667 = 3.300 related to generator 9 as the actual driving force. As shown in Fig. 16, Generator 9, as the correct source location, behaves clearly apart from the other generators. Table 8 also illustrates the results of the method presented in [30]. Based on what is seen in Table 8, the largest value is 13.072, which belongs to the ninth generator, the one who experiences the most energy increase.

E. GENERATOR 2 TURBINE GOVERNING REFERENCE SIGNAL MODULATION
The turbine governing reference signal of generator 2 is modulated as τ m 2 (t) = τ m 2 o (t) + 0.1sin(2π(1.333)t). This scenario further explains the resonance challenge with a predominant natural frequency of 1.333. The main component (i.e., 1.333 Hz) and the second to sixth harmonics (i.e., 2.67, 4.0, 5.33, 6.67, and 8.0 Hz) of the voltage amplitude signal of generator's 2 terminal bus are shown in Fig. 17. The results obtained through simulations are provided in Table 9.    Generator 2 has minimal participation in this local oscillatory mode, which was the primary purpose of checking this unique condition. During this natural mode, generator 4 oscillates against generators 5, 6 and 7 [29]. The generators' angular velocities and active powers envelope estimated by PMUs are shown in Fig. 18.
This scenario has its own charm, which we will get to in the following. Let's start with the implementation of the first stage. We get ψ 1 Min ω f = 1.333 = 294.04 and   The value of ψ 2 Min ω f = 1.333 = 524.73 is allocated for all generators 1, 3, 5, 8, and 9; So we ignore it. As the achievement of this stage, we focus on generators 2 and 9 as potential candidates. It should be noted that the simultaneous reference of two of the decision functions to a generator in no way indicates that the generator is the source (For instance, assigning ψ 1 Min ω f = 1.333 and ψ 2 Max ω f = 1.333 to generator 2), and the issue is more complicated than that. Here, since there are only two candidate generators, the algorithm jumps to the third stage while ignoring the second step. Now the third step comes into play, we obtain ξ Min As a result, the value of 2.077 is allocated to η ω f = 1.333 , which targets generator 2 as the external disturbance source. The simulation results of the energy-based algorithm in Fig. 19 show that since the generator 2 goes through a different oscillation path, it is the source of the FOs, which confirms the results of the proposed methodology. Table 10 also illustrates the results of the method presented in [30]. Based on what is seen in Table 10, The infinite value in the table belongs to the second generator, the one with the most energy increase.
As can be clearly seen, the proposed algorithm very accurately and quickly identifies the location of the source of forced oscillations. It is worth mentioning that, this algorithm can only pinpoint the forced oscillations sources with the generator side origin. Extending this methodology to load-side oscillating sources (cyclic loads such as arc furnaces, nuclear accelerators) may require further investigation, which is beyond the scope of this paper and is one of the authors' future research interests.

V. CONCLUSION
This paper has discovered a novel phase-driven horizon on the issue of locating sources of forced oscillations (FOs). The proposed perspective can reliably respond to any sources on the generation side. For the first time, three decision functions have been established. Accordingly, the whole work is divided into three main stages. First of all, the power system is reduced to its generator buses. The communication links between each pair of nodes are then weighted based on the phase difference in angular velocities of the rotors and the generators' active powers. The network details is then entered into the first decision function. No matter how large the network, as a result of this analysis, up to four options will remain as candidates. The second decision function then takes effect, which offers precisely two generators as possible choices. In the final step, the third function, based on the previous steps' outcomes, reveals the source's exact location. However, depending on the nature of the source, the second and third stages may not be required. In the entire procedure, no simplification assumptions are made. All desirable real-time measurements have been obtained using PMUs with a reporting rate of 60 samples per second installed at all generators terminal. The proposed methodology has been evaluated on the New England 10-machine 39-bus benchmark power systems modelled in the RTDS/RSCAD Simulator. The results obtained through the simulations performed in the MATLAB environment demonstrate the presented framework's capability to address the FOs source location problem considering real-world constraints. The results also indicate that the proposed approach meticulously releases the sources' location with various kinds of origin, including limit cycle and turbine governing reference signal modulation under the most challenging full resonant circumstances.