A power system oscillation source positioning method based on extended dynamic modal decoupling
Patent Information
- Application Number
- CN202611451941.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-09-18
- Publication Date
- 2026-10-23
AI Technical Summary
当多个振荡模式的频率较为接近时,不同滤波频带容易发生重叠,造成模式分量泄漏,并可能改变不同监测位置之间原有的幅值和相位关系,使预分离后的振荡分量难以完整保持同步量测数据中的时空关联
1、本发明在同一目标频带内联合保留多个频率相近的目标振荡模式,通过高维提升观测空间及Koopman模态分析确定各目标振荡模式对应的有效共轭模态对,减少逐模态窄带预分离引起的模式泄漏及跨监测位置相位关系失真。
Smart Images

Figure CN122962333A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of power system stability analysis technology, and specifically to a method for locating power system oscillation sources based on extended dynamic mode decoupling. Background Technology
[0002] With the large-scale integration of new energy power generation, energy storage devices, flexible DC transmission, and other power electronic interface equipment, the dynamic characteristics of power systems are becoming increasingly complex. Under the influence of different control links, continuous disturbances, and dynamic coupling between devices, multiple oscillation modes may exist simultaneously in the system, and the frequencies of different oscillation modes may be quite similar. These oscillation modes are superimposed on synchronous measurement data from multiple monitoring locations and propagate along the power network, causing monitoring locations far from the actual oscillation excitation location to exhibit large oscillation amplitudes or strong modal responses.
[0003] Existing oscillation analysis methods based on synchronous measurement data typically perform narrowband filtering or frequency band extraction separately according to different target frequencies. When the frequencies of multiple oscillation modes are close, different filter frequency bands are prone to overlap, causing mode component leakage and potentially altering the original amplitude and phase relationships between different monitoring locations. This makes it difficult for the pre-separated oscillation components to fully maintain the spatiotemporal correlation in the synchronous measurement data.
[0004] Dynamic mode decomposition and its extended methods can obtain information such as oscillation frequency, damping and spatial modes from measurement data, but the modes obtained are usually in complex form, making it difficult to directly use them as physical time series for comparing time-leading relationships between different monitoring locations; if cross-correlation or phase comparison is performed directly on the original multimodal measurement data, it is easily affected by other oscillation modes and periodic multi-peaks.
[0005] Furthermore, existing oscillation source localization methods, which determine the oscillation source based on oscillation amplitude, energy, or modal participation, can only reflect the coupling degree between the corresponding monitoring location and the target oscillation mode, and cannot directly indicate that the target oscillation mode propagates outward from that monitoring location. In cases of network resonance, branch coupling, and far-end response amplification, non-source locations may have larger oscillation responses, leading to misidentification of the oscillation source.
[0006] Meanwhile, during continuous analysis within a sliding time window, the output order of similar oscillation modes may change in different time windows. Relying solely on output sequence number or frequency for mode correlation can easily lead to mismatches between different physical oscillation modes. Summary of the Invention
[0007] The purpose of this invention is to provide a power system oscillation source localization method based on extended dynamic mode decoupling to solve the above-mentioned technical problems, thereby enabling reliable differentiation of similar oscillation modes and oscillation source localization by combining the modal participation characteristics of the target oscillation mode and the propagation characteristics along the power system topology.
[0008] The objective of this invention can be achieved through the following technical solutions: A method for locating power system oscillation sources based on extended dynamic mode decoupling includes: S1. Acquire synchronous measurement data from multiple monitoring locations of the power system under a unified time reference, determine a target frequency band containing at least two target oscillation modes with similar frequencies, suppress out-of-target frequency band components of the synchronous measurement data to obtain target frequency band measurement data, and jointly retain the target oscillation mode within the target frequency band. S2. Construct an observation function based on the target frequency band measurement data at each monitoring location, map the target frequency band measurement data to a high-dimensional enhanced observation space, and without performing narrowband pre-separation for different target oscillation modes, construct a finite-dimensional Koopman operator and a reduced-order Koopman operator corresponding to the extended dynamic mode decomposition based on the high-dimensional enhanced observation space, perform eigenvalue decomposition on the reduced-order Koopman operator, and map the mode information in the reduced-order space to the high-dimensional enhanced observation space to obtain the eigenvalues, right Koopman mode vectors, and left Koopman eigenvectors corresponding to each target oscillation mode, and determine the effective conjugate mode pairs corresponding to each target oscillation mode based on the eigenvalues; S3. For each target oscillation mode, the corresponding effective conjugate mode pair is conjugated and merged on the same type of physical quantity observation coordinates pre-selected at each monitoring location to reconstruct the real-valued single-mode time series of the target oscillation mode at each monitoring location. S4. Based on the right Koopman mode vector and the left Koopman feature vector, generate a mode contribution index for each monitoring location corresponding to the target oscillation mode; based on the power system topology, perform relative lead relationship analysis on the real-valued single-mode time series corresponding to the topologically adjacent monitoring locations, and generate an outward propagation direction index for each monitoring location corresponding to the target oscillation mode based on the relative lead relationship; S5. The monitoring position where the modal contribution index meets the preset contribution condition and the outward propagation direction index meets the preset direction condition is determined as the candidate position of the oscillation source corresponding to the target oscillation mode.
[0009] Furthermore, the synchronous measurement data includes basic measurement vectors corresponding to each monitoring location at multiple sampling times, and the basic measurement vectors include active power and reactive power; The synchronous measurement data is sequentially subjected to time alignment, outlier processing, detrending, and noise suppression, and the target frequency band is determined based on the spectral peak value. For the j-th monitoring location, a local boosting observation vector is constructed based on the observation function: in, This indicates that the j-th monitoring location is at the sampling time. The corresponding basic measurement vector, Let s represent the s-th observation function. The observation function represents the local boosted observation dimension corresponding to the j-th monitoring location. The observation function includes the active power, the reactive power, the square of the active power, the square of the reactive power, and the product of the active power and the reactive power. The local boosted observation vectors corresponding to each of the monitoring locations are concatenated in a fixed order to obtain a global high-dimensional boosted observation vector: Where T represents transpose, J represents the number of monitoring locations, and each monitoring location uses a local observation dictionary composed of the same type of observation functions.
[0010] Furthermore, the construction of the finite-dimensional Koopman operator and the reduced-order Koopman operator corresponding to the extended dynamic mode decomposition includes: Based on the global high-dimensional boosted observation vector corresponding to the consecutive sampling times, M sets of adjacent snapshot pairs are formed, and the autocorrelation empirical matrix G and the time-shifted empirical matrix H are constructed: Where M represents the number of adjacent snapshot pairs, and * represents the conjugate transpose; Construct the finite-dimensional Koopman operator K based on the autocorrelation empirical matrix G and the time-shift empirical matrix H: in, This represents the Moore-Penrose generalized inverse of the autocorrelation empirical matrix G; The autocorrelation empirical matrix G is subjected to truncated singular value decomposition, and the sum of the first r left singular vectors, right singular vectors, and singular values is obtained. And construct the reduced-order Koopman operator: Where r represents the truncation rank, which is determined based on the singular value spectrum transition and cumulative energy of the autocorrelation empirical matrix G.
[0011] Furthermore, the step of performing eigenvalue decomposition on the reduced-order Koopman operator and mapping the modal information in the reduced-order space to the high-dimensional enhanced observation space includes: Eigenvalue decomposition is performed on the reduced-order Koopman operator: Where W represents the right eigenvector matrix of the reduced-order space, Λ represents the discrete-time eigenvalue matrix, and the m-th discrete-time eigenvalue of the discrete-time eigenvalue matrix is... ; Based on the right eigenvector matrix W of the reduced-order space and the matrices corresponding to the first r left singular vectors The right Koopman mode matrix Φ is obtained, and the left Koopman eigenvector matrix Ξ is obtained from the Moore-Penrose generalized inverse of the right Koopman mode matrix Φ. The discrete-time feature values Mapping to the continuous-time domain yields continuous-time eigenvalues. Modal frequencies Damping ratio : Where Δt represents the sampling period, Im(·) represents taking the imaginary part, and Re(·) represents taking the real part; From the candidate modes whose modal frequencies are located within the target frequency band, candidate modes that simultaneously satisfy the preset damping condition, numerical stability condition, and conjugate matching condition are selected, and two candidate modes with approximately conjugate eigenvalues and conjugate matching errors below a preset threshold are determined as an effective conjugate mode pair.
[0012] Furthermore, the step of conjugating and merging the corresponding effective conjugate mode pairs on the pre-selected observation coordinates of the same type of physical quantity at each of the monitoring locations includes: The target mode matrix ΦM is formed by combining the target modes and their conjugate modes corresponding to each target oscillation mode, and the mode amplitude vector is calculated based on the global high-dimensional boosted observation vector corresponding to the initial time. Where b represents the modal amplitude vector, Represents the target mode matrix Moore-Penrose generalized inverse; For each effective conjugate mode pair, the discrete-time eigenvalue with a positive imaginary part is selected as the representative eigenvalue. The representative feature value The corresponding right Koopman mode vector is determined as the representative mode vector. And extract the corresponding modal amplitude from the modal amplitude vector b. ; The pre-selected coordinates of the same type of physical quantity at each monitoring location are defined as the active power observation coordinates. This represents the selection vector for extracting the active power observation coordinates corresponding to the j-th monitoring location from the global high-dimensional enhanced observation vector. The real-valued single-mode time series of the m-th target oscillation mode corresponding to the j-th monitoring location is: in, denoted by t_k, where k represents the sampling time number. T represents the transpose, t_k represents the k-th sampling time, j represents the monitoring location number, m represents the target oscillation mode number, and P represents the active power observation coordinate identifier.
[0013] Furthermore, generating the modal contribution index for each monitoring location corresponding to the target oscillation mode includes: set up For the m-th target oscillation mode, the element of the right Koopman mode vector at the s-th lifting observation coordinate is... For the element of the left Koopman eigenvector corresponding to the m-th target oscillation mode at the s-th lift observation coordinate, calculate the data-driven participation factor corresponding to the m-th target oscillation mode at the s-th lift observation coordinate: set up Given the set of lift observation coordinates for the j-th monitoring location in the high-dimensional lift observation space, the data-driven participation factors belonging to the same monitoring location are weighted and summarized according to the observation variable weights, and then normalized across all monitoring locations to obtain the modal contribution index for the j-th monitoring location corresponding to the m-th target oscillation mode: in, Let J represent the weight of the observation variable corresponding to the s-th improved observation coordinate, J represent the number of monitoring locations, q represent the traversal sequence number of the monitoring locations, and ε represent a positive number to prevent the denominator from being zero.
[0014] Furthermore, the relative leading relationship analysis of the real-valued single-modal time series corresponding to the topologically adjacent monitoring locations includes: For the topologically adjacent i-th and j-th monitoring locations, obtain the real-valued single-mode time series corresponding to the m-th target oscillation mode. and Calculate the signed normalized cross-correlation: in, and Representing the real-valued single-modal time series respectively and Time average, Indicates relative time delay; Based on the modal frequency of the m-th target oscillation mode Determine the period of the target oscillation mode The maximum search delay is determined based on the period of the target oscillation mode. ,exist The maximum positive correlation peak value within the time delay range: in, This indicates the peak delay corresponding to the maximum positive correlation peak. This represents the maximum positive correlation peak value; When the peak delay When the peak delay is determined, the i-th monitoring location is determined to be ahead of the j-th monitoring location; At that time, it is determined that the j-th monitoring position is ahead of the i-th monitoring position.
[0015] Furthermore, generating the outward propagation direction index of each monitoring location corresponding to the target oscillation mode includes: After excluding the peak neighborhood corresponding to the largest positive correlation peak, the second largest positive correlation peak is determined. Based on the largest positive correlation peak, the second largest positive correlation peak, and the minimum distinguishable delay, the valid event of directional evidence is determined. in, For the relevant peak threshold, The peak separation threshold is... The minimum resolvable time delay is determined based on the sampling period and the target oscillation mode period. It is an indicator function; Topological adjacency relationships with a valid directional evidence event value of 1 are determined as valid directional relationships. set up Given the set of topologically adjacent locations corresponding to the j-th monitoring location, and based on the maximum positive correlation peak value, the peak delay, and the topological adjacency weight corresponding to each effective directional relationship, calculate the outward propagation direction index corresponding to the m-th target oscillation mode for the j-th monitoring location: in, q represents the topological adjacency weight and the topological adjacent monitoring location index. When there is no valid directional relationship in the set of topological adjacent locations, the propagation direction state corresponding to the j-th monitoring location is determined to be direction uncertain.
[0016] Furthermore, determining the monitoring location where the modal contribution index satisfies a preset contribution condition and the outward propagation direction index satisfies a preset direction condition as the candidate oscillation source location corresponding to the target oscillation mode includes: Let the contribution threshold be The direction threshold is The set of candidate oscillation source locations corresponding to the m-th target oscillation mode is determined as follows: Wherein, the preset contribution condition corresponding to the modal contribution index and the preset direction condition corresponding to the outward propagation direction index are parallel admission conditions, and the modal contribution index and the outward propagation direction index are not weighted and merged; The monitoring location where the modal contribution index meets the preset contribution condition and the outward propagation direction index does not meet the preset direction condition is determined as the main coupling contribution location. The monitoring location where the modal contribution index satisfies the preset contribution condition and the propagation direction state is uncertain is determined as the main coupling contribution location; The monitoring locations that were not identified as candidate oscillation source locations and major coupling contribution locations are determined as non-source response locations.
[0017] Furthermore, the target oscillation mode is extracted in multiple sliding time windows respectively, and cross-time window target mode matching is performed on the candidate modes in adjacent sliding time windows; Let candidate mode a and candidate mode b be two adjacent candidate modes in the sliding time window. Based on the right Koopman mode vector φa corresponding to candidate mode a and the right Koopman mode vector φa corresponding to candidate mode b... Calculate spatial modal similarity: Based on the modal frequency corresponding to the candidate mode a The modal frequency corresponding to the candidate mode b The matching cost is constructed based on the spatial modal similarity: Where α and β represent non-negative weights; The candidate modes with the lowest matching cost and below the preset matching threshold are determined as the same physical target oscillation modes in adjacent sliding time windows; set up For the set of sliding time windows where the m-th target oscillation mode is effectively identified and the cross-time window target mode matching is completed, the contribution-direction joint persistence index of the j-th monitoring position corresponding to the m-th target oscillation mode is calculated based on the case where the modal contribution index and the outward propagation direction index in each sliding time window simultaneously meet the corresponding admission conditions: in, This represents the modal contribution index corresponding to the w-th sliding time window. This represents the outward propagation direction index corresponding to the w-th sliding time window. Represents a set The number of sliding time windows in the time frame, where w represents the index of the sliding time window; The monitoring location where the contribution-direction joint persistence index meets the preset persistence condition is determined as a candidate location for a stable oscillation source.
[0018] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention jointly retains multiple target oscillation modes with similar frequencies within the same target frequency band. By enhancing the observation space in high dimensions and using Koopman mode analysis, it determines the effective conjugate mode pairs corresponding to each target oscillation mode, thereby reducing mode leakage and phase relationship distortion across monitoring positions caused by narrow-band pre-separation of each mode.
[0019] 2. This invention merges the effective conjugate mode pairs corresponding to each target oscillation mode on the same type of physical quantity observation coordinates uniformly selected at each monitoring location, and reconstructs a real-valued single-mode time series specific to the mode, so that the propagation relationship analysis between different monitoring locations is based on the same target oscillation mode, the same type of physical quantity and the same time reference.
[0020] 3. This invention generates modal contribution indices based on the left and right Koopman modal information, and generates outward propagation direction indices based on the relative leading relationship between topologically adjacent monitoring positions. The modal contribution condition and propagation direction condition are used as parallel admission conditions to determine the candidate positions of oscillation sources, thereby reducing the possibility that high-response, high-coupling non-source positions are misjudged as oscillation sources. Attached Figure Description
[0021] Figure 1 This is a schematic diagram of the overall technical process of the power system oscillation source localization method based on extended dynamic mode decoupling of the present invention; Figure 2 This is a schematic diagram illustrating the joint preservation and extended dynamic mode decoupling of multiple target oscillation modes within the target frequency band of the present invention; Figure 3 This is a schematic diagram illustrating the mapping of reduced-order modal information to a higher-dimensional enhanced observation space and the reconstruction of real-valued single-modal time series according to the present invention. Figure 4 This is a schematic diagram illustrating the determination of oscillation source candidate locations based on modal contribution index and outward propagation direction index according to the present invention. Figure 5 This is a schematic diagram of the output of the target oscillation mode matching and oscillation source localization results across time windows in this invention. Detailed Implementation
[0022] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. These embodiments are based on the technical solution of the present invention and provide detailed implementation methods and specific operating procedures; however, the scope of protection of the present invention is not limited to the following embodiments.
[0023] like Figures 1 to 5 The method for locating power system oscillation sources based on extended dynamic mode decoupling, as shown, includes: S1. Acquire synchronous measurement data from multiple monitoring locations in the power system under a unified time reference. These monitoring locations include synchronous power plant grid connection points, inverter-type resource substation grid connection points, energy storage substation grid connection points, flexible DC converter stations, flexible AC transmission line connection points, large load connection points, substations, and critical line endpoints. The synchronous measurement data uses dynamic measurement data with a unified timestamp, including voltage and current phasors acquired by a synchronous phasor measurement unit, and active and reactive power calculated based on the voltage and current phasors.
[0024] The synchronous measurement data corresponding to each monitoring location are time-aligned, outlier sample processed, detrended, and noise suppressed. A target frequency band containing at least two similar target oscillation modes is determined based on the spectral peak values. Only out-of-band suppression is applied to the synchronous measurement data; no separate narrowband filters are used for different target oscillation modes within the target frequency band, allowing multiple target oscillation modes to be jointly preserved within the same time window and under the same preprocessing conditions.
[0025] S2. Construct observation functions based on target frequency band measurement data from each monitoring location, and map the target frequency band measurement data from each monitoring location to a unified high-dimensional enhanced observation space. In this high-dimensional enhanced observation space, construct a finite-dimensional Koopman operator corresponding to extended dynamic mode decomposition based on data from continuous sampling times. Reduce the order of the finite-dimensional Koopman operator to obtain a reduced-order Koopman operator, and then perform eigenvalue decomposition on the reduced-order Koopman operator. Subsequently, remap the modal information obtained from the reduced-order space to the high-dimensional enhanced observation space to obtain eigenvalues, right Koopman mode vectors, and left Koopman eigenvectors corresponding to different target oscillation modes. Determine the effective conjugate mode pairs corresponding to each target oscillation mode based on modal frequency, damping, numerical stability, and conjugate matching relationships. Multiple similar target oscillation modes in the target frequency band are distinguished by different effective conjugate mode pairs generated during the same Koopman analysis process.
[0026] S3. For each target oscillation mode, pre-select the same type of physical quantity observation coordinates at each monitoring location, such as using active power observation coordinates. Then, conjugate and merge the representative mode and its conjugate mode corresponding to the target oscillation mode to reconstruct the real-valued single-mode time series of the target oscillation mode at each monitoring location. Thus, each real-valued single-mode time series has the same target oscillation mode, physical quantity type, and time reference.
[0027] S4. Calculate the modal contribution index of different monitoring locations to the corresponding target oscillation mode based on the right Koopman mode vector and the left Koopman eigenvector; at the same time, based on the power system topology, perform relative lead relationship analysis on the real-valued single-mode time series corresponding to topologically adjacent monitoring locations, and form an outward propagation direction index based on the effective lead of each monitoring location relative to its topologically adjacent monitoring locations.
[0028] S5. The monitoring locations where the modal contribution index meets the preset contribution conditions and the outward propagation direction index meets the preset direction conditions are determined as candidate oscillation source locations for the corresponding target oscillation mode.
[0029] In this embodiment, by forming a continuous data processing chain, the mode decoupling results of the target oscillation mode are simultaneously used to generate modal contribution evidence and propagation direction evidence. The modal contribution index describes the degree of coupling between the monitoring location and the target oscillation mode, while the outward propagation direction index describes whether the target oscillation mode propagates from the corresponding monitoring location to topologically adjacent monitoring locations. By using both types of evidence to jointly determine oscillation source candidates, the location of the oscillation source can be further distinguished from non-source locations with high response and high coupling.
[0030] Furthermore, to verify the processing method of jointly preserving similar target oscillation modes within the same target frequency band, comparative experiments were conducted under the conditions of nominal frequencies of 2.9 Hz and 3.1 Hz, signal-to-noise ratio of 20 dB, and effective data window of 10 s. The results of the unified target frequency band joint EDMD scheme and the mode-by-mode narrowband pre-filtering scheme are as follows: Table 1. Comparison of Unified Target Frequency Band Joint EDMD and Mode-by-Mode Narrowband Pre-Filtering
[0031] The above comparison results only correspond to the 2.9 Hz and 3.1 Hz modes, 20 dB signal-to-noise ratio, 10 s effective data window and [0.9f, 1.1f] filtering rules used in this embodiment. It does not mean that all narrowband filtering methods cannot be used under all parameter conditions. This comparison is used to illustrate that when the preset filter passbands of similar modes highly overlap, retaining the mixed mode in the same target frequency band and completing the mode differentiation by subsequent Koopman feature decomposition can reduce mode leakage and changes in relative relationships caused by mode-by-mode pre-separation.
[0032] Specifically, suppose the power system contains J monitoring locations, and the j-th monitoring location at the sampling time... The corresponding basic measurement vector is represented as Each basic measurement vector includes at least active power P and reactive power Q.
[0033] The power system stability analysis sequentially performs time alignment, outlier processing, detrending, and noise suppression on the synchronous measurement data of each monitoring location, so that each monitoring location has a corresponding basic measurement vector at the same sampling time.
[0034] Construct a locally boosted observation vector for the j-th monitoring location: in, Let s represent the s-th observation function. The observation function represents the local boosted observation dimension corresponding to the j-th monitoring location. The observation function includes the active power, the reactive power, the square of the active power, the square of the reactive power, and the product of the active power and the reactive power. The local boosted observation vectors corresponding to each of the monitoring locations are concatenated in a fixed order to obtain a global high-dimensional boosted observation vector: Where T represents transpose, J represents the number of monitoring locations, and each monitoring location uses a local observation dictionary composed of the same type of observation function. The use of a structurally consistent local observation dictionary at each monitoring location ensures that the same physical quantity and its nonlinear combinations have corresponding lift-up observation coordinates at different monitoring locations.
[0035] By mapping and stitching synchronous measurement data from multiple monitoring locations according to a fixed structure, the dynamic information of multiple monitoring locations is brought into a unified high-dimensional enhanced observation space, and the correspondence between enhanced observation coordinates and specific monitoring locations and physical quantities is preserved. This provides a data foundation for extracting multiple target oscillation modes, calculating the modal contribution at the monitoring location level, and extracting unified physical quantity observation coordinates in the same Koopman operator.
[0036] Specifically, continuously sampled data are divided into M pairs of adjacent snapshots, and an autocorrelation empirical matrix G and a time-shifted empirical matrix H are constructed based on the global high-dimensional boosted observation vectors corresponding to adjacent sampling times. The superscript H indicates conjugate transpose.
[0037] Construct a finite-dimensional Koopman operator based on the autocorrelation empirical matrix G and the time-shift empirical matrix H: Where G† is the Moore-Penrose generalized inverse of the autocorrelation empirical matrix G.
[0038] Further, a truncated singular value decomposition is performed on the autocorrelation empirical matrix G, and a structure is constructed based on the first r retained singular values and their corresponding singular vectors. Construct a reduced-order Koopman operator: The truncation rank r can be determined based on the singular value spectrum transition, cumulative energy, and the retention of target oscillation modes within the target frequency band. For example, in an embodiment of dual-target oscillation modes, the effective data length is 10 s and the sampling frequency is 60 Hz. After performing truncation singular value decomposition on adjacent snapshot matrices, the first four singular values are significantly higher than the subsequent singular values. Therefore, r=4 is chosen to retain the cosine and sine directions corresponding to the two real-value oscillation modes, respectively.
[0039] The finite-dimensional Koopman operator is used to characterize the data-driven evolution relationship between adjacent sampled states in the high-dimensional enhanced observation space. By truncating the singular value decomposition, a reduced-order Koopman operator is further formed, which can form a lower-dimensional eigenvalue decomposition object while preserving the main dynamic information of the target oscillation mode, and provide a computational basis for subsequently determining the effective conjugate mode pairs corresponding to different target oscillation modes.
[0040] Specifically, eigenvalue decomposition is performed on the reduced-order Koopman operator: Where W represents the right eigenvector matrix of the reduced-order space, Λ represents the discrete-time eigenvalue matrix, and its m-th discrete-time eigenvalue is denoted as . .
[0041] Mapping the right eigenvector in the reduced-order space back to the high-dimensional enhanced observation space yields the right Koopman mode matrix: The left Koopman eigenvector matrix is then determined based on the generalized inverse of the right Koopman mode matrix. The discrete-time feature values Mapping to the continuous-time domain yields continuous-time eigenvalues. Modal frequencies Damping ratio : Where Δt represents the sampling period, Im(·) represents taking the imaginary part, and Re(·) represents taking the real part; This allows two candidate modes with approximately conjugate eigenvalues and conjugate matching errors below a threshold to form an effective conjugate mode pair.
[0042] The mathematical eigenvalue decomposition results of the reduced-order Koopman operator are converted into modal information corresponding to the actual target oscillation modes, and non-target modes are eliminated using frequency, damping, numerical stability, and conjugate relations. When multiple similar oscillation modes are contained in the same target frequency band, it is not necessary to perform pre-separation through multiple independent narrowband filters. Instead, the target oscillation modes are distinguished by different effective conjugate mode pairs in the same reduced-order Koopman operator.
[0043] Specifically, real-valued single-mode time series reconstruction is performed for each effective conjugate mode pair.
[0044] In power system stability analysis, the selected target modes and their conjugate modes are combined to form the target mode matrix ΦM. The mode amplitude vector is then calculated based on the global high-dimensional boosted observation vector at the initial time. Where b represents the modal amplitude vector, Represents the target mode matrix The Moore-Penrose generalized inverse; for each effective conjugate mode pair, the discrete-time eigenvalue with a positive imaginary part is selected as the representative eigenvalue. The representative feature value The corresponding right Koopman mode vector is determined as the representative mode vector. And extract the corresponding modal amplitude from the modal amplitude vector b. .
[0045] For the same target oscillation mode, all monitoring locations pre-select the same type of physical quantity observation coordinates. In a preferred embodiment, each monitoring location selects the lift observation coordinates corresponding to the active power P. Let... Given the selection vector used to extract the active power observation coordinates of the j-th monitoring location from the global high-dimensional observation space, the real-valued single-mode time series representation of the m-th target oscillation mode at the j-th monitoring location is as follows: By merging the representative mode with its conjugate mode, the complex Koopman mode is reconverted into a real-valued time-domain signal. The same type of physical quantity is used at each monitoring location for the same target oscillation mode; different physical quantities are not temporarily mixed between different monitoring locations.
[0046] This makes the input for subsequent topology propagation analysis no longer the raw mixed measurement data containing multiple target oscillation modes, nor the multidimensional complex Koopman modes, but the mode-specific real-valued time series of a specific effective conjugate mode pair. This makes the comparison between different monitoring locations correspond to the same target oscillation mode, the same physical quantity, and the same time reference.
[0047] In the similar dual-mode embodiment, the unified target frequency band is set to [2.5, 3.5] Hz, which simultaneously contains two target oscillation modes with nominal frequencies of 2.8 Hz and 3.2 Hz. The truncation rank is set to r=4, and two pairs of effective conjugate modes are obtained using the same reduced-order Koopman operator. Active power coordinates are uniformly selected at each monitoring location for conjugate merging. The main identification and reconstruction results are as follows: Table 2. Results of joint extraction and real-value reconstruction of similar dual modes
[0048] Specifically, let's set This represents the element of the right Koopman mode vector corresponding to the m-th target oscillation mode at the s-th lift observation coordinate. Let represent the element corresponding to the left Koopman eigenvector at the s-th lift observation coordinate. Then, the data-driven participation factor corresponding to the m-th target oscillation mode at the s-th lift observation coordinate is: Since a monitoring location typically corresponds to multiple boosted observation coordinates, data-driven participation factors belonging to the same monitoring location are further aggregated to the monitoring location level.
[0049] set up Let represent the set of lift observation coordinates for the j-th monitoring location in the high-dimensional lift observation space, and ws represent the weights of the observation variables corresponding to the lift observation coordinates. Then, the monitoring location-level modal contribution index is: Where q represents the traversal sequence number of the monitoring location, and ε is a positive number to prevent the denominator from being zero.
[0050] Modal information, originally defined at the high-dimensional elevated observation coordinate level, is aggregated to the actual power system monitoring location level. This identifies which monitoring locations have a significant modal contribution to the corresponding target oscillation modes, providing the first type of evidence for further identifying candidate oscillation source locations from these high-contribution locations. However, this indicator itself is not directly used as the sole criterion for determining oscillation sources.
[0051] Specifically, relative leading relationship analysis is performed based on the real-valued single-modal time series corresponding to topologically adjacent monitoring locations.
[0052] For the topologically adjacent i-th and j-th monitoring locations, obtain the real-valued single-mode time series corresponding to the same target oscillation mode m. and Calculate the signed normalized cross-correlation between the two: in, and Representing the real-valued single-modal time series respectively and The time mean; Indicates relative time delay. Based on the modal frequency of the target oscillation mode. Determine the oscillation period and set the maximum positive and negative search delay accordingly, then search for the maximum positive correlation peak within the preset search range.
[0053] The relative lead relationship between two topologically adjacent monitoring locations is determined based on the sign of the peak delay. In this process, the cross-correlation function retains its sign rather than its absolute value, ensuring that large negative correlation peaks appearing in the half-cycle out-of-phase case are not used as leading relationships in the same direction.
[0054] The temporal sequence of the same target oscillation mode between two adjacent monitoring locations in the power system topology is converted into computable propagation direction information. This allows the oscillation source localization to consider not only the degree of participation of a certain location in the target oscillation mode, but also the propagation order of the target oscillation mode on the power system topology.
[0055] Specifically, after obtaining the maximum positive correlation peak value ij,m and the peak delay ij,m, the peak neighborhood near the maximum positive correlation peak value is excluded, and then the second largest positive correlation peak value is determined.
[0056] in, For the relevant peak threshold, The peak separation threshold is... The minimum resolvable time delay is determined based on the sampling period and the target oscillation mode period. It is an indicator function; Topological adjacency relationships with a valid directional evidence event value of 1 are determined as valid directional relationships. set up Given the set of topologically adjacent locations corresponding to the j-th monitoring location, and based on the maximum positive correlation peak value, the peak delay, and the topological adjacency weight corresponding to each effective directional relationship, calculate the outward propagation direction index corresponding to the m-th target oscillation mode for the j-th monitoring location: in, q represents the topological adjacency weight and the topological adjacent monitoring location index. When there is no valid directional relationship in the set of topological adjacent locations, the propagation direction state corresponding to the j-th monitoring location is determined to be direction uncertain.
[0057] In one specific implementation, mode 1 propagates along ABCD, with each adjacent edge obtaining a peak delay of +2 sampling points; mode 2 propagates along DCBA, with each adjacent edge obtaining a peak delay of -2 sampling points. The maximum positive correlation peak value of each edge is greater than 0.999, the difference between the maximum peak and the second largest positive peak is greater than 0.14, and the adjacent delay of 33.33 ms exceeds the minimum resolvable delay. Finally, the directional index of mode 1 is [1.000, 0.500, 0.500, 0.000], and the directional index of mode 2 is [0.000, 0.500, 0.500, 1.000].
[0058] This eliminates unreliable propagation direction evidence formed when periodic signals generate multiple similar correlation peaks, the sampling period is insufficient to distinguish fine hour delays, and the correlation degree is insufficient, thus avoiding forcibly giving a direction judgment when the propagation direction evidence is insufficient.
[0059] Specifically, modal contribution conditions and propagation direction conditions are treated as parallel admission conditions that cannot be mutually compensated.
[0060] Let the contribution threshold be The direction threshold is The set of candidate oscillation source locations corresponding to the m-th target oscillation mode is determined as follows: That is, a monitoring location is added to the set of candidate oscillation source locations only when it simultaneously meets the modal contribution condition and the outward propagation direction condition.
[0061] When the modal contribution index of the monitoring location meets the preset contribution condition, but the outward propagation direction index does not meet the preset direction condition, or the propagation direction state corresponding to the monitoring location is uncertain, the monitoring location is determined as the main coupling contribution location; the monitoring location that is not determined as a candidate location for oscillation source and the main coupling contribution location is determined as a non-source response location.
[0062] In one specific implementation, the contribution threshold is set to 0.30, and the direction threshold is set to 0.75. For mode 1, the response amplitude at position C is 1.10, and the modal contribution index is 0.3869, both higher than the 0.90 and 0.3338 corresponding to position A, but the direction index at position C is only 0.500. Therefore, position C is determined as the main coupling contribution position. The direction index at position A is 1.000, satisfying both admission conditions. Therefore, position A is determined as a candidate oscillation source position for mode 1. Mode 2 uses the same judgment process. Position D is determined as a candidate oscillation source position, while B, which has a higher response, is determined as the main coupling contribution position.
[0063] Significant coupling with the target oscillation mode and the outward propagation of the target oscillation mode from that location constitute two independent pieces of evidence. By requiring both pieces of evidence to be true simultaneously, we can avoid misidentifying network resonance, far-end amplification, and highly coupled response locations as oscillation sources when determining the oscillation source solely based on the maximum response amplitude or the maximum modal contribution.
[0064] In the dual admissions implementation, the contribution threshold is set to 0.30, and the direction threshold is set to 0.75. For mode 1, propagation is set from A to B, C, and D; for mode 2, propagation is set from D to C, B, and A, with non-source locations C and B having higher active responses. The response amplitude, modal contribution index, direction index, and classification results for each location are as follows: Table 3. Classification results based on dual admission criteria of modal contribution and propagation direction.
[0065] Further comparison of the dual admission method with the determination methods based on the maximum response amplitude and the maximum modal contribution yields the following results: Table 4. Comparison of results from different oscillation source determination methods
[0066] Specifically, the target oscillation mode extraction, modal contribution index calculation, outward propagation direction index calculation, and oscillation source candidate position determination are repeatedly performed within multiple sliding time windows, and the target oscillation modes in adjacent time windows are matched for identity.
[0067] Let candidate modes a and candidate mode b be in two adjacent sliding time windows, and calculate the spatial mode similarity based on their right Koopman mode vectors: Based on the modal frequency corresponding to the candidate mode a The modal frequency corresponding to the candidate mode b The matching cost is constructed based on the spatial modal similarity: Here, α and β are non-negative weights. Matching is performed between candidate modes in adjacent sliding time windows, and the candidate mode with the lowest matching cost and below a preset matching threshold is identified as the same physical target oscillation mode. Thus, the identity of the target oscillation mode does not depend on the fixed index of the Koopman eigende decomposition output, but is confirmed based on both the frequency similarity and the consistency of the spatial modal structure.
[0068] In one specific implementation, the 10 s of valid data is divided into three sliding time windows of 6 s each. Adjacent time windows shift by 2 s, and the output order of the two candidate modes is swapped in the second time window. After matching using frequency distance and spatial modal similarity, the 2.8 Hz mode and the 3.2 Hz mode still maintain a correct correspondence in the three time windows.
[0069] set up For the set of sliding time windows where the m-th target oscillation mode is effectively identified and cross-time window matching is completed, calculate the contribution-direction joint persistence index for the j-th monitoring location: in, This represents the modal contribution index corresponding to the w-th sliding time window. This represents the outward propagation direction index corresponding to the w-th sliding time window. Represents a set The number of sliding time windows is defined, where w represents the sliding time window index. A sliding time window is only included in the joint persistence statistics if both the modal contribution condition and the propagation direction condition are met simultaneously within the same sliding time window. Based on the comparison between the contribution-direction joint persistence index and the preset persistence conditions, candidate locations of stable oscillation sources are determined.
[0070] On the one hand, cross-time window matching is performed by combining frequency distance and spatial modal similarity to avoid identity mismatch due to output order swapping of target oscillation modes with similar frequencies. On the other hand, by statistically analyzing the persistence ratio of the two types of admission conditions that are simultaneously met within the same time window, the high-modal contribution evidence and outward propagation direction evidence appearing in different time windows are avoided from being spliced into persistent evidence of the same oscillation source, thereby improving the consistency of continuous oscillation source localization results.
[0071] Furthermore, based on the scenario of Example 2, repeated random experiments were conducted to verify the stability and applicable boundaries under different measurement noise and target mode frequency intervals. Each experiment only provided the positioning program with the superimposed P and Q mixed measurements, sampling period, unified target frequency band, and topological adjacency relationship, without providing the actual source location, actual propagation direction, and separated single-mode components.
[0072] Table 5. Source location identification results under different signal-to-noise ratio conditions
[0073] By further altering the frequency interval between the two target oscillation modes at a signal-to-noise ratio of 10 dB, the following results were obtained: Table 6. Source location identification results under different target mode frequency intervals
[0074] The above results are reproducible simulation results of synthetic synchronous measurements. The actual resolvability of similar modes is affected by the data window length, sampling rate, signal-to-noise ratio, spatial mode differences, boosting of the observation dictionary, and truncation rank. When the effective mode, propagation direction, or cross-time window matching does not meet the credibility conditions, the corresponding results are marked as invalid or uncertain, rather than forcing the conclusion of a unique oscillation source by lowering the threshold.
[0075] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments according to this application. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, devices, components, and / or combinations thereof.
[0076] It should be noted that the terms "first," "second," etc., used in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this application described herein can be implemented in sequences other than those illustrated or described herein.
[0077] The present invention has been further described above with reference to specific embodiments. However, it should be understood that the specific description herein should not be construed as limiting the nature and scope of the present invention. Various modifications made to the above embodiments by those skilled in the art after reading this specification are all within the scope of protection of the present invention.
Claims
1. A method for locating oscillation sources in power systems based on extended dynamic mode decoupling, characterized in that, include: S1. Acquire synchronous measurement data from multiple monitoring locations of the power system under a unified time reference, determine a target frequency band containing at least two target oscillation modes with similar frequencies, suppress out-of-target frequency band components of the synchronous measurement data to obtain target frequency band measurement data, and jointly retain the target oscillation mode within the target frequency band. S2. Construct an observation function based on the target frequency band measurement data at each monitoring location, map the target frequency band measurement data to a high-dimensional enhanced observation space, and without performing narrowband pre-separation for different target oscillation modes, construct a finite-dimensional Koopman operator and a reduced-order Koopman operator corresponding to the extended dynamic mode decomposition based on the high-dimensional enhanced observation space, perform eigenvalue decomposition on the reduced-order Koopman operator, and map the mode information in the reduced-order space to the high-dimensional enhanced observation space to obtain the eigenvalues, right Koopman mode vectors, and left Koopman eigenvectors corresponding to each target oscillation mode, and determine the effective conjugate mode pairs corresponding to each target oscillation mode based on the eigenvalues; S3. For each target oscillation mode, the corresponding effective conjugate mode pair is conjugated and merged on the same type of physical quantity observation coordinates pre-selected at each monitoring location to reconstruct the real-valued single-mode time series of the target oscillation mode at each monitoring location. S4. Based on the right Koopman mode vector and the left Koopman feature vector, generate a mode contribution index for each monitoring location corresponding to the target oscillation mode; based on the power system topology, perform relative lead relationship analysis on the real-valued single-mode time series corresponding to the topologically adjacent monitoring locations, and generate an outward propagation direction index for each monitoring location corresponding to the target oscillation mode based on the relative lead relationship; S5. The monitoring position where the modal contribution index meets the preset contribution condition and the outward propagation direction index meets the preset direction condition is determined as the candidate position of the oscillation source corresponding to the target oscillation mode.
2. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 1, characterized in that, The synchronous measurement data includes basic measurement vectors corresponding to each monitoring location at multiple sampling times, and the basic measurement vectors include active power and reactive power. The synchronous measurement data is sequentially subjected to time alignment, outlier processing, detrending, and noise suppression, and the target frequency band is determined based on the spectral peak value. For the j-th monitoring location, a local boosting observation vector is constructed based on the observation function: in, This indicates that the j-th monitoring location is at the sampling time. The corresponding basic measurement vector, Let s represent the s-th observation function. The observation function represents the local boosted observation dimension corresponding to the j-th monitoring location. The observation function includes the active power, the reactive power, the square of the active power, the square of the reactive power, and the product of the active power and the reactive power. The local boosted observation vectors corresponding to each of the monitoring locations are concatenated in a fixed order to obtain a global high-dimensional boosted observation vector: Where T represents transpose, J represents the number of monitoring locations, and each monitoring location uses a local observation dictionary composed of the same type of observation functions.
3. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 2, characterized in that, The construction of the finite-dimensional Koopman operator and the reduced-order Koopman operator corresponding to the extended dynamic mode decomposition includes: Based on the global high-dimensional boosted observation vector corresponding to the consecutive sampling times, M sets of adjacent snapshot pairs are formed, and the autocorrelation empirical matrix G and the time-shifted empirical matrix H are constructed: Where M represents the number of adjacent snapshot pairs, and * represents the conjugate transpose; Construct the finite-dimensional Koopman operator K based on the autocorrelation empirical matrix G and the time-shift empirical matrix H: in, This represents the Moore-Penrose generalized inverse of the autocorrelation empirical matrix G; The autocorrelation empirical matrix G is subjected to truncated singular value decomposition, and the sum of the first r left singular vectors, right singular vectors, and singular values is obtained. And construct the reduced-order Koopman operator: Where r represents the truncation rank, which is determined based on the singular value spectrum transition and cumulative energy of the autocorrelation empirical matrix G.
4. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 3, characterized in that, The step of performing eigenvalue decomposition on the reduced-order Koopman operator and mapping the modal information in the reduced-order space to the high-dimensional enhanced observation space includes: Eigenvalue decomposition is performed on the reduced-order Koopman operator: Where W represents the right eigenvector matrix of the reduced-order space, Λ represents the discrete-time eigenvalue matrix, and the m-th discrete-time eigenvalue of the discrete-time eigenvalue matrix is... ; Based on the right eigenvector matrix W of the reduced-order space and the matrices corresponding to the first r left singular vectors The right Koopman mode matrix Φ is obtained, and the left Koopman eigenvector matrix Ξ is obtained from the Moore-Penrose generalized inverse of the right Koopman mode matrix Φ. The discrete-time feature values Mapping to the continuous-time domain yields continuous-time eigenvalues. Modal frequencies Damping ratio : Where Δt represents the sampling period, This indicates taking the imaginary part. Indicates taking the real part; From the candidate modes whose modal frequencies are located within the target frequency band, candidate modes that simultaneously satisfy the preset damping condition, numerical stability condition, and conjugate matching condition are selected, and two candidate modes with approximately conjugate eigenvalues and conjugate matching errors below a preset threshold are determined as an effective conjugate mode pair.
5. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 4, characterized in that, The step of conjugate merging of the corresponding effective conjugate mode pairs on the pre-selected observation coordinates of the same type of physical quantity at each of the monitoring locations includes: The target mode matrix ΦM is formed by combining the target modes and their conjugate modes corresponding to each target oscillation mode, and the mode amplitude vector is calculated based on the global high-dimensional boosted observation vector corresponding to the initial time. Where b represents the modal amplitude vector, Represents the target mode matrix Moore-Penrose generalized inverse; For each effective conjugate mode pair, the discrete-time eigenvalue with a positive imaginary part is selected as the representative eigenvalue. The representative feature value The corresponding right Koopman mode vector is determined as the representative mode vector. And extract the corresponding modal amplitude from the modal amplitude vector b. ; The pre-selected coordinates of the same type of physical quantity at each monitoring location are defined as the active power observation coordinates. This represents the selection vector for extracting the active power observation coordinates corresponding to the j-th monitoring location from the global high-dimensional enhanced observation vector. The real-valued single-mode time series of the m-th target oscillation mode corresponding to the j-th monitoring location is: in, denoted by , where k represents the sampling time number; T represents the transpose. denoted by k, j represents the monitoring location number, m represents the target oscillation mode number, and P represents the active power observation coordinate identifier.
6. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 5, characterized in that, The generation of modal contribution indices for each monitoring location corresponding to the target oscillation mode includes: set up For the m-th target oscillation mode, the element of the right Koopman mode vector at the s-th lifting observation coordinate is... For the element of the left Koopman eigenvector corresponding to the m-th target oscillation mode at the s-th lift observation coordinate, calculate the data-driven participation factor corresponding to the m-th target oscillation mode at the s-th lift observation coordinate: set up Given the set of lift observation coordinates for the j-th monitoring location in the high-dimensional lift observation space, the data-driven participation factors belonging to the same monitoring location are weighted and summarized according to the observation variable weights, and then normalized across all monitoring locations to obtain the modal contribution index for the j-th monitoring location corresponding to the m-th target oscillation mode: in, Let J represent the weight of the observation variable corresponding to the s-th improved observation coordinate, J represent the number of monitoring locations, q represent the traversal sequence number of the monitoring locations, and ε represent a positive number to prevent the denominator from being zero.
7. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 6, characterized in that, The relative leading relationship analysis of the real-valued single-modal time series corresponding to the topologically adjacent monitoring locations includes: For the topologically adjacent i-th and j-th monitoring locations, obtain the real-valued single-mode time series corresponding to the m-th target oscillation mode. and Calculate the signed normalized cross-correlation: in, and Representing the real-valued single-modal time series respectively and Time average, Indicates relative time delay; Based on the modal frequency of the m-th target oscillation mode Determine the period of the target oscillation mode The maximum search delay is determined based on the period of the target oscillation mode. ,exist The maximum positive correlation peak value within the time delay range: in, This indicates the peak delay corresponding to the maximum positive correlation peak. This represents the maximum positive correlation peak value; When the peak delay When the peak delay is determined, the i-th monitoring location is determined to be ahead of the j-th monitoring location; At that time, it is determined that the j-th monitoring position is ahead of the i-th monitoring position.
8. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 7, characterized in that, The generation of the outward propagation direction index corresponding to each of the monitoring locations for the target oscillation mode includes: After excluding the peak neighborhood corresponding to the largest positive correlation peak, the second largest positive correlation peak is determined. Based on the largest positive correlation peak, the second largest positive correlation peak, and the minimum distinguishable delay, the valid event of directional evidence is determined. in, For the relevant peak threshold, The peak separation threshold is... The minimum resolvable time delay is determined based on the sampling period and the target oscillation mode period. It is an indicator function; The topological adjacency relationship with a valid directional evidence event value of 1 is determined as a valid directional relationship; set up Given the set of topologically adjacent locations corresponding to the j-th monitoring location, and based on the maximum positive correlation peak value, the peak delay, and the topological adjacency weight corresponding to each effective directional relationship, calculate the outward propagation direction index corresponding to the m-th target oscillation mode for the j-th monitoring location: in, q represents the topological adjacency weight and the topological adjacent monitoring location index. When there is no valid directional relationship in the set of topological adjacent locations, the propagation direction state corresponding to the j-th monitoring location is determined to be direction uncertain.
9. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 8, characterized in that, The step of determining the monitoring location where the modal contribution index meets a preset contribution condition and the outward propagation direction index meets a preset direction condition as the candidate location of the oscillation source corresponding to the target oscillation mode includes: Let the contribution threshold be The direction threshold is The set of candidate oscillation source locations corresponding to the m-th target oscillation mode is determined as follows: Wherein, the preset contribution condition corresponding to the modal contribution index and the preset direction condition corresponding to the outward propagation direction index are parallel admission conditions, and the modal contribution index and the outward propagation direction index are not weighted and merged; The monitoring location where the modal contribution index meets the preset contribution condition and the outward propagation direction index does not meet the preset direction condition is determined as the main coupling contribution location. The monitoring location where the modal contribution index satisfies the preset contribution condition and the propagation direction state is uncertain is determined as the main coupling contribution location; The monitoring locations that were not identified as candidate oscillation source locations and major coupling contribution locations are determined as non-source response locations.
10. The power system oscillation source localization method based on extended dynamic mode decoupling according to claim 9, characterized in that, Also includes: The target oscillation mode is extracted in multiple sliding time windows, and cross-time window target mode matching is performed on the candidate modes in adjacent sliding time windows. Let candidate mode a and candidate mode b be two adjacent candidate modes in the sliding time window. Based on the right Koopman mode vector φa corresponding to candidate mode a and the right Koopman mode vector φa corresponding to candidate mode b... Calculate spatial modal similarity: Based on the modal frequency corresponding to the candidate mode a The modal frequency corresponding to the candidate mode b The matching cost is constructed based on the spatial modal similarity: Where α and β represent non-negative weights; The candidate modes with the lowest matching cost and below the preset matching threshold are determined as the same physical target oscillation modes in adjacent sliding time windows; set up For the set of sliding time windows where the m-th target oscillation mode is effectively identified and the cross-time window target mode matching is completed, the contribution-direction joint persistence index of the j-th monitoring position corresponding to the m-th target oscillation mode is calculated based on the case where the modal contribution index and the outward propagation direction index in each sliding time window simultaneously meet the corresponding admission conditions: in, This represents the modal contribution index corresponding to the w-th sliding time window. This represents the outward propagation direction index corresponding to the w-th sliding time window. Represents a set The number of sliding time windows in the time frame, where w represents the index of the sliding time window; The monitoring location where the contribution-direction joint persistence index meets the preset persistence condition is determined as a candidate location for a stable oscillation source.