Fault propagation blocking and transient stability coordination control method for chain microgrid cluster

Abstract

The application discloses a kind of chain microgrid group-oriented fault propagation blocking and transient stability collaborative control method.The method collects the high-frequency electrical quantity data of each node of chain microgrid group, generates synchronous transient data set;Detect transient power oscillation event and extract oscillation characteristic parameter.On this basis, the chain oscillation propagation model reflecting the coupling relationship of adjacent nodes is constructed, especially the energy flow conservation mechanism is quantified The along-chain propagation characteristics of oscillation energy;With the maximum global oscillation suppression index as the target, a virtual impedance optimization model is established combining with the cascade dissipation criterion, and the optimal parameters meeting the energy step-by-step attenuation requirement are calculated.On this basis, the virtual impedance dynamic reconstruction is executed on each node inverter by adaptive control law.The application solves the problem of oscillation energy accumulation and step-by-step amplification in chain topology, and realizes the rapid blocking of fault propagation.

Description

Technical Field

[0001] This invention relates to the field of power system stability control technology, and in particular to a method for coordinated control of fault propagation blocking and transient stability for chain-type microgrid groups. Background Technology

[0002] Chain-type microgrids, as topologies distributed along linear geographical spaces, are widely found in areas far from power sources, such as long-distance transmission corridors, narrow industrial parks, plateau hinterlands, and some inland mountainous regions. In such systems, random disturbances at the source and load ends can easily induce transient power oscillations. If not suppressed in time, the oscillation energy will continue to accumulate and propagate along the interconnects, threatening the synchronization stability of the inverters and the safe operation of the system. Therefore, researching efficient oscillation blocking and suppression technologies is beneficial to ensuring the resilience and reliability of high-proportion power electronic systems.

[0003] Currently, microgrid group oscillation suppression strategies are mostly based on virtual impedance control technology, simulating the inertia and damping characteristics of synchronous generators to improve system stability. Mainstream solutions employ fixed parameter configurations or use the Prony algorithm to identify the dominant modes (frequency, damping ratio) of the system, then utilize consensus algorithms or centralized optimization methods to uniformly adjust the virtual impedance parameters. These methods are well-established in radial or ring-shaped microgrid groups, focusing on solving global low-frequency power distribution and frequency recovery problems, and typically assume that the damping contribution of each node to the system is similar or uniform.

[0004] However, existing technologies, when applied to chain topologies, neglect the physical characteristics of unidirectional energy transfer and cascading along the chain structure, leading to the system's inability to establish an effective energy dissipation barrier. Therefore, further research and innovation are needed to address these issues in existing technologies. Summary of the Invention

[0005] The purpose of this invention is to provide a method for coordinated control of fault propagation blocking and transient stability in chain-type microgrid groups.

[0006] According to one aspect of this application, a method for coordinated control of fault propagation blocking and transient stability in chain-type microgrid groups includes:

[0007] High-frequency electrical quantity data of each node in the chain microgrid group are collected, and multi-node spatiotemporal alignment processing is performed to generate a synchronous transient dataset.

[0008] Based on synchronous transient datasets, transient power oscillation events are detected, and oscillation feature parameters are extracted.

[0009] Based on the oscillation characteristic parameters, a chain oscillation propagation model reflecting the coupling relationship between adjacent nodes is constructed to quantify the propagation characteristics of oscillation along the link;

[0010] Based on the chain oscillation propagation model, a virtual impedance optimization model with chain coordination constraints is established. With the goal of maximizing the global oscillation suppression index, the optimal virtual impedance parameters of each node are calculated.

[0011] The optimal virtual impedance parameters are converted into control commands, and the virtual impedance is dynamically reconstructed for each node inverter to block the propagation of oscillations.

[0012] Beneficial effects: This invention solves the problem of oscillatory energy accumulation and step-by-step amplification in chain topologies, achieving rapid blocking of fault propagation. The related technical effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0013] Figure 1 The flowchart illustrates a collaborative control method for fault propagation blocking and transient stability in a chain-type microgrid group, as provided in this application embodiment.

[0014] Figure 2 This is a flowchart illustrating the detection of transient power oscillation events and the extraction of oscillation characteristic parameters provided in an embodiment of this application.

[0015] Figure 3 The flowchart for constructing a chain oscillation propagation model that reflects the coupling relationship between adjacent nodes is provided in the embodiments of this application.

[0016] Figure 4 A flowchart for oscillation source identification based on Granger causality test provided for embodiments of this application.

[0017] Figure 5 This is a flowchart illustrating the establishment of a virtual impedance optimization model containing chained coordination constraints, provided for embodiments of this application. Detailed Implementation

[0018] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0019] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention 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 embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0020] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0021] Specifically, traditional uniform or evenly distributed virtual impedance configuration strategies cannot match the energy gradient requirements of chain propagation. The residual oscillating energy upstream will be transmitted to downstream like waves and amplified, causing the tail node to often accumulate energy or even become unstable due to insufficient damping.

[0022] Furthermore, static parameter adjustments based on modal characteristics lack the ability to respond in real time to instantaneous energy flow, making it difficult to adaptively adjust according to the deviation of energy dissipation rate during millisecond-level transient processes. This results in long oscillation smoothing periods, making it difficult to meet the requirements for rapid fault interruption.

[0023] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments.

[0024] In the following text, the sensitivity matrix is ​​also called the sensitivity matrix.

[0025] Firstly, an exemplary scheme for a collaborative control method of fault propagation blocking and transient stability for chain-type microgrid groups is provided, particularly the control framework, which supports both signal analysis-based processing paths and optional paths based on energy physics mechanisms, achieving flexible collaborative control. Specifically, it includes:

[0026] Step 101: Collect high-frequency electrical quantity data of each node in the chain microgrid group, perform multi-node spatiotemporal alignment processing, and generate a synchronous transient dataset.

[0027] In this embodiment, a chain-type microgrid cluster specifically refers to a topology formed by connecting multiple distributed generation units, energy storage units, and loads sequentially in series via power lines. This structure is commonly found in island microgrid clusters distributed along coastlines, long-distance railway traction power supply systems, or linearly distributed industrial parks. Unlike traditional radial or ring networks, disturbances in a chain-type topology exhibit a step-by-step propagation characteristic.

[0028] At each node of the chain-like microgrid cluster, typically located at the point of common connection (PCC), a high-frequency data acquisition terminal is deployed. High-frequency electrical quantity data includes, but is not limited to, three-phase instantaneous voltage, three-phase instantaneous current, and calculated instantaneous active and reactive power. The sampling frequency is typically set to 10kHz or higher to meet the analysis requirements of subsynchronous oscillations and high-frequency switching harmonics.

[0029] Considering the geographically dispersed nature of the nodes, each data acquisition terminal is equipped with a GPS or BeiDou high-precision timing module to timestamp each frame of data at the microsecond level. After data aggregation, interpolation algorithms (such as cubic spline interpolation) are used to resample the data from all nodes onto a unified time axis, suppressing minor deviations in sampling time and generating a synchronized transient dataset. This step provides a time-synchronized basis for subsequent analysis of phase differences and energy flow between adjacent nodes.

[0030] Step 102: Based on the synchronous transient dataset, detect transient power oscillation events and extract oscillation feature parameters.

[0031] The system monitors power changes in the synchronous transient data set in real time. When the system is subjected to disturbances such as sudden load changes, short-circuit faults, or drastic fluctuations in the output of new energy sources, transient power oscillations may occur. The detection process is usually based on time-frequency analysis methods to identify oscillation signals that exceed the normal fluctuation range.

[0032] Once an oscillation event is detected, the system immediately extracts oscillation characteristic parameters from the data. These parameters mainly include: the dominant oscillation frequency (or angular frequency), the time-varying envelope of the oscillation amplitude, the oscillation attenuation factor, and the phase distribution of each node relative to the reference node. These characteristic parameters are the direct basis for subsequently establishing the propagation model and calculating control parameters.

[0033] Step 103: Based on the oscillation characteristic parameters, construct a chain oscillation propagation model that reflects the coupling relationship between adjacent nodes and quantify the propagation characteristics of oscillation along the link.

[0034] In practice, the chain oscillation propagation model can employ different modeling strategies. One option is to use a transfer function model based on signal processing to identify the input-output relationship between adjacent nodes to describe the coupling. Another option is to use a chain energy flow conservation model based on physical mechanisms, using energy density and energy flow rate to describe the accumulation and dissipation of oscillating energy.

[0035] Regardless of the modeling method used, the goal is to quantify the propagation characteristics of oscillations along the link. This includes determining whether the oscillation propagates from the beginning to the end of the link or in the opposite direction, and calculating the degree of attenuation or amplification of the oscillation during propagation, i.e., the propagation attenuation factor. This allows the controller to distinguish between the location of the oscillation source and the disturbed node, and to adopt targeted control strategies.

[0036] Step 104: Based on the chain oscillation propagation model, establish a virtual impedance optimization model containing chain coordination constraints, and calculate the optimal virtual impedance parameters of each node with the goal of maximizing the global oscillation suppression index.

[0037] The optimization problem involves solving for the virtual impedance parameters required by each inverter node, typically including virtual resistance and virtual inductance. The virtual impedance optimization model comprises two elements: an objective function and constraints. The objective function maximizes the global oscillation suppression metric. This metric can be the damping ratio of the dominant mode, corresponding to the signal processing path, or the attenuation rate of the end-to-end oscillation energy, corresponding to the energy physical path. Maximizing this metric ensures that oscillations decay and subside relatively quickly.

[0038] The constraints specifically include chain coordination constraints to accommodate the unique characteristics of chain topologies. For example, the virtual impedance of downstream nodes must be greater than that of upstream nodes to form a tiered damping distribution, preventing energy reflection or local accumulation. Furthermore, the optimization model must also consider the inverter's own capacity limitations and small-signal stability constraints. By solving this optimization model, the optimal virtual impedance parameters for each node are obtained.

[0039] Step 105: Convert the optimal virtual impedance parameters into control commands and perform virtual impedance dynamic reconstruction on each node inverter to block oscillation propagation.

[0040] Based on this, the calculated optimal parameters are translated into specific control commands and sent to the local controllers of each node via the communication network. Upon receiving the commands, each node inverter adjusts the virtual impedance value in its internal control loop. This process is called virtual impedance dynamic reconfiguration.

[0041] Reconstruction can be achieved through static parameter switching, where parameters are adjusted to their optimal values ​​all at once; or through dynamic adaptive adjustment, where parameters are continuously adjusted based on real-time feedback deviations. Employing coordinated impedance reconstruction increases the system's damping for specific oscillation modes or constructs a dissipative barrier that hinders the propagation of oscillation energy, thus preventing the escalation of fault effects, smoothing transient fluctuations, and restoring system stability.

[0042] In conjunction with the first aspect, some implementations of the first aspect describe optional implementation methods for high-frequency data sensing and oscillation event detection, providing a data foundation and triggering mechanism for subsequent modeling and control.

[0043] Step 201: Process the synchronous transient dataset using the sliding window short-time Fourier transform and calculate the time-frequency energy spectrum of the transient power.

[0044] Accordingly, the system employs a sliding window short-time Fourier transform (STFT) to process the power signal in the synchronous transient dataset. Specifically, the system sets a time window of fixed length, such as 200 ms, which slides forward along the time axis in fixed steps (e.g., 20 ms).

[0045] For each data segment within a time window, a window function (such as the Hanning window) is applied for weighting to reduce spectral leakage, followed by a Fast Fourier Transform (FFT). This process converts the time-domain power signal into a time-frequency domain representation, yielding the time-frequency energy spectrum of the transient power. This energy spectrum reflects the energy intensity distribution at different frequency components at different times.

[0046] Step 202: Calculate the energy ratio of the time-frequency energy spectrum within the preset oscillation frequency band. When the energy ratio continuously exceeds the trigger threshold, determine that a transient power oscillation event has occurred.

[0047] The system pre-defines the oscillation frequency band of interest, for example, [0.1Hz, 100Hz]. This band covers the range of low-frequency electromechanical oscillations and subsynchronous oscillations commonly found in microgrids. For each instant, the total energy E within this frequency band is calculated. _osc (t) and the total energy E of the entire frequency band _total The ratio of (t) to the energy percentage R _osc (t)=E _osc (t) / E _total (t).

[0048] Simultaneously, a preset trigger threshold is set, such as 0.3% or 30%. The system continuously monitors the energy percentage. To avoid false alarms caused by noise, the state where the energy percentage exceeds the threshold is usually required to last for a certain duration, such as 100ms. When this duration condition is met, the system determines that a transient power oscillation event has occurred, triggering subsequent parameter identification and control procedures. The detection logic based on energy percentage can distinguish between normal load fluctuations (usually energy distributed across the entire frequency band or DC component) and oscillation faults with specific frequency characteristics.

[0049] Step 203: Extract the dominant oscillation mode frequency, oscillation amplitude envelope, and node phase distribution of the transient power oscillation event as oscillation characteristic parameters.

[0050] After confirming the occurrence of oscillations, the system further analyzes the time-frequency energy spectrum. Accordingly, the frequency corresponding to the energy peak is searched within the oscillation band and determined as the dominant oscillation mode frequency. If multiple peaks exist, multimodal oscillations may exist; the mode with the higher energy is selected as the dominant mode.

[0051] Next, the amplitude envelope of the oscillation signal is extracted using Hilbert transform or directly based on the STFT amplitude spectrum, reflecting the time-varying characteristics of the oscillation intensity. Simultaneously, the phase angle of each node at the dominant frequency is calculated, yielding the node phase distribution. By comparing the phases of different nodes, the phase difference between adjacent nodes can be calculated, which helps determine the propagation direction and coupling properties of the oscillation. The extracted dominant oscillation frequency, amplitude envelope, and phase distribution data together constitute the oscillation characteristic parameters, providing input for subsequent model construction.

[0052] In conjunction with the first aspect, among other implementations of the first aspect, an optional scheme based on the physical mechanism of energy flow is provided. This embodiment reveals the accumulation, transmission, and dissipation mechanism of oscillations in chain topology, which has a clearer physical meaning and stronger robustness; correspondingly, it includes:

[0053] Step 301: The chain oscillation propagation model is specifically constructed as a chain energy flow conservation model.

[0054] In this step, the electromechanical oscillation process in the microgrid is mapped as an energy flow process within the network. The chain-like energy flow conservation model no longer focuses on the waveform details of voltage or current, but rather on the oscillating energy, including electromagnetic energy and electromechanical potential energy, its storage at each node, and its transmission between nodes. This model consists of three parts: the energy density equation, the energy flow rate equation, and the state-space equation.

[0055] Step 302: Calculate the instantaneous oscillation energy density of each node based on the node voltage oscillation components, current oscillation components, and dominant oscillation angular frequency in the oscillation characteristic parameters of the synchronous transient dataset.

[0056] Furthermore, the instantaneous oscillation energy density E is defined. _osc_i (t). This energy density comprises three physical components: electric field energy, magnetic field energy, and the equivalent electromechanical potential energy caused by power oscillations. The specific calculation formula is as follows:

[0057] E _osc_i (t)=(1 / 2)×C _node_i ×Δu _i (t) 2 +(1 / 2)×L _node_i ×Δi _osc_i (t) 2 +(1 / (2×(ω _dom ) 2×S _base_i ))×Δp _i (t) 2 ;

[0058] The above formulas are expressed in a per-unit system, E _osc_i (t) represents the dimensionless per-unit energy density, C _node_i The equivalent capacitance parameter of node i is the equivalent synthesis of the inverter output filter capacitor and the distributed capacitance of the tie line to ground at that node. _node_i The equivalent inductance parameter of node i is derived from the equivalent synthesis of the inverter output filter inductance and the series inductance of the tie line at that node. C can be determined through conventional equivalent circuit analysis based on the actual circuit topology parameters of the microgrid. _node_ i L _node_i ;

[0059] Δu _i (t) represents the voltage oscillation component at node i; Δi _osc_i (t) represents the current oscillation component at node i; ω _dom The angular frequency of the dominant oscillation mode; Δp _i (t) represents the active power oscillation component at node i; S _base_i Let be the baseline capacity of node i.

[0060] The third term in the above formula (1 / (2×(ω)) _dom ) 2 ×S _base_i ))×Δp _i (t) 2 By mapping the amplitude of power oscillations to the form of potential energy, the dynamic response of a generator or inverter can be incorporated into the framework of energy conservation.

[0061] Step 303: Based on the tie line impedance parameters and oscillation phase difference between adjacent nodes, establish an energy flow rate equation describing the transfer of oscillation energy between adjacent nodes.

[0062] In a chain structure, energy is primarily transferred between adjacent nodes via interconnects. Let P be the energy flow rate from node i to node i+1. _flow_i (t). This flow rate is driven by two parts: one part is the diffusion flow driven by the energy density gradient, and the other part is the coupling flow driven by the phase difference. The specific formula is as follows:

[0063] P _flow_i (t)=K _transfer_i ×(E _osc_i (t)-E _osc_i+1 (t))+G _couple_i ×sqrt(E _osc_i (t)×E _osc_i+1 (t))×sin(Δφ_i (t));

[0064] Among them, P _flow_i (t) represents the per-unit energy flux; K _transfer_i E is the per-unit energy transfer coefficient determined by the tie line impedance; _osc_i (t) and E _osc_i+1 (t) represents the energy density of adjacent nodes; G _couple_i Δφ is the per-unit oscillatory coupling conductivity between adjacent nodes. _i (t) represents the oscillation phase difference between adjacent nodes; sqrt represents the square root operation.

[0065] In other embodiments, K _transfer_i The specific value is determined by the tie line impedance parameters and the node rated voltage; G _couple_i The specific values ​​are determined by the tie-line admittance parameters and the node voltage amplitude. These coefficients can be derived from the power balance equation based on the specific line parameters of the microgrid.

[0066] In this equation, when there is an energy density difference between adjacent nodes, energy will diffuse from the high-density region to the low-density region; at the same time, the phase difference between nodes will also forcefully pull the energy exchange like a synchronous torque.

[0067] Step 304: Based on the chain topology, integrate the instantaneous oscillation energy density and energy flow rate equations of each node into a state space equation, and construct a tridiagonal state matrix that reflects the chain-like adjacent coupling characteristics.

[0068] Based on the law of conservation of energy, the rate of energy change at any node equals the energy flowing into that node minus the energy flowing out, minus the node's own dissipation. Solving the energy balance equations of all nodes simultaneously forms the state-space equations of the entire system, specifically:

[0069] dE _osc / dt=A _chain ×E _osc +B _diss ×G _damp ;

[0070] Among them, E _osc G is the state vector composed of the energy densities of each node; _damp A is the equivalent damping term generated by the virtual impedance of each node. _chain Corresponding state matrix, B _diss The corresponding damping term gain matrix.

[0071] In some scenarios, this equation is a small-signal linearization model, suitable for the small-signal fluctuation range corresponding to the operating point, G _damp Specifically, it refers to the linearized equivalent damping coefficient vector.

[0072] In particular, since each node in a chain topology is only directly connected to its immediate and adjacent nodes, the state matrix A _chain It exhibits a tridiagonal structure. That is, only the elements on the main diagonal and the two secondary diagonals above and below it are non-zero, while the remaining elements are all zero. The sparse tridiagonal structure is not only a mathematical fingerprint of the chain topology, but also reduces the complexity of subsequent optimization calculations, making real-time control of large microgrid groups possible.

[0073] In conjunction with the first aspect, and in some other implementations of the first aspect, a specific technical solution for parameter decoupling and optimization configuration based on the cascade dissipation criterion is described. This embodiment proposes a control concept for cascade dissipation, introduces a dimensionality reduction method for optimal impedance angle decoupling, solves the problem of energy accumulation at the end of the chain caused by impedance configuration homogenization in traditional methods, and reduces the computational burden of online optimization.

[0074] Step 401: The global oscillation suppression index is specifically set as the full-link oscillation energy attenuation rate, which is determined by the energy dissipation power generated by the virtual impedance of each node.

[0075] The optimization objective is not the modal damping ratio, but rather the energy decay rate, which reflects the physical nature of the system. The total energy of the entire link oscillation is defined as the sum of the instantaneous oscillation energy densities of all nodes. The energy decay rate of the entire link oscillation reflects how quickly the system dissipates oscillation energy.

[0076] According to the law of conservation of energy, the rate of energy decay equals the negative of the total power dissipation of the system. In microgrid inverters, virtual impedance (mainly virtual resistance components) is the primary dissipation element. Therefore, the optimization objective function J... _opt This can be expressed as maximizing the total power dissipation:

[0077] J _opt =max∑P _diss_i =max∑(3×R _v_i ×(I _osc_i ) 2 );

[0078] Among them, P _diss_i R represents the power dissipation generated by the virtual impedance of node i; _v_i I is the virtual resistance of node i; _osc_i Let be the effective value of the oscillation current at node i; ∑ represents the summation over all nodes i in the chained microgrid group, and max corresponds to the maximum value operator. By maximizing this index, the oscillation energy can be forced to be converted into heat energy at a faster speed, which is reflected in the numerical decay in the control algorithm, thus quickly suppressing the oscillation.

[0079] Step 402: Establish a cascade dissipation criterion as a chain coordination constraint. The cascade dissipation criterion limits the ratio of the virtual impedance of the downstream node to the virtual impedance of the upstream node along the oscillation propagation direction.

[0080] In other words, the cascade dissipation criterion is confirmed, and the chain coordination is constrained accordingly; the direction in which the oscillation propagates is determined, and the virtual impedance of the downstream and upstream nodes in that direction is defined, limiting the proportional relationship between the two.

[0081] Correspondingly, traditional virtual impedance configurations typically use the same or similar parameters at each node, which can lead to insufficient damping in a chain topology. That is, energy transmitted from upstream cannot be absorbed downstream, causing energy reflection at the end and forming standing waves. Therefore, this embodiment proposes a cascade dissipation criterion.

[0082] According to this criterion, downstream nodes must have a stronger energy dissipation capacity than upstream nodes to absorb the remaining energy transmitted from upstream along the direction of oscillation energy propagation. This criterion requires that the virtual impedance (or virtual resistance) of downstream nodes must be greater than that of upstream nodes, forming a spatial impedance gradient.

[0083] It should be understood that gradient configuration is similar to cascade energy dissipation in flood control systems, where energy is reduced layer by layer during the process of energy transfer.

[0084] Step 403: Based on the propagation attenuation factor obtained from the propagation characteristics of the oscillation along the link, determine the minimum impedance increment gradient required to maintain the cascade dissipation criterion, and constrain the search space of the optimal virtual impedance parameter accordingly (minimum impedance increment gradient).

[0085] In this step, the propagation attenuation factor β is introduced. _prop The value is 0 < β _prop <1 describes the natural attenuation ratio of oscillating energy after it has been transmitted through a line.

[0086] If the line impedance is small, β _prop When the value is close to 1, the energy is transferred downstream with almost no loss. In this case, the downstream needs to be configured with a larger impedance to dissipate the energy.

[0087] Based on the energy flow continuity equation, the minimum impedance increasing gradient required to maintain the cascade dissipation criterion satisfies the following analytical inequality:

[0088] R _v_i+1 / R _v_i ≥sqrt(1 / β _prop_i );

[0089] Among them, R _v_i+1 R is the virtual resistance of the downstream node i+1 along the propagation direction; _v_i β is the virtual resistance of upstream node i;_prop_i is the energy propagation attenuation factor between node i and i+1.

[0090] This formula shows that the smaller the propagation attenuation, the greater the corresponding β. _prop_i The closer to 1, the more impedance multiples are required (sqrt(1 / β)). _prop_i The larger it is.

[0091] In other embodiments, the propagation attenuation factor can be specifically defined as the ratio of the instantaneous oscillation energy density of node i+1 to that of node i under oscillatory steady-state propagation conditions;

[0092] This ratio can be obtained by calculating the average ratio of the energy densities of adjacent nodes based on the initial stage data after oscillations are detected in the synchronous transient dataset.

[0093] The initial stage is the natural oscillation stage before the intervention of virtual impedance control, reflecting the propagation attenuation characteristics of the tie line itself.

[0094] Based on this, the inequality constitutes the optimization variable R. _v The lower bound of the search space ensures that the solved parameters must satisfy the cascade dissipation requirements.

[0095] Step 404: Decouple the virtual impedance optimization model and fix the virtual impedance angle parameter of each node to the analytical optimal angle that maximizes the equivalent damping conductance. The analytical optimal angle is 45 degrees at the dominant oscillation frequency.

[0096] Accordingly, this embodiment uses an analytical method to determine the virtual resistance R. _v With virtual inductance L _v The optimal proportional relationship is as follows. The equivalent damping conductance provided by the virtual impedance does not increase monotonically with the virtual resistance, but is limited by the line inductive reactance. Furthermore, by finding the extremum of the equivalent damping conductance function, it can be proven that when the virtual resistance equals the virtual reactance, i.e., the virtual resistance R... _v =Virtual impedance X _v =ω _dom ×L _v At that time, the damping effect is relatively strong.

[0097] In the complex plane, the virtual impedance angle θ _v It should be set to 45 degrees (or π / 4 radians). Therefore, in this step, the virtual impedance angle parameter of all nodes is directly fixed to 45 degrees, i.e., constraint L is set. _v_i =R _v_i / ω _dom , where i is the node index. This operation will reduce the time required to optimize R. _v and L _v The bivariate problem is decoupled into requiring only the optimization of R. _vFor single-variable problems (or impedance magnitudes), the computational workload is reduced by more than half.

[0098] Step 405: Based on the fixed virtual impedance angle parameters, the virtual impedance optimization model is transformed into a univariate distribution optimization problem for the virtual impedance magnitude.

[0099] Based on this, the optimization model was rewritten to include only the virtual impedance magnitude |Z _v_i The objective function becomes finding the sequence of moduli {|Z}. _v_1 |,|Z _v_2 |,...,|Z _v_N The |} structure results in a relatively large overall energy decay rate for the system. In other words, N represents the total number of nodes in the chained microgrid group, and i is the node index.

[0100] The optimization problem then becomes a univariate distributed optimization problem. Since the constraints (cascade dissipation criterion) are linear and the objective function (dissipated power) is a convex function with respect to the modulus after decoupling, the problem can be solved using a highly efficient algorithm, and even has an analytical solution under certain conditions (such as ignoring inverter capacity limitations).

[0101] Step 406: Solve the univariate distribution optimization problem to obtain the optimal virtual impedance magnitude of each node that satisfies the cascade dissipation criterion, and generate the optimal virtual impedance parameters by combining the analytical optimal angle.

[0102] Alternatively, the above problem can be solved using convex optimization algorithms (such as the interior-point method) or fast gradient descent to obtain the optimal virtual impedance magnitude |Z| for each node. _v_i | _opt .

[0103] Based on this, and using a fixed 45-degree angle, the modulus is converted back into specific resistance and inductance parameters, using the following formula:

[0104] R _v_i_opt =|Z _v_i | _opt ×cos(π / 4);

[0105] L _v_i_opt =(|Z _v_i | _opt ×sin(π / 4) / ω _dom ;

[0106] Among them, R _v_i_opt L represents the optimal virtual resistance of node i; _v_i_opt Let ω be the optimal virtual inductance of node i; cos and sin are trigonometric functions; ω _dom The dominant oscillation angular frequency is given above. This set of parameters represents the optimal virtual impedance parameters, satisfying both the physical requirements of cascade dissipation and ensuring real-time computation.

[0107] As an example, an alternative implementation of adaptive dynamic execution based on energy feedback is described. A continuous adaptive control law based on energy dissipation rate feedback is proposed, which allows the virtual impedance to automatically adjust according to the real-time strength of the oscillation, achieving both flexibility and robustness in control. Specifically, this embodiment can be implemented using the following scheme:

[0108] Step 501: Calculate in real time the dissipation deviation between the current energy dissipation rate of each node and the target energy dissipation rate determined based on the cascade dissipation criterion.

[0109] During the execution phase, the controller monitors the actual status of each node in real time. Accordingly, it calculates the actual energy dissipation rate P of node i. _diss_actual_i (t). Optionally, the current I can be measured in real time. _osc_i (t) and the current virtual resistance value R _v_i (t) is calculated as follows:

[0110] P _diss_actual_i (t)=3×R _v_i (t)×I _osc_i (t) 2 .

[0111] Simultaneously, based on the cascade dissipation criterion and the current oscillation energy density, the target energy dissipation rate P that this node should provide is calculated. _diss_target_i (t). The target dissipation rate reflects the level of power required to prevent energy buildup. The difference between the two is calculated to obtain the dissipation deviation e. _diss_i (t)=P _diss_target_i (t)-P _diss_actual_i (t); when e _diss_i When (t)>0, the current virtual impedance is insufficient to provide the required dissipation, and energy is accumulating; conversely, it means that the damping is too large, which is beneficial to stability but may affect the accuracy of steady-state voltage.

[0112] Step 502: Based on the sensitivity matrix (also known as the sensitivity matrix) of the dissipation deviation and energy dissipation rate to the virtual impedance parameter, construct an adaptive control law based on energy dissipation rate feedback.

[0113] In this embodiment, an adaptive control law is designed to dynamically adjust the virtual impedance parameter Θ. _i =[R _v_i L _v_i ] T This control law is designed based on Lyapunov stability theory to ensure that the defined energy function (usually the square of the deviation) decreases monotonically with time. The constructed adaptive control law is shown below:

[0114] dΘ _i (t) / dt=Γ _gain ×S _diss_i T (t)×e _diss_i (t);

[0115] Where, dΘ _i (t) / dt represents the time derivative vector of the virtual impedance parameter; Γ _gain The positive definite adaptive gain matrix determines the adjustment speed; in other words, it can be determined through simulation and debugging based on actual control bandwidth requirements and system response speed requirements; S _diss_i (t) represents the energy dissipation rate P. _diss For the virtual impedance parameter Θ _i The sensitivity (gradient) vector, i.e., S _diss_i =ЯP _diss / ЯΘ _i ; T Я represents the matrix transpose, and Я corresponds to the partial derivative.

[0116] The control law adjusts the impedance parameter along the steepest descent direction (gradient direction) that minimizes dissipation deviation. Calculate S. _diss_i Normalization can be used in such cases. That is, in practical engineering implementation, a small regularization term is usually added to the denominator, i.e., an adaptive normalization law is used.

[0117] Step 503: Calculate the time derivative of the virtual impedance parameter in real time using the adaptive control law, and dynamically adjust the virtual impedance parameter of each node inverter based on this (time derivative) until the dissipation deviation converges to zero.

[0118] The local controller executes the above adaptive law in each control cycle to calculate the current required rate of change of parameters dΘ. _i / dt. Next, the integrator is used to update the current virtual impedance value Θ. _i (t+Δt), the specific formula can be expressed as:

[0119] Θ _i (t+Δt)=Θ _i (t)+(dΘ _i (t) / dt)×Δt;

[0120] Where Δt is the control period. The updated parameters are immediately sent to the inverter's underlying control loop (such as a voltage and current dual closed loop) to take effect.

[0121] Through a closed-loop feedback mechanism, when the oscillation suddenly intensifies, the dissipation deviation increases, driving the virtual resistance to increase rapidly, forming strong damping. When the oscillation subsides, the dissipation deviation tends to zero, and the virtual resistance automatically stabilizes at the minimum value required to maintain system stability. This mechanism solves the problem of insufficient damping under strong disturbances and overdamping affecting voltage under weak disturbances with fixed parameters.

[0122] As another example, an optional implementation process for signal modeling based on transfer function identification and causal analysis is provided. This embodiment employs signal processing and statistical analysis methods to construct a system model based on the correlation between input and output data, suitable for scenarios where accurate physical parameters are difficult to obtain or where the black-box characteristics of the system are of concern. Specifically, it includes:

[0123] Step 601: The chain oscillation propagation model is specifically constructed as a transfer function model that includes transfer gain, transfer damping ratio and transfer delay parameters.

[0124] In this alternative embodiment, the system does not focus on internal energy, but rather treats the relationship between adjacent nodes as a signal transmission system. Assume the power fluctuation ΔP of node i... _i The input is the power fluctuation ΔP at node i+1. _i+1 If it is the output, then the relationship between the two can be expressed using the transfer function H. _i (s) Description. Transfer function model H _i (s) is specifically formalized as a series connection of a second-order oscillatory element and a time-delay element, that is:

[0125] H _i (s)=K _gain_i ×(ω _n_i 2 / (s 2 +2×ζ _trans_i ×ω _n_i ×s+ω _n_i 2 ))×exp(-τ _delay_i ×s);

[0126] Among them, K _gain_i The transfer gain reflects the amplification factor of the oscillation amplitude; ζ _trans_i The transmission damping ratio reflects the attenuation characteristics during transmission; ω _n_i τ is the natural oscillation frequency; _delay_i The delay parameter is used to pass the time; s is the Laplace operator.

[0127] Step 602: Calculate the cross-power spectral density and self-power spectral density between adjacent nodes using the node power oscillation components in the synchronous transient dataset.

[0128] In other embodiments, this step may involve constructing a chain-like oscillation propagation model that reflects the coupling relationship between adjacent nodes based on the synchronous transient dataset and oscillation characteristic parameters, and quantifying the propagation characteristics of oscillation along the link.

[0129] Optionally, non-parametric spectral analysis methods are used to process the data. For the power oscillation time series of synchronously acquired nodes i and i+1, their self-power spectral density S is calculated. _xx (f) and cross-power spectral density S _xy (f).

[0130] S _xy (f)=E[X(f)×Y * (f)];

[0131] Where X(f) and Y(f) are the Fourier transforms of the input and output signals, respectively; * E[...] represents conjugate; E[...] represents expectation operation, which is usually achieved through multi-segment averaging.

[0132] Step 603: Based on the cross-power spectral density and the self-power spectral density, the frequency response function between adjacent nodes is calculated using the Hall first-order H1 estimation method. The Hall first-order H1 estimation method can also be simply referred to as the H1 estimation method.

[0133] Alternatively, using the H1 estimation method, assuming that the noise mainly exists at the output, the frequency response function FRF can be calculated. _i (f)=S _xy (f) / S _xx (f); This frequency response function describes the amplitude and phase frequency characteristics of the oscillation at different frequencies in a non-parametric form.

[0134] Step 604: Near the dominant oscillation frequency determined by the oscillation characteristic parameters, fit the frequency response function to a second-order oscillation system model to obtain the transmission gain, transmission damping ratio and transmission delay parameters of each adjacent link.

[0135] Optionally, a least-squares curve fitting algorithm is used to calculate the non-parametric FRF within the frequency band near the dominant oscillation frequency. _i (f) Approximating the parameterized model. K is identified by minimizing the fitting error. _gain_i ζ _trans_i and τ _delay_i The specific value.

[0136] The parameters above quantify the dynamic behavior of oscillations propagating along the link.

[0137] Step 605: Construct a vector autoregressive model for the power oscillation components of adjacent nodes and calculate the hysteresis coefficient matrix.

[0138] After establishing the model, Granger causality analysis was employed. A bivariate vector autoregressive (VAR) model was constructed, which can be described by the following formula:

[0139] y _t =∑A _k ×y _t-k +ε _t ;

[0140] Among them, y _t =[ΔP _i (t), ΔP _i+1 (t)] T Let ΔP be a vector containing power data from two nodes. _i (t) represents the power oscillation component of node i at time t, ΔP _i+1 (t) represents the power oscillation component of node i+1 adjacent to node i at time t; A _k The coefficient matrix is ​​a k-order lag matrix; ε _t Let y be the residual vector. _t-k The corresponding power oscillation component vector with lag order k.

[0141] Step 606: Perform an F-test on the hysteresis coefficient matrix to determine the Granger causality and significance level of oscillation propagation between adjacent nodes.

[0142] In other words, if statistical testing proves that the lag coefficient matrix is ​​non-zero and not due to random fluctuations in the data, it is called significantly non-zero; otherwise, it is not significant.

[0143] Analyze the coefficient matrix A using the F-test statistic. _k The significance of non-diagonal elements is determined by the following: if the coefficient from node i to node i+1 is significantly non-zero, and conversely, it is not significant, then a one-way Granger causal relationship from i to i+1 exists, meaning that node i is the cause of the oscillation of node i+1.

[0144] Step 607: Based on the Granger causality relationship and the differences in the oscillation start time in the synchronous transient dataset, identify the location of the dominant oscillation source node in the chain microgrid group.

[0145] Connecting the causal directions of all adjacent node pairs along the entire link forms a causal chain diagram. Typically, the node with the causal arrows pointing in opposite directions (i.e., all adjacent nodes are affected by it) is the dominant oscillation source node. Simultaneously, by using the starting time of the auxiliary reference oscillation waveform, the node where waveform distortion first occurs further corroborates the location of the source.

[0146] Step 608: Calculate the dominant oscillation frequency, waveform kurtosis coefficient, and autocorrelation decay rate of the power oscillation component.

[0147] To further refine the identification of oscillation sources, the system calculates statistical characteristics to distinguish oscillation types. The dominant frequency f is calculated. _dom Calculate the waveform kurtosis coefficient K. _kurtosis =E[(x-μ) 4 ] / σ 4 , used to measure the sharpness of a waveform; where K _kurtosis σ represents the kurtosis of the data sequence, a statistic describing the steepness and thickness of the data distribution's tails; x represents a single sample value in the data sequence; μ represents the population mean of the data sequence; and σ represents the population standard deviation of the data sequence. Calculate the autocorrelation decay rate R. _decay This reflects how quickly the correlation of a signal decreases over time.

[0148] Step 609: When the dominant oscillation frequency is in the preset low-frequency range and the autocorrelation decay rate is less than the preset randomness threshold, the oscillation source type is determined to be new energy power output fluctuation; when the dominant oscillation frequency is in the preset medium-frequency range and the waveform kurtosis coefficient satisfies the sinusoidal distribution characteristics, the oscillation source type is determined to be electromechanical oscillation.

[0149] Furthermore, the system has a pre-defined criterion database, as follows:

[0150] Fluctuations in new energy output (such as wind gusts) are characterized by low frequency (e.g., f). _dom <0.3Hz), and exhibits strong randomness, manifested by a rapid decay of the autocorrelation function (R0). _decay Less than the randomness threshold, for example, 0.5).

[0151] Electromechanical oscillations (such as the oscillation of a synchronous machine) are characterized by intermediate frequencies (e.g., 0.3 Hz). <f _dom <2.5Hz), the waveform is close to a standard sine wave, the kurtosis coefficient is close to 1.5 (the kurtosis of a sine wave), and the autocorrelation is well maintained.

[0152] Through multi-dimensional feature matching, the system not only tells maintenance personnel where the problem is, but also indicates what caused the problem and takes different coping strategies, such as smoothing control for new energy fluctuations and damping for electromechanical oscillations.

[0153] As another example, an alternative implementation method for decision-making based on modal sensitivity analysis and distributed optimization is provided. This method follows the modal analysis approach in control theory and is suitable for scenarios with high computational resource requirements but needing accurate control of the damping ratio of specific oscillation modes. Specifically, it demonstrates how to extract modes using an improved Prony algorithm, construct an optimization model based on the sensitivity matrix, and solve for the optimal parameters using the Distributed Alternating Direction Multiplier Method (ADMM). Further, this embodiment includes:

[0154] Step 701: Optimize the Prony algorithm using singular value decomposition to extract the dominant oscillation mode;

[0155] Construct a Hankel matrix based on a synchronous transient dataset, and perform singular value decomposition on the Hankel matrix;

[0156] The effective mode order is determined based on the cumulative energy contribution of singular values, and the linear prediction coefficients are solved using the least squares method.

[0157] The eigenvalues ​​are calculated based on linear prediction coefficients, and the frequency and damping factor of the dominant oscillation mode are extracted.

[0158] Before performing sensitivity calculations, accurate information about the dominant oscillation mode of the current system must be obtained. Although the STFT can provide frequency and amplitude, its accuracy in identifying damping factors is limited. Therefore, this embodiment employs an improved Prony algorithm.

[0159] Specifically, based on the power oscillation sequence x(n) in the synchronous transient dataset, the Hankel matrix H is constructed. _ankel Let the matrix have L rows and M columns. Typically, L > M and L + M = N. _sample (Total number of samples)

[0160] Accordingly, H _ankel =[x(0),x(1),...,x(M-1);x(1),x(2),...,x(M);...,;x(L-1),x(L),...,x(N _sample -1)];where, H _ankel is the constructed Hankel matrix; x(n) represents the oscillation data at discrete time points; the semicolon indicates a newline in the matrix.

[0161] Next, singular value decomposition (SVD) is performed on the Hankel matrix. That is:

[0162] H _ankel =U×Σ×V T ;

[0163] In the formula, U and V are orthogonal matrices; Σ represents the matrix containing singular values ​​σ. _i A diagonal matrix.

[0164] The effective mode order is determined based on the cumulative energy contribution of singular values. The energy proportion R of singular values ​​is calculated. _energy_i =(σ _i ) 2 / ∑(σ _k ) 2 ;σ _i Let be the i-th singular value, which is the eigenvalue obtained by singular value decomposition of the data matrix; ∑(σ _k ) 2The summation of the squares of all singular values ​​represents the total energy of the singular values; k is the index of the singular value, and all singular values ​​obtained from the decomposition are traversed. The orders corresponding to the top P singular values ​​whose cumulative proportion exceeds a preset threshold (e.g., 0.95) are selected as the effective mode order. This step eliminates the interference of noisy modes. Next, the least squares method is used to solve for the linear prediction coefficients, and the eigenvalues ​​λ are calculated. _i =α _i +j×ω _i Accurately extract the frequency ω of the dominant oscillation mode. _i With damping factor α _i Where j is the imaginary unit.

[0165] Step 702: Specifically set the global oscillation suppression index as the damping ratio increment of the dominant oscillation mode of the chain microgrid group.

[0166] In this embodiment, the optimization objective is set to maximize the damping ratio of the dominant mode. The damping ratio ζ of the current dominant mode is defined. _curr =-α _i / sqrt((α _i ) 2 +(ω _i ) 2 The optimization objective is to achieve a target damping ratio ζ after control is applied. _target As large as possible, or the damping ratio increment Δζ=ζ _target -ζ _curr maximize.

[0167] Step 703: Based on the oscillation characteristic parameters, calculate the sensitivity matrix of the damping ratio increment to the virtual resistance and virtual inductance parameters of each node.

[0168] Accordingly, the sensitivity matrix S needs to be calculated. _zeta This matrix describes the extent to which small changes in the virtual impedance parameters of each node affect the dominant mode damping ratio.

[0169] Based on the eigenvalue perturbation theory, the sensitivity calculation formula is as follows:

[0170] Яλ _i / ЯX _v_k =(w _left_i T ×(ЯA _sys / ЯX _v_k )×v _right_i ) / (w _left_i T ×v _right_i );

[0171] In the formula, Яλ _i / ЯX _v_k For the eigenvalue λ _iFor the virtual impedance parameter X of the k-th node _v_k (Including R) _v_k and L _v_k The partial derivative of ); w _left_i and v _right_i The eigenvalues ​​are λ. _i The corresponding left and right eigenvectors; A _sys This is the system state matrix.

[0172] Based on this, the damping ratio sensitivity S is obtained using the chain rule. _zeta_k ,Right now:

[0173] S _zeta_k =[Яζ _i / ЯR _v_k ,Яζ _i / ЯL _v_k ];ζ _i Let S be the damping ratio of the i-th oscillation mode of the system; the matrix S is... _zeta_k This quantifies how much the damping ratio of the dominant mode will increase for every unit increase in the virtual resistance and virtual inductance of the k-th node. Typically, nodes located at the antinodes of oscillations have high sensitivity and are the main nodes for suppressing oscillations.

[0174] Step 704: Using the sensitivity matrix, construct a virtual impedance optimization model with the objective function of maximizing the weighted damping ratio increment, and satisfying the small-signal stability margin constraint of the inverter and the impedance coordination constraint of adjacent nodes.

[0175] Based on the sensitivity matrix, the nonlinear damping control problem is linearized into the following optimization model:

[0176] The objective function is maxJ = ∑(W) _k ×S _zeta_k ×ΔΘ _k Its constraints include:

[0177] Parameter limiting constraint, Θ _min ≤Θ _base +ΔΘ _k ≤Θ _max ;

[0178] Chain coordination constraints, ΔR _v_k+1 ≥ΔR _v_k This is similar to cascaded dissipation, preventing energy from accumulating.

[0179] Small-signal stability margin constraints ensure that parameter adjustments do not disrupt the stability of the underlying voltage and current loops.

[0180] In the formula, W _k ΔΘ is the weight coefficient for node k, which is usually set according to the node's capacity or location; _kLet Θ be the virtual impedance increment vector to be solved. _base Θ is the original fundamental value of the virtual impedance. _min Θ _max These are the lower and upper threshold values ​​for the virtual impedance parameter, respectively, limiting the range of values ​​that can be adjusted for the parameter, ΔR. _v_k+1 This represents the virtual resistance increment of the adjacent node k+1.

[0181] Step 705: Decompose the virtual impedance optimization model into local optimization subproblems corresponding to each node, and define boundary coordination variables between adjacent nodes.

[0182] Optionally, this embodiment employs the Alternating Direction Multiplier Method (ADMM). The global optimization problem is decomposed into N decoupled local subproblems. For node k, its local state variable x is defined. _k (i.e., local virtual impedance increment) and boundary coordination variable z that interacts with neighboring nodes. _k Constructing the augmented Lagrangian function L _aug_k , can be represented as:

[0183] L _aug_k =f _k (x _k )+y _k T ×(A _k ×x _k -z _k )+(ρ / 2)×||A _k ×x _k -z _k || 2 ;

[0184] In the formula, f _k (x _k y is the local objective function of the node; _k For dual multipliers; ρ is the penalty factor; A _k Let be the coupling matrix, describing the constraint relationship between node k and k-1, k+1, where ||…|| represents the Euclidean norm and L2 norm.

[0185] Step 706: Control each node to solve the local optimization subproblem in parallel to update the local virtual impedance parameters, and exchange boundary coordination variables with adjacent nodes to update the dual variables.

[0186] The ADMM algorithm is executed in three steps, which are repeated in each iteration step t:

[0187] x-update step (parallel), each node has fixed z and y, solve for minL _aug_k Get x _k_t+1 Since the problem has been decomposed, it is usually a low-dimensional quadratic programming problem, and the solution speed is relatively fast; x_k_t+1 This represents the optimal solution of the x variable for node k at iteration step t+1;

[0188] Communication exchange, each node will update x _k_t+1 Send to physically adjacent nodes;

[0189] z-Update Step (Coordination): Each node collects neighbor information and calculates the boundary coordination variable z. _k_t+1 This step involves calculating the average consistency value that satisfies the chain coordination constraints.

[0190] y-Update step (dual), update dual multiplier y _k_t+1 =y _k_t +ρ×(A _k ×x _k_t+1 -z _k_t+1 ).

[0191] Step 707: Iteratively execute the steps of solving the local optimization subproblem and updating the boundary coordination variables until the local virtual impedance parameters converge and the optimal virtual impedance parameters are obtained.

[0192] The system continuously executes an iterative process. In each iteration, the original residual r is calculated. _prim =||A×xz|| and dual residual r _dual =||ρ×A T ×(z _t+1 -z _t )||.

[0193] In the above, A is the global coupling matrix; x is the set of original optimization variables for each node, which is the solution result of the x-update step of the ADMM algorithm; z is the global coupling variable of the ADMM algorithm, used to associate the coupling constraints of each node; ρ is the penalty factor, used to adjust the penalty intensity for deviations in coupling constraints; z _t z _t+1 These are the values ​​of the global coupling variables at iteration steps t and t+1.

[0194] The algorithm is considered to have converged when both the original residual and the dual residual are less than the preset convergence threshold.

[0195] At this point, each node obtains x. _k This is the globally optimal virtual impedance parameter.

[0196] Among them, the ADMM algorithm utilizes the sparsity of chain topology and can usually converge within a few dozen iterations, meeting the requirements of near real-time control.

[0197] According to one aspect of this application, an exemplary scheme for offline stability boundary construction and smooth switching control is described, wherein a stability boundary lookup table solves the problem of complex constraint handling in online optimization. Specifically, the inverter small-signal stability margin constraint is limited by querying a pre-built stability boundary lookup table, which is pre-built through the following steps:

[0198] Step 801: For the typical operating conditions of the chain microgrid group, calculate the gain crossover frequency and phase crossover frequency of the open-loop transfer function of the inverter of each node.

[0199] The small-signal stability of an inverter is typically measured using the phase margin (PM) and gain margin (GM) in a Bode plot. However, PM and GM have complex nonlinear relationships with virtual impedance parameters, making them difficult to directly incorporate into online optimization constraints. Therefore, this embodiment employs an offline pre-calculation strategy.

[0200] Accordingly, typical operating conditions of the microgrid cluster are traversed (e.g., different load levels, different grid impedances). For each operating point, an open-loop transfer function of the inverter including virtual impedance is established. Through numerical analysis, the gain crossover frequency ω when the amplitude of the open-loop transfer function is 1 is calculated. _gc And the phase crossover frequency ω when the phase is -180 degrees. _pc .

[0201] Step 802: Based on the gain crossover frequency, the preset phase margin index is converted into a linear upper bound constraint that limits the ratio of virtual inductance to virtual resistance.

[0202] To ensure that the system has sufficient damping, the phase margin PM is typically required to be greater than or equal to the lower limit of the phase margin PM. _min For example, 45 degrees. Studies have shown that the virtual inductance L _v This introduces phase lag, reducing PM. At ω _gc At this point, to meet the phase margin requirement, the ratio of virtual inductance to virtual resistance must be limited. Through a first-order Taylor expansion or geometric approximation, this constraint can be transformed into R... _v -L _v Linear constraints on a plane, namely:

[0203] L _v ≤k _slope ×R _v +b _intercept ;

[0204] In the formula, k _slope and b _intercept Is with ω _gc and PM _minThe relevant coefficients. This inequality shows that if the virtual inductance is to be increased, the virtual resistance must be increased accordingly to maintain the phase margin.

[0205] Step 803: Based on the phase crossover frequency, the preset gain margin index is converted into a secondary lower bound constraint that limits the weighted magnitude of the virtual resistance and virtual inductance.

[0206] Among them, the gain margin GM measures the robustness of the system to changes in gain, and it is generally required that GM ≥ the lower limit of the gain margin GM. _min For example, 3dB. At ω _pc At this point, the open-loop gain of the system must be less than 1 / GM. _min The virtual impedance must have a sufficiently large magnitude to attenuate gain in the high-frequency range. This requirement translates into a quadratic constraint (outside-circle constraint) on the parametric plane, namely:

[0207] (R _v ) 2 +(ω _pc ) 2 ×(L _v ) 2 ≥(K _gm ) 2 ;

[0208] In the formula, K _gm It is by GM _min A defined lower bound for the minimum impedance magnitude. This constraint requires that the virtual impedance parameter must lie outside an elliptical region centered at the origin.

[0209] Step 804: Determine the feasible region of the virtual impedance parameter based on the linear upper bound constraint and the quadratic lower bound constraint, and discretize the feasible region to generate a stability boundary lookup table.

[0210] Based on this, a safe and feasible region for the virtual impedance parameter is formed, typically a sector-shaped area or a region with a corner cut off. To facilitate online table lookup, operating parameters (such as output power) are used as indexes, and the corresponding constraint boundary parameters (k...) are... _slope b _intercept K _gm It is stored in the microprocessor's non-volatile memory.

[0211] During online operation, the controller reads the boundary parameters based on the current operating conditions and adds them to the optimization model as linear or quadratic inequality constraints, ensuring that the calculated parameters will never compromise the stability of the inverter itself.

[0212] Step 805: Predict the zero-crossing point of the phase of the oscillation waveform at each node.

[0213] This step describes a smooth switching strategy based on discrete events. Direct parameter step switching may induce new transient shocks, especially when the current amplitude is large. The optimal switching time is at the zero-crossing point of the oscillating current, when the inductor energy storage is minimal and the magnetic flux disturbance caused by the switching is minimized.

[0214] The system utilizes a phase-locked loop (PLL) or a linear predictor, based on the current angular frequency ω. _dom and phase φ _i Predict the nearest one or more current zero-crossing moments in the future, i.e.:

[0215] t _zero_i =t _curr +(n×π-φ _i ) / ω _dom ;

[0216] In the formula, n is an integer used to select future zero-crossing points, and π is the constant of pi.

[0217] Step 806: Based on the propagation delay of the oscillation along the link, arrange the switching sequence of each node in the order of the oscillation propagation direction.

[0218] Considering the propagation characteristics of oscillations in a chain-like microgrid, if all nodes switch simultaneously, overvoltage may occur in the middle of the link due to the superposition of control actions. Therefore, this embodiment adopts a delayed switching strategy.

[0219] Assume the oscillation propagates from node 1 to node N, and the propagation delay between adjacent nodes is τ. _delay Then the switching time t of each node _switch_i It should meet the following requirements:

[0220] t _switch_i+1 ≥t _switch_i +τ _delay ;

[0221] Based on this timing constraint and the zero-crossing point, the system selects the moment closest to the target timing from the candidate zero-crossing points to generate the final switching sequence.

[0222] Step 807: At the zero-crossing point of each node phase determined by the switching sequence, the virtual impedance parameters are updated to minimize the transient disturbances introduced by the switching process.

[0223] When the system clock reaches the predetermined t _switch_i At this time, the controller of node i immediately updates the virtual impedance coefficient. Because the current i at this time... _osc ≈0, although the voltage jump caused by parameter jumps cannot be eliminated, the jump caused by the resistive component is minimized, and the continuity of the inductor current is best maintained.

[0224] In summary, the timing control strategy smooths out the control intervention process and avoids the risk of secondary disturbances caused by treating the disease.

[0225] According to another aspect of this application, an exemplary scheme for online correction and multi-round correction strategy of sensitivity matrix is ​​described, especially the closed-loop correction mechanism when the control effect is not ideal. It can be used to solve the problem of inaccurate sensitivity calculation caused by model parameter drift. By correcting the sensitivity matrix online and performing multi-round correction, the system eventually converges to a stable state.

[0226] Step 901: Monitor the control effect in real time and calculate the deviation between the actual damping ratio increment and the expected damping ratio increment.

[0227] After the virtual impedance parameters are issued and executed, the system enters the monitoring and evaluation phase. An evaluation time window is set, for example, T_eval = 500ms. Within this window, the Prony algorithm is used to re-identify the current dominant mode damping ratio ζ. _actual .

[0228] At the same time, the expected target damping ratio ζ calculated by the optimization model is read. _expected Calculate the deviation ratio or correction factor κ between the two. _correct The corresponding formula is:

[0229] κ _correct =Δζ _actual / Δζ _expected ;

[0230] Where, Δζ _actual =ζ _actual -ζ _initial ;ζ _initial The initial damping ratio before control; Δζ _expected To optimize the increment of the model's predictions. If κ _correct A value less than 1, such as <0.8, indicates that the actual system response is less sensitive than the model prediction, i.e., the sensitivity matrix S... _zeta There is a discrepancy.

[0231] Step 902: Correct the sensitivity matrix online based on the deviation ratio to reflect the difference between the actual system and the mathematical model.

[0232] Accordingly, the system uses the aforementioned correction factor to perform a linear correction on the sensitivity matrix, that is:

[0233] S _zeta_corrected =κ _correct ×S _zeta ;

[0234] Among them, S _zeta_corrected This is the corrected sensitivity matrix. This step is equivalent to calibrating the controller's feel in real time based on the actual input-output ratio.

[0235] Step 903: Calculate the incremental correction amount of the virtual impedance based on the corrected sensitivity matrix, and solve it using the pseudo-inverse method.

[0236] Furthermore, using the corrected sensitivity matrix, the parameter fine-tuning required to suppress the residual damping deviation is calculated. The target residual damping ratio increment Δζ is defined. _residual =ζ _target -ζ _actual Construct the following incremental equation and solve it:

[0237] ΔX _v_correct =(S _zeta_corrected ) + ×Δζ _residual ;

[0238] Where, ΔX _v_correct This is the vector of additional corrections for the virtual impedance; + This represents the pseudo-inverse operation of a matrix, used to handle cases where a system of equations is overdetermined or underdetermined.

[0239] To prevent excessive correction from causing new instability, the calculated ΔX _v_correct Amplitude limiting: |ΔX _v_correct |≤ΔX _max_correct ;ΔX _max_correct The upper limit of the vector of additional corrections corresponding to virtual impedance.

[0240] Step 904: Generate a correction instruction containing the fast ramp time, and perform multiple rounds of iterative correction until the oscillation suppression meets the requirements.

[0241] The calculated additional correction amount is packaged into a command and issued. Accordingly, this correction uses a shorter ramp buffer time T. _ramp_correct For example, setting it to 20ms is better than the usual 50ms.

[0242] Among them, Cmd _vi_correct_i ={ΔR _v_i_correct ΔL _v_i_correct , t _switch_correct T _ramp_correct};

[0243] Or rather, Cmd _vi_correct_i The virtual impedance supplementary correction amount, ΔR, is packaged and sent to the i-th node of the system. _v_i_correct The virtual resistance correction amount ΔL represents at node i. _v_i_correct The virtual inductance addition correction amount represents node i, t _switch_correct T represents the correction trigger time. _ramp_correct This represents the correction slope buffer time.

[0244] After executing the instruction, the system returns to the monitoring and evaluation phase. If the oscillation persists, meaning the evaluation level remains "moderate," the above correction and adjustment process is repeated. A maximum number of iterations is set, for example, three, to prevent the system from entering an infinite loop. This multi-round iterative mechanism of monitoring-correction-re-execution constitutes a robust defense for the system.

[0245] In some embodiments, parts of the method of the present invention may also be:

[0246] Optionally, the analysis window length is set to 200ms, containing 2000 sampling points, the window sliding step is 20ms, and the Hanning window function w is used. _hann (n) Weighted, calculate the short-time Fourier transform (STFT) for the power oscillation component of the i-th node. _i (t, f), the specific formula is:

[0247] STFT _i (t, f) = Σ n=0 N_win-1 Δp _i (t+n×T _s )×w _hann (n)×exp(-j2πfnT _s );

[0248] Obtain the time-frequency spectrum matrix STFT _i (t, f), frequency resolution Δf = f _s / N _win =5Hz;

[0249] The above represents the short-time Fourier transform value of the i-th node at time t and frequency f, N _win T is the length of the window function. _s For the sampling period, exp(-j2πfnT) _s ) represents the complex exponential term of the Fourier transform, j is the imaginary unit, and Δp _i (t+n×T _s ) represents the i-th node at t+n×T _s The power oscillation component at time f _s The sampling frequency.

[0250] Furthermore, the window function length for transient power oscillations is defined to focus on the frequency band, which can be f. _osc ∈[0.5Hz, 50Hz], calculate the instantaneous energy E' within this frequency band. _osc_i (t), is:

[0251] E' _osc_i (t)=Σ f∈[0.5，50] |STFT _i (t, f)| 2 ×Δf;

[0252] Simultaneously calculate the total energy E across the entire frequency band. _total_i (t), to obtain the energy ratio R of the oscillation frequency band. _osc_i (t), that is:

[0253] R _osc_i (t)=E' _osc_i (t) / E _total_i (t);

[0254] Optionally, an oscillation detection threshold system can be set, as follows:

[0255] Energy threshold E _th =0.01×P _rated_i 2 ;P _rated_i Let i be the rated power of node i;

[0256] Energy percentage threshold R _th =0.3;

[0257] Duration threshold T _duration_th =100ms;

[0258] Specifically, the oscillation event at time t is determined when the following conditions are met:

[0259] Condition C1, E' _osc_i (t)>E _th ;

[0260] Conditions C2, R _osc_i (t)>R _th ;

[0261] The duration for which conditions C3, C1, and C2 are simultaneously satisfied exceeds T. _duration_th .

[0262] Optionally, traverse all nodes and time windows to generate an oscillation event flag matrix. _osc (i, t); if Flag _osc (i, t) = 1, indicating that oscillation has been detected. If Flag _osc (i, t) = 0 indicates that no oscillation was detected.

[0263] Determine the start time t of the oscillation event _start_osc With end time t _end_osc Output the oscillation event detection results.

[0264] Optionally, the analysis period [t] after the onset of the oscillation event can be selected. _start_osc , t _start_osc +T _prony ], where T _prony=500ms, for the power oscillation component Δp at the i-th node _i (t) Applying the Prony algorithm for mode decomposition, the signal is modeled as the sum of multiple decaying sinusoidal components, i.e.:

[0265] Δp _i (t)≈Σ k=1 M A _ik ×exp(σ _ik ×t)×cos(2πf _ik ×t+φ _ik );

[0266] Where M is the modal order, we can set M=6, A _ik Let σ be the amplitude of the k-th mode. _ik f is the attenuation factor. _ik φ is the oscillation frequency. _ik This is the initial phase.

[0267] Based on this, a data matrix is ​​constructed, the eigenvalue problem is solved, and least squares fitting is performed to obtain the modal parameters.

[0268] Furthermore, the energy contribution E of each mode is calculated. _mode_ik It can be described by the following formula:

[0269] E _mode_ik =(A _ik ) 2 / Σ k=1 M (A _ik ) 2 ;

[0270] Among them, the energy contribution of the screening is greater than the threshold E. _mode_th The mode with a value of 0.1 is selected as the dominant mode, and a maximum of 3 dominant modes are retained.

[0271] For the k-th dominant mode, calculate its damping ratio ζ. _ik =-σ _ik / sqrt((σ _ik ) 2 +(2πf _ik ) 2 );

[0272] When ζ _ik A value less than 0 indicates negative damping, corresponding to oscillatory divergence. _ik A value greater than 0 indicates positive damping, corresponding to oscillation decay.

[0273] Compare the dominant mode frequencies identified at different nodes, when |f _ik -f _jk |<Δf_tol (Δf) _tol When the frequency is 0.5Hz, it is assumed that node i and node j have the same oscillation mode; the set of participating nodes for each oscillation mode in the chained microgrid is counted. The set of oscillation mode parameters is output. Among them, f _ik f _jk The oscillation frequencies corresponding to nodes i and j are Δf, respectively. _tol Corresponding frequency tolerance threshold.

[0274] Optionally, based on the oscillation event detection results and the oscillation mode parameter set, the spatiotemporal distribution characteristics of the oscillation in the chain microgrid group are analyzed, and the power oscillation component Δp of each node is determined. _i Perform a Hilbert transform on (t) to obtain the analytic signal Δp. _i_analytic (t)=Δp _i (t)+j×H[Δp _i (t)]; H[…] corresponds to the Hilbert transform operation.

[0275] Extracting the oscillation amplitude envelope A _env_i (t)=|Δp _i_analytic (t)|;

[0276] Extracting the instantaneous phase φ _inst_i (t)=arg[Δp _i_analytic [(t)]; arg[…] represents the phase angle function.

[0277] In summary, this scheme employs chain-like energy flow conservation modeling and a tiered dissipation criterion. A physical energy flow model is established, quantifying the propagation attenuation factor of oscillation energy along the chain, from which the minimum impedance increase factor of downstream nodes relative to upstream nodes (analytical inequality) is derived. This spatial gradient configuration strategy forces downstream nodes to possess stronger energy throughput and dissipation capabilities than upstream nodes, constructing a progressively increasing energy dissipation barrier and suppressing the risk of standing wave superposition and tail-end divergence during oscillation propagation along the chain.

[0278] Building upon this foundation, an adaptive control law based on energy dissipation rate feedback and an optimal impedance angle decoupling and dimensionality reduction technique are introduced. The bivariate nonlinear optimization is decoupled into a single-variable modulus optimization. A continuous-time dynamic feedback loop is constructed using Lyapunov stability theory. The controller can monitor the deviation between the actual dissipated power and the target value in real time and automatically adjust the virtual impedance within milliseconds. This mechanism enables the control strategy to instantaneously follow the intensity of fault disturbances, allowing sufficient damping torque to be output shortly after a fault occurs.

[0279] The optional embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details of the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solution of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for coordinated control of fault propagation blocking and transient stability in chain-type microgrid groups, characterized in that, include: High-frequency electrical quantity data of each node in the chain microgrid group are collected, and multi-node spatiotemporal alignment processing is performed to generate a synchronous transient dataset. Based on synchronous transient datasets, transient power oscillation events are detected, and oscillation feature parameters are extracted. Based on the oscillation characteristic parameters, a chain oscillation propagation model reflecting the coupling relationship between adjacent nodes is constructed to quantify the propagation characteristics of oscillation along the link; Based on the chain oscillation propagation model, a virtual impedance optimization model with chain coordination constraints is established. With the goal of maximizing the global oscillation suppression index, the optimal virtual impedance parameters of each node are calculated. The optimal virtual impedance parameters are converted into control commands, and the virtual impedance is dynamically reconstructed for each node inverter to block the propagation of oscillations.

2. The method according to claim 1, characterized in that, Detect transient power oscillation events and extract oscillation characteristic parameters, including: The time-frequency energy spectrum of transient power is calculated by processing synchronous transient datasets using sliding window short-time Fourier transform. Calculate the energy proportion of the time-frequency energy spectrum within the oscillation frequency band, and determine that a transient power oscillation event has occurred when the energy proportion continuously exceeds the trigger threshold; The dominant oscillation mode frequency, oscillation amplitude envelope, and node phase distribution of transient power oscillation events are extracted as oscillation characteristic parameters.

3. The method according to claim 1, characterized in that, Construct a chain-like oscillation propagation model that reflects the coupling relationship between adjacent nodes, including: The chain oscillation propagation model is specifically constructed as a transfer function model that includes parameters such as transfer gain, transfer damping ratio, and transfer delay. Using the node power oscillation components in the synchronous transient dataset, the cross power spectral density and self power spectral density between adjacent nodes are calculated; Based on the cross-power spectral density and the self-power spectral density, the frequency response function between adjacent nodes is calculated using the Hall first-order estimation method; Near the dominant oscillation frequency determined by the oscillation characteristic parameters, the frequency response function is fitted to a second-order oscillation system model to obtain the transmission gain, transmission damping ratio, and transmission delay parameters of each adjacent link.

4. The method according to claim 3, characterized in that, The quantification of the propagation characteristics of oscillations along the link also includes oscillation source identification based on Granger causality tests, specifically: A vector autoregression model is constructed for the power oscillation components of adjacent nodes, and the hysteresis coefficient matrix is ​​calculated. An F-test was performed on the hysteresis coefficient matrix to determine the Granger causality and significance level of oscillation propagation between adjacent nodes; By combining Granger causality and the differences in oscillation start times in the synchronous transient dataset, the location of the dominant oscillation source node in the chain microgrid group is identified.

5. The method according to claim 1, characterized in that, A virtual impedance optimization model with chained coordination constraints is established to maximize the global oscillation suppression index, including: The global oscillation suppression index is specifically set as the damping ratio increment of the dominant oscillation mode of the chain microgrid group; Based on the oscillation characteristic parameters, the sensitivity matrix of the damping ratio increment to the virtual resistance and virtual inductance parameters of each node is calculated. Using the sensitivity matrix, a virtual impedance optimization model is constructed with the objective function of maximizing the weighted damping ratio increment, while satisfying the small-signal stability margin constraint of the inverter and the impedance coordination constraint of adjacent nodes.

6. The method according to claim 5, characterized in that, The optimal virtual impedance parameters for each node are calculated, specifically by using the distributed alternating directional multiplier method to solve the virtual impedance optimization model: The virtual impedance optimization model is decomposed into local optimization subproblems corresponding to each node, and boundary coordination variables between adjacent nodes are defined. Control each node to solve local optimization subproblems in parallel to update local virtual impedance parameters, and exchange boundary coordination variables with adjacent nodes to update dual variables; The local optimization subproblem solution and boundary coordination variable update steps are executed iteratively until the local virtual impedance parameters converge, thus obtaining the optimal virtual impedance parameters.

7. The method according to claim 1, characterized in that, Construct a chain-like oscillation propagation model that reflects the coupling relationship between adjacent nodes, including: The chain oscillation propagation model is specifically constructed as a chain energy flow conservation model; The instantaneous oscillation energy density of each node is calculated based on the node voltage oscillation components, current oscillation components, and dominant oscillation angular frequency in the oscillation characteristic parameters of the synchronous transient dataset. Based on the tie line impedance parameters and oscillation phase difference between adjacent nodes, an energy flow rate equation describing the transfer of oscillation energy between adjacent nodes is established. Based on the chain topology, the instantaneous oscillation energy density and energy flow rate equations of each node are integrated into a state space equation, and a tridiagonal state matrix reflecting the chain-like adjacent coupling characteristics is constructed.

8. The method according to claim 7, characterized in that, A virtual impedance optimization model with chained coordination constraints is established to maximize the global oscillation suppression index, including: The global oscillation suppression index is specifically set as the full-link oscillation energy decay rate, which is determined by the energy dissipation power generated by the virtual impedance of each node. A cascade dissipation criterion is established as a chain coordination constraint, which limits the ratio of the virtual impedance of the downstream node to the virtual impedance of the upstream node along the oscillation propagation direction. Based on the propagation attenuation factor obtained from the propagation characteristics of quantized oscillations along the link, the minimum impedance increment gradient required to maintain the cascade dissipation criterion is determined, thereby constraining the search space for the optimal virtual impedance parameters.

9. The method according to claim 8, characterized in that, Calculate the optimal virtual impedance parameters for each node, including: The virtual impedance optimization model is decoupled, and the virtual impedance angle parameter of each node is fixed to the analytical optimal angle that maximizes the equivalent damping conductance. The analytical optimal angle is 45 degrees at the dominant oscillation frequency. Based on the fixed virtual impedance angle parameter, the virtual impedance optimization model is transformed into a single-variable distribution optimization problem for the virtual impedance magnitude. Solve the univariate distribution optimization problem to obtain the optimal virtual impedance magnitude of each node that satisfies the cascade dissipation criterion, and generate the optimal virtual impedance parameters by combining the analytical optimal angle.

10. The method according to claim 8, characterized in that, Perform virtual impedance dynamic reconfiguration on each node inverter, including: Real-time calculation of the dissipation deviation between the current energy dissipation rate of each node and the target energy dissipation rate determined based on the cascade dissipation criterion; Based on the sensitivity matrix of dissipation deviation and energy dissipation rate to virtual impedance parameters, an adaptive control law based on energy dissipation rate feedback is constructed. The time derivative of the virtual impedance parameter is calculated in real time using an adaptive control law, and the virtual impedance parameter of each node inverter is dynamically adjusted accordingly until the dissipation deviation converges to zero.

Images