A method and system for tracing multi-modal oscillation of offshore new energy access direct current grid

By constructing a weighted and directed system topology graph and performing edge betweenness centrality analysis, the problem of inaccurate oscillation source location in the flexible direct transmission system of new energy power plants was solved, enabling accurate identification and tracing of the oscillation source.

CN120879565BActive Publication Date: 2026-06-26SHANGHAI UNIVERSITY OF ELECTRIC POWER +3

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
SHANGHAI UNIVERSITY OF ELECTRIC POWER
Filing Date
2025-08-27
Publication Date
2026-06-26

AI Technical Summary

Technical Problem

Existing technologies cannot effectively handle the complex multimodal oscillation problem of new energy power plants transmitted through flexible DC transmission systems, especially the low positioning accuracy or erroneous positioning of oscillation sources under the characteristics of high power electronics and multi-timescale control links.

Method used

The new energy transmission system is divided into several energy subsystems. A weighted and directed system topology is constructed. Through edge betweenness centrality analysis and modular partitioning, the dynamic interaction between energy subsystems is quantified, the dominant oscillation propagation path is identified, and the oscillation source is determined.

Benefits of technology

It enables accurate location of the oscillation source, improves the accuracy and reliability of oscillation source tracing, simplifies the analysis of complex systems, and avoids errors.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN120879565B_ABST
    Figure CN120879565B_ABST
Patent Text Reader

Abstract

The present application relates to a kind of offshore new energy access direct current grid multimodal oscillation tracing method and system, belong to wind power generation system technical field, solve the problem of the low positioning accuracy of the existing new energy station through flexible transmission system oscillation source positioning error even positioning problem.The method comprises: the new energy is divided into several energy subsystems by control link and constructs system topology diagram through flexible transmission system;Wherein, the system topology diagram is weighted directed graph, node is each energy subsystem, dynamic energy interaction between each energy subsystem is directed edge, and the edge weight is based on the state variable of each energy subsystem, and the obtained interactive energy weight factor;Based on the interactive energy weight factor between each energy subsystem, the edge betweenness centrality value between each energy subsystem is obtained;Based on the edge betweenness centrality value between each energy subsystem, the system topology diagram is modularized and partitioned, and the energy subsystem of oscillation occurrence is obtained.The accurate positioning of oscillation source is realized.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of wind power generation system technology, and in particular to a method and system for tracing the source of multi-mode oscillations when offshore new energy is connected to a DC power grid. Background Technology

[0002] With the rapid development of renewable energy sources such as wind power and photovoltaics, the proportion of renewable energy connected to the power system has been increasing year by year. Large-scale new energy sources such as offshore wind power and photovoltaics have become the mainstream technical solution for connecting to the grid through flexible direct current transmission (VSC-HVDC) systems. However, new energy power plants transmitted through flexible direct current transmission systems exhibit a high proportion of power electronics characteristics. The interaction of multiple time-scale control links leads to complex multimodal oscillation problems. These oscillations are characterized by strong nonlinearity and wideband coupling.

[0003] Existing impedance analysis methods based on impedance models are only applicable to simple single-input single-output systems and cannot handle complex coupling and modal analysis methods involving multiple converters. Modal decomposition methods rely on linearized models and preset modes, making it difficult to correlate with actual physical subsystems. Traditional complex network methods do not consider the directionality and time-varying nature of energy interaction, resulting in difficulties in quantifying oscillation energy transfer paths and a lack of effective oscillation source partitioning and localization mechanisms. These shortcomings severely restrict the oscillation analysis and suppression capabilities of complex new energy DC power grids. Summary of the Invention

[0004] Based on the above analysis, the embodiments of the present invention aim to provide a method and system for tracing the source of multi-mode oscillations when marine new energy is connected to a DC power grid, in order to solve the problem of low positioning accuracy or even incorrect positioning of oscillation sources in existing new energy power plants via flexible DC transmission systems.

[0005] The objective of this invention is mainly achieved through the following technical solutions:

[0006] This invention provides a method for tracing the source of multi-mode oscillations when connecting offshore renewable energy sources to a DC power grid, comprising the following steps:

[0007] The new energy transmission system is divided into several energy subsystems according to the control links, and a system topology diagram is constructed. The system topology diagram is a weighted and directed graph, where the nodes of the system topology diagram are each energy subsystem, the directed edges are the dynamic energy interaction relationships between each energy subsystem, and the edge weights are the interaction energy weight factors obtained based on the state variables of each energy subsystem.

[0008] Based on the interaction energy weighting factors among the energy subsystems, the edge betweenness centrality values ​​among the energy subsystems are obtained.

[0009] Based on the edge betweenness centrality values ​​between each energy subsystem, the system topology graph is modularly partitioned to obtain the energy subsystems where oscillations occur.

[0010] Furthermore, based on the state variables of each energy subsystem, the interaction energy weighting factors between the energy subsystems are obtained, including:

[0011] Based on the state variables of each energy subsystem, the aperiodic components of the interactive energy derivatives received by each energy subsystem are obtained;

[0012] Integrating the aperiodic components of the interaction energy derivatives received by each energy subsystem yields the aperiodic components of the interaction energy received by each energy subsystem; wherein, the aperiodic components of the interaction energy received by each energy subsystem are the sum of the aperiodic components of the interaction energy transferred from other energy subsystems to each energy subsystem.

[0013] Based on the aperiodic components of the interactive energy received and emitted by each energy subsystem, the interactive energy weighting factor between each energy subsystem is obtained.

[0014] Furthermore, based on the state variables of each energy subsystem, the aperiodic component of the interactive energy derivative received by each energy subsystem is obtained using the following formula:

[0015]

[0016] in, H represents the aperiodic component of the interactive energy derivative received by energy subsystem i; α represents the oscillation decay coefficient; a,i H represents the current-voltage coupling coefficient of energy subsystem i; b,i A represents the voltage-current drive coefficient of energy subsystem i; Fx,i and Let A represent the magnitude and phase of the generalized voltage interaction quantity Fx of energy subsystem i, respectively; Fy,i and Let A represent the magnitude and phase of the generalized current interaction quantity Fy of energy subsystem i, respectively; x,i and Let A represent the magnitude and phase of the generalized voltage state variable Δx of energy subsystem i, respectively; y,j and Let Δy represent the magnitude and phase of the generalized current state variable Δy of the energy subsystem j, respectively.

[0017] Furthermore, the method for obtaining the interaction energy weighting factor between energy subsystems based on the aperiodic components of the interaction energy received and emitted by each energy subsystem includes:

[0018] The difference between the aperiodic component of the interactive energy received by each energy subsystem and the aperiodic component of the interactive energy emitted by each energy subsystem is taken as the effective interactive energy between each energy subsystem.

[0019] Based on the effective interaction energy between each energy subsystem, the interaction energy weighting factor between each energy subsystem is obtained using the following formula:

[0020]

[0021] Where, ω ij The interaction energy weighting factor between energy subsystems i and j is represented by k1, k2, and k3, which represent the weights of the interaction energy, the rate of change of interaction energy, and the energy density of interaction energy between energy subsystems, respectively. This represents the normalized value of the effective interactive energy between energy subsystems i and j. This represents the normalized value of the effective rate of change of interactive energy between energy subsystems i and j. The energy density normalized value represents the effective interaction energy between energy subsystems i and j.

[0022] Furthermore, based on the interaction energy weighting factors among the energy subsystems, the edge betweenness centrality values ​​among the energy subsystems are obtained using the following formula:

[0023]

[0024] in, Let represent the edge betweenness centrality value between energy subsystems i and j; wl(i,j) represents the sum of the interaction energy weight factors of all paths that minimizes the sum of the interaction energy weight factors between energy subsystems i and j; σ ij σ represents the interaction energy weight factor and the minimum number of paths between energy subsystems i and j; ij (e(s,t)) represents the number of paths in the path with the minimum interaction energy weight factor between energy subsystems i and j, which includes a path between energy subsystems s and t; ω st This represents the interaction energy weighting factor between energy subsystems s and t.

[0025] Furthermore, based on the edge betweenness centrality values ​​between each energy subsystem, the system topology graph is iteratively modularly partitioned to obtain the energy subsystems where oscillations occur; wherein, in one iteration:

[0026] Obtain the initial module region for this iteration;

[0027] Calculate the edge betweenness centrality values ​​between nodes in the initial module region;

[0028] Remove the path edge with the highest edge betweenness centrality value in the initial module region;

[0029] When the initial module region is still a module region after removing the path edge, continue to remove the path edge with the highest edge betweenness centrality value among the remaining path edges until the initial module region is divided into two independent regions, and then calculate the modularity of the system topology graph after removing the path edge;

[0030] When the modularity is 0, the region where the energy subsystem from which the energy flows out of the removed path edge is located is used as the initial modular region for the next iteration;

[0031] When the modularity is not 0, all energy subsystems in the region where the energy subsystem from the removed path edge is located are energy subsystems where oscillations occur, and the iteration ends.

[0032] Furthermore, the modularity Q of the system topology is obtained using the following formula:

[0033]

[0034] Where m represents the sum of the absolute values ​​of the interaction energy weights of all path edges after removing the path edges in the system topology graph; c represents a module region after removing the path edges. This represents the set of all module regions after removing path edges; ω represents the set of energy subsystems in module region c; ij The interaction energy weighting factor between energy subsystem i and energy subsystem j; The sum of the out-degree interactive energy weighting factors of energy subsystem i; I represents the sum of the in-degree interactive energy weighting factors of energy subsystem j; c This indicates the validity of module region c.

[0035] Furthermore, when the total number of path edges between the energy subsystems in the module region is at least 1, the validity I of the module region is... c =1;

[0036] Otherwise, the validity of the module region I c =0.

[0037] Furthermore, during the first iteration, the system topology diagram serves as the initial module region.

[0038] On the other hand, the present invention provides a multi-mode oscillation source tracing system for marine new energy access to DC grid, the system comprising: a data acquisition module, an interactive energy calculation module, and an oscillation source determination module;

[0039] The acquisition module M1 is used to acquire the instantaneous changes in the output voltage and current of the new energy through the flexible direct transmission system at the sampling interval of the dominant oscillation mode period, so as to obtain the instantaneous changes in the electrical quantities output by the ports of each energy subsystem in the new energy through the flexible direct transmission system, and use them as the state variables of each energy subsystem.

[0040] The interactive energy calculation module M2 is used to calculate the interactive energy weighting factor between each energy subsystem based on the state variables of each energy subsystem of the new energy transmission system, so as to construct a weighted and directed system topology graph; wherein, the nodes in the system topology graph are each energy subsystem; the edges are the energy interaction relationships between each energy subsystem; the weight of the edge is the interactive energy weighting factor between each energy subsystem; and the direction of the edge is the direction of net energy flow between each energy subsystem.

[0041] The oscillation source determination module M3 is used to perform graph theory deduction based on the interaction energy weight factors between the energy subsystems to determine the energy subsystem in which the oscillation occurs.

[0042] Compared with the prior art, the present invention can achieve at least one of the following beneficial effects:

[0043] 1. This invention constructs a weighted and directed system topology diagram of the new energy transmission system via flexible direct transmission, quantifies the dynamic interaction relationship between energy subsystems, and combines improved edge betweenness centrality analysis to accurately identify the propagation path of the dominant oscillation and achieve accurate positioning of the oscillation source.

[0044] 2. When calculating the interaction energy weighting factor, this invention comprehensively considers multiple factors such as effective interaction energy, effective interaction energy change rate, and effective interaction energy energy density, so that the weighting factor can more comprehensively reflect the actual interaction relationship between energy subsystems, thereby improving the accuracy of oscillation source tracing.

[0045] 3. This invention modularizes the system topology by iteratively partitioning it and calculating the modularity of the system. Through multiple iterations, the range of the oscillation source can be gradually narrowed down, and the specific energy subsystem can be determined. This avoids the difficulties and errors caused by analyzing the entire complex system at once, and improves the reliability of oscillation source localization.

[0046] In this invention, the above-described technical solutions can be combined with each other to achieve more preferred combinations. Other features and advantages of this invention will be set forth in the following description, and some advantages may become apparent from the description or be learned by practicing the invention. The objects and other advantages of this invention can be realized and obtained from what is particularly pointed out in the description and drawings. Attached Figure Description

[0047] The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. Throughout the drawings, the same reference numerals denote the same parts.

[0048] Figure 1 This is a flowchart illustrating a method for tracing the source of multi-mode oscillations when connecting marine new energy sources to a DC power grid, as described in an embodiment of the present invention.

[0049] Figure 2 This is a schematic diagram of the structure of the new energy power station's flexible direct transmission system in an embodiment of the present invention;

[0050] Figure 3 This is a schematic diagram of the topology of the new energy power station's flexible direct transmission system in an embodiment of the present invention;

[0051] Figure 4 This is a schematic diagram of a multi-mode oscillation tracing system for marine new energy access to a DC power grid, as described in an embodiment of the present invention. Detailed Implementation

[0052] Preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings, which form part of this application and are used together with the embodiments of the present invention to illustrate the principles of the present invention, but are not intended to limit the scope of the present invention.

[0053] A specific embodiment of the present invention discloses a method for tracing the source of multi-mode oscillations when connecting marine new energy sources to a DC power grid, such as... Figure 1 As shown, it includes the following steps S1-S3:

[0054] Step S1: Divide the new energy transmission system into several energy subsystems according to the control links and construct a system topology graph; wherein, the system topology graph is a weighted and directed graph, the nodes of the system topology graph are each energy subsystem, the directed edges are the dynamic energy interaction relationships between each energy subsystem, and the edge weights are the interaction energy weight factors obtained based on the state variables of each energy subsystem.

[0055] Specifically, the aforementioned new energy transmission system via flexible direct current transmission is a novel power transmission and distribution method that directly transmits new energy sources (such as wind power and photovoltaics) to the user side or the main power grid through flexible direct current transmission technology (VSC-HVDC). For example... Figure 2As shown, this embodiment takes a four-terminal flexible DC transmission system including a direct-drive wind farm and a photovoltaic power station as an example. The sending-end station includes a direct-drive wind farm PMSG, a photovoltaic power station PV, and sending-end flexible DC converter stations MMC1 and MMC2. Flexible DC converter stations MMC1 and MMC2 receive the electrical energy generated by the direct-drive wind farm and the photovoltaic power station, respectively, and convert the wind power and photovoltaic power into DC power, which is then connected to the DC bus. The receiving-end station includes flexible DC converter stations MMC3 and MMC4, as well as two AC power grids G1 and G2. Receiving-end converters MMC3 and MMC4 invert the DC bus power into AC power and connect it to the AC power grid to realize the external transmission of electrical energy. The DC transmission lines of the sending / receiving flexible DC converter stations are used for power exchange between the sending-end and receiving-end flexible DC converter stations.

[0056] Among them, L f and R f L represents the filter inductance and parasitic resistance of the grid-side converter in a direct-drive wind farm, respectively. s and R s These represent the filter inductance and parasitic resistance of the grid-connected inverter for a photovoltaic power station, respectively; L l1 L l2 L g1 L g2 and R l1 R l2 R g1 R g2 These represent the inductance and resistance of an AC transmission line, respectively; C f C s C t1 C t2 C t3 C t4 These represent the capacitances to ground at the common point of the converter and converter station outlets, respectively; L r1 L r2 L r3 L r4 and R r1 R r2 R r3 R r4 These represent the output line inductance and resistance of converters MMC1, MMC2, MMC3, and MMC4, respectively; r1 i r2 i r3 i r4 These represent the AC output currents of converters MMC1, MMC2, MMC3, and MMC4, respectively; U dc1 U dc2 U dc3 U dc4 These represent the DC-side output voltages of converters MMC1, MMC2, MMC3, and MMC4, respectively; I 12 I13 I 24 I 34 These represent the currents of each DC line; i g1 i g2 U g1 U g2 These represent the current flowing into the power grid and the port voltage, respectively; u f u s u r1 u r2 u pcc3 u pcc4 This indicates the voltage at the common point.

[0057] More specifically, based on the control components of the new energy transmission system via flexible direct transmission, energy models that can characterize the lower-order subsystems of a high-dimensional, multi-order system are constructed, such as... Figure 3 As shown, the structure of the new energy transmission system via flexible direct current transmission is represented in graphical topology. In the energy subsystems of the wind-driven flexible direct current converter station (PMSG-MMC1), A1 represents the outer loop d-axis energy subsystem of the grid-side converter in the direct-drive wind farm; A2 represents the inner loop d-axis energy subsystem of the grid-side converter in the direct-drive wind farm; A3 represents the inner loop q-axis energy subsystem of the grid-side converter in the direct-drive wind farm; A4 represents the d-axis energy subsystem of the AC transmission line on the grid-side converter side of the direct-drive wind farm; and A5 represents the q-axis energy subsystem of the grid-side converter in the direct-drive wind farm. A6 represents the q-axis energy subsystem of the AC transmission line on the MMC1 side of the flexible DC converter station; A7 represents the q-axis energy subsystem of the AC transmission line on the MMC1 side of the flexible DC converter station; A8 represents the inner ring d-axis energy subsystem of the MMC1 flexible DC converter station; A9 represents the inner ring q-axis energy subsystem of the MMC1 flexible DC converter station; A10 represents the outer ring d-axis energy subsystem of the MMC1 flexible DC converter station; A11 represents the outer ring q-axis energy subsystem of the MMC1 flexible DC converter station.

[0058] In the various energy subsystems of the photovoltaic-flexible DC converter station (PV-MMC2), B1 represents the outer ring d-axis energy subsystem of the photovoltaic grid-connected inverter; B2 represents the inner ring d-axis energy subsystem of the photovoltaic grid-connected inverter; B3 represents the inner ring q-axis energy subsystem of the photovoltaic grid-connected inverter; B4 represents the d-axis energy subsystem of the AC transmission line on the grid-connected inverter side of the photovoltaic station; B5 represents the q-axis energy subsystem of the AC transmission line on the grid-connected inverter side of the photovoltaic station; B6 represents the d-axis energy subsystem of the AC transmission line on the MMC2 side of the flexible DC converter station; B7 represents the q-axis energy subsystem of the AC transmission line on the MMC2 side of the flexible DC converter station; B8 represents the inner ring d-axis energy subsystem of the MMC2 flexible DC converter station; B9 represents the inner ring q-axis energy subsystem of the MMC2 flexible DC converter station; B10 represents the outer ring d-axis energy subsystem of the MMC2 flexible DC converter station; and B11 represents the outer ring q-axis energy subsystem of the MMC2 flexible DC converter station.

[0059] In the various energy subsystems of the AC power grid G1, C1 represents the outer loop d-axis energy subsystem of the flexible DC MMC3; C2 represents the outer loop q-axis energy subsystem of the flexible DC MMC3; C3 represents the inner loop d-axis energy subsystem of the flexible DC MMC3; C4 represents the inner loop q-axis energy subsystem of the flexible DC MMC3; C5 represents the phase-locked loop 3 energy subsystem; C6 represents the power grid 1 d-axis energy subsystem; and C7 represents the power grid 1 q-axis energy subsystem.

[0060] In the energy subsystems of AC power grid G2, D1 represents the outer loop d-axis energy subsystem of flexible DC MMC4; D2 represents the outer loop q-axis energy subsystem of flexible DC MMC4; D3 represents the inner loop d-axis energy subsystem of flexible DC MMC4; D4 represents the inner loop q-axis energy subsystem of flexible DC MMC4; D5 represents the phase-locked loop 4 energy subsystem; D6 represents the 2nd d-axis energy subsystem of the power grid; and D7 represents the 2nd q-axis energy subsystem of the power grid.

[0061] E1 represents the flexible DC MMC1 DC energy subsystem; E2 represents the flexible DC MMC2 DC energy subsystem; E3 represents the flexible DC MMC3 DC energy subsystem; and E4 represents the flexible DC MMC4 DC energy subsystem.

[0062] Furthermore, the direction of the directed edges in the system topology graph represents the net energy flow direction in the dynamic energy interaction relationship between each energy subsystem. For example, for the outer ring d-axis energy subsystem A1 of the direct-drive wind farm grid-side converter and the inner ring d-axis energy subsystem A2 of the direct-drive wind farm grid-side converter, there exists an interactive energy E transferred from energy subsystem A1 to energy subsystem A2. A12 And the interactive energy E transferred from energy subsystem A2 to energy subsystem A1 A21 ; when E A12 >E A21 When the energy flow is such that the net energy flow between energy subsystem A1 and energy subsystem A2 is energy subsystem A1 → energy subsystem A2.

[0063] Based on the state variables of each energy subsystem, the weighting factors of the interactive energy between each energy subsystem are obtained and used as the edge weights of the system topology graph, including steps S11-S13:

[0064] Step S11: Based on the state variables of each energy subsystem, obtain the aperiodic component of the interactive energy derivative received by each energy subsystem.

[0065] Specifically, the general mathematical model of each energy subsystem of the new energy power plant's flexible direct transmission system is as follows:

[0066]

[0067] Among them, C k Indicates the equivalent capacitance; L kIndicates equivalent inductance; H a This represents the current-voltage coupling coefficient, i.e., the dynamic feedback strength of current change on voltage. For example, it can take a value of 0.1 to 1.0 S; H b This represents the voltage-current drive coefficient, i.e., the dynamic driving strength of voltage changes on current. For example, it can take a value of 0.5 to 1.5Ω⁻¹; R k Δx represents the dynamic damping coefficient, i.e., the energy dissipation capacity of the energy subsystem; for example, it can be taken as 0.05 to 0.2Ω; Δx represents the generalized voltage state quantity, i.e., voltage change; Δy represents the generalized current state quantity, i.e., current change; F x F represents a generalized voltage interaction quantity, i.e., an external input in terms of voltage; y It represents a generalized current interaction quantity, that is, an external input in terms of current.

[0068] The mathematical models of each energy subsystem are rearranged to obtain the following equation:

[0069]

[0070] The energy conservation equation is obtained by first integrating equation (2). By analyzing the dynamic energy interaction relationship between energy subsystems, we can obtain:

[0071]

[0072] in, It represents voltage potential energy, reflecting the electric field energy stored in the energy subsystem; Represents the kinetic energy of an electric current, reflecting the magnetic field energy of the energy subsystem; R k H a ∫Δy 2 dt represents the energy consumed by the resistor, reflecting the amount of current flowing through the resistor R. k The cumulative time-domain loss; -H a ∫F y Δydt-H b ∫F x Δxdt represents the work done by the external disturbance, which is the energy injected or absorbed by the system by the external disturbance, i.e., the interactive energy; Const represents the total energy reference of the energy subsystem, which is a constant determined by the generalized voltage state and generalized current state of each energy subsystem at the initial moment.

[0073] Define the dynamic energy function E of the energy subsystem as:

[0074] E = E S -E D -E T (4)

[0075]

[0076] Among them, E S Indicates stored energy; E D E represents dissipated energy. T Represents interactive energy.

[0077] More specifically, by calculating the derivative of formula (4) with respect to time, we can obtain:

[0078]

[0079] From formula (6), it can be seen that the derivative of the system's dynamic energy satisfies energy conservation:

[0080]

[0081] in, and These represent the time derivatives of stored energy, dissipated energy, and interactive energy, respectively. The sign and magnitude of the derivatives reflect the trend and rate of change of the corresponding energy terms over time.

[0082] When the new energy power station is disturbed by the flexible direct transmission system and the induced mode is α+jω i During oscillation, where α represents the damping coefficient of the oscillation; ω i ω represents the angular velocity of the oscillation, i.e., the oscillation frequency; j represents the imaginary unit in mathematics. Based on this, all electrical quantities of the energy subsystem are represented in the following form:

[0083]

[0084] Among them, A x and Let A represent the magnitude and phase of the generalized voltage state quantity Δx, respectively; y and Let A represent the magnitude and phase of the generalized voltage state quantity Δy, respectively; Fx and Represent the magnitude and phase of the generalized voltage interaction quantity Fx, respectively; A Fy and These represent the magnitude and phase of the generalized current interaction quantity Fy, respectively.

[0085] Substituting formula (8) into formula (7) yields:

[0086]

[0087]

[0088] More specifically, by extracting the non-periodic components of the derivatives of each energy term, we can obtain:

[0089]

[0090] Therefore, based on the state variables of each energy subsystem, the aperiodic component of the interactive energy derivative received by each energy subsystem can be obtained using the following formula:

[0091]

[0092] in, H represents the aperiodic component of the interactive energy derivative received by energy subsystem i; α represents the oscillation decay coefficient; a,i H represents the current-voltage coupling coefficient of energy subsystem i; b,i A represents the voltage-current drive coefficient of energy subsystem i; Fx,i and Let A represent the magnitude and phase of the generalized voltage interaction quantity Fx of energy subsystem i, respectively; Fy,i and Let A represent the magnitude and phase of the generalized current interaction quantity Fy of energy subsystem i, respectively; x,i and Let A represent the magnitude and phase of the generalized voltage state variable Δx of energy subsystem i, respectively; y,j and Let Δy represent the magnitude and phase of the generalized current state variable Δy of the energy subsystem j, respectively.

[0093] Specifically, the aperiodic component is the decaying component separated from the total interaction energy derivative, characterizing the non-oscillatory energy transfer between subsystems.

[0094] For example, the above general calculation method will be explained using the calculation of energy exchange between sending-end stations as an example:

[0095] For calculating the interactive energy at the sending-end power station, the non-periodic component of the derivative of the interactive energy term of the outer ring d-axis energy subsystem of the grid-side converter in the direct-drive wind farm is considered. for:

[0096]

[0097] in, This represents the rate of change of interactive energy transferred from the q-axis energy subsystem of the AC transmission line in the direct-drive wind farm to the d-axis energy subsystem of the outer ring of the grid-side converter in the direct-drive wind farm; i fdc0 ω1 represents the steady-state value of the DC-side current of the converter; K represents the output angular frequency of the phase-locked loop of the direct-drive wind farm converter; f The gain of the equivalent pulse width modulation for the grid-side converter in a direct-drive wind farm; This indicates the gain of the equivalent power calculation module for the grid-side converter in a direct-drive wind farm. and These represent the amplitude and phase of the q-axis component of the converter output current, respectively. and These represent the amplitude and phase of the DC-side voltage of the converter, respectively. This represents the rate of change of interactive energy transferred from the inner ring d-axis energy subsystem of the grid-side converter in the direct-drive wind farm to the outer ring d-axis energy subsystem of the grid-side converter in the direct-drive wind farm. and These represent the magnitude and phase of the equivalent voltage on the d-axis of the inner loop of the converter, respectively. This represents the rate of change of interactive energy transferred from the d-axis energy subsystem of the grid-side converter-side transmission line in a direct-drive wind farm to the d-axis energy subsystem of the outer ring of the converter. and These represent the magnitude and phase of the d-axis component of the common point voltage at the converter outlet, respectively.

[0098] The non-periodic component of the derivative of the interaction energy term of the outer ring d-axis energy subsystem of the flexible DC converter station MMC1 for:

[0099]

[0100] in, This represents the rate of change of interactive energy transferred from the outer ring q-axis energy subsystem of the MMC1 of the flexible DC converter station to the outer ring d-axis energy subsystem of the MMC1 of the flexible DC converter station; and These represent the magnitude and phase of the q-axis component of the common point voltage at the MMC1 outlet of the converter station, respectively. and These represent the magnitude and phase of the d-axis component of the common point voltage at the MMC1 outlet of the converter station, respectively. This represents the rate of change of interactive energy transferred from the d-axis energy subsystem of the AC transmission line on the MMC1 side of the flexible DC converter station to the d-axis energy subsystem of the outer ring of the MMC1 of the flexible DC converter station; and These represent the amplitude and phase of the d-axis component of the output current of the converter station MMC1, respectively. This represents the rate of energy change within the outer ring d-axis energy subsystem of the MMC1 in the flexible DC converter station; and These represent the amplitude and phase of the d-axis command value of the MMC1 output current of the converter station, respectively. This represents the rate of change of interactive energy transferred from the d-axis energy subsystem of the AC transmission line on the grid-side converter side of the direct-drive wind farm to the d-axis energy subsystem of the MMC1 outer ring of the flexible DC converter station; and These represent the amplitude and phase of the d-axis component of the current in the AC transmission line between the direct-drive wind farm and the flexible DC converter station MMC1, respectively.

[0101] The non-periodic component of the derivative of the interaction energy term of the outer loop q-energy subsystem of the flexible DC converter station MMC1 for:

[0102]

[0103] in, This represents the rate of change of interactive energy transferred from the outer ring d-axis energy subsystem of the MMC1 of the flexible DC converter station to the outer ring q-axis energy subsystem of the MMC1 of the flexible DC converter station; This represents the rate of change of interactive energy transferred from the q-axis energy subsystem of the AC transmission line on the grid-side converter side of the direct-drive wind farm to the outer ring q-axis energy subsystem of the MMC1 flexible DC converter station; and These represent the amplitude and phase of the q-axis component of the output current of the converter station MMC1, respectively. This represents the rate of energy change within the outer ring q-axis energy subsystem of the MMC1 in the flexible DC converter station; and These represent the amplitude and phase of the q-axis command value of the MMC1 output current of the converter station, respectively. This represents the rate of change of interactive energy transferred from the q-axis energy subsystem of the AC transmission line on the MMC1 side of the flexible DC converter station to the q-axis energy subsystem of the outer ring of the MMC1 of the flexible DC converter station; and These represent the amplitude and phase of the q-axis component of the current in the AC transmission line between the direct-drive wind farm and the flexible DC converter station MMC1, respectively.

[0104] Step S12: Integrate the aperiodic component of the interaction energy derivative received by each energy subsystem to obtain the aperiodic component of the interaction energy received by each energy subsystem; wherein, the aperiodic component of the interaction energy received by each energy subsystem is the sum of the aperiodic components of the interaction energy transferred from other energy subsystems to each energy subsystem.

[0105] Specifically, the aperiodic component of the interaction energy derivative is integrated to obtain the cumulative effect of energy changes in each energy subsystem. The aperiodic component is usually related to the damping characteristics and stability level of the system and can reflect the dissipation characteristics of the components in the system. By integrating the aperiodic component of the interaction energy derivative, the amount of energy change within a certain time range can be calculated, thereby obtaining the aperiodic component of the interaction energy received by each energy subsystem.

[0106] For example, the aperiodic component of the interactive energy received by each energy subsystem is represented as follows:

[0107]

[0108] in, E represents the aperiodic component of the interactive energy received by the i-th energy subsystem; T_dc,ji This represents the aperiodic component of the interactive energy received by the i-th energy subsystem from the j-th energy subsystem.

[0109] Step S13: Based on the aperiodic components of the interactive energy received and emitted by each energy subsystem, obtain the interactive energy weighting factor between each energy subsystem.

[0110] Specifically, by obtaining the interaction energy weighting factors between each energy subsystem, the importance and degree of influence of energy interactions between different energy subsystems can be assessed and quantified.

[0111] Furthermore, the method for obtaining the interaction energy weighting factor between energy subsystems based on the aperiodic components of the interaction energy received and emitted by each energy subsystem includes:

[0112] The difference between the aperiodic component of the interactive energy received by each energy subsystem and the aperiodic component of the interactive energy emitted by each energy subsystem is taken as the effective interactive energy between each energy subsystem.

[0113] Specifically, aperiodic components are usually related to the basic operating state and long-term trend of energy subsystems. Therefore, they can reflect the energy interaction of energy subsystems without the influence of periodic fluctuations. By calculating the effective interaction energy between energy subsystems, the actual energy exchange between energy subsystems can be quantified, which can serve as a basis for evaluating the energy flow and dynamic behavior within the system.

[0114] Based on the effective interaction energy between each energy subsystem, the interaction energy weighting factor between each energy subsystem is obtained using the following formula:

[0115]

[0116] Where, ω ij The interaction energy weighting factor between energy subsystems i and j is represented by k1, k2, and k3, which represent the weights of the interaction energy, the rate of change of interaction energy, and the energy density of interaction energy between energy subsystems, respectively. For example, the weights can be set to k1 = 0.3, k2 = 0.5, and k3 = 0.2. This represents the normalized value of the effective interactive energy between energy subsystems i and j. This represents the normalized value of the effective rate of change of interactive energy between energy subsystems i and j. The energy density normalized value represents the effective interaction energy between energy subsystems i and j.

[0117] It should be noted that the normalized value of the effective interaction energy between energy subsystem i and energy subsystem j is... The following formula is used to calculate:

[0118]

[0119] Where, q ij It represents the effective interactive energy between energy subsystems i and j, and is the net energy exchange between energy subsystems i and j.

[0120] The normalized value of the effective interactive energy change rate between energy subsystems i and j is calculated using the following formula:

[0121]

[0122] Among them, v ij It represents the rate of change of effective interactive energy between energy subsystems i and j, and is the derivative of the effective interactive energy, reflecting the acceleration / deceleration trend of energy interaction.

[0123] The normalized energy density value of the effective interaction energy between energy subsystems i and j is calculated using the following formula:

[0124]

[0125] Where, p ij The energy density represents the effective interaction energy between energy subsystems i and j, reflecting the intensity of energy interaction.

[0126] Step S2: Based on the interaction energy weighting factors between each energy subsystem, obtain the edge betweenness centrality values ​​between each energy subsystem.

[0127] Specifically, the edge betweenness centrality value is used to quantify the pivotal nature of the interaction path in the entire system's energy transmission network, representing the importance of the connection between two energy subsystems in the overall energy flow of the system.

[0128] Furthermore, based on the interaction energy weighting factors among the energy subsystems, the edge betweenness centrality values ​​among the energy subsystems are obtained using the following formula:

[0129]

[0130] in, Let represent the edge betweenness centrality value between energy subsystems i and j; wl(i,j) represents the sum of the interaction energy weight factors of all paths that minimizes the sum of the interaction energy weight factors between energy subsystems i and j; σ ij σ represents the interaction energy weight factor and the minimum number of paths between energy subsystems i and j; ij (e(s,t)) represents the number of paths in the path with the minimum interaction energy weight factor between energy subsystems i and j, which includes paths between energy subsystems s and t; ω stThis represents the interaction energy weighting factor between energy subsystems s and t.

[0131] Specifically, two directly connected energy subsystems may also be connected through several intermediate energy subsystems, forming several connection paths. The weighting factor is a value that reflects the intensity of energy interaction. At the same time, due to the rapid response of power electronic devices, oscillating energy preferentially propagates along the path with the minimum time constant. This path usually coincides with the path with the minimum weight. Therefore, the smaller the sum of weights, the smaller the energy transmission delay, and the higher the probability of it becoming the dominant oscillation channel.

[0132] Step S3: Based on the edge betweenness centrality values ​​between each energy subsystem, the system topology graph is modularly partitioned to obtain the energy subsystems where oscillations occur.

[0133] Specifically, in power electronic systems, oscillations do not spread uniformly but propagate preferentially along paths with high betweenness centrality. By using the edge betweenness centrality value, connections that play a key role in the energy flow of the system can be identified. These connections may be the main channels of energy flow or key points of interaction between energy subsystems. Modular partitioning can be performed using the edge betweenness centrality value, that is, dividing the system into relatively independent modules. Since the connections within these modules are relatively tight, while the connections between modules are relatively weak, the complexity of the system can be simplified.

[0134] In power grid systems, oscillations are usually associated with imbalances or instabilities in energy flow. By modular partitioning, regions in the system that may have imbalances or instabilities in energy flow can be identified. These regions are caused by overly complex interactions between energy subsystems or overly concentrated energy flow, which are the energy subsystems where oscillations occur.

[0135] Furthermore, based on the edge betweenness centrality values ​​between each energy subsystem, the system topology graph is iteratively modularly partitioned to obtain the energy subsystems where oscillations occur; wherein, in one iteration, steps S31-S35 are included:

[0136] Step S31: Obtain the initial module region for this iteration.

[0137] Furthermore, in the first iteration, the system topology diagram serves as the initial module region.

[0138] Specifically, the initial module region for this iteration is used to determine the topology that needs to be analyzed in this iteration. In the first iteration, the entire system topology graph is used as the initial module region, which means that the interaction edges of all energy subsystems are used as the analysis objects.

[0139] Step S32: Calculate the edge betweenness centrality values ​​between each node in the initial module region.

[0140] Specifically, using the formula in step S2, the edge betweenness centrality value is calculated for all connected edges within the current initial module region. This is used to identify the most critical energy interaction paths within the current region. It should be noted that the structure of the topology graph of the initial module region changes each time a path edge is removed. The edge betweenness centrality value is calculated based on the shortest path between all node pairs in the topology of the initial module region. Therefore, when the topology changes, the shortest path edges also change. Thus, it is necessary to recalculate the edge betweenness centrality value between each node in the current initial module region to accurately reflect the situation of the initial module region.

[0141] Step S33: Remove the path edge with the highest edge betweenness centrality value in the initial module region.

[0142] Specifically, the path edge with the highest edge betweenness centrality value is removed to cut off the main channel for oscillation energy transmission, thereby disrupting the oscillation propagation path.

[0143] Step S34: When the initial module region is still a module region after removing the path edge, continue to remove the path edge with the highest edge betweenness centrality value among the remaining path edges until the initial module region is divided into two independent regions, and calculate the modularity of the system topology graph after removing the path edge.

[0144] Specifically, since oscillations in power systems usually originate from specific equipment or control links, their energy interactions are characterized by strong local coupling and weak global correlation. Therefore, by dividing the region into two independent regions, the core oscillation region and the non-oscillation region can be forcibly separated, providing a basis for subsequent precise location of the energy subsystem that causes the oscillation.

[0145] Furthermore, the modularity Q of the system topology is obtained using the following formula:

[0146]

[0147] Where m represents the sum of the absolute values ​​of the interaction energy weights of all path edges after removing the path edges in the system topology graph; c represents a module region after removing the path edges. This represents the set of all module regions after removing path edges; ω represents the set of energy subsystems in module region c; ij The interaction energy weighting factor between energy subsystem i and energy subsystem j; The sum of the out-degree interactive energy weighting factors of energy subsystem i; I represents the sum of the in-degree interactive energy weighting factors of energy subsystem j; c This indicates the validity of module region c.

[0148] Specifically, the modularity is used to determine whether the partitioned module region truly represents an independent oscillation energy region; more specifically, the validity of the module region c is determined by I. c , is a binary criterion used to determine whether the currently segmented module region constitutes a physically meaningful oscillating energy module region.

[0149] Furthermore, the effectiveness I of the module region is determined when the total number of path edges between each energy subsystem in the module region is at least 1. c =1;

[0150] Otherwise, the validity of the module region I c =0.

[0151] Specifically, when the segmented module region is only an isolated energy subsystem, the segmented module region c is invalid, and I... c Set to 0; and when at least one actual energy interaction path exists in the segmented module region, then module region c is valid, i.e., I c Set to 1.

[0152] Step S35: When the modularity is 0, the region where the energy subsystem from which the energy flowed out of the removed path edge is located is taken as the initial modular region for the next iteration.

[0153] Specifically, when the modularity is 0, it indicates that the currently segmented module region is still in a globally strongly coupled state, and the oscillation energy has not been localized. Therefore, it needs to be segmented again. At the same time, since the oscillation energy usually flows out from the energy subsystem where the oscillation occurs, it is necessary to continue to segment the region where the energy subsystem where the energy flows out is located.

[0154] When the modularity is not 0, all energy subsystems in the region where the energy subsystem from the removed path edge is located are energy subsystems where oscillations occur, and the iteration ends.

[0155] Specifically, when the modularity is not 0, it indicates that there is significant interactive energy accumulation between energy subsystems in the region where the energy subsystem is located, which can be identified as an energy subsystem where oscillation occurs.

[0156] Another embodiment of the present invention discloses a multi-mode oscillation tracing system for offshore new energy access to DC power grids, such as... Figure 4 As shown, the system includes: a data acquisition module, an interactive energy calculation module, and an oscillation source determination module;

[0157] The acquisition module M1 is used to acquire the instantaneous changes in the output voltage and current of the new energy through the flexible direct transmission system at the sampling interval of the dominant oscillation mode period, so as to obtain the instantaneous changes in the electrical quantities output by the ports of each energy subsystem in the new energy through the flexible direct transmission system, and use them as the state variables of each energy subsystem.

[0158] The interactive energy calculation module M2 is used to calculate the interactive energy weighting factor between each energy subsystem based on the state variables of each energy subsystem of the new energy transmission system, so as to construct a weighted and directed system topology graph; wherein, the nodes in the system topology graph are each energy subsystem; the edges are the energy interaction relationships between each energy subsystem; the weight of the edge is the interactive energy weighting factor between each energy subsystem; and the direction of the edge is the direction of net energy flow between each energy subsystem.

[0159] The oscillation source determination module M3 is used to perform graph theory deduction based on the interaction energy weight factors between the energy subsystems to determine the energy subsystem in which the oscillation occurs.

[0160] The above system and method embodiments are based on the same principles, and their related aspects can be referenced from each other to achieve the same technical effects. For specific implementation processes, please refer to the foregoing embodiments, which will not be repeated here.

[0161] In summary, the multi-mode oscillation tracing method and system for offshore new energy access to DC grids according to embodiments of the present invention have the following beneficial effects:

[0162] 1. This invention constructs a weighted and directed system topology diagram of the new energy transmission system via flexible direct transmission, quantifies the dynamic interaction relationship between energy subsystems, and combines improved edge betweenness centrality analysis to accurately identify the propagation path of the dominant oscillation and achieve accurate positioning of the oscillation source.

[0163] 2. When calculating the interaction energy weighting factor, this invention comprehensively considers multiple factors such as effective interaction energy, effective interaction energy change rate, and effective interaction energy energy density, so that the weighting factor can more comprehensively reflect the actual interaction relationship between energy subsystems, thereby improving the accuracy of oscillation source tracing.

[0164] 3. This invention modularizes the system topology by iteratively partitioning it and calculating the modularity of the system. Through multiple iterations, the range of the oscillation source can be gradually narrowed down, and the specific energy subsystem can be determined. This avoids the difficulties and errors caused by analyzing the entire complex system at once, and improves the reliability of oscillation source localization.

[0165] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for tracing the source of multi-mode oscillations when connecting marine new energy sources to a DC power grid, characterized in that, Includes the following steps: The new energy transmission system is divided into several energy subsystems according to the control links, and a system topology diagram is constructed. The system topology diagram is a weighted and directed graph, where the nodes of the system topology diagram are each energy subsystem, the directed edges are the dynamic energy interaction relationships between each energy subsystem, and the edge weights are the interaction energy weight factors obtained based on the state variables of each energy subsystem. Based on the interaction energy weighting factors among the energy subsystems, the edge betweenness centrality values ​​among the energy subsystems are obtained. Based on the edge betweenness centrality values ​​between each energy subsystem, the system topology graph is modularly partitioned to obtain the energy subsystems in which oscillations occur. Based on the state variables of each energy subsystem, the interaction energy weighting factors between the energy subsystems are obtained, including: Based on the state variables of each energy subsystem, the aperiodic components of the interactive energy derivatives received by each energy subsystem are obtained; Integrating the aperiodic components of the interaction energy derivatives received by each energy subsystem yields the aperiodic components of the interaction energy received by each energy subsystem; wherein, the aperiodic components of the interaction energy received by each energy subsystem are the sum of the aperiodic components of the interaction energy transferred from other energy subsystems to each energy subsystem. Based on the aperiodic components of the interactive energy received and emitted by each energy subsystem, the interaction energy weighting factor between each energy subsystem is obtained; wherein, the difference between the aperiodic components of the interactive energy received and emitted by each energy subsystem is taken as the effective interaction energy between each energy subsystem; based on the effective interaction energy between each energy subsystem, the interaction energy weighting factor between each energy subsystem is obtained using the following formula: in, The interaction energy weighting factor between energy subsystem i and energy subsystem j; These represent the weights of the interaction energy, the rate of change of interaction energy, and the energy density of interaction energy between energy subsystems, respectively. This represents the normalized value of the effective interactive energy between energy subsystems i and j. This represents the normalized value of the effective rate of change of interactive energy between energy subsystems i and j. The energy density normalized value represents the effective interaction energy between energy subsystems i and j.

2. The method according to claim 1, characterized in that, Based on the state variables of each energy subsystem, the aperiodic component of the interactive energy derivative received by each energy subsystem is obtained using the following formula: in, This represents the aperiodic component of the interactive energy derivative received by energy subsystem i; Indicates the oscillation damping coefficient; This represents the current-voltage coupling coefficient of energy subsystem i; Represents the voltage-current drive coefficient of energy subsystem i; and Let i represent the generalized voltage interaction quantities of the energy subsystem i. The amplitude and phase; and Let i represent the generalized current interaction quantities of the energy subsystem i. The amplitude and phase; and Let i represent the generalized voltage state quantities of the energy subsystem i. The amplitude and phase; and Let each represent a generalized current state quantity of the energy subsystem j. The amplitude and phase.

3. The method according to claim 1, characterized in that, Based on the interaction energy weighting factors among the energy subsystems, the edge betweenness centrality values ​​among the energy subsystems are obtained using the following formula: in, This represents the edge betweenness centrality value between energy subsystems i and j; Let represent the sum of the interaction energy weight factors between energy subsystems i and j, and the minimum sum of the interaction energy weight factors of each path. This represents the interaction energy weight factor and the minimum number of paths between energy subsystems i and j; This represents the number of paths that include the path between energy subsystems s and t, and the minimum interaction energy weight factor between energy subsystems i and j. This represents the interaction energy weighting factor between energy subsystems s and t.

4. The method according to claim 3, characterized in that, Based on the edge betweenness centrality values ​​between each energy subsystem, the system topology graph is iteratively modularly partitioned to obtain the energy subsystems where oscillations occur; wherein, in one iteration: Obtain the initial module region for this iteration; Calculate the edge betweenness centrality values ​​between nodes in the initial module region; Remove the path edge with the highest edge betweenness centrality value in the initial module region; When the initial module region is still a module region after removing the path edge, continue to remove the path edge with the highest edge betweenness centrality value among the remaining path edges until the initial module region is divided into two independent regions, and then calculate the modularity of the system topology graph after removing the path edge; When the modularity is 0, the region where the energy subsystem from which the energy flows out of the removed path edge is located is used as the initial modular region for the next iteration; When the modularity is not 0, all energy subsystems in the region where the energy subsystem from the removed path edge is located are energy subsystems where oscillations occur, and the iteration ends.

5. The method according to claim 4, characterized in that, The modularity Q of the system topology graph is obtained using the following formula: Where m represents the sum of the absolute values ​​of the interaction energy weights of all path edges after removing the path edges in the system topology graph; c represents a module region after removing the path edges. This represents the set of all module regions after removing path edges; This represents the set of energy subsystems in module region c; The interaction energy weighting factor between energy subsystem i and energy subsystem j; The sum of the out-degree interactive energy weighting factors of energy subsystem i; The sum of the in-degree interactive energy weight factors of energy subsystem j; This indicates the validity of module region c.

6. The method according to claim 5, characterized in that, The effectiveness of a module region is achieved when the total number of path edges between energy subsystems in the module region is at least 1. ; Otherwise, the validity of the module region .

7. The method according to any one of claims 4-6, characterized in that, In the first iteration, the system topology diagram is used as the initial module region.

8. A multi-mode oscillation tracing system for marine new energy access to a DC power grid, characterized in that, The system includes: a data acquisition module, an interactive energy calculation module, and an oscillation source determination module; The acquisition module M1 is used to acquire the instantaneous changes in the output voltage and current of the new energy through the flexible direct transmission system at the sampling interval of the dominant oscillation mode period, so as to obtain the instantaneous changes in the electrical quantities output by the ports of each energy subsystem in the new energy through the flexible direct transmission system, and use them as the state variables of each energy subsystem. The interactive energy calculation module M2 is used to calculate the interactive energy weighting factor between each energy subsystem based on the state variables of each energy subsystem of the new energy transmission system, so as to construct a weighted and directed system topology graph. In the system topology graph, nodes represent each energy subsystem; edges represent the energy interaction relationships between each energy subsystem; the weight of the edge is the interactive energy weighting factor between each energy subsystem; and the direction of the edge is the direction of net energy flow between each energy subsystem. The calculation of the interactive energy weighting factor between each energy subsystem based on the state variables of each energy subsystem includes: obtaining the aperiodic component of the interactive energy derivative received by each energy subsystem based on the state variables of each energy subsystem; integrating the aperiodic component of the interactive energy derivative received by each energy subsystem to obtain the aperiodic component of the interactive energy received by each energy subsystem; obtaining the effective interactive energy between each energy subsystem based on the difference between the aperiodic component of the interactive energy received and the aperiodic component of the interactive energy emitted by each energy subsystem; and using the following formula to obtain the interactive energy weighting factor between each energy subsystem: ;in, The interaction energy weighting factor between energy subsystem i and energy subsystem j; These represent the weights of the interaction energy, the rate of change of interaction energy, and the energy density of interaction energy between energy subsystems, respectively. This represents the normalized value of the effective interactive energy between energy subsystems i and j. This represents the normalized value of the effective rate of change of interactive energy between energy subsystems i and j. The normalized value of the energy density representing the effective interaction energy between energy subsystems i and j; The oscillation source determination module M3 is used to perform graph theory deduction based on the interaction energy weight factors between the energy subsystems to determine the energy subsystem in which the oscillation occurs.