Thermal power coupling compressed air energy storage power station collaborative voltage regulation method

Abstract

The application discloses a kind of thermal power coupling compressed air energy storage power station coordination voltage regulating method, comprising: response switching request, by accurate calculation and fusion static, transient and dynamic three kinds of stable boundary, determine a comprehensive stable boundary, according to the boundary planning generates a piecewise optimization target loss of excitation trajectory of guiding CAES motor excitation current smooth decline;In the process that CAES motor executes the loss of excitation trajectory, based on the system model of online identification, the excitation current of thermal power unit is adjusted in coordination, to compensate the reactive power reduced by CAES motor due to loss of excitation in real time, dynamically, realize the constant of system total reactive power output;When monitoring that the reactive power output of CAES motor meets the preset near-zero off-grid condition, accurately control its circuit breaker action.The application changes the hard switching of CAES unit into active, controllable soft landing process, eliminates reactive power impact, improves the smoothness of switching process and the transient stability of grid voltage.

Description

Technical Field

[0001] This invention relates to the field of power system automatic control technology, and in particular to a method for coordinated voltage regulation in a thermal power coupled compressed air energy storage power station. Background Technology

[0002] As the global energy structure transitions towards a higher proportion of renewable energy, the demand for highly flexible regulation resources in the power system is becoming increasingly urgent. Coal-fired power plants coupled with compressed air energy storage (FC-CAES), as a novel hybrid energy storage technology integrating the stable power support capabilities of traditional thermal power units with the rapid start-up, shutdown, and energy time-shifting characteristics of compressed air energy storage systems (CAES), have demonstrated enormous application potential in multiple fields such as grid peak shaving, frequency regulation, and voltage support. This technology, through deep coupling of thermal power and energy storage units at the thermodynamic and electrical levels, improves the operational flexibility and economy of traditional thermal power plants, enabling them to better adapt to the volatility and intermittency of renewable energy generation. Therefore, researching and developing advanced control strategies for such complex coupled systems, especially control technologies that enable rapid and smooth switching between multiple operating modes (such as energy storage, energy release, and independent operation), is of crucial research significance and technological value for ensuring the safe and stable operation of high-proportion renewable energy power systems.

[0003] Currently, research on the operation and control of thermal power coupled compressed air energy storage (CAES) power plants mainly focuses on the overall energy management, efficiency optimization, and collaborative operation strategies under steady-state or quasi-steady-state conditions. At the voltage control level, existing technical solutions typically employ the automatic voltage regulation (AVR) logic of traditional synchronous generator sets, passively responding to and correcting voltage deviations at the generator terminals through the excitation system. When mode switching of the CAES unit is involved, such as when the compressor motor needs to be disconnected from the grid upon completion of energy storage, the standard operating procedure is usually to execute a direct circuit breaker tripping command. Changes in reactive power on the grid side rely on the thermal power unit's own AVR system or other dynamic reactive power compensation devices (such as SVC and SVG) for post-event compensation. At the collaborative control level, some solutions propose coordination strategies based on fixed-sequence logic, where after the CAES unit operates, a preset delay is allowed before instructing the thermal power unit to perform the corresponding compensation operation. In addition, some studies have applied advanced control algorithms such as Model Predictive Control (MPC) to optimize and regulate voltage fluctuations during normal system operation. By continuously optimizing the excitation commands of thermal power units and CAES units, they aim to achieve a comprehensive optimization of voltage deviation and regulation time. These methods together constitute the basic technical framework for ensuring the operation of such coupled power plants.

[0004] However, existing technologies still face profound technical challenges when dealing with the specific and drastic transient process of CAES unit mode switching. These challenges mainly stem from the passive response nature of their control philosophy, leading to an inherent contradiction between the instantaneous performance of the system during switching and the stability of the process. Specifically, this manifests in the following interrelated technical problems: Existing technologies cannot eliminate the instantaneous reactive power surges and voltage transient stability risks caused by switching operations. Furthermore, the lack of precise quantification of the risk of loss of synchronization makes it difficult for control strategies during switching to balance speed and safety. Summary of the Invention

[0005] Objective: This invention provides a method for coordinated voltage regulation of thermal power coupled compressed air energy storage (CAES) power plants, aiming to solve the grid voltage stability problem caused by the switching of compressed air energy storage (CAES) unit modes.

[0006] According to one aspect of the present invention, a method for coordinated voltage regulation in a thermal power coupled compressed air energy storage power station includes:

[0007] In response to the switching request, the target demagnetization trajectory of the compressed air energy storage motor is planned and generated.

[0008] The target demagnetization trajectory is executed, and the excitation current of the thermal power unit is adjusted in coordination to compensate for the reactive power changes during the demagnetization process of the compressed air energy storage motor.

[0009] When the reactive power output of the compressed air energy storage motor meets the preset disconnection conditions, it will be disconnected from the grid.

[0010] The planned target demagnetization trajectory for the compressed air energy storage motor includes:

[0011] Collect initial operating state data including initial excitation current, and obtain active power and initial power angle from it;

[0012] Calculate the initial internal potential based on the initial excitation current;

[0013] Based on active power, initial power angle and initial internal potential, three stability boundaries are calculated, including static, transient and dynamic stability boundaries.

[0014] The three stability boundaries are fused to determine the comprehensive stability boundary, and the target demagnetization trajectory is generated based on the comprehensive stability boundary.

[0015] The calculation of the transient stability boundary includes:

[0016] Determine the power angle stability margin based on the preset critical power angle and initial power angle;

[0017] By combining the power angle stability margin, the preset unit inertia constant, and the transient time constant, the transient stability boundary is calculated.

[0018] Among them, the three stability boundaries are integrated to determine the comprehensive stability boundary, including:

[0019] Preset confidence weights are assigned to static, transient, and dynamic stable boundaries, and then weighted and fused to obtain a weighted average boundary.

[0020] Choose the minimum value among the three stable boundaries as the minimum boundary;

[0021] The comprehensive stability boundary is determined based on the weighted average boundary and the minimum boundary.

[0022] The planning and generation of the target demagnetization trajectory for the compressed air energy storage motor is a piecewise function constructed by sequentially connecting the following parts:

[0023] The main descent stage is set to execute the main excitation reduction at a basically constant rate.

[0024] A transition section used to mitigate rate changes;

[0025] And a terminal approximation section used to make the rate of change of the excitation current approach zero when it reaches the target value.

[0026] The main decreasing segment is composed of a linear function i f (t)=i f0 -k fast *t defines the boundary, where i f (t) represents the target excitation current at time t, i f0 k is the initial excitation current obtained from the initial operating state data. fast The demagnetization rate is determined based on the comprehensive stability boundary.

[0027] The transition segment is composed of the exponential function i f (t)=i f1 *exp(-(t-t1) / τ smooth )+i f_target Define, where i f1 t1 and t2 represent the excitation current value and time point at the end of the main drop phase, respectively, and τ smooth i is the smoothing time constant. f_target The target excitation current is the transition current;

[0028] The terminal approximation segment is composed of the cosine function i f (t)=i f_min +(i f2 -i f_min )*cos(π*(t-t2) / (2*(t final -t2))) delimits, where i f2 t2 and t2 are the excitation current value and time point at the end of the transition section, respectively. f_mint is the preset minimum excitation current. final This is the point in time when the demagnetization process ends.

[0029] The process, after planning and generating the target demagnetization trajectory for the compressed air energy storage motor, also includes:

[0030] Based on the target demagnetization trajectory, the power angle trajectory and reactive power prediction trajectory during the demagnetization process are calculated in advance;

[0031] Verify whether the power angle trajectory and the reactive power prediction trajectory meet the preset power angle constraints and reactive power change rate constraints;

[0032] If any constraint is not satisfied, the gradient projection method is used to iteratively optimize the parameters of the target demagnetization trajectory until all constraints are satisfied.

[0033] The coordinated adjustment is based on the control gain determined through online system identification, which includes:

[0034] Pseudo-random binary sequence disturbances are injected into the excitation systems of compressed air energy storage motors and thermal power units, respectively.

[0035] Simultaneously record their respective reactive power responses;

[0036] Based on the disturbance and response, the transfer function models of the compressed air energy storage motor and the thermal power unit are identified using the recursive least squares method to determine the control gain.

[0037] After identifying the transfer function model, the process also includes:

[0038] Based on the transfer function model, a coupled state-space model is constructed to characterize the dynamic interaction between the compressed air energy storage motor and the thermal power unit.

[0039] By applying the Routh-Hurwitz criterion to analyze the coupled state-space model, the stability region of the control gain that ensures the stability of the closed-loop system is determined, and the control gain is identified.

[0040] Determining the control gain includes:

[0041] A comprehensive performance index is constructed, consisting of a weighted sum of total reactive power deviation, excitation current variation, and adjustment time.

[0042] Under the constraint of the control gain stability region, the particle swarm optimization algorithm is used to optimize the comprehensive performance index in order to find the optimal control gain.

[0043] Beneficial effects: Through the above technical solutions, this invention transforms the hard switching of the CAES unit into an active and controllable soft landing process, eliminating reactive power impact and improving the smoothness of the switching process and the transient stability of the grid voltage. Attached Figure Description

[0044] Figure 1 This is a flowchart of a method for coordinated voltage regulation in a thermal power coupled compressed air energy storage power station.

[0045] Figure 2 This is a flowchart of the target demagnetization trajectory for planning and generating a compressed air energy storage motor.

[0046] Figure 3 This is a flowchart for calculating the transient stability boundary.

[0047] Figure 4 It is a flowchart for integrating three stability boundaries to determine the comprehensive stability boundary. Detailed Implementation

[0048] Example 1: This example provides the technical background and system environment for the application of the present invention, and describes the basic data acquisition steps required to implement the present invention.

[0049] In a specific application scenario, the thermal power coupled compressed air energy storage (CAES) power plant system used in this invention mainly includes a boiler, high-pressure cylinder, low-pressure cylinder, and generator G on the thermal power unit side, and a motor M, compressor, air storage tank, turbine, and generator G on the compressed air energy storage side. The feedwater system of the thermal power unit is deeply coupled with the compression and expansion heat exchange processes of the CAES system to improve energy utilization efficiency. However, when the CAES unit needs to be disconnected from the grid (e.g., when energy storage is completed or during fault switching) or connected to the grid, its characteristics as a high-power inductive / capacitive load / power source will change abruptly, resulting in a step change in the reactive power exchanged to the grid, which in turn causes severe fluctuations in the grid bus voltage or even transient voltage instability problems.

[0050] To address this issue, one approach is to employ advanced control strategies for rapid voltage deviation compensation. For example, an automatic voltage control method based on Model Predictive Control (MPC) establishes a transient model of the FC-CAES system, sets minimum voltage deviation and shortest response time as the combined objectives, and uses MPC to continuously optimize control variables such as the excitation regulator and turbine valve openings. Furthermore, based on the difference in response speed between the thermal power unit and the CAES system, reactive power compensation tasks are dynamically allocated: the fast-responding CAES system handles rapidly fluctuating voltage deviations, while the thermal power unit with a wider adjustment range handles sustained voltage deviations.

[0051] In this embodiment, it can be understood that the aforementioned MPC control method is essentially a responsive or compensatory control strategy after a voltage deviation occurs. Although it can optimize the compensation process, it does not solve the root problem of reactive power step jump caused by the sudden change in the state of the CAES unit. Accordingly, at the moment of switching, the power grid will still inevitably suffer from the initial reactive power impact and voltage drop; only the subsequent recovery process is optimized.

[0052] In light of the above embodiments, and to address the problems existing in the prior art, the applicant conducted in-depth research and discovered that when the CAES compressor motor, which serves as a high-power inductive load, is disconnected using a direct tripping method, the large amount of reactive power it consumes disappears instantaneously within a timescale of nanoseconds to microseconds. This near-vertical dQ / dt step change represents a significant reactive power surge for the power grid, causing a momentary rise in the bus voltage at the connection point. Existing responsive compensation methods, whether AVR or SVG for thermal power units, have response times on the order of tens to hundreds of milliseconds, and cannot trace back and eliminate this initial shock originating from the physical action itself. This forces the system to withstand the shock before compensation, and the risk of voltage transient instability always exists.

[0053] Due to the aforementioned risks of instantaneous impact and loss of synchronization, there is a lack of scientific theoretical basis for planning a switching process that is both fast and relatively safe. To mitigate these risks, engineering practice often employs extremely conservative strategies, such as gradually reducing the motor load over several minutes before disconnection. However, this severely sacrifices the rapid response capability that the CAES system should possess. Conversely, if speed is prioritized, simplified, experience-based control logic must be relied upon, making it impossible to accurately predict the safe switching rate under specific operating conditions. This is because the stability boundary of the switching process is not a fixed value, but a complex multidimensional constraint surface determined by static power angle stability, transient initial swing stability, and dynamic damping characteristics. Current technology lacks the means to model and quantify this comprehensive stability boundary online and accurately. Therefore, it is impossible to plan a control trajectory with optimal dynamic performance and theoretically guaranteed safety for the switching process (e.g., the de-excitation process of the motor), making the switching control method overly conservative or risky. To address this, the following implementation example is provided:

[0054] According to one aspect of this application, a method for coordinated voltage regulation in a thermal power coupled compressed air energy storage power station includes:

[0055] In response to the switching request signal, the initial operating status data of the compressed air energy storage motor is collected. The initial operating status data includes the initial excitation current and the initial reactive power.

[0056] Based on the initial operating status data, a target demagnetization trajectory is planned and generated to actively reduce the initial excitation current, and the total reactive power demand of the system is determined.

[0057] Based on the target demagnetization trajectory, the compressed air energy storage motor is demagnetized, and based on the total reactive power demand of the system and the reactive power collected in real time, the excitation current of the thermal power unit is adjusted in coordination to generate and issue excitation adjustment commands for the compressed air energy storage motor and the thermal power unit.

[0058] The instantaneous excitation current and instantaneous reactive power of the compressed air energy storage motor are monitored. When both meet the preset disconnection conditions, a circuit breaker tripping command is generated and sent to disconnect the compressed air energy storage motor from the grid.

[0059] To address this technical challenge, this invention proposes an active control method based on pre-planning and process coordination. The method begins with switching trigger detection and initial state acquisition. Specifically, it includes the following steps:

[0060] S1.1 Switching Request Verification and Timing Calibration

[0061] When a switching request signal CMD is received from the CAES control system switch Upon receiving the signal, full-state synchronous data acquisition is immediately initiated. To ensure the validity and reliability of the handover request signal, a triple redundancy check logic is preferably used for confirmation. For example, the received original handover request signal CMD... switch_raw The process involves three checks: duration (e.g., requiring a signal duration greater than 100ms), signal amplitude (e.g., conforming to standard TTL levels), and checksum verification (e.g., CRC-16 check). Only after all three checks pass will a confirmation switching command (CMD) be generated. switch_confirmed And mark a switching trigger time t with precision down to the millisecond level or higher. trigger Simultaneously, a high-speed data logger is activated, increasing the sampling rate from the usual 100Hz to 1kHz, providing a time reference for subsequent precise control.

[0062] S1.2, CAES synchronous acquisition of motor electrical quantities

[0063] In t trigger At any given moment, the three-phase instantaneous current i is acquired from the current transformer at the CAES synchronous motor port using synchronous sampling technology. a (t), i b (t), i c (t), and collect the three-phase instantaneous voltage u from the voltage transformer. a (t), u b (t), u c(t). For subsequent calculations based on phasors and flux linkage, the acquired three-phase time-domain variables need to be transformed into a synchronously rotating dq coordinate system. Specifically, a Clarke transformation is performed to convert the current quantities {i} in the ABC three-phase stationary coordinate system. a i b i c} and voltage quantity {u a ,u b ,u c Transforming to the αβ two-phase stationary coordinate system, we obtain the αβ axis component i. α i β u α u β Subsequently, the electrical angle θ is calculated in real time based on the phase-locked loop (PLL) circuit. pll Perform the Park transformation to further rotate the αβ axis components to the dq synchronous rotating coordinate system, and obtain the d-axis current i. d q-axis current i q d-axis voltage u d and q-axis voltage u q The synchronous sampling technique here is based on the same θ. pll The transformation allows all electrical quantities to be snapshotted at the same time and described in the same reference frame, eliminating calculation errors caused by asynchronous sampling or inconsistent phase references.

[0064] S1.3, Precise Calculation of Power and Excitation State

[0065] Based on this, the current power state of the motor is calculated. Preferably, to obtain a more accurate instantaneous response, instantaneous power theory is used for calculation, rather than an average value over a power frequency cycle. Instantaneous active power P e and instantaneous reactive power Q CAES_init The calculation formula is: P e =(3 / 2)*(u d *i d +u q *i q );Q CAES_init =(3 / 2)*(u q *i d -u d *i q Simultaneously, the current excitation current i is read from the excitation controller. f0 The initial power angle δ0 is obtained through a power angle measuring device or by calculation based on voltage and current phasors. For example, one accurate way to calculate the power angle is: δ0 = arctan(u q / (u d +R s *i d )); where R sThe value is the pre-determined stator resistance.

[0066] The accurate calculation of the power angle employs a modified two-step method. The internal potential phasor, E, is calculated using voltage and current phasors. vector =U vector +jX q ·I vector ;where U vector =u d +j·u q I is the terminal voltage phasor. vector =i d +j·i q X is the stator current phasor. q This is the q-axis synchronous reactance.

[0067] The arctangent function in the four quadrants is used to calculate the work angle: δ0 = atan2(E q E d ); where atan2 is the arctangent function that can correctly handle all quadrants, E d and E q These are the d-axis and q-axis components of the internal potential, respectively. This method avoids the problem of the denominator being zero, in u d +R s ·i d It can still give the correct result when the value is close to zero.

[0068] In the actual implementation, a numerical stability check is also added: if |u d +R s ·i d If |<ε (where ε is a preset small positive number, such as 0.001pu), then switch to the power angle calculation method based on power balance: δ0=asin(P e ·X d / (V·E0)), where P e V represents the active power, V represents the bus voltage amplitude, and E0 represents the internal electromotive force amplitude. This dual protection mechanism ensures the robustness of power angle calculation under various operating conditions.

[0069] S1.4 Acquisition of Coordinated Status of Thermal Power Units

[0070] At the same time, the initial reactive power output Q of the thermal power unit is acquired from the DCS system of the thermal power unit via an industrial bus protocol (such as Modbus protocol). thermal_init Adding the initial reactive power of CAES to the initial reactive power of thermal power, we obtain the total reactive power baseline or target value Q that the system needs to maintain at the moment of switching. target =Q CAES_init +Q thermal_init .

[0071] Example 2: In an exemplary embodiment, a method for coordinated voltage regulation of a thermal power coupled compressed air energy storage power station is provided, such as... Figure 1 As shown, it includes:

[0072] Step 1: Respond to the switching request and plan and generate the target demagnetization trajectory of the compressed air energy storage motor.

[0073] In this embodiment, the target demagnetization trajectory is defined as a pre-calculated target excitation current function with time as the independent variable, denoted as i. f *(t). This function describes the state from the switching trigger time t. trigger Initially, the excitation current of the CAES synchronous motor should decrease smoothly and controllably over time until it reaches a preset minimum value. Planning and generating this trajectory is a prerequisite for achieving active control and smooth transition in this invention. Specifically, the rapid, step-like changes in the reactive power of the CAES motor are transformed into a controlled, slope-constrained gradual process lasting from hundreds of milliseconds to several seconds.

[0074] In other words, by actively controlling the excitation current—an internal factor—precise management of the reactive power—an external characteristic—is achieved, thereby avoiding voltage surges. Specific planning methods will be detailed in subsequent embodiments, and typically consider various system stability constraints to ensure that the demagnetization process itself does not cause motor step loss or system oscillation.

[0075] Step 2: Execute the target demagnetization trajectory and coordinate the adjustment of the excitation current of the thermal power unit to compensate for the reactive power changes during the demagnetization process of the compressed air energy storage motor.

[0076] On the one hand, the excitation system of the CAES motor, as a high-precision actuator, needs to strictly track and execute the target demagnetization trajectory i generated in step 1. f *(t). That is, at any time t, the actual excitation current i f_CAES_actual (t) should be as close as possible to the target value i. f *(t). On the other hand, a two-way coupled coordinated control mechanism is initiated between the thermal power unit and the CAES generator. This coordinated regulation refers to establishing a closed-loop feedback control system that monitors in real time the total reactive power Q provided by the CAES and thermal power units in the power grid. total (k)=Q CAES (k)+Q thermal (k), and compare it with the total reactive power demand Q of the system determined before the switch. target By comparing the results, a reactive power deviation ΔQ(k) = Q is formed. total (k)-Q target Based on this deviation, the cooperative control law will simultaneously calculate the regulation rate di of the CAES excitation current. f_CAES / dt and the regulation rate di of the excitation current of the thermal power unit f_thermal / dt. The principle is that when CAES reduces its reactive power output Q due to executing a demagnetization trajectory... CAES When this occurs, ΔQ(k) becomes negative. At this time, the control law will instruct the excitation system of the thermal power unit to increase excitation in order to improve its reactive power output Q. thermal This brings ΔQ(k) back to near zero. In this way, the increased reactive power generated by the thermal power units precisely compensates for the decreased reactive power generated by CAES, ensuring that the sum of the two remains dynamically stable at the target value Q. target Through real-time dynamic compensation, the total reactive power exchange remains stable from the perspective of the power grid bus, thus achieving voltage stability.

[0077] Step 3: When the reactive power output of the compressed air energy storage motor meets the preset disconnection conditions, disconnect it from the grid.

[0078] As step 2 continues, the excitation current i of the CAES motor... f_CAES And its reactive power output Q CAES It will descend continuously along a predetermined trajectory. This step aims to precisely determine the optimal time to physically disconnect the CAES motor from the power grid. The preset disconnection condition is usually a multi-dimensional logical judgment condition designed to minimize the impact at the moment of disconnection.

[0079] Specifically, the condition includes at least the following:

[0080] (1) Excitation current i of CAES motor f_CAES_actual It has dropped to a preset minimum excitation threshold i f_min The following threshold is typically set to the rated excitation current i. f_rated The minimum value, for example, 5% (i f_min =0.05*i f_rated Maintaining a minimum excitation level is to preserve basic controllability of the motor before disconnecting it from the grid, thus preventing it from entering an asynchronous operating state.

[0081] (2) Residual reactive power output Q of CAES motor residual The absolute value is less than a preset reactive power threshold, such as the rated reactive power Q. rated 2% (|Q) residual |<0.02*Q rated This is the criterion for achieving a reactive soft landing.

[0082] Optionally, to further improve the robustness of the judgment, a stability criterion can be added, for example, requiring the standard deviation σ of the residual reactive power fluctuation within the most recent time window (e.g., 100ms). Q_residual The value should be less than a minimum (e.g., 0.01 pu) to confirm that the reactive power output has entered a steady state.

[0083] When the above conditions are met simultaneously, the system generates an offline enable signal ENABLE. disconnect Preferably, at the instant the motor phase current naturally crosses zero, a tripping command (CMD) is sent to the output circuit breaker. open Disconnecting the current at zero point minimizes the generation of electric arcs, further reducing electrical shock. Through this step, the CAES motor is smoothly and safely disconnected with almost no reactive power exchanged with the grid, thus completing the entire coordinated switching process.

[0084] Example 3: This example describes how to plan and generate the target demagnetization trajectory of the compressed air energy storage motor after responding to a switching request, such as... Figure 2 As shown, it specifically includes:

[0085] Step 2.1: Collect initial operating state data including initial excitation current, and obtain active power and initial power angle from it.

[0086] It can provide accurate initial conditions for subsequent calculations, including the initial excitation current i. f0 Active power P e Initial work angle δ0, d-axis / q-axis current i d / i q A series of state variables.

[0087] Step 2.2: Calculate the initial internal electromotive force (EMF) based on the initial excitation current. This yields the initial internal EMF E0, which accurately reflects the electromagnetic state inside the motor. In traditional simplified models, a linear relationship E0 = K is typically used. f *i f0 Calculations are performed, but in actual operation, due to nonlinear effects such as magnetic circuit saturation, this linear relationship will produce large errors. To overcome this defect, this embodiment preferably employs a refined calculation method that considers multiple nonlinear effects.

[0088] Specifically, based on the preset magnetic circuit saturation characteristics, a nonlinear mapping relationship between the initial excitation current and the excitation flux is established.

[0089] Magnetic circuit saturation characteristics refer to a curve describing the nonlinear relationship between the motor's excitation current and the generated magnetic flux linkage. This curve is obtained in advance through no-load tests on the motor and stored in the controller's database.

[0090] In practical implementation, multiple sets of data points {i} obtained from the experiment can be used. f_test (i),Φ gap (i)}, through numerical fitting methods, a continuous nonlinear mapping function f is generated. sat For example, in a linear segment (such as i) f<0.8*i f_rated Using a linear function Ψ f =L ad0 *i f In the saturation region, a modified Frolich equation is used for fitting, which takes the form: Ψ f =(a*i f ) / (b+i f ); where Ψ f The air gap flux linkage generated by the excitation winding, measured in Weber (Wb); f L is the excitation current, measured in amperes (A). ad0 The d-axis mutual inductance is in the unsaturated state, expressed in Henry (H); a and b are fitting coefficients identified using optimization algorithms such as least squares. This is achieved through the function f. sat This can be determined based on the initial excitation current i f0 The excitation flux linkage Ψ, taking into account the basic saturation effect, was calculated. f0 =f sat (i f0 ).

[0091] Furthermore, taking into account the cross saturation effect and the dynamic influence of real-time temperature on motor parameters, the excitation flux obtained from the nonlinear mapping relationship is corrected.

[0092] The cross-saturation effect refers to the phenomenon where the magnetic fields of the d-axis and q-axis magnetic circuits influence each other; that is, the current on one axis affects the permeability on the other axis. To account for this effect, an equivalent magnetomotive force F can be introduced. eq The concept of F is calculated using the formula: F eq =sqrt((i f0 +i d ) 2 +(k cross *i q ) 2 ); where i d and i q The initial d-axis and q-axis currents obtained from step 2.1; k cross The cross-saturation coefficient is a dimensionless parameter with typical values ​​ranging from 0.6 to 0.8, reflecting the contribution of the q-axis magnetomotive force to the saturation of the d-axis magnetic circuit. At this point, the total effective flux linkage should be determined by F... eq Through the saturation curve f sat Obtain, i.e., Ψ total =f sat (F eq The flux linkage component along the d-axis is distributed proportionally to the magnetomotive force.

[0093] The dynamic impact of real-time temperature on motor parameters refers to how the temperature of the motor windings and core changes their resistance and permeability, thereby affecting electromagnetic relationships. Specifically, the excitation winding temperature T can be read in real time from a preset temperature sensor. f and stator winding temperature T s Based on the temperature value, relevant motor parameters are adjusted. For example, the correction formula for the excitation winding resistance is: R f_corrected =R f_20℃ *(1+0.004*(T f -20)); where R f_20℃ This is the reference resistance value at 20℃, and 0.004 is the temperature coefficient of resistance of copper.

[0094] Similarly, the permeability of the iron core also varies with temperature, and its effect can be reflected as a correction factor k for the magnetic flux linkage. temp =1+c temp *(T s -20); where c temp This is the temperature effect coefficient of magnetic permeability.

[0095] Furthermore, the initial internal potential is derived from the corrected excitation flux linkage. Combining the above corrections, the final and accurate total d-axis flux linkage Ψ is obtained. d_final The initial internal potential E0 is then calculated from the precise magnetic flux and the synchronous angular velocity ω (for a 50Hz system, ω = 2π * 50 rad / s): E0 = ω * Ψ d_final Compared to the result calculated by the linear model, this E0 value more realistically reflects the internal operating conditions of the motor under specific loads and temperatures, providing a high-precision input for subsequent stability calculations.

[0096] Step 2.3: Based on the active power, initial power angle and initial internal potential, calculate the static stability boundary, transient stability boundary and dynamic stability boundary respectively.

[0097] This step is used to determine the upper limit of the demagnetization rate. The demagnetization rate is not a fixed empirical value, but must simultaneously satisfy stability constraints at three different time scales, as follows:

[0098] Calculate the static stability boundary B static This boundary condition describes the limiting condition that the rate of change of excitation current must satisfy to maintain synchronization between the motor and the power grid during an extremely slow demagnetization process. It is primarily related to the system's power angle stability margin. Its calculation formula is: B static =|di f / dt| static =(P e *X d ) / (3*V*E0*sin(δ0)); where, P e X represents the initial active power (W);d Here, d is the synchronous reactance (Ω); V is the voltage amplitude of the grid bus (V); E0 is the precise initial internal potential (V) calculated in step 2.2; and δ0 is the initial power angle (rad). The physical meaning of this formula is that the rate of decrease in internal potential caused by demagnetization cannot exceed the rate at which the system compensates for its power transmission capacity by increasing the power angle; otherwise, static instability will result.

[0099] Calculate the transient stability boundary B transient This boundary focuses on the synchronous stability of the motor rotor during the first swing or the first few swing cycles in a relatively rapid demagnetization process, i.e., preventing the motor from losing synchronism due to power imbalance during transient processes. In this embodiment, the evaluation method is as follows: based on the preset critical power angle and initial power angle, the power angle stability margin is determined; combining the power angle stability margin, the preset unit inertia constant, and the transient time constant, the transient stability boundary is calculated, such as... Figure 3 As shown. Specifically, the calculation formula is:

[0100] B transient =|di f / dt| transient =(2*H*ω s *(δ critical -δ0)) / (T d0 '*cos(δ0)); where H is the inertial time constant (s) of the unit, reflecting the magnitude of the unit's moment of inertia; ω s δ is the synchronous angular velocity (rad / s); critical The transient stability critical power angle (rad) of the system is related to the network structure and is usually taken as a conservative preset value, such as 75 degrees (approximately 1.309 rad); δ0 is the initial power angle (rad); T d0 ' is the d-axis transient open-circuit time constant (s), reflecting the rate of change of the excitation winding flux linkage. This formula is derived based on the equal area criterion in transient stability theory, which ensures that the acceleration energy accumulated by the rotor due to power imbalance during demagnetization can be absorbed by the subsequent deceleration process, preventing the power angle from exceeding the critical value and causing loss of synchronism.

[0101] Calculate the dynamic stability boundary B dynamic This boundary primarily focuses on whether power oscillations can be effectively damped and eventually attenuated after the system is disturbed, i.e., small-signal stability. It ensures that the demagnetization process does not excite or exacerbate the inherent low-frequency oscillation modes in the system. Its calculation formula is: B dynamic =|di f / dt| dynamic =(ξ*ω n *E0) / K fWhere ξ is the desired minimum damping ratio of the system, a dimensionless parameter, which is usually set to be greater than 0.3 to obtain good dynamic performance, and preferably set to 0.707 to obtain the best response characteristics; ω n The natural oscillation frequency (rad / s) of the system's electromechanical oscillation mode is related to the system's equivalent reactance and moment of inertia; E0 is the initial internal potential (V); K f This is the equivalent gain coefficient of the excitation system. This boundary ensures that the system has sufficient positive damping for even small disturbances during demagnetization, enabling it to quickly suppress power oscillations and avoid dynamic instability.

[0102] Through the above calculations, three upper limits for the demagnetization rate, B, representing different physical constraints, were obtained. static B transient and B dynamic This provides a comprehensive theoretical basis for subsequently determining the final, safe demagnetization trajectory.

[0103] Example 4, in another specific embodiment, describes how to further fuse the three stability boundaries—static, transient, and dynamic—to determine the comprehensive stability boundary after calculating them separately. The fusion process is a prudent decision-making process that integrates multiple factors and considers both safety redundancy and dynamic adaptability. Its final output is the comprehensive stability boundary k. max This will serve as the highest rate constraint for subsequent demagnetization trajectory planning. The fusion decision-making process can be decomposed into two progressive stages: weighted fusion and adaptive adjustment.

[0104] Preliminary boundary fusion is performed: Preset confidence weights are assigned to the static, transient, and dynamic stable boundaries, and then weighted and fused to obtain a weighted average boundary; the minimum value among the three stable boundaries is selected as the minimum boundary; based on the weighted average boundary and the minimum boundary, the comprehensive stable boundary is determined. Figure 4 As shown.

[0105] Through the above steps, the three boundaries B static B transient B dynamic The stability limits of the system are characterized from different physical dimensions, and their importance may vary under different operating conditions. For example, in a switching scenario where a transient process dominates the contradiction, B... transient The constraints should be given greater attention. Therefore, a preferred implementation is to introduce a weighted fusion mechanism.

[0106] Specifically, for the statically stable boundary B static Transient stability boundary B transient With dynamic stable boundary B dynamic Each is assigned a preset credibility weight wstatic w transient with w dynamic These weights are all dimensionless parameters and satisfy w static +w transient +w dynamic =1. The value of the weight reflects the designer's emphasis on different stability problems. For example, w can be set to... static =0.3, w transient =0.5, w dynamic =0.2. The transient stability boundary is assigned the highest weight here because active demagnetization is a drastic transient process, and preventing initial swing loss is the primary task; static stability is a fundamental constraint, while dynamic stability mainly focuses on subsequent small-signal oscillations, hence its relatively lower weight. Furthermore, the weighted average boundary B is calculated using the following formula. weighted =w static *B static +w transient *B transient +w dynamic *B dynamic .

[0107] At the same time, to ensure the relative safety of the system, the most stringent single constraint condition also needs to be considered. Therefore, the minimum value among the three stability boundaries is chosen as the minimum boundary B. min The calculation method is as follows:

[0108] B min =min{B static B transient B dynamic}. B min This represents the most conservative rate limit that cannot be exceeded under any circumstances.

[0109] Based on the weighted average boundary B weighted With minimum boundary B min Determine the preliminary comprehensive stability boundary k max_raw A conservative fusion strategy that balances average performance and extreme security is: k max_raw =min{B weighted 1.2*B min The significance of this formula lies in using a weighted average as a benchmark, while simultaneously imposing a constraint that does not exceed the strictest limit (B). min A hard cap of 1.2 times is used to prevent the weighted average process from masking an extremely stringent short-board constraint due to the other two boundary values ​​being too large, thus leading to unsafe results.

[0110] In some alternative implementations, the weights can be assigned not as fixed values, but dynamically adjusted based on the initial operating conditions. For example, when the initial power angle δ0 of the system is large and close to the static stability limit, w can be dynamically increased. static The weight of w can be dynamically increased when the system's inertial constant H is small and transient stability is poor. transient The weight.

[0111] After obtaining the initial comprehensive boundary k max_raw Subsequently, this invention further proposes an adaptive adjustment mechanism to achieve a prudent determination of this boundary. Specifically, the steps are as follows: evaluating the performance indicators of historical switching processes and obtaining the reactive power margin of the current thermal power unit; dynamically adjusting the safety factor based on historical performance indicators and reactive power margin; and applying the adjusted safety factor to finally determine the comprehensive stability boundary.

[0112] It is important to note that a fixed safety margin cannot adapt to changes in system operating conditions and the accumulation of historical experience. By introducing an adaptive safety factor based on historical data and current operating conditions, the decision-making regarding the demagnetization rate can be made more intelligent and precise.

[0113] Preferably, it is necessary to evaluate the performance metrics of historical handover processes. The controller maintains a historical database storing key performance data from the most recent N successful handover processes (e.g., N=100). From this database, a comprehensive historical performance metric P can be calculated. history .

[0114] For example, P history P can be calculated using the following formula: history =avg(δ max _ history / δ critical )+2*avg(ΔV history / 0.05);

[0115] Where, δ max_history δ represents the maximum power angle offset during the historical handover process. critical The critical work angle for transient stability; ΔV history P represents the maximum bus voltage deviation during the historical switching process (in per-unit value); avg() represents the average value of the most recent N records. history It is a dimensionless comprehensive index. The smaller the value, the more stable and better the historical switching performance.

[0116] At the same time, it is necessary to obtain the reactive power margin Q of the current thermal power unit. margin This value can be provided in real time by the DCS system of the thermal power unit, and its calculation method is Q. margin =Q thermal_max -Q thermal_init ; where Qthermal_max Q represents the maximum reactive power that a thermal power unit is allowed to generate under its current active power conditions. thermal_init This represents the current reactive power output. Q margin This represents the maximum reactive power compensation capacity that the thermal power units, as backup units, can provide.

[0117] Based on this, the safety factor α is dynamically adjusted according to historical performance indicators and reactive power margin. safe α safe It is a dimensionless coefficient between 0 and 1. Its adjustment logic can be a piecewise function or a set of fuzzy logic rules. For example, the following piecewise function rule can be used: If P history If the value is less than 0.5 (indicating excellent historical performance), then the basic safety factor is set to 0.85; if 0.5 ≤ P history If P < 0.8 (good historical performance), then the basic safety factor is set to 0.70; if P history If the historical performance is ≥0.8 (average or poor), then the basic safety factor is set to 0.55.

[0118] Furthermore, the basic safety factor can be fine-tuned based on the reactive power margin: if Q margin Very plentiful (e.g., Q) margin >0.3*Q rated_thermal If the margin is tight, the calculated safety factor can be increased by 5% (i.e., multiplied by 1.05). Conversely, if the margin is tight, the safety factor can be adjusted downward accordingly.

[0119] Furthermore, by applying the adjusted safety factor, the final comprehensive stability boundary k is determined. max The calculation formula is: k max =α safe *k max_raw .

[0120] It can be seen that, through the above weighted fusion and adaptive adjustment, the final comprehensive stable boundary k is obtained. max It is no longer a rigid, conservative value designed based on worst-case scenarios. It takes into account various physical constraints, incorporates the system's historical operating experience and the support capabilities of current cooperating units, and improves the dynamic performance and efficiency of the demagnetization process while ensuring safety.

[0121] Example 5: This example is based on the comprehensive stability boundary k determined in Example 4. max The target demagnetization trajectory of the compressed air energy storage motor is planned and generated by constructing a piecewise function consisting of the following sequentially connected parts.

[0122] In this embodiment, using a single function (such as a pure exponential function) throughout the entire demagnetization process is not optimal. This is because the control objectives differ at different stages of the demagnetization process: initially, speed is prioritized, i.e., reducing the excitation current as quickly as possible to shorten the total switching time while satisfying stability constraints; in the final stage, stability and accuracy are prioritized, i.e., smoothly approaching the target zero point to avoid overshoot and oscillation. Therefore, a preferred implementation is to use a piecewise function to construct the target demagnetization trajectory i. f *(t) represents the function segment that specifically optimizes the dynamic performance of a particular stage.

[0123] Specifically, the piecewise function is composed of the following parts connected in sequence: a main descent segment that performs the main excitation current reduction at a basically constant rate. In this invention, this segment is also referred to as the rapid demagnetization segment or the linear descent segment, the purpose of which is to quickly complete most of the excitation current reduction task at a constant slope close to the maximum safe rate in the early stage of demagnetization.

[0124] The transition section is used to mitigate the rate change. This section is also known as the smooth transition section. Its function is to establish a buffer between the high-speed main descent section and the low-speed approximation section, so as to avoid the excitation control system from being impacted or oscillating due to sudden rate changes (i.e., excessive acceleration), and to make the first derivative of the trajectory (i.e., the rate) smooth and continuous.

[0125] The final approximation stage is used to ensure that the rate of change of the excitation current approaches zero when it reaches the target value. This stage is also known as the precise approximation stage. It enables a smooth soft landing at the final minimum excitation current value at the end of the demagnetization process, so that the rate of change when reaching the target point is exactly zero, thereby eliminating overshoot or steady-state error.

[0126] To further clarify, this embodiment explains the specific mathematical definitions of the above-mentioned piecewise functions:

[0127] The main decreasing segment is composed of the linear function i f (t)=i f0 -k fast *t defines the boundary, where i f (t) represents the target excitation current at time t, i f0 k is the initial excitation current obtained from the initial operating state data. fast The demagnetization rate is determined based on the comprehensive stability boundary.

[0128] Specifically, k fast The value of is directly derived from the comprehensive stability boundary k calculated in Example 4. max To allow for a certain safety margin, it can be set to k. fast =α margin *k max ;where α marginA safety factor less than 1, such as 0.9, is used. The duration of this segment is determined by its target reduction rate. For example, the target of the main reduction segment can be set to reduce the excitation current to 30% of its initial value. Then, the time point t1 of this segment can be calculated by the following formula: t1=(i f0 -0.3*i f0 ) / k fast .

[0129] The transition segment is composed of the exponential function i f (t)=i f1 *exp(-(t-t1) / τ smooth )+i f_target Define, where i f1 t1 and t2 represent the excitation current value and time point at the end of the main drop phase, respectively, and τ smooth i is the smoothing time constant. f_target The target excitation current is the transition current.

[0130] In this function, i f1 The value is equal to i f0 -k fast *t1. Smoothing time constant τ smooth This is the key parameter that determines the decay rate of this segment, and its value needs to be related to the rate k of the main decay segment. fast The initial velocity is matched with that of the next segment (approximation segment) to ensure the slope continuity at the connection point. f_target This is the asymptotic target value of the exponential function; it is not necessarily the final target excitation current, but rather an intermediate parameter used to construct a suitable decay curve. The end time point t2 of this segment can be set as t2 = t1 + k. τ *τ smooth , where k τ It is a coefficient, for example, taking 2, which means that the transition process is basically completed after twice the time constant.

[0131] In some alternative implementations, the transition segment can also employ other functions that can provide a smooth rate change, such as a cubic polynomial interpolation function. By setting the function values ​​and first derivative values ​​at the two endpoints t1 and t2 to uniquely determine the polynomial coefficients, a smooth transition effect can also be achieved.

[0132] The terminal approximation segment is composed of the cosine function i f (t)=i f_min +(i f2 -i f_min )*cos(π*(t-t2) / (2*(t final -t2))) delimits, where i f2 t2 and t2 are the excitation current value and time point at the end of the transition section, respectively. f_min t is the preset minimum excitation current.final This is the point in time when the demagnetization process ends.

[0133] The i here f2 This is the function value of the transition segment at time t2. f_min This is the ultimate target value for demagnetization. To maintain controllability of the motor, this value is not zero, but a small positive value, such as 5% of the rated excitation current, i.e., i.e. f_min =0.05*i f_rated .

[0134] t final This is the endpoint of the entire demagnetization trajectory. A key characteristic of this cosine function is that its time variable t changes from t2 to t... final At t=t2, the phase of the cosine function changes from 0 to π / 2. This causes the function value to be i. f2 ; at t=t final When the function value is i f_min Wherein, since the derivative of the cosine function is zero at π / 2, this function causes the excitation current to reach i f_min At that instant, its rate of change di f / dt is exactly equal to zero, achieving a shock-free soft landing.

[0135] By sequentially connecting the three function segments mentioned above, and ensuring the continuity of function values ​​and first derivatives at connection points t1 and t2, a globally smooth (C1 continuous) and performance-optimized target demagnetization trajectory i can be constructed. f *(t).

[0136] Furthermore, to address the actual execution deviations caused by inaccurate model parameters or external disturbances, this invention also proposes a real-time trajectory correction mechanism.

[0137] The process of executing the target demagnetization trajectory also includes: real-time monitoring of the tracking error between the actual excitation current and the target demagnetization trajectory; online estimation of the constant deviation and drift rate of the tracking error model using the recursive least squares method; and generation of bias compensation and slope compensation based on the constant deviation and drift rate to construct a corrected trajectory under the comprehensive stability boundary constraints to guide subsequent control.

[0138] Specifically, the mechanism is implemented as follows: the actual excitation current i is monitored in real time at a relatively high sampling frequency (e.g., every 5ms). factual (t), and calculate its relationship with the current target trajectory i. f The tracking error e between *(t) track (t)=i factual (t)-i f *(t).

[0139] Establish an online model to describe the dynamic characteristics of this error, such as a linear model e. track (t) = a0 + a1*t + n(t); where a0 is a constant deviation or bias; a1 is the linear drift rate of the error; and n(t) is random noise. The Recursive Least Squares (RLS) algorithm is used, based on the continuous e track (t) data, online and in real time, to estimate the values ​​of parameters a0 and a1.

[0140] A compensation signal is generated based on the estimated error model parameters. When the constant deviation |a0| exceeds a preset threshold (e.g., 0.02*i), a compensation signal is generated. f0 When ), a bias compensation amount Δi is generated. f_bias =-a0. When the drift rate |a1| exceeds a preset threshold (e.g., 0.001*k). max When ), a slope compensation quantity Δk is generated. comp =-a1.

[0141] Furthermore, the compensation amount is superimposed on the original target trajectory to construct the real-time corrected trajectory i. f_corrected *(t), its expression is: i f_corrected *(t)=i f *(t)+Δi f_bias +Δk comp *(tt current ); where t current This is the current moment. This corrected trajectory will serve as the actual tracking target for the excitation controller. When generating the corrected trajectory, it is also necessary to verify whether its slope remains within the comprehensive stability boundary k. max Within the constraints, to ensure the safety of the correction operation.

[0142] Through a closed-loop correction mechanism of online estimation and compensation, uncertainties in the execution process can be effectively overcome, enabling the actual demagnetization process to reproduce the safe and efficient trajectory planned in theory.

[0143] According to one aspect of this application, an online trajectory correction mechanism based on real-time tracking error is also provided, which solves the trajectory tracking deviation problem caused by model uncertainty and external disturbances.

[0144] In practical implementation, a dynamic evolution model of the tracking error is established. The tracking error e is defined as follows: track (t)=i f_actual (t)-i f *(t), where i f_actual (t) represents the actual measured value of the excitation current, i f *(t) represents the target trajectory. The dynamic characteristics of this error are described by the following state-space model:

[0145] Error state equation: x e (k+1)=A e ·x e (k)+w(k); Observation equation: e track (k)=C e ·x e (k)+v(k); where x e =[e bias ;e drift [ ] represents the error state vector, containing the constant deviation e bias and drift rate e drift A e =[1,Ts;0,1] is the state transition matrix, Ts is the sampling period; C e =[1,0] is the observation matrix; w(k) and v(k) are the process noise and measurement noise, respectively.

[0146] The recursive least squares (RLS) method is used to estimate the error model parameters online. The recursive formula for the RLS algorithm is: K(k) = P(k-1)·φ(k) / (λ+φ) T (k)·P(k-1)·φ(k)); θ*(k)=θ*(k-1)+K(k)·[e track (k)-φ T (k)·θ*(k-1)]; P(k)=(IK(k)·φ T (k))·P(k-1) / λ; where θ*=[a*0;a*1] is the parameter vector to be estimated, corresponding to constant bias and drift rate; φ(k)=[1;k·Ts] is the regression vector; K(k) is the Kalman gain; P(k) is the estimation error covariance matrix; λ is the forgetting factor, typically taking values ​​of 0.95~0.99.

[0147] A corrected trajectory is generated based on the estimation results. When |a*0|>0.02·i is detected... f0 Or |a*1|>0.001·k max At that time, the corrected trajectory is generated: i f_corrected *(t)=i f *(t)-a*0-a*1·(tt current ); where t current For the current moment. The adjustment is applied gradually to avoid abrupt changes: i f_final *(t)=(1-α(t))·i f *(t)+α(t·i f_corrected *(t); where α(t) is a time-varying weighting coefficient that smoothly transitions from 0 to 1 in approximately 100ms.

[0148] This correction process is constrained by the comprehensive stability boundary. Before generating the corrected trajectory, the corrected demagnetization rate needs to be verified: |di f_corrected / dt|≤0.9·k max If the constraints are exceeded, the correction amount will be subject to saturation limits to prioritize stability.

[0149] Through this online correction mechanism, the system can cope with uncertainties such as excitation system parameter drift and load disturbance, so that the actual demagnetization process always closely follows the theoretical optimal trajectory, and the tracking error can be controlled within ±3%.

[0150] Example 6: Description of generating the initial target demagnetization trajectory i f Following *(t), to ensure the safety and feasibility of this trajectory in practical applications, an offline verification and optimization step is preferably added. Specifically,

[0151] After planning and generating the target demagnetization trajectory of the compressed air energy storage motor, the process also includes: based on the target demagnetization trajectory, pre-calculating the power angle trajectory and reactive power prediction trajectory during the demagnetization process; verifying whether the power angle trajectory and reactive power prediction trajectory meet the preset power angle constraints and reactive power change rate constraints; if any constraint is not met, the gradient projection method is used to iteratively optimize the parameters of the target demagnetization trajectory until all constraints are met.

[0152] The first step is to pre-calculate the power angle trajectory and reactive power prediction trajectory during the demagnetization process based on the target demagnetization trajectory.

[0153] Pre-calculation involves using a high-precision dynamic simulation model to calculate the planned i... f Using *(t) as input, the dynamic response of the generator system is calculated. Specifically, the power angle trajectory δ(t) is solved numerically. The rotor dynamics of the generator are described by the classical swing equation: d(Δω) / dt=(P m -P e (t)) / (2*H); dδ / dt=ω s *Δω; where Δω is the deviation of the rotor angular velocity from the synchronous angular velocity (rad / s); P m ω represents the mechanical power input to the motor (W), which can be assumed to be constant during the short period of demagnetization; H is the inertial time constant of the unit (s); ω s P is the synchronous angular velocity (rad / s); e (t) represents the electromagnetic power (W) of the motor. It should be noted that the electromagnetic power P here... e (t) is a quantity that changes with time because it depends on the internal potential E(t), which in turn is directly controlled by the demagnetization trajectory i being executed. f *(t), that is, E(t)≈K f *if *(t). Therefore, a high-precision numerical integration method, such as the fourth-order Runge-Kutta method, is required to iteratively solve the problem with a very small time step (e.g., 1 ms) to obtain the entire demagnetization time [0, t]. final The complete trajectory of the work angle δ(t) and angular velocity deviation Δω(t) within the range.

[0154] After obtaining the power angle trajectory δ(t), the reactive power prediction trajectory Q of the CAES motor can be further calculated. CAES (t) and the trajectory of collaborative compensation demand for thermal power units Q thermal_required (t). Its calculation formula is: Q CAES (t)=(V 2 / X q )-(V*E(i f *(t)) / X d )*cos(δ(t)); Q thermal_required (t)=Q target -Q CAES (t); where X q The q-axis synchronous reactance is Ω.

[0155] The second step is to check whether the power angle trajectory and the reactive power prediction trajectory meet the preset power angle constraints and reactive power change rate constraints.

[0156] The power angle constraint is the predicted peak power angle δ throughout the entire dynamic process. peak =max(δ(t)) must be kept within the safety limit, i.e., δ peak <δ limit δ here limit The critical work angle δ for transient stability is usually taken as the value. critical The discount value, such as δ limit =0.9*δ critical To retain sufficient stability margin.

[0157] The reactive power change rate constraint refers to the finite reactive power response rate of the thermal power unit, which acts as a co-compensator. Therefore, it is necessary to verify whether the rate of change in compensation demand caused by CAES demagnetization is within the capacity range of the thermal power unit. Specifically, Q is calculated. thermal_required The time derivative of (t) is used to obtain its maximum rate of change k. Q_max =max(|dQ thermal_required / dt|), and check whether it satisfies k Q_max <k Q_limit .

[0158] The limiting rate of change k here Q_limit It can be determined based on the technical parameters of the thermal power unit's excitation system, for example, by taking its rated reactive power Q. thermal_ratedDivide by the principal time constant τ of its excitation system thermal .

[0159] The third step is to use the gradient projection method to iteratively optimize the parameters of the target demagnetization trajectory if any constraint is not satisfied, until all constraints are satisfied.

[0160] If the above tests find that any constraint is violated (e.g., δ), peak ≥δ limit or k Q_max ≥k Q_limit If ), then it indicates the trajectory i of the initial plan. f *(t) is too aggressive and needs adjustment. At this point, an automatic optimization program is started. A comprehensive constraint violation function V is defined. total For example, V total =w angle *max(0,δ peak -δ limit ) 2 +w rate *max(0,k Q_max -k Q_limit ) 2 ;where w angle and w rate V represents the weighting coefficient. The value of this function is proportional to the severity of the constraint violation; V increases if and only if all constraints are satisfied. total It is zero. This will constitute i. f Key parameters of *(t), such as the demagnetization rate k in the main falling section. fast The smoothing time constant τ of the transition section smooth The duration Δt of the approximation segment approach These, along with others, form a parameter vector θ to be optimized.

[0161] Iterative optimization is performed using the gradient projection method. In each iteration, V is numerically calculated using the finite difference method. total The gradient ΔV(θ) is relative to the parameter vector θ. This gradient indicates the direction of parameter adjustment that most quickly reduces the violation. According to θ... new =Proj c (θ old -α*ΔV(θ old The rule updates the parameter vector.

[0162] Where, θ old and θ new These are the parameter vectors before and after the update, respectively; α is the learning rate or step size, which determines the adjustment magnitude in each iteration; Proj c It is a projection operator used to satisfy the updated parameter θ new It will not exceed its own reasonable physical constraints C (e.g., k) fastIt cannot be negative, and it cannot exceed k. max ).

[0163] This cycle of previewing, verifying, and optimizing will continue until V is reached. total The convergence occurs within a sufficiently small tolerance range (e.g., 1e-6), or after reaching a preset maximum number of iterations. The final output is an optimized target demagnetization trajectory i that has undergone complete dynamic process verification and meets safety and feasibility requirements. f_optimal *(t).

[0164] Example 7, in another specific embodiment, describes how to design the key parameters of the coordinated control controller. It is important to note that this approach does not rely on fixed, offline-tuned model parameters, but rather obtains the system's most realistic current dynamic characteristics through online identification, and then performs rigorous stability analysis and controller design based on these characteristics.

[0165] Specifically, coordinated regulation is based on the control gain determined by online system identification. Online system identification includes: injecting pseudo-random binary sequence disturbances into the excitation systems of the compressed air energy storage motor and the thermal power unit respectively; synchronously recording their respective reactive power responses; and using the recursive least squares method to identify the transfer function models of the compressed air energy storage motor and the thermal power unit based on the disturbances and responses, so as to determine the control gain.

[0166] Specifically, the online identification process is performed when the system is in a relatively steady state and the conditions for performing the identification operation are met (e.g., in the preparation phase before a planned task switchover). Pseudo-random binary sequence (PRBS) perturbations are injected into the excitation systems of the compressed air energy storage motor and the thermal power unit, respectively. PRBS signals are deterministic signals with spectral characteristics similar to white noise, but with a defined amplitude and repeatable generation, making them ideal excitation sources for system identification. During the design phase, key parameters need to be determined; for example, a 7th-order M-sequence with a sequence length of 2 can be selected. 7 -1=127; symbol width T bit The selection needs to consider the system's response time, for example, choosing 100ms; the disturbance amplitude A test The value needs to be small enough not to cause significant disturbance to the system, but large enough to ensure a sufficient signal-to-noise ratio for the response signal, for example, 2% of the rated excitation value. The reactive power response of each component is recorded synchronously. While injecting the PRBS disturbance signal, the actual values ​​of the excitation current and reactive power output of both the CAES and the thermal power unit are recorded simultaneously at a high sampling rate (e.g., 100Hz).

[0167] The transfer function models of the compressed air energy storage motor and the thermal power unit were identified using the recursive least squares (RLS) method. The collected input (excitation current disturbance) and output (reactive power response) data sequences were preprocessed (e.g., DC component removal, window function application), and an autoregressive exogenous (ARX) model was established to describe the dynamic relationship between the two, for example, Q(k) + a1*Q(k-1) = b0*i f (kd)+b1*i f (kd-1). The model parameters {a1,b0,b1} can be estimated online using the RLS algorithm. Subsequently, this discrete model can be transformed into a more physically meaningful continuous-domain first-order inertial element transfer function model G(s)=K / (1+s*τ) using methods such as bilinear transformation. Through this step, the real-time gain K on the CAES side can be obtained. CAES and time constant τ CAES And K on the thermal power side thermal and τ thermal .

[0168] After identifying the transfer function model, the process also includes: constructing a coupled state-space model based on the transfer function model to characterize the dynamic interaction between the compressed air energy storage motor and the thermal power unit; applying the Routh-Hurwitz criterion to analyze the coupled state-space model, analyzing the control gain stability domain that ensures the stability of the closed-loop system, and determining the control gain.

[0169] This step forms the theoretical foundation for controller parameter design. Specifically, it involves constructing a four-dimensional state-space model with state vectors...

[0170] X can be chosen as X=[Q] CAES Q thermal i f_CAES i f_thermal ] T Based on the identified transfer function and the designed cooperative control law (di f_CAES / dt=-k1*(Q CAES +Q thermal -Q target ),di f_thermal / dt=k2*(-(Q CAES +Q thermal -Q target From this, we can write the system's state matrix A. The elements of this matrix contain the system's physical parameters {K, τ} and the controller parameters {k1, k2} to be designed.

[0171] The stability of the system is determined by the eigenvalues ​​of the state matrix A. To analytically obtain the stability condition, the closed-loop characteristic polynomial det(λI-A)=0 can be calculated, and the Routh-Hurwitz criterion can be applied. By analyzing the condition that the elements in the first column of the Routh table are greater than zero, the stability inequality that the controller gains k1 and k2 must satisfy can be derived, for example, in the form k1*k2. <C critical ;

[0172] Where C critical It is a system physical parameter {K CAES ,τ CAES ,K thermal ,τ thermal The critical constant is determined by}. The set of all (k1,k2) parameter pairs that satisfy this inequality constitutes the control gain stability region that guarantees the stability of the closed-loop system.

[0173] Furthermore, to evaluate the robustness of the designed stability region to parameter uncertainties,

[0174] After the control gain stability domain is determined, the following steps are also included: applying random perturbations to the system parameters in the coupled state-space model through Monte Carlo simulation; statistically analyzing the distribution of the stability probability and stability margin of the system under the perturbations; and calculating robustness indices based on the distribution of stability probability and stability margin to evaluate the sensitivity of the control gain stability domain to parameter changes.

[0175] Specifically, in the identified nominal values ​​of parameters (such as K) CAES ,τ CAES Near the target parameter, a random perturbation with a specific statistical distribution (e.g., mean 0, standard deviation 10% of the parameter value) is applied, generating thousands of possible system parameter samples. For each parameter sample, the stability boundary C is recalculated. critical Furthermore, for a given pair of controller parameters (k1, k2), we can statistically analyze how many of these thousands of simulations still satisfy the stability condition k1*k2. <C critical This ratio is the stability probability P. stable The stability margin SM can also be calculated as 1 - (k1 * k2) / C. critical The mean and variance of the variance are calculated, and a robustness metric, such as R0, is defined. robust =μ SM -2*σ SM This step allows for the selection of controller parameters that, even with significant uncertainties in system parameters, can still guarantee system stability with a high probability and possess sufficient stability margin, thereby enhancing the robustness of the control system to some extent.

[0176] According to one aspect of this application, in this embodiment, the control gains k1 and k2 are defined as the transfer coefficients from reactive power deviation to excitation regulation rate, where k1 is used for the de-excitation control of the CAES motor and k2 is used for the compensation control of the thermal power unit. Since their control objectives are opposite (one reduces reactive power, the other increases reactive power), in the stability analysis of the closed-loop system, the actual stability condition can also be expressed as: Stability criterion: |k1·K CAES ·k2·K thermal | <C critical ;where K CAES and K thermal These are the static gains from excitation to reactive power for the two subsystems, C. critical =1 / (τ CAES +τ thermal ) 2 , τ CAES and τ thermal Let be their respective time constants. The physical meaning of this criterion is that the total open-loop gain of the cooperative control loop cannot exceed the critical value, otherwise oscillations will occur.

[0177] Furthermore, considering that in actual control, k1 takes a positive value (reducing excitation) and k2 takes a negative value (increasing excitation), the stability region can be more accurately expressed as: {(k1,k2)|k1>0,k2<0,k1·|k2| <C critical / (K CAES ·K thermal )}.

[0178] Example 8: This example, based on the control gain stability region analyzed in Example 7, describes an optimal method for determining and applying the optimal control gain. It is understood that meeting stability requirements is a fundamental requirement for controller design; however, within the stability region, different combinations of control gains (k1, k2) will result in drastically different dynamic response qualities (such as response speed, control accuracy, and stability). This example addresses how to find a set of control gains that optimizes the overall system performance while maintaining stability, and how to enable it to intelligently adapt to real-time operating conditions.

[0179] This includes: constructing a comprehensive performance index consisting of a weighted sum of total reactive power deviation, excitation current variation, and adjustment time; and using a particle swarm optimization algorithm to optimize the comprehensive performance index under the constraint of the control gain stability domain, so as to find the optimal control gain.

[0180] Specifically, the first step is to construct a comprehensive performance index J. This index function J aims to unify multiple, sometimes even conflicting, control objectives (e.g., fast response but smooth regulation) into a single scalar function to facilitate optimization.

[0181] A preferred form of the J function is as follows:

[0182] J=w1*∫[0,T](Q total (t)-Q target ) 2 dt+w2*∫[0,T](di f / dt) 2 dt+w3*t settling ; where the first integral term ∫(Q) total (t)-Q target ) 2 dt is the integral of the square of the total reactive power deviation (ISE), used to measure the accuracy of control. Q total (t) represents the total reactive power provided by CAES and thermal power units at any time t; Q target The total reactive power is the target value. The smaller this term is, the more accurately the total reactive power tracks the target value, and the smaller the fluctuation. The second integral term is ∫(di f / dt) 2 dt is the square integral of the rate of change of the control variable, representing the control energy or regulation stability. f / dt can represent the weighted norm of the CAES and the rate of change of excitation current on both sides of the thermal power plant. The smaller the value of this term, the smoother the adjustment action of the excitation system, which is beneficial to extending equipment life and reducing disturbance to the system. The third term t settling It is an indicator of the system's settling time or response speed. It is defined as the time from the occurrence of the disturbance to the total reactive power deviation |Q|. total (t)-Q target |Initially entered and permanently remained within a very small error band (e.g., ±2%*Q) target The time required for each item is calculated as follows: w1, w2, and w3 are the weight coefficients for each item, all of which are positive real numbers and are usually normalized. Designers can adjust these three weights to express their preferences for the three objectives of accuracy, stability, and speed. For example, in scenarios requiring rapid response, the weight of w3 can be increased. In some more refined implementations, these weights can be systematically determined using scientific decision-making methods such as the Analytic Hierarchy Process (AHP) to reduce subjectivity and trial-and-error costs.

[0183] The second step is to optimize the overall performance index using the Particle Swarm Optimization (PSO) algorithm. PSO is a global stochastic search algorithm that simulates the foraging behavior of bird flocks, and is particularly suitable for solving optimization problems with complex and nonlinear objective functions, such as in this example. Specifically: Search space: The search dimension of the algorithm is 2, that is, the cooperative control gain pair (k1, k2) to be optimized. Fitness function: Each particle represents a candidate (k1, k2) solution. By applying this set of gains in a simplified closed-loop system simulation model, the value of its corresponding overall performance index J is calculated. The fitness function F of the particle is defined as F = 1 / (J + ε); where ε is a small positive number to avoid zero in the denominator. The smaller the value of J, the higher the fitness F. Constraints: The optimization process must be carried out within the control gain stability domain determined in Example 7. This is a hard constraint. Specifically, the position (k1, k2) of each particle must meet the stability margin requirement, for example, SM = 1 - (k1 * k2) / C. critical >0.3, meaning the system must always maintain a stability margin of over 30%. Particles exceeding this constraint will have their fitness set to a minimum (or zero), thus being naturally eliminated during evolution. Following the standard PSO algorithm (initializing the particle swarm, iteratively updating velocity and position, updating individual optima and global optima), after a certain number of iterations, the algorithm finally converges to the global optimal particle position, which is the solved optimal control gain (k1). opt k2 opt ).

[0184] After obtaining the optimal control gain, this invention further proposes an adaptive application strategy to enhance its robustness and performance under varying operating conditions. Specifically, when applying the optimal control gain in coordinated regulation, the strategy further includes: dynamically adjusting the control dead zone based on the standard deviation of the real-time reactive power fluctuation; and segmenting and scaling the optimal control gain according to the amplitude of the real-time reactive power deviation to form the actual control gain applied under different deviations.

[0185] The method for dynamically adjusting the control dead zone is: control dead zone ε q It is no longer a fixed empirical value, but is adjusted based on the real-time cleanliness of the system. The adjustment rule can be expressed as: ε q (t)=ε q 0*(1+β*σ q (t)); where ε q (t) represents the actual control dead zone at the current moment; ε q0 A base dead zone value (e.g., 0.01 pu); σ q (t) represents the standard deviation of total reactive power fluctuation within the most recent time window, reflecting the current noise or disturbance level of the system; β is a positive adaptation coefficient. When the system operates smoothly and the noise is low (σ...q (t) small), the dead zone is automatically narrowed, improving the sensitivity and accuracy of control; when the system has large noise or high-frequency disturbances (σ q (t) large), the dead zone is automatically widened, avoiding unnecessary frequent adjustment actions (i.e., control chattering) caused by the controller's excessive response to noise.

[0186] Segmented scaling of the optimal control gain, also known as variable gain or gain scheduling strategy, aims to match the controller's response strength to the magnitude of the deviation. A preferred segmented scaling rule is as follows: when the absolute value of the total reactive power deviation |ΔQ| is very small (e.g., |ΔQ| < 0.05 pu), it indicates that the system is in a fine-tuning stage, at which point a smaller gain is used, for example, a gain k in practical applications. actual =0.5*k opt This allows for small deviations and weak control, avoiding overshoot and improving steady-state accuracy. When |ΔQ| is within the normal range (e.g., 0.05pu≤|ΔQ|<0.15pu), the optimal gain itself is used, i.e., k actual =k opt To achieve optimal overall performance, when |ΔQ| is very large (e.g., |ΔQ|≥0.15pu), it indicates that the system has suffered a significant impact and requires a rapid response. In this case, a larger gain, such as k, can be used. actual =1.5*k opt This allows for large deviations and strong control, bringing the system state back to the vicinity of the target as quickly as possible.

[0187] It can be seen that through this nonlinear adaptive gain strategy, the controller can intelligently adjust its control strength according to the severity of the situation, just like an experienced operator, thereby achieving a better control effect than a fixed gain controller under various operating conditions.

[0188] Example 9: This example further describes the coordinated adjustment steps. Compared to the feedback adjustment strategy based on real-time deviation in Example 8, this example introduces the concept of model predictive control (MPC), which uses the system model to predict future dynamic responses and pre-plans the optimal control layout, thereby achieving faster and smoother control effects.

[0189] In this embodiment, the coordinated adjustment includes a predictive control loop that continuously provides the system with forward-looking control commands through a rolling optimization mechanism. Specifically, the workflow of the predictive control loop is as follows:

[0190] The first step is to predict the total reactive power trajectory over multiple future time steps based on the system model.

[0191] The controller uses mathematical models that can accurately describe the dynamic behavior of the system to predict the future state of the system.

[0192] The model is preferably a discrete-time state-space model, and its state transition equation is:

[0193] X(k+1)=A d *X(k)+B d *U(k); Y(k)=C d *X(k); where k is the current discrete time; X(k) is the system's state vector, including parameters such as Q. CAES Q thermal i f_CAES i f_thermal Key dynamic variables include: U(k) is the control input vector, i.e., the adjustment command applied to CAES and the excitation system of the thermal power unit; Y(k) is the system output, i.e., the total reactive power Q. total (k); A d B d C d This is the system matrix obtained through system identification (as described in Example 7) and discretization. Based on the state X(k) measured at the current time k, the controller performs recursive calculations using this model to predict the state from time k+1 to the future k+N. p The response trajectory of the total reactive power of the system at time t = {Y(k+1|k), Y(k+2|k), ..., Y(k+N)} p |k)}。 Here, N p Known as the prediction time domain, it defines the length of future time that the controller can see.

[0194] The second step involves using rolling optimization to find an optimal control sequence that minimizes the deviation between the total reactive power trajectory and the system's total reactive power demand. Having obtained future predictions, the controller needs to plan an optimal sequence of future control actions {U(k|k), U(k+1|k), ..., U(k+N)} at the current time k. c -1|k)}, where N c Known as the control time domain (N c ≤N p The optimal definition is achieved through a quadratic cost function J. mpc To quantify, its typical form is: J mpc =Σ[j=1,N p ]||Y(k+j|k)-Q target || 2 Qy +Σ[j=0,N c -1]||ΔU(k+j|k)|| 2 RuThe cost function consists of two parts: the first part is Σ||Y(k+j|k)-Q _target || 2 qy The penalty is applied to the predicted output trajectory and the target value Q over the entire prediction time domain. target The deviation between them. The weight matrix Qy reflects the required tracking accuracy at different times. Part Two Σ||ΔU(k+j|k)|| 2 Ru The penalty applies to the magnitude of the change in control action ΔU within the control time domain. The weight matrix Ru is used to suppress overly drastic control actions and ensure the smoothness of the regulation process. The controller needs to solve for J under various physical constraints (e.g., rate of change limit of excitation current, amplitude limit, etc.). mpc Minimize the optimal control sequence U opt ={U*(k|k),U*(k+1|k),...}. This solution process is repeated in each control cycle, i.e., rolling optimization.

[0195] The third step is to extract the first element of the optimal control sequence to form the feedforward control command. This is done after calculating the complete optimal control sequence U. opt Subsequently, the controller does not execute all of them, but only extracts and executes the first control action U*(k|k) in the sequence. This U*(k|k) is then used as the feedforward control command u at the current moment. predict In the next control cycle k+1, the controller remeasures the actual state of the system and, starting from this new point, repeats the entire process of prediction-optimization-first element extraction. This mechanism enables the control to not only be forward-looking but also to continuously utilize the latest actual measurement information to revise its subsequent plans, thus exhibiting strong robustness to disturbances and model uncertainties.

[0196] To further improve the performance and adaptability of the MPC controller, rolling optimization involves dynamically updating the output constraints, control constraints, and control weight matrices used for rolling optimization based on the different stages of the de-excitation process, the prediction uncertainty of the total reactive power trajectory, and the real-time stability margin of the system. This means that the optimization problem of MPC itself is not static, but rather changes with the circumstances. Specifically, updates are made according to the different stages of the de-excitation process: in the initial main descent stage of de-excitation, speed is more important than accuracy, so the constraint boundary for the predicted output Y(k+j|k) can be appropriately relaxed; while in the final terminal approximation stage, accuracy is the primary objective, so the constraint boundary needs to be tightened.

[0197] Updates based on prediction uncertainty: The controller can estimate its own prediction uncertainty, i.e., prediction covariance, using methods such as Kalman filtering. When the prediction uncertainty is large, to ensure robustness, the output constraints need to be tightened accordingly, leaving more safety margin.

[0198] Update based on the system's real-time stability margin: the stability margin index SM calculated in Example 7 can be used. rel Introduced into the cost function. When SM rel When the system is close to the stability boundary, the element values ​​in the control weight matrix Ru should be dynamically increased to punish drastic control actions and make the control behavior more conservative and gentle.

[0199] To suppress disturbances and errors not covered by the model, the collaborative tuning also includes a feedback correction loop.

[0200] The structure and function of the feedback correction loop are as follows: Based on the deviation between the real-time collected reactive power and the total reactive power demand of the system, a feedback correction command is generated; and the feedback correction command is combined with the feedforward control command to form an excitation regulation command issued to the compressed air energy storage motor and the thermal power unit.

[0201] In practice, a conventional feedback controller (e.g., the PI controller detailed in Example 8) operates in parallel with an MPC controller. This feedback controller is based on the actual measured total reactive power Q. total_actual (k) and objective value Q target The deviation ΔQ between actual (k) generates a feedback correction instruction u feedback The final actual control command u issued to the excitation system final It is a combination of feedforward and feedback instructions, i.e., u final =u predict +u feedback .

[0202] This dual-loop structure of predictive feedforward + feedback correction combines the ability of MPC to actively respond to known dynamics and predictable disturbances with the ability of feedback control to passively eliminate unknown disturbances and model mismatch errors, thus achieving a high degree of unity between the speed, stability and robustness of the control system.

[0203] Example 10 describes a preferred engineering implementation detail for improving the overall performance and reliability of the coordinated voltage regulation method of the present invention, mainly involving a high-fidelity data processing method at the input end of the control system and a smooth control command generation method at the output end. It is understood that advanced control algorithms rely on high-quality input data and output commands that can be accurately executed by physical devices.

[0204] In a specific embodiment, to achieve accurate and delay-free acquisition of state information by the control system, this invention proposes a high-frequency data acquisition and filtering preprocessing method. Specifically, the control system acquires data with a sampling period higher than that of conventional power system monitoring and control (SCADA) systems. For example, it acquires the CAES and the raw reactive power signal Q of the thermal power unit with a sampling period of 5ms (i.e., a sampling rate of 200Hz). CAES_raw and Q thermal_raw Since the raw signal acquired by high-frequency sampling inevitably contains measurement noise and random disturbances, directly using this signal for feedback control can easily cause jitter in the control output. Therefore, this invention preferably uses a Kalman filter to perform online filtering processing on the raw signal.

[0205] The Kalman filter is an optimal linear state estimator that can make an optimal estimate of the system state in a dynamic system with uncertainties, based on a series of incomplete and noisy measurements. In this embodiment, its application process is as follows: A simple state-space model describing the dynamic behavior of the reactive power signal is established; at each sampling time k, the filter performs two steps: prediction and update. Its update equation is:

[0206] Q'(k)=Q'(k|k-1)+K f *[Q raw [(k)-H*Q'(k|k-1)]; where Q'(k) is the optimal estimated value output after filtering at time k; Q'(k|k-1) is the predicted value at time k based on the state at time k-1; Q raw (k) represents the original measurement value at time k; H is the observation matrix; K f Kalman gain. Kalman gain K f It dynamically adjusts based on the prediction error covariance and measurement noise covariance. Compared to traditional moving averages or low-pass filters, Kalman filters can effectively filter out noise while reducing signal phase delay, providing smooth, real-time, and high-fidelity system state feedback for subsequent collaborative control algorithms (such as Examples 8 and 9).

[0207] On the other hand, to ensure that the regulation commands calculated by the control algorithm can be executed safely and smoothly by the excitation system, this invention proposes a multi-level limiting and smooth output method for control commands. The cooperative control algorithm (such as in Embodiment 8 or 9) outputs an ideal excitation current regulation rate di. f / dt or target excitation current i f_cmdSubsequently, the command is not directly sent to the excitation power unit, but first passes through a command conditioning module. This module performs multi-level limiting operations to protect the physical equipment. Rate of change limit: The rate of change of the command must not exceed the maximum dynamic response rate of the excitation system itself, i.e., |di|. f / dt| <k max_dynamic Absolute value limit: The target value of the instruction must not exceed the safe operating range of the excitation current, for example, 0. f <1.2*i f_rated , where 1.2*i f_rated The maximum allowable short-time excitation current. Acceleration limit (optional): To further protect the equipment, the rate of change of the command's rate of change (i.e., acceleration) can also be limited, |d 2 i f / dt 2 | max .

[0208] After clipping, the command sequence may produce non-smooth inflection points at the clipping boundaries. To address this issue, the module further performs smoothing output processing. A preferred approach is to use cubic spline interpolation. Specifically, using a series of clipped discrete command points as interpolation nodes, the controller calculates a set of piecewise cubic polynomials i. f_cmd (t)=a3*t 3 +a2*t 2 The nodes are connected using +a1*t+a0. The properties of cubic spline interpolation ensure that at the nodes, not only are the function values ​​continuous, but their first and second derivatives are also continuous. This continuous curve, composed of a set of smooth polynomials, will serve as the final command signal sent to the excitation controller. In this way, the commands sent to the power electronic devices are continuous and have smooth first derivatives. This reduces the impact on the power electronic devices to some extent, avoids the generation of high-frequency harmonics, and makes the final excitation current response more accurate and smooth.

[0209] Example 11 describes the online monitoring and anomaly protection mechanisms established by the present invention during the demagnetization process to ensure system safety and stability, as well as the precise timing determination logic designed at the end of the process to achieve seamless disconnection from the grid. These mechanisms are key aspects of the present invention's transition from theory to engineering practice, ensuring robustness and achieving the desired final effect.

[0210] Specifically, the present invention executes the target demagnetization trajectory i f Throughout the entire process of *(t), the demagnetization process monitoring and anomaly protection module runs in parallel. The functions of this module include:

[0211] Track tracking error monitoring: Real-time calculation of actual excitation current i f_CAES_actual ​​With target trajectory i f The normalization deviation e of *(t) track =|if_ CAES_actual -i f *(t)| / i f0 If this deviation continues to exceed a preset threshold (e.g., 15%) for a period of time, it indicates that the excitation system may have a fault or slow response, and cannot effectively track the planned trajectory. At this point, the protection mechanism is triggered, and a preferred course of action is to proactively reduce the risk, for example, by adjusting the original demagnetization rate k. demag (or k) fast The demagnetization trajectory is reduced by 20% and replanned based on this corrected rate.

[0212] Power angle change rate monitoring: Real-time calculation of the generator power angle change rate dδ / dt. The power angle change rate is a very sensitive indicator for measuring transient stability. If dδ / dt exceeds an emergency threshold (e.g., 10° / s), this usually indicates that the system is rapidly sliding towards the loss-of-synchronization boundary. At this time, an emergency stability control procedure should be initiated immediately, for example, by adjusting the current loss-of-excitation rate k. demag Immediately halve or even zero the excitation, thus halting the demagnetization process and prioritizing system synchronization. Simultaneously, reduce the maximum power angle offset δ during this event. max And the recovery time t recover Key metrics are recorded for subsequent performance evaluation and model calibration.

[0213] Through the above-mentioned online monitoring and protection, the present invention adds a safety lock to the actual execution process of the theoretically safe and feasible demagnetization planning (such as in Example 6), so as to maintain the system stability to the greatest extent even when faced with unexpected disturbances or system anomalies.

[0214] As the demagnetization process nears its end, this invention activates a zero-excitation grid disconnection timing precision determination module to minimize the disturbance to the power grid caused by the final physical disconnection action. The logic of this module is as follows: when the actual excitation current i of CAES is detected... f_CAES_actual First drop to minimum excitation threshold i f_min (e.g., 0.05*i) f_rated When the following condition is met, the system enters the grid disconnection preparation phase. During this phase, the controller continuously monitors the remaining reactive power Q of the CAES motor at a higher frequency and with greater accuracy. residual .

[0215] To achieve accurate determination, the following three conditions must be met simultaneously: Current amplitude condition: i f_CAES_actual f_min This results in the motor being in a state of deep demagnetization. Reactive power amplitude condition: the absolute value of the remaining reactive power |Q residual ​| Less than a very small threshold, such as 2% of the rated reactive power (|Q) residual |<0.02*Q rated This is the core of achieving zero reactive power disconnection.

[0216] Reactive power stability condition: To prevent erroneous operation during oscillations when reactive power crosses zero, a stability criterion is added. Specifically, the remaining reactive power Q is calculated within the most recent short time window (e.g., 100ms). residual standard deviation of fluctuation σ Q_residual σ is required. Q_residual The value is less than a preset minimum (e.g., 0.01 pu) to confirm that the remaining reactive power has stabilized rather than instantaneously crossing zero.

[0217] The controller will generate the offline enable signal ENABLE only if all three conditions mentioned above are met simultaneously at a certain moment. disconnect In the optional preferred solution, when generating ENABLE... disconnect After receiving the signal, the controller does not immediately issue a tripping command, but waits for an optimal electrical timing. Specifically, it monitors the instantaneous values ​​of the three-phase stator current of the motor, and only issues a tripping command (CMD) to the corresponding circuit breaker the instant it detects that the current in any one phase naturally crosses zero. open This zero-current breaking technology can suppress the electric arc generated when the circuit breaker contacts are disconnected, which not only extends the life of the circuit breaker but also further reduces the electromagnetic transient process generated by the disconnection operation. Through this set of precise timing and condition judgment logic, this invention enables the final disconnection action of the CAES motor to be completed at an electromagnetic static point that is both grid-friendly and equipment-friendly.

[0218] Example 12: This example provides a specific, non-limiting numerical calculation case to illustrate in detail the complete calculation process of the multi-constraint integrated stability boundary in Examples 3 and 4. Those skilled in the art can clearly understand and reproduce the implementation details of the algorithm of this invention based on the steps shown in this example.

[0219] Assuming a specific switching scenario, the following initial operating state parameters and preset parameters of the thermal power coupled compressed air energy storage (CAES) power plant and grid system were obtained using the method in Example 1: Initial operating state data: Active power P output by the CAES motor e =0.8 (per unit, pu). Grid bus voltage V = 1.0 (pu). CAES motor initial power angle δ0 = 30° (i.e., π / 6 rad). Initial internal electromotive force E0 calculated using the refined method in Example 3 = 1.1 (pu). Motor and system parameters: d-axis synchronous reactance X d=1.2(pu). The unit's inertial time constant H = 4.0(s).

[0220] d-axis transient open-circuit time constant T d 0' = 1.0 (s). Synchronous angular velocity ω s =314.159 (rad / s) (corresponding to a 50Hz system). Equivalent gain coefficient K of the excitation system. f =1.5(pu / pu). The natural frequency ω of the electromechanical oscillation. n =5.0 (rad / s). Control and constraint parameters: Transient stability critical work angle δ critical =75° (i.e., 75*π / 180rad≈1.309rad). The desired minimum damping ratio of the system is ξ=0.707.

[0221] Based on the parameters given above, calculate the comprehensive stability boundary k. max The steps are as follows:

[0222] Step 1: Calculate the three independent stability boundaries (corresponding to Example 3) and calculate the static stability boundary B. static According to formula B static =(P e *X d Substitute the values ​​into the equation: ) / (3*V*E0*sin(δ0)).

[0223] B static =(0.8*1.2) / (3*1.0*1.1*sin(30°))=0.96 / (3.3*0.5)=0.96 / 1.65≈0.5818(pu / s); This value indicates that, to maintain static stability, the rate of decrease of the excitation current should theoretically not exceed 0.5818 per second. Calculate the transient stability boundary B. transient Unify the angle unit to radians: δ0 = π / 6 rad ≈ 0.5236 rad; δ critical ≈1.309 rad. According to formula B transient =(2*H*ω s *(δ critical -δ0)) / (T d Substitute the value into 0'*cos(δ0)).

[0224] B transient=(2*4.0*314.159*(1.309-0.5236)) / (1.0*cos(30°))=(8.0*314.159*0.7854) / 0.866=1973.92 / 0.866≈2279.35 (A / s, here we assume the relationship between the excitation current base value and pu value, for simplification, we still use pu / s as the unit); It should be noted that the calculated result here is relatively large, indicating that the transient stability constraint is relatively loose under this specific working condition. This may be due to the large inertia constant H and sufficient power angle margin. To maintain unit consistency and rationality, we re-examine the dimensions of this boundary, which should be consistent with the excitation current change rate. Assuming the excitation current base value is 1000A, then B transient ≈2.279 (pu / s). Calculate the dynamic stability boundary B. dynamic According to formula B dynamic =(ξ*ω n *E0) / K f Substitute the values:

[0225] B dynamic =(0.707*5.0*1.1) / 1.5=3.8885 / 1.5≈2.5923(pu / s);

[0226] Thus, three independent boundary values ​​were obtained: B static ≈0.5818,B transient ≈2.279,B dynamic ≈2.5923 (all units are pu / s).

[0227] The second step is to merge the three boundaries to determine the preliminary integrated boundary (corresponding to Example 4) and calculate the weighted average boundary B. weighted The exemplary weights given in Example 4 are used: w static =0.3,w transient =0.5,w dynamic =0.2. B weighted =0.3*0.5818+0.5*2.279+0.2*2.5923=0.1745+1.1395+0.5185=1.8325(pu / s); Determine the minimum boundary B. min =min{0.5818,2.279,2.5923}=0.5818(pu / s); It can be seen that under this operating condition, static stability is the most important limiting factor for the system. Determine the preliminary synthetic boundary k. max_raw : Adopting a conservative fusion strategy k max_raw =min{B weighted 1.2*B min}=min{1.8325,1.2*0.5818}=min{1.8325,0.6982}=0.6982(pu / s);

[0228] Step 3: Perform adaptive adjustments to determine the final synthesis boundary and evaluate historical performance and current margins:

[0229] Data was retrieved from the historical database to calculate the comprehensive historical performance index P. history =0.9. According to the rules of Example 4, P history A value ≥0.8 falls within the range of average or poor historical performance. Meanwhile, it is assumed that the current reactive power margin Q obtained from the thermal power unit's DCS system... margin The value is too low, failing to meet the conditions for increasing the safety factor. The safety factor α is dynamically adjusted. safe According to P history =0.9, according to the piecewise function rules, the basic safety factor should be 0.55. Due to insufficient reactive power margin, this factor will not be adjusted upwards. Therefore, the final determined dynamic safety factor α is... safe =0.55. Calculate the final integrated stability boundary k. max :k max =α safe *k max_raw k max =0.55*0.6982≈0.384(pu / s).

[0230] Through the above comprehensive multi-level calculations, the final integrated stability boundary k under this specific working condition is obtained. max ≈0.384 pu / s. This value will be used as the demagnetization rate k in the main descending segment when planning the demagnetization trajectory in Example 5. fast The upper limit benchmark (e.g., k) fast 0.9*k can be taken max (≈0.3456 pu / s). This embodiment demonstrates how the present invention transforms complex theoretical models and multi-dimensional constraints into specific control parameters that can guide engineering practice, fully reflecting the advanced nature, rigor, and practicality of the present invention.

[0231] Example 13 describes an important preferred extension of the cooperative voltage regulation method of the present invention, embodying the system's intelligence and self-learning capabilities. This extension is executed after a complete FC-CAES cooperative switching task is completed. It quantitatively evaluates the control performance of this switching process and uses the evaluation results to iteratively optimize the control parameters, thereby enabling the system to perform better in future switching tasks. This constitutes a complete closed-loop adaptive optimization process of execution-evaluation-learning.

[0232] In a specific embodiment, the performance evaluation and parameter self-optimization method after the switchover includes the following steps: First, the collection and calculation of key performance indicators (KPIs) during the switchover process. After the CAES motor successfully disconnects from the grid and the switchover task is completed, the control system retrieves the entire switchover process (from t...) from the historical recorder. trigger to t disconnect The system records high-frequency data. Based on this data, the system automatically calculates a series of predefined KPIs to comprehensively and quantitatively evaluate the effectiveness of the control.

[0233] Preferred KPIs may include: Voltage stability index: maximum voltage deviation ΔV max ΔV max =max(|V bus (t)-V nominal |); where V bus (t) is the time series of the bus voltage, V nominal This is the rated voltage. This indicator directly reflects the effectiveness in suppressing voltage surges in the power grid. Duration of voltage fluctuation t fluctuation This refers to the first deviation of the voltage from the normal range (e.g., ±1% * V). nominal The time it takes for the system to eventually stabilize within this range. This metric measures how quickly the system recovers to stability.

[0234] Reactive power control accuracy index: Total reactive power deviation integral E Q_integral E Q_integral =∫[0,t disconnect ]|Q total (t)-Q target |dt. This index measures the accuracy of total reactive power in tracking the target value throughout the entire coordinated control process; the smaller the value, the more precise the coordinated control. Coordinated smoothness index: Smoothness of reactive power transfer from thermal power plants. smooth :S smooth =max(|dQ thermal / dt|) / Q thermal_rated This indicator measures the maximum rate of change in the response process of a thermal power unit when undertaking reactive power compensation tasks, reflecting the smoothness of the coordination process and avoiding excessive impact on the thermal power unit.

[0235] The second step involves parameter optimization decisions based on the KPI evaluation results. The system compares the calculated KPIs with a set of preset performance benchmarks. For example, the performance benchmark might be: ΔV max_limit =1%,t fluctuation_limit =500ms. If all KPIs in this switch are better than the performance baseline, the current control parameters (such as the cooperative control gains k1, k2) are considered to be performing well and no adjustment is needed. If any KPI fails to meet the baseline requirement (for example, ΔV in this switch), the switch will fail. maxIf the KPI reaches 1.2%, the online parameter optimization program will be activated. This program will fine-tune the core parameters of the collaborative controller based on the current KPI performance.

[0236] A preferred optimization algorithm is gradient descent. Define a performance cost function J. perf It is a weighted sum of the excesses of each KPI over the benchmark. J is calculated through sensitivity analysis of the control system model. perf The gradient ΔJ with respect to the parameters to be optimized (e.g., k1, k2) perf (k i Furthermore, update the parameter k according to the following rules: i_new =k i_old -η*ΔJ perf (k i ); where k i_new and k i_old These are the parameter values ​​before and after the update, respectively; η is a small learning rate, such as 0.01.

[0237] The third step involves maintaining the historical database and iteratively updating the optimal parameter set. To achieve long-term, cross-task learning and evolution, this invention also includes a historical database (DB). switching The system's maintenance mechanism. After each task switch, the system stores the feature parameter set of the current task in this database. This feature parameter set is a vector containing the input, process, and result of the current task, for example, {τ demag ,k1,k2,ΔV max ,t fluctuation E Q_integral ,S smooth To prevent data from growing indefinitely, the database preferably employs a sliding window mechanism, for example, retaining only the most recent 100 switch records. Based on this database, the system periodically, or after each update, uses a weighted average method to update a set of globally optimal control parameters. When calculating the weighted average, each historical record can be assigned a weight, which is related to the performance of that record (e.g., 1 / J). perf The parameters used in a historical handover case are weighted proportionally to the actual performance of the handover. This means that the better-performing historical handover cases are, the greater the weight given to their parameters in calculating the globally optimal parameter set. This set of optimal control parameters, updated iteratively using weighted historical data, will serve as the controller's initial default parameters for the next handover task.

[0238] According to one aspect of this application, a real-time trajectory correction process based on error compensation is described. Specifically, the tracking error between the actual excitation current and the target demagnetization trajectory is monitored in real time.

[0239] Tracking error is defined as the difference between the actual measured excitation current value and the target value of the planned trajectory at the same moment. The sampling period is set to 5ms, and the actual excitation current i is obtained through a high-precision current transformer. f_actual (t), and the target trajectory i stored in the controller memory. f *(t) performs real-time comparisons and calculates the error sequence e. track (k)=if actual (kTs)-i f *(kTs); where k is the discrete-time index and Ts=5ms is the sampling period. The error sequence stores the most recent 200 sampling points in a circular buffer for subsequent model estimation.

[0240] Furthermore, the constant bias and drift rate of the tracking error model are estimated online using the recursive least squares method.

[0241] Establish a parameterized model of the error e track (k) = a0 + a1·k·Ts + n(k); where a0 represents the constant deviation component, in amperes (A) or per unit (pu); a1 represents the linear drift rate of the error, in A / s or pu / s; n(k) is zero-mean white noise. The initialization parameters of the RLS algorithm are set as follows: θ*(0) = [0;0], that is, the initial assumption is no bias; P(0) = 100·I2, where I2 is a 2×2 identity matrix, indicating that the initial estimate has a large uncertainty; the forgetting factor λ = 0.98, which enables the algorithm to track time-varying parameters. The parameter is updated once in each sampling period. When the parameter change is less than 1% in 10 consecutive periods, the estimate is considered to have converged.

[0242] Based on this, bias compensation and slope compensation are generated according to constant deviation and drift rate to construct a corrected trajectory under integrated stability boundary constraints to guide subsequent control.

[0243] The compensation amount is generated following a hierarchical strategy. When |a0|>0.02·i f0 At that time, bias compensation Δi is generated. f_bias =-0.8·a0, where 0.8 is the compensation coefficient to avoid overcompensation; when |a1|>0.001·k max At that time, slope compensation Δk is generated. comp =-0.7·a1. The corrected trajectory is constructed using a smooth transition:

[0244] i f_corrected *(t)=i f *(t)+Δi f_bias ·S(tt detect )+Δk comp ·(tt current )·S(ttdetect ); where S(τ)=(1-exp(-τ / τ) smooth )) is a smoothing function, τ smooth =50ms, t detect The moment requiring correction was detected. Verification is required after correction. |di f_corrected / dt|≤0.9·k max If the conditions are not met, the compensation amount will be reduced proportionally.

[0245] This correction mechanism can effectively compensate for the static error and dynamic drift of the excitation system. Experimental data shows that the trajectory tracking accuracy improved from ±5% to ±2% after the correction was enabled, while maintaining the stability of the system.

[0246] Through the complete closed loop of execution-evaluation-learning described above, the cooperative voltage regulation method of the present invention is no longer a static, one-time tuned control system, but an intelligent control system that can learn from experience, continuously optimize itself, and adapt to changes in power grid operating conditions.

[0247] This invention solves the technical challenge of scientifically balancing speed and safety in switching control strategies by accurately quantifying the comprehensive stability boundary of the switching process online. This is achieved by uniformly modeling and calculating stability constraints at three different time scales: static, transient, and dynamic, thus realizing the determined boundary k. max_raw This comprehensiveness avoids the potential instability risks caused by a single consideration. Furthermore, by introducing P based on historical performance data... history And the current reactive power margin Q of thermal power units margin Adaptive safety factor α safe This approach transforms the determination of the stability boundary from a static, conservative, one-size-fits-all value into a dynamic, intelligent decision that reflects the system's historical performance and current support capabilities. It converts the abstract risk of out-of-sync into a concrete, actionable rate ceiling k. max This provides a solid theoretical basis for subsequent demagnetization trajectory planning, enabling the switching process to be executed efficiently while ensuring safety, thus overcoming the dilemma of traditional methods that rely on experience and cannot balance speed and safety.

[0248] This invention transforms the instantaneous hard disconnection of the CAES unit into a fully controllable reactive power soft landing process lasting hundreds of milliseconds by constructing a mathematically smooth demagnetization trajectory pre-verified by the dynamic process. This solves the problem of transient bus voltage stability caused by instantaneous reactive power impact. One aspect of this is the utilization of a piecewise trajectory i composed of linear, exponential, and cosine functions. f*(t), each segment optimizes the speed of the initial switching phase, the stability of the middle phase, and the accuracy of the final phase. In particular, the terminal approach segment achieves zero change rate of excitation current when reaching the target point, physically realizing a smooth transition. On the other hand, by pre-calculating the power angle and reactive power before trajectory execution, and using the gradient projection method to iteratively optimize the trajectory that violates the constraints, the planned route map i is improved. f_optimal *(t) is theoretically safe and feasible. This technology, which actively creates and manages the transition state, eliminates the step abrupt change in dQ / dt, enabling the grid voltage to remain highly stable throughout the switching process. Its maximum deviation can be reduced from 3-5% in traditional methods to less than 0.5%.

[0249] This invention achieves precise, rapid, and stable reactive power compensation of thermal power units during de-excitation by acquiring the most realistic dynamic characteristics of the system and designing an optimal and adaptive controller based on these characteristics. This ensures the dynamic balance of the system's total reactive power. This effect is achieved by injecting a PRBS signal into the system and using the RLS algorithm to identify the transfer function model G(s) online, overcoming the performance degradation problem caused by traditional controllers relying on offline, fixed model parameters. Furthermore, a stable control gain domain is analytically derived using the Routh-Hurwitz criterion, and a comprehensive performance index J that balances accuracy, stability, and speed is optimized using the Particle Swarm Optimization (PSO) algorithm. The resulting control gain (k1) opt k2 opt Theoretically, it achieves Pareto optimality. By dynamically adjusting the control dead zone and employing a piecewise scaling gain strategy, the controller can intelligently adapt to different disturbance levels and deviation amplitudes. These combined technical features enable the cooperative control system to achieve high performance and strong robustness.

[0250] This invention integrates the forward-looking nature of Model Predictive Control (MPC) with the robustness of traditional feedback control to construct a high-performance cooperative control system capable of anticipating the future and making timely corrections, further improving the response speed and control accuracy of reactive power compensation. The optimal feedforward control command u is pre-calculated through rolling optimization. predict This control method, which anticipates deviations before they occur, shortens the system's response delay. Meanwhile, the parallel feedback correction loop handles unknown disturbances not covered by the model and parameter mismatch errors, generating feedback commands u. feedback By u predict and u feedback The control system synthesizes the signals to form the final control command. This dual-loop structure enables the control system to actively respond to predictable dynamic processes while also passively and robustly eliminating unpredictable errors. Its cooperative control performance (such as the integral E of reactive power tracking error) is excellent. Q_integral This is an improvement over a single feedback control strategy.

[0251] This invention establishes a mathematical relationship based on the total reactive power deviation for bidirectional, real-time, closed-loop regulation of the CAES (Chemical Active Power System) and thermal power unit excitation. This achieves dynamic balance of total reactive power during the active de-excitation transition state, resolving secondary voltage fluctuations caused by compensation lag or incoordination. The de-excitation process of the CAES and the compensation process of the thermal power unit are transformed from two independent, time-sequential open-loop actions into a single, interdependent closed-loop system driven by a common deviation ΔQ. Under this mechanism, the reactive power compensation response of the thermal power unit is no longer a blind action based on a preset delay, but rather its real-time response speed and amplitude precisely and inversely match the rate of decrease of reactive power in the CAES. This dynamic balance of one decreasing and one increasing, with a constant total amount, makes the two units appear from the perspective of the power grid as a virtual unit with a constant total reactive power output. This maintains high bus voltage stability throughout the switching process, avoiding the under-compensation or over-compensation problems that may occur in traditional methods.

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

Claims

1. A method for coordinated voltage regulation in a thermal power coupled compressed air energy storage power station, characterized in that, include: In response to the switching request, the target demagnetization trajectory of the compressed air energy storage motor is planned and generated. The target demagnetization trajectory is executed, and the excitation current of the thermal power unit is adjusted in coordination to compensate for the reactive power changes during the demagnetization process of the compressed air energy storage motor. When the reactive power output of the compressed air energy storage motor meets the preset disconnection conditions, it will be disconnected from the grid. The planned generation of the target demagnetization trajectory for the compressed air energy storage motor includes: Collect initial operating state data including initial excitation current, and obtain active power and initial power angle from it; Calculate the initial internal potential based on the initial excitation current; Based on active power, initial power angle and initial internal potential, three stability boundaries are calculated, including static, transient and dynamic stability boundaries. The three stability boundaries are fused to determine the comprehensive stability boundary, and the target demagnetization trajectory is generated based on the comprehensive stability boundary; The planned target demagnetization trajectory for the compressed air energy storage motor is constructed by creating a piecewise function consisting of the following sequentially connected parts: The main descent stage is set to execute the main excitation reduction at a basically constant rate. A transition section used to mitigate rate changes; And a terminal approximation section used to make the rate of change of the excitation current approach zero when it reaches the target value; The main decreasing segment is composed of the linear function i f (t)=i f0 -k fast *t defines the boundary, where i f (t) represents the target excitation current at time t, i f0 k is the initial excitation current obtained from the initial operating state data. fast The demagnetization rate is determined based on the comprehensive stability boundary. The transition segment is composed of the exponential function i f (t)=i f1 *exp(-(t-t1) / τ smooth )+i f_target Define, where i f1 t1 and t2 represent the excitation current value and time point at the end of the main drop phase, respectively, and τ smooth i is the smoothing time constant. f_target The target excitation current is the transition current; The terminal approximation segment is composed of the cosine function i f (t)=i f_min +(i f2 -i f_min )*cos(π*(t-t2) / (2*(t final -t2))) delimits, where i f2 t2 and t2 are the excitation current value and time point at the end of the transition section, respectively. f_min The preset minimum excitation current, t final This is the point in time when the demagnetization process ends.

2. The method according to claim 1, characterized in that, Calculating the transient stability boundary includes: Determine the power angle stability margin based on the preset critical power angle and initial power angle; By combining the power angle stability margin, the preset unit inertia constant, and the transient time constant, the transient stability boundary is calculated.

3. The method according to claim 1, characterized in that, The three stability boundaries are combined to determine the integrated stability boundary, including: Preset confidence weights are assigned to static, transient, and dynamic stable boundaries, and then weighted and fused to obtain a weighted average boundary. Choose the minimum value among the three stable boundaries as the minimum boundary; The comprehensive stability boundary is determined based on the weighted average boundary and the minimum boundary.

4. The method according to claim 1, characterized in that, After planning the target demagnetization trajectory for generating the compressed air energy storage motor, the following steps are also included: Based on the target demagnetization trajectory, the power angle trajectory and reactive power prediction trajectory during the demagnetization process are calculated in advance; Verify whether the power angle trajectory and the reactive power prediction trajectory meet the preset power angle constraints and reactive power change rate constraints; If any constraint is not satisfied, the gradient projection method is used to iteratively optimize the parameters of the target demagnetization trajectory until all constraints are satisfied.

5. The method according to claim 1, characterized in that, Coordinated adjustment is based on control gains determined through online system identification, which includes: Pseudo-random binary sequence disturbances are injected into the excitation systems of compressed air energy storage motors and thermal power units, respectively. Simultaneously record their respective reactive power responses; Based on the disturbance and response, the transfer function models of the compressed air energy storage motor and the thermal power unit are identified using the recursive least squares method to determine the control gain.

6. The method according to claim 5, characterized in that, After identifying the transfer function model, the following is also included: Based on the transfer function model, a coupled state-space model is constructed to characterize the dynamic interaction between the compressed air energy storage motor and the thermal power unit. By applying the Routh-Hurwitz criterion to analyze the coupled state-space model, the stability region of the control gain that ensures the stability of the closed-loop system is determined, and the control gain is identified.

7. The method according to claim 6, characterized in that, Determining the control gain includes: A comprehensive performance index is constructed, consisting of a weighted sum of total reactive power deviation, excitation current variation, and adjustment time. Under the constraint of the control gain stability region, the particle swarm optimization algorithm is used to optimize the comprehensive performance index in order to find the optimal control gain.

Images