Cooperative voltage regulation method for thermal power coupling compressed air energy storage power station

Abstract

The invention discloses a thermal power coupling compressed air energy storage power station cooperation voltage regulation method comprising the following steps: responding to a switching request, accurately calculating and fusing static, transient and dynamic stable boundaries, determining a comprehensive stable boundary, and adjusting the voltage of the thermal power coupling compressed air energy storage power station; according to the boundary planning, generating a segmented optimization target excitation loss track for guiding the CAES motor excitation current to smoothly descend; in the process that the CAES motor executes the field loss track, the excitation current of the thermal power generating unit is cooperatively adjusted based on an on-line identification system model so as to dynamically compensate reactive power, reduced due to field loss, of the CAES motor in real time, and constant total reactive power output of the system is achieved; and when it is monitored that the reactive power output of the CAES motor meets the preset near-zero off-network condition, the circuit breaker is accurately controlled to act. Hard switching of the CAES unit is converted into an active and controllable soft landing process, reactive impact is eliminated, and the smoothness of the switching process and the transient stability of power grid voltage are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of power system automatic control, and in particular to a collaborative voltage regulation method for a thermal power-coupled compressed air energy storage power station. Background Art

[0002] As the global energy mix shifts toward a high proportion of renewable energy, power systems are increasingly in need of highly flexible regulation resources. As a novel hybrid energy storage technology that combines the stable power support capabilities of traditional thermal power units with the rapid start-up and shutdown and energy time-shifting characteristics of compressed air energy storage systems (CAES), the FC-CAES power plant demonstrates significant potential for application in multiple areas, including grid peak regulation, frequency regulation, and voltage support. By deeply coupling thermal power and energy storage units at the thermodynamic and electrical levels, this technology enhances the operational flexibility and economic efficiency 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, particularly those enabling rapid and smooth switching between multiple operating modes (such as energy storage, release, and independent operation), is of vital research significance and technological value for ensuring the safe and stable operation of power systems with a high proportion of renewable energy.

[0003] Currently, research on the operational control of thermal power-coupled compressed air energy storage (CAES) power plants focuses on overall system energy management, efficiency optimization, and coordinated operation strategies under steady-state or quasi-steady-state conditions. At the voltage control level, existing technical solutions typically rely on the automatic voltage regulation (AVR) logic of traditional synchronous generators, passively responding to and providing feedback correction to generator-side voltage deviations through the excitation system. When it comes to switching CAES units, for example, when energy storage is completed and the compressor motor needs to be disconnected from the grid, the standard operating procedure is typically to directly open the circuit breaker. Grid-side reactive power changes rely on the thermal power unit's own AVR system or other dynamic reactive power compensation devices (such as SVCs and SVGs) for post-compensation. At the coordinated control level, some solutions propose coordination strategies based on fixed sequential logic. After the CAES unit activates, a preset delay passes before instructing the thermal power unit to perform the corresponding compensation operation. In addition, research is also applying advanced control algorithms such as Model Predictive Control (MPC) to optimize and regulate voltage fluctuations during normal system operation. This involves rolling optimization of the excitation commands for thermal power units and CAES units to achieve a comprehensive optimization of voltage deviation and regulation time. These methods collectively constitute the fundamental technical framework currently used to ensure the operation of such coupled power plants.

[0004] However, existing technologies still face significant technical challenges in handling the specific and drastic transient process of CAES unit mode switching. These challenges stem primarily from the passive response nature of their control philosophy, which leads to an inherent contradiction between the system's transient switching performance and process stability. Specifically, these challenges manifest themselves in the following interrelated technical issues: Existing technologies are unable to eliminate the transient reactive power surges and voltage transient stability risks caused by switching operations. Furthermore, due to the lack of precise quantification of the risk of loss of step, control strategies for the switching process struggle to balance speed and safety. Summary of the Invention

[0005] Purpose of the invention: The present invention provides a collaborative voltage regulation method for a thermal power-coupled compressed air energy storage power station, aiming to solve the problem of grid voltage stability caused by mode switching of compressed air energy storage (CAES) units.

[0006] Technical solution: According to one aspect of the present invention, a method for coordinated voltage regulation of a thermal power plant coupled with a 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] Execute the target demagnetization trajectory and coordinately adjust 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;

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

[0010] The target demagnetization trajectory of the compressed air energy storage motor is planned and generated, including:

[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 transient stability boundary includes:

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

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

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

[0019] Assign preset credibility weights to the static, transient and dynamic stability boundaries respectively, and perform weighted fusion to obtain the weighted average boundary;

[0020] The minimum value among the three stable boundaries is selected as the minimum boundary;

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

[0022] Among them, planning and generating the target demagnetization trajectory of the compressed air energy storage motor is to construct a piecewise function consisting of the following parts connected in sequence, including:

[0023] Set the main step-down section to execute the main excitation reduction at a basically constant rate;

[0024] The transition section used to ease the rate change;

[0025] and a terminal approach section for making the rate of change of the excitation current approach zero when the excitation current reaches the target value.

[0026] Among them, the main descending section is composed of the linear function i f (t)=i f0 -k fast *t defines, where i f (t) is the target excitation current corresponding to time t, i f0 is the initial excitation current obtained from the initial operating state data, k fast is the demagnetization rate determined based on the comprehensive stability boundary;

[0027] The transition section is composed of the exponential function i f (t)=i f1 *exp(-(t-t1) / τ smooth )+i f_target Define, where i f1 and t1 are the excitation current value and time point at the end of the main drop section, τ smooth is the smoothing time constant, i f_target is the transition target excitation 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))) defined, where i f2 t and t2 are the excitation current value and time point at the end of the transition period, respectively. f_minis the preset minimum excitation current, t final is the end time point of the demagnetization process.

[0029] Among them, after planning and generating the target demagnetization trajectory of the compressed air energy storage motor, it also includes:

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

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

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

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

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

[0035] Synchronously record the respective reactive power responses;

[0036] Based on the disturbance and response, the recursive least squares method is used to identify the transfer function model of the compressed air energy storage motor and the thermal power unit to determine the control gain.

[0037] After the transfer function model is identified, the following steps are also included:

[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] The Routh-Hurwitz criterion is applied to analyze the coupled state space model, and the control gain stability domain that ensures the stability of the closed-loop system is analyzed to determine the control gain.

[0040] Determining the control gain includes:

[0041] Construct a comprehensive performance index consisting of the weighted sum of total reactive power deviation, excitation current variation, and regulation 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 solve the optimal control gain.

[0043] Beneficial effect: Through the above technical solution, the present invention transforms the hard switching of the CAES unit into an active and controllable soft landing process, eliminating reactive impact, improving the smoothness of the switching process and the transient stability of the grid voltage. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 The present invention is a flow chart of a collaborative voltage regulation method for a thermal power plant coupled with a compressed air energy storage power station.

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

[0046] Figure 3 It is the flow chart for calculating transient stability boundary.

[0047] Figure 4 It is a flow chart for fusing three stability boundaries to determine the comprehensive stability boundary. DETAILED DESCRIPTION

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

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

[0050] To address this issue, one approach is to employ advanced control strategies to rapidly compensate for voltage deviations. 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 combined objectives, and uses MPC to continuously optimize control variables such as the excitation regulator and turbine valve opening. Reactive power compensation tasks are dynamically allocated based on the response speed differences between the thermal power units and the CAES system. This means that the fast-responding CAES system handles rapidly fluctuating voltage deviations, while the thermal power units with a wide adjustment range handle persistent voltage deviations.

[0051] In this embodiment, it is understood that the aforementioned MPC control method is essentially a reactive or compensatory control strategy after a voltage deviation occurs. While it can optimize the compensation process, it does not address the root cause of reactive power surges caused by sudden changes in CAES unit states. Accordingly, at the switching moment, the grid will still inevitably experience an initial reactive power surge and voltage sag, though the subsequent recovery process is optimized.

[0052] In combination with the above embodiments, in order to solve the problems existing in the prior art, the applicant conducted in-depth research and found that: when the CAES compressor motor, which is a high-power inductive load, is decoupled by directly disconnecting the switch, the large amount of reactive power it consumes will disappear instantaneously within a time scale of nanoseconds to microseconds. This nearly vertical dQ / dt step change is a large reactive impact on the power grid, resulting in an instantaneous increase in the bus voltage at the connection point. The existing responsive compensation means, whether it is the AVR or SVG of the thermal power unit, have a response time of tens to hundreds of milliseconds, and are unable to trace and eliminate this initial impact caused by the physical action itself. The system must first withstand the impact and then compensate, and the risk of voltage transient instability always exists.

[0053] Due to the above-mentioned transient shock and loss of step risks, there is a lack of scientific theoretical basis for planning a switching process that is both fast and relatively safe. In order to avoid risks, extremely conservative strategies are often adopted in engineering practice, such as spending several minutes to gradually reduce the motor load before decoupling, but this seriously sacrifices the rapid response capability that the CAES system should have. On the contrary, if speed is to be pursued, one can only rely on simplified, experience-based control logic, and it is impossible to accurately predict how fast the switching rate is safe under specific operating conditions. The reason is that the stability boundary of the switching process is not a fixed value, but a complex multi-dimensional constraint surface determined by static power angle stability, transient first swing stability and dynamic damping characteristics. The existing technology lacks the means to accurately model and quantify this comprehensive stability boundary online, so it is impossible to plan a control trajectory with optimal dynamic performance and theoretically guaranteed safety for the switching process (for example, the de-excitation process of the motor), making the switching control method too conservative or risky. To this end, the following embodiments are provided:

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

[0055] In response to the switching request signal, collecting initial operating state data of the compressed air energy storage motor, the initial operating state data including initial excitation current and initial reactive power;

[0056] Based on the initial operating status data, a target deexcitation trajectory is generated for proactively reducing the initial excitation current and determining the total reactive power demand of the system.

[0057] The compressed air energy storage motor is demagnetized according to the target demagnetization trajectory, and the excitation current of the thermal power unit is coordinated and adjusted according to the total reactive power demand of the system and the reactive power collected in real time, so as to generate and issue the excitation adjustment instructions of the compressed air energy storage motor and the thermal power unit;

[0058] Monitor the instantaneous excitation current and instantaneous reactive power of the compressed air energy storage motor. When both meet the preset off-grid conditions, generate and send a circuit breaker opening command to disconnect the compressed air energy storage motor from the grid.

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

[0060] S1.1. Switching request verification and timing calibration

[0061] When the switching request signal CMD is obtained from the CAES control system switch When the full state synchronous data collection is started immediately. The switching request signal here is preferably confirmed by triple redundancy check logic to meet its validity and reliability. For example, the received original switching request signal CMD switch_raw Perform duration check (e.g., signal duration must be greater than 100ms), signal amplitude check (e.g., TTL level conforming to the standard), and checksum verification (e.g., CRC-16 check). Only after all three checks are passed will the confirmation switching command CMD be generated. switch_confirmed And mark a switching trigger moment t with millisecond or higher precision trigger At the same time, the high-speed data recorder is started, and the sampling rate is increased from the conventional 100Hz to 1kHz, providing a time reference for subsequent precise control.

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

[0063] In t trigger At this moment, the three-phase instantaneous current i is collected from the current transformer of the CAES synchronous motor port through 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). In order to perform subsequent calculations based on phasors and flux linkages, the collected three-phase time domain variables need to be converted to a synchronously rotating dq coordinate system. Specifically, Clarke transformation is performed to convert the current quantities {i a ,i b ,i c} and voltage {u a ,u b ,u c}Transformed into the αβ two-phase stationary coordinate system, we get the αβ axis component i α 、i β 、u α 、u β Then, the electrical angle θ is outputted in real time based on the phase-locked loop (PLL) circuit. pll , perform Park transformation, 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 technology here is based on the same θ pll The transformation enables all electrical quantities to be snapshotted at the same moment and described in the same reference frame, eliminating the calculation errors introduced by asynchronous sampling or inconsistent phase reference.

[0064] S1.3. Accurate calculation of power and excitation state

[0065] On this basis, the current power state of the motor is calculated. Preferably, in order to obtain a more accurate instantaneous response, the instantaneous power theory is used for calculation rather than the average value within a power frequency cycle. Instantaneous active power P e and instantaneous reactive power Q CAES_init The calculation formula is: e =(3 / 2)*(u d *i d +u q *i q );Q CAES_init =(3 / 2)*(u q *i d -u d *i q ). At the same time, the current excitation current i is read from the excitation controller f0 The initial power angle δ0 is obtained by a power angle measurement device or by calculation based on voltage and current phasors. For example, a method for accurately calculating the power angle is: δ0=arctan(u q / (u d +R s *i d )); where R sis the pre-measured stator resistance value.

[0066] The improved two-step method is used to accurately calculate the power angle. The internal potential phasor is calculated by the voltage phasor and the current phasor: E vector =U vector +jX q I vector ; where U vector =u d +j·u q is the terminal voltage phasor, I vector =i d +j·i q is the stator current phasor, X q is the q-axis synchronous reactance.

[0067] The power angle is calculated using the four-quadrant inverse tangent function: δ0=atan2(E q ,E d ); where atan2 is the inverse tangent function that can correctly handle all quadrants, E d and E q are the d-axis and q-axis components of the internal potential respectively. This method avoids the problem of the denominator being zero. d +R s ·i d It can still give correct results 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 |<ε (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 is the active power, V is the bus voltage amplitude, and E0 is the internal potential amplitude. This dual protection mechanism ensures the robustness of the power angle calculation under various operating conditions.

[0069] S1.4. Acquisition of the coordinated status of thermal power units

[0070] At the same time, the initial reactive output Q of the thermal power side is collected from the DCS system of the thermal power unit through the industrial bus protocol (such as Modbus protocol). thermal_init Add the initial reactive power of CAES and the initial reactive power of thermal power to obtain the total reactive power benchmark or target value Q that the system needs to maintain at the switching moment. target =Q CAES_init +Q thermal_init .

[0071] Example 2: In an exemplary embodiment, a method for collaborative voltage regulation of a thermal power plant coupled with a compressed air energy storage power station is provided, such as Figure 1 Shown, including:

[0072] Step 1: In response to a switching request, a target demagnetization trajectory of the compressed air energy storage motor is planned and generated.

[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 time from the switching trigger time t trigger Initially, the excitation current of the CAES synchronous motor should be smoothly and controllably reduced over time to a preset minimum value. Planning and generating this trajectory is a prerequisite for the present invention to achieve active control and smooth transitions. Specifically, the rapid, step-like change in the reactive power of the CAES motor is converted into a controlled, slope-constrained gradual change lasting hundreds of milliseconds to several seconds.

[0074] In other words, by actively controlling the internal factor of excitation current, we can precisely manage the external characteristic of reactive power, thereby avoiding voltage surges. The specific planning method will be detailed in subsequent examples. It generally takes into account various stability constraints of the system to ensure that the demagnetization process itself does not cause motor loss of steps or system oscillation.

[0075] Step 2: Execute the target demagnetization trajectory and coordinately adjust the excitation current of the thermal power unit to compensate for the reactive power change 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 to the target value i as possible f *(t). On the other hand, a bidirectional coupling coordinated control mechanism is initiated between the thermal power unit and the CAES motor. The coordinated regulation here refers to the establishment of a closed-loop feedback control system that monitors the total reactive power Q provided by the CAES and thermal power units in the power grid in real time. total (k)=Q CAES (k)+Q thermal (k), and compare it with the total reactive power demand Q determined before switching target Compare and form a reactive deviation ΔQ(k)=Q 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 reactive output Q due to the execution of demagnetization trajectory CAES When ΔQ(k) is negative, the control law will instruct the excitation system of the thermal power unit to increase the excitation to increase its reactive output Q thermal , thereby pulling ΔQ(k) back to near zero. In this way, the reactive power generated by the thermal power unit just makes up for the reactive power reduced by the CAES, so that the sum of the two is always 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 grid bus, achieving voltage stability.

[0077] Step 3: When the reactive output of the compressed air energy storage motor meets the preset off-grid condition, the compressed air energy storage motor is disconnected from the grid.

[0078] As step 2 continues, the excitation current i of the CAES motor f_CAES and its reactive output Q CAES The energy level will continue to decrease along a predetermined trajectory. This step aims to accurately determine the optimal time to physically disconnect the CAES motor from the grid. The preset disconnection condition is typically a multi-dimensional logical judgment condition designed to minimize the impact of the disconnection.

[0079] Specifically, the conditions include at least:

[0080] (1) CAES motor excitation current i f_CAES_actual has dropped to a preset minimum excitation threshold i f_min Below. This threshold is usually set to the rated excitation current i f_rated The minimum value, such as 5% (i f_min =0.05*i f_rated The purpose of retaining a minimum excitation is to maintain basic controllability of the motor before disconnection from the grid and to prevent it from entering an asynchronous operation state.

[0081] (2) The residual reactive output Q of the CAES motor residual The absolute value of the reactive power is less than a preset reactive power threshold, such as the rated reactive power Q rated 2% of (|Q residual |<0.02*Q rated ). This is the criterion for achieving 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 in the most recent time window (such as 100ms) to be Q_residual Less than a minimum value (such as 0.01pu) to confirm that the reactive power output has entered a steady state.

[0083] When the above conditions are met at the same time, the system generates the off-grid enable signal ENABLE disconnect , and preferably, when the motor phase current is detected to cross the zero point naturally, a tripping command CMD is sent to the output circuit breaker open Disconnecting at the zero current point minimizes arcing and further reduces electrical shock. This step allows the CAES motor to be smoothly and safely disconnected while exchanging virtually no reactive power to the grid, completing the 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 the switching request, such as Figure 2 As shown, specifically including:

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

[0086] Able to provide accurate initial conditions for subsequent calculations, including the initial excitation current i f0 , active power P e , initial power angle δ0, d-axis / q-axis current i d / i q And a series of state quantities.

[0087] Step 2.2, calculate the initial internal potential based on the initial excitation current. The initial internal potential E0 that can accurately reflect the internal electromagnetic state of the motor can be obtained. In the traditional simplified model, the linear relationship E0=K f *i f0 However, 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 adopts a refined calculation method that takes into account 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 linkage is established.

[0089] The magnetic circuit saturation characteristic is a curve that describes the nonlinear relationship between the motor's excitation current and the resulting magnetic flux. This curve is obtained in advance through a no-load test of the motor and stored in the controller's database.

[0090] In the specific implementation, multiple groups of data points {i f_test (i),Φ gap (i)}, through the numerical fitting method, 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 the linear function Ψ f =L ad0 *i f In the saturation section, the modified Frolich equation is used for fitting, which is in the form of: f =(a*i f ) / (b+i f ); where Ψ f is the air gap flux generated by the excitation winding, in Weber (Wb); i f is the excitation current, in ampere (A); L ad0 is the d-axis mutual inductance in the unsaturated state, in Henry (H); a and b are fitting coefficients identified by optimization algorithms such as the least squares method. sat , which can be calculated based on the initial excitation current i f0 , calculate the excitation flux Ψ considering the basic saturation effect f0 =f sat (i f0 ).

[0091] Furthermore, the cross-saturation effect and the dynamic influence of real-time temperature on motor parameters are comprehensively considered, and the excitation flux obtained by the nonlinear mapping relationship is corrected.

[0092] The cross-saturation effect refers to the phenomenon that the magnetic fields of the d-axis and q-axis magnetic circuits affect each other, that is, the current on one axis affects the magnetic permeability on the other axis. To account for this effect, the equivalent magnetomotive force F can be introduced. eq The concept of F is calculated as follows: eq =sqrt((i f0 +i d ) 2 +(k cross *i q ) 2 ); where i d and i q are the initial d-axis and q-axis currents obtained from step 2.1; k cross is the cross saturation coefficient, a dimensionless parameter with a typical value range of 0.6 to 0.8, reflecting the contribution of the q-axis magnetic potential to the saturation of the d-axis magnetic circuit. At this time, the total effective magnetic flux should be given by F eq Through the saturation curve f sat Get, that is, Ψ total =f sat (F eq ), and the magnetic flux component in the d-axis direction is obtained by distributing it in proportion to the magnetomotive force.

[0093] The dynamic effect of real-time temperature on motor parameters means that the temperature of the motor winding and core will change its resistance and magnetic permeability, thereby affecting the electromagnetic relationship. 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, the relevant motor parameters are corrected. 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℃ is the reference resistance value at 20°C, and 0.004 is the resistance temperature coefficient of copper.

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

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

[0096] In step 2.3, the static stability boundary, transient stability boundary, and dynamic stability boundary are calculated based on the active power, initial power angle, and initial internal potential.

[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 meet the stability constraints under three different time scales, as follows:

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

[0099] Calculate transient stability boundary B transient This boundary focuses on the synchronous stability of the motor rotor in the first swing or the first few swing cycles during the rapid demagnetization process, that is, to prevent the motor from losing step due to power imbalance in the transient process. In this embodiment, the evaluation method is: based on the preset critical power angle and the initial power angle, the power angle stability margin is determined; combined with 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 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 inertia time constant of the unit (s), reflecting the magnitude of the unit's rotational inertia; ω s is the synchronous angular velocity (rad / s); δ critical is the transient stability critical power angle of the system (rad), which is related to the network structure and is usually set to a conservative preset value, such as 75 degrees (about 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 in the excitation winding flux. This formula is derived based on the equal-area criterion in transient stability theory. This ensures that during demagnetization, the accumulated acceleration energy from rotor power imbalance can be absorbed during the subsequent deceleration process, preventing the power angle from exceeding the critical value and causing a loss of step.

[0101] Calculate the dynamic stability boundary B dynamic This boundary focuses on whether the power oscillation can be effectively damped and eventually decayed after the system is disturbed, that is, small signal stability. It ensures that the demagnetization process will not excite or aggravate the inherent low-frequency oscillation mode in the system. Its calculation formula is: B dynamic =|di f / dt| dynamic =(ξ*ω n *E0) / K f; Wherein, ξ is the minimum damping ratio that the system is expected to achieve, which is a dimensionless parameter. To obtain good dynamic performance, it is usually set to be greater than 0.3, and preferably set to 0.707 to obtain the best response characteristics; ω n is the natural oscillation frequency of the system's electromechanical oscillation mode (rad / s), which is related to the system's equivalent reactance and moment of inertia; E0 is the initial internal potential (V); K f 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, quickly suppressing power oscillations and avoiding dynamic instability.

[0102] Through the above calculations, three upper limits of the demagnetization rate B representing different physical constraints are obtained: static 、B transient and B dynamic , providing a comprehensive theoretical basis for the subsequent determination of the final and safe demagnetization trajectory.

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

[0104] Perform preliminary boundary fusion: assign preset credibility weights to the static stability boundary, transient stability boundary and dynamic stability boundary respectively, and perform weighted fusion to obtain the weighted average boundary; select the minimum value of the three stability boundaries as the minimum boundary; determine the comprehensive stability boundary based on the weighted average boundary and the minimum boundary. Figure 4 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 working conditions. For example, for a switching scenario where a transient process occupies the dominant contradiction, B transient The constraints of should be given higher attention. Therefore, a preferred implementation is to introduce a weighted fusion mechanism.

[0106] Specifically, the static stability boundary B static 、Transient stability boundary B transient and the dynamic stability boundary B dynamic Assign the preset credibility weight wstatic 、w transient With w dynamic These weights are dimensionless parameters and satisfy w static +w transient +w dynamic = 1. The value of the weight reflects the importance the designer attaches to different stability issues. For example, w static =0.3, w transient =0.5, w dynamic =0.2. The transient stability boundary is given the highest weight here because active demagnetization itself is a violent transient process, and preventing the first swing out of step is the primary task; static stability is a basic constraint, while dynamic stability focuses on subsequent small signal oscillations, so the weight is relatively low. Further, the weighted average boundary B is calculated by the following formula weighted =w static *B static +w transient *B transient +w dynamic *B dynamic .

[0107] At the same time, in order to ensure the relative safety of the system, the most stringent single constraint condition needs to be considered. Therefore, the minimum value among the three stable boundaries is selected as the minimum boundary B min , which is calculated as follows:

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

[0109] According to the weighted average boundary B weighted With the minimum boundary B min , determine the preliminary comprehensive stability boundary k max_raw A conservative fusion strategy that takes into account both average performance and extreme safety is: k max_raw =min{B weighted ,1.2*B min The significance of this formula is to take the weighted average as the benchmark, but at the same time impose a constraint of not exceeding the most stringent constraint (B min ) A hard cap of 1.2x is used to prevent the weighted average process from masking an extremely stringent short board constraint due to the excessive size of the other two boundary values, thus leading to unsafe results.

[0110] In some other optional implementations, the weighting may not be fixed, but may be dynamically adjusted according to the initial operating conditions. For example, when the system initial power angle δ0 is large and close to the static stability limit, w may be dynamically increased. static The weight of w can be dynamically increased when the system inertia constant H is small and the transient stability is poor. transient The weight of .

[0111] After obtaining the initial integrated boundary k max_raw The present invention further proposes an adaptive adjustment mechanism to achieve prudent determination of this boundary. Specifically, this step involves evaluating historical switching performance indicators and obtaining the reactive power margin of the current thermal power unit; dynamically adjusting the safety factor based on the historical performance indicators and reactive power margin; and applying the adjusted safety factor to ultimately determine the comprehensive stability boundary.

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

[0113] Preferably, the performance index of the historical switching process needs to be evaluated. The controller maintains a historical database that stores the key performance data of the most recent N (e.g., N=100) successful switching processes. From this, a comprehensive historical performance index P can be calculated. history .

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

[0115] Among them, δ max_history is the maximum power angle offset in the historical switching process; δ critical is the transient stability critical power angle; ΔV history is the maximum bus voltage deviation during the historical switching process (in per unit); avg() represents the average value of the most recent N records. history It is a dimensionless comprehensive indicator. The smaller its value, the smoother 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 is calculated as Q margin =Q thermal_max -Q thermal_init ; where Qthermal_max Q is the maximum reactive power allowed to be generated by the thermal power unit under the current active power, thermal_init is the current reactive power output. margin It represents the maximum reactive power compensation capability that a thermal power unit serving as a reserve can provide.

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

[0118] Furthermore, the basic safety factor can be fine-tuned according to the reactive margin: if Q margin Very sufficient (such as Q margin >0.3*Q rated_thermal ), the safety factor calculated above 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, the adjusted safety factor is applied to finally determine the comprehensive stability boundary k max The calculation formula is: max =α safe *k max_raw .

[0120] It can be seen that through the above-mentioned weighted fusion and adaptive adjustment, the final comprehensive stable boundary k max It is no longer a rigid, conservative value based on worst-case design. It takes into account multiple physical constraints, incorporates the system's historical operating experience and the support capabilities of current collaborative 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 , planning and generating the target demagnetization trajectory of the compressed air energy storage motor, specifically constructing a piecewise function consisting of the following parts connected in sequence.

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

[0123] Specifically, the piecewise function consists of the following sequentially connected parts: A main ramp-down section is set to execute the main excitation current reduction at a substantially constant rate. In the present invention, this section is also referred to as the rapid demagnetization section or the linear ramp-down section. Its purpose is to quickly complete the majority of the excitation current reduction task at a constant slope close to the maximum safe rate in the early stages of demagnetization.

[0124] The transition section is used to ease the rate change. This section is also called the smooth transition section. Its function is to create a buffer between the high-speed main descent section and the low-speed approach section. It prevents the excitation control system from being impacted or oscillated by sudden rate changes (i.e., excessive acceleration), and ensures that the first-order derivative of the trajectory (i.e., rate) is smooth and continuous.

[0125] The terminal approach stage is used to make the rate of change of the excitation current approach zero when it reaches the target value. This stage, also known as the precise approach stage, can smoothly land at the final minimum excitation current value at the end of the demagnetization process, making the rate of change when reaching the target point exactly zero, thereby eliminating overshoot or steady-state error.

[0126] For further clarification, this embodiment describes the specific mathematical definitions of the above-mentioned piecewise functions:

[0127] The main descending segment is composed of the linear function i f (t)=i f0 -k fast *t defines, where i f (t) is the target excitation current corresponding to time t, i f0 is the initial excitation current obtained from the initial operating state data, k fast is the demagnetization rate 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 , in order to leave a certain safety margin, it can be set to k fast =α margin *k max ; where α marginis a safety factor less than 1, such as 0.9. The duration of this section is determined by its target reduction. For example, the target of the main reduction section can be set to reduce the excitation current to 30% of the initial value. Then the time point t1 of this section can be calculated by the following formula: t1=(i f0 -0.3*i f0 ) / k fast .

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

[0130] In this function, i f1 The value is equal to i f0 -k fast *t1. Smoothing time constant τ smooth It is the key parameter that determines the attenuation speed of this section, and its value needs to be consistent with the rate k of the main drop section. fast Match the initial rate of the next segment (approaching segment) to ensure the slope of the connection point is continuous. f_target is the asymptotic target value of the exponential function. It is not necessarily the final target excitation current, but an intermediate parameter used to construct a suitable attenuation 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, 2, which means that the transition process is basically completed after twice the time constant.

[0131] In some optional embodiments, the transition section may also adopt other functions that can provide smooth rate changes, such as a cubic polynomial interpolation function. By setting the function values ​​and first-order derivative values ​​of 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))) defined, where i f2 t and t2 are the excitation current value and time point at the end of the transition period, respectively. f_min is the preset minimum excitation current, tfinal is the end time point of the demagnetization process.

[0133] The i here f2 is the function value of the transition section at time t2. f_min It is the final target value of demagnetization. In order to maintain the controllability of the motor, this value is not zero, but a small positive value, such as 5% of the rated excitation current, i f_min =0.05*i f_rated .

[0134] t final The key feature of the cosine function is that its time variable t changes from t2 to t final When t=t2, the phase of the cosine function changes from 0 to π / 2. This makes the function value i f2 ; at t = t final When the function value is i f_min Since the derivative of the cosine function at π / 2 is zero, the function makes the excitation current reach i f_min At the moment of f / dt is exactly equal to zero, achieving a soft landing without shock.

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

[0136] Furthermore, in order to deal with actual execution deviations caused by inaccurate model parameters or external disturbances, the present 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 generating bias compensation and slope compensation based on the constant deviation and drift rate to construct a corrected trajectory under the comprehensive stability boundary constraint to guide subsequent control.

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

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

[0140] The compensation signal is generated based on the estimated error model parameters. When the constant deviation |a0| exceeds the preset threshold (for example, 0.02*i f0 ), generates an offset compensation amount Δi f_bias =-a0. When the drift rate |a1| exceeds the preset threshold (for example 0.001*k max ), a slope compensation value Δk is generated comp =-a1.

[0141] Furthermore, the compensation amount is added to 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 is the current moment. This corrected trajectory will be used as the actual tracking target of the excitation controller. When generating the corrected trajectory, it is also necessary to check whether its slope is still within the comprehensive stability boundary k max within the constraints of the system to ensure the safety of the correction operation.

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

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

[0144] In the specific implementation, a dynamic evolution model of tracking error is established. The tracking error e is defined as track (t)=i f_actual (t)-i f *(t), where i f_actual (t) is the actual excitation current measurement value, i f *(t) is the target trajectory. The dynamic characteristics of the 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 ] is the error state vector, which contains the constant deviation e bias and drift rate e drift ; A e =[1,Ts;0,1] is the state transfer matrix, Ts is the sampling period; C e =[1,0] is the observation matrix; w(k) and v(k) are 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 of 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 the 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, which is typically 0.95~0.99.

[0147] Generate a corrected trajectory based on the estimation results. f0 or |a*1|>0.001·k max When , the corrected trajectory is generated: i f_corrected *(t)=i f *(t)-a*0-a*1·(tt current );where t current is the current moment. The correction is applied gradually to avoid sudden changes: i f_final *(t)=(1-α(t))·i f *(t)+α(t)·i f_corrected *(t); where α(t) is the time-varying weighting coefficient, which smoothly transitions from 0 to 1 with a transition time of about 100ms.

[0148] The correction process is subject to the constraints of the comprehensive stability boundary. Before generating the correction trajectory, the corrected demagnetization rate needs to be verified: f_corrected * / dt|≤0.9·k max If the constraint is exceeded, the correction amount is saturated and limited, giving priority to ensuring stability.

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

[0150] Example 6: Describing the generation of the initial target demagnetization trajectory i f After *(t), in order to ensure the safety and feasibility of the 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 method also includes: pre-calculating the power angle trajectory and reactive power prediction trajectory during the demagnetization process based on the target demagnetization trajectory; checking 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, using the gradient projection method 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 is to convert the planned i f *(t) is used as input to solve the dynamic response of the generator system. Specifically, the power angle trajectory δ(t) is numerically integrated and solved. The rotor dynamics of the generator is described by the classic 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 is the mechanical power input to the motor (W), which can be assumed to be a constant value during the short period of demagnetization; H is the inertia time constant of the unit (s); ω s is the synchronous angular velocity (rad / s); P e (t) is the electromagnetic power of the motor (W). It should be noted that the electromagnetic power P here is e (t) is a time-varying quantity 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, it is necessary to use a high-precision numerical integration method, such as the fourth-order Runge-Kutta method, to iteratively solve with a very small time step (for example, 1ms) to obtain the entire demagnetization time [0,t final ]The complete trajectory of the power angle δ(t) and angular velocity deviation Δω(t) within .

[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 collaborative compensation demand trajectory Q for thermal power units 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 is the q-axis synchronous reactance Ω.

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

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

[0157] The reactive power change rate constraint means that the reactive power response rate of the thermal power unit as the cooperative compensation party is limited. Therefore, it is necessary to verify whether the compensation demand change rate caused by CAES demagnetization is within the capacity of the thermal power unit. Specifically, calculate Q 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] Here the limiting rate of change k Q_limit It can be determined based on the technical parameters of the excitation system of the thermal power unit, for example, the rated reactive power Q thermal_ratedDivide by the main time constant τ of its excitation system thermal .

[0159] In the third step, if any constraint is not satisfied, the parameters of the target demagnetization trajectory are iteratively optimized using the gradient projection method until all constraints are satisfied.

[0160] If the above test finds that any constraint is violated (for example, δ peak ≥δ limit or k Q_max ≥k Q_limit ), it indicates that the initial planned trajectory i f *(t) is too aggressive and needs to be adjusted. At this point, an automatic optimization program is started. Define a comprehensive constraint violation function V 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 is the weight coefficient. The value of this function is proportional to the severity of the constraint violation. If and only if all constraints are satisfied, V total Is zero. It will constitute i f *(t) key parameters, such as the demagnetization rate k of the main drop section fast , smoothing time constant τ of the transition section smooth , duration of the approach segment Δt approach Etc., constitute a parameter vector θ to be optimized.

[0161] The gradient projection method is used for iterative optimization. In each iteration, V is numerically calculated by the finite difference method. total The gradient ΔV(θ) relative to the parameter vector θ. This gradient indicates the parameter adjustment direction that can reduce the violation the fastest. new =Proj c (θ old -α*ΔV(θ old ))’s rule updates the parameter vector.

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

[0163] This rehearsal-test-optimization cycle will continue until V total Converge to a sufficiently small tolerance range (such as 1e-6), or reach the preset maximum number of iterations. The final output is a target demagnetization trajectory i that has been fully dynamic process verified and meets the requirements of safe and feasible optimization. f_optimal *(t).

[0164] In another specific embodiment, the design of key parameters for a collaborative control controller is described. It is important to note that this approach does not rely on fixed, offline tuned model parameters, but rather uses online identification to obtain the most realistic dynamic characteristics of the system, based on which rigorous stability analysis and controller design are performed.

[0165] Specifically, the coordinated regulation is based on the control gain determined through online system identification. The 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 based on the disturbance and response, using the recursive least squares method to identify the transfer function model of the compressed air energy storage motor and the thermal power unit to determine the control gain.

[0166] Specifically, the online identification process is carried out when the system is in a relatively steady state and has the conditions for performing identification operations (for example, in the preparation stage before the planned switching task). Pseudo-random binary sequence (PRBS) disturbances are injected into the excitation systems of the compressed air energy storage motor and the thermal power unit respectively. The PRBS signal is a deterministic signal with a spectral characteristic similar to white noise, but with a fixed amplitude and repeatable generation. It is an ideal excitation source for system identification. During design, its 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; code element width T bit The choice of needs to consider the response time of the system, for example, 100ms; the disturbance amplitude A test The PRBS perturbation signal should be small enough to avoid noticeable system disturbances, but large enough to ensure a sufficient signal-to-noise ratio for the response signal. For example, 2% of the rated excitation value should be used. Simultaneously record the reactive power responses of each. While injecting the PRBS perturbation signal, simultaneously record the actual excitation current and reactive power output values ​​of the CAES and thermal power units at a high sampling rate (e.g., 100 Hz).

[0167] The recursive least squares (RLS) method is used to identify the transfer function model between the compressed air energy storage motor and the thermal power unit. The collected input (excitation current disturbance) and output (reactive power response) data sequences are preprocessed (such as removing the DC component and applying a window function), and an autoregressive exogenous (ARX) model is 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 by online iteration of the RLS algorithm. Subsequently, this discrete model can be converted into a more physically meaningful continuous domain first-order inertial link transfer function model G(s)=K / (1+s*τ) through 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, it 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, resolving the control gain stability domain to ensure the stability of the closed-loop system, and determining the control gain.

[0169] This step is the theoretical basis for controller parameter design. Specifically, a four-dimensional state space model is constructed, whose state vector

[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 ))), the state matrix A of the system can be written. The elements of this matrix contain the physical parameters of the system {K, τ} and the controller parameters to be designed {k1, k2}.

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

[0172] Among them C critical is a system physical parameter {K CAES ,τ CAES ,K thermal ,τ thermal The set of all (k1, k2) parameter pairs that satisfy this inequality constitutes the control gain stability region that ensures the stability of the closed-loop system.

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

[0174] After parsing the control gain stability domain, it also includes: applying random perturbations to the system parameters in the coupled state-space model through Monte Carlo simulation; statistically analyzing the distribution of the system's stability probability and stability margin under disturbance; and calculating the robustness index based on the distribution of the stability probability and stability margin to evaluate the sensitivity of the control gain stability domain to parameter changes.

[0175] Specifically, after the identified parameter nominal value (such as K CAES ,τ CAES ), a random perturbation with a specific statistical distribution (for example, mean 0 and standard deviation 10% of the parameter value) is applied to generate thousands of possible system parameter samples. For each set of parameter samples, the stability boundary C is recalculated. critical Furthermore, for a given controller parameter pair (k1, k2), we can count how many times in these thousands of simulations still meet the stability condition k1*k2 <C critical , this ratio is the stability probability P stable The stability margin SM=1-(k1*k2) / C can also be calculated. critical The mean and variance of the data, and define a robustness indicator, such as R robust =μ SM -2*σ SM Through this step, controller parameters that can guarantee system stability with high probability and sufficient stability margin even when there is large uncertainty in the system parameters can be selected, thereby enhancing the robustness of the control system to a certain extent.

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

[0177] Furthermore, considering that k1 takes a positive value (decrease excitation) and k2 takes a negative value (increase excitation) in actual control, the stability domain 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 describes a preferred method for determining and applying optimal control gains based on the control gain stability domain analyzed in Example 7. It is understood that ensuring stability is a fundamental requirement for controller design, but within the stability domain, different control gain combinations (k1, k2) can result in vastly different dynamic response qualities (e.g., response speed, control accuracy, and smoothness). This example addresses the question of how to find a set of control gains that optimizes overall system performance while ensuring stability and intelligently adapts to real-time operating conditions.

[0179] It includes: constructing a comprehensive performance index consisting of the weighted sum of total reactive power deviation, excitation current change and adjustment time; under the constraint of the control gain stability domain, using the particle swarm optimization algorithm to optimize the comprehensive performance index to solve 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 (for example, fast response but smooth regulation) into a scalar function to facilitate optimization.

[0181] A preferred J function form 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 ; Among them, the first integral term ∫(Q total (t)-Q target ) 2 dt is the integral of the square of the total reactive power error (ISE), which is used to measure the accuracy of control. total (t) is the total reactive power provided by CAES and thermal power generation units at any time t; Q target is the target total reactive power. The smaller the value of this term is, the more accurately the total reactive power tracks the target value and the smaller the fluctuation is. The second integral term ∫(di f / dt) 2 dt is the square integral of the rate of change of the controlled variable, representing the control energy or regulation stability. f / dt can represent the weighted norm of the rate of change of the excitation current on both sides of CAES and thermal power. The smaller the value of this term, the smoother the adjustment action of the excitation system, which is beneficial to extend the life of the equipment and reduce the disturbance to the system. settling It is an indicator of the system's adjustment time or response speed. It is defined as the time from the occurrence of the disturbance to the total reactive deviation |Q total (t)-Q target |First entry and permanent maintenance within a very small error band (e.g. ±2%*Q target ). w1, w2, and w3 are weight coefficients for each item, all positive real numbers and typically normalized. Designers can adjust these weights to express their preference for accuracy, stability, and speed. For example, in scenarios requiring rapid response, the weight of w3 can be increased. In more sophisticated 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 use the particle swarm optimization (PSO) algorithm to optimize the comprehensive performance index. PSO is a global random search algorithm that simulates the foraging behavior of bird flocks. It is particularly suitable for solving optimization problems with complex and nonlinear objective functions like this example. Specifically: Search space: The search dimension of the algorithm is 2, that is, the collaborative control gain pair (k1, k2) to be optimized. Fitness function: Each particle represents a candidate (k1, k2) solution. By applying this group of gains in a simplified closed-loop system simulation model, the value of the corresponding comprehensive 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 the denominator being zero. The smaller the J value, 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, it is required that the position (k1, k2) of each particle must meet the stability margin requirements, for example SM=1-(k1*k2) / C critical >0.3, which means that the system is required to always maintain a stability margin of more than 30%. For particles that exceed this constraint, their fitness will be set to a minimum value (or zero), so that they will be naturally eliminated during the evolution process. Through the standard PSO algorithm process (initializing the particle swarm, iteratively updating the speed and position, updating the individual optimum and the global optimum), after a certain number of iterations, the algorithm finally converges to the global optimal particle position, which is the optimal control gain (k1 opt ,k2 opt ).

[0184] After obtaining the optimal control gain, the present invention further proposes an adaptive application strategy to enhance its robustness and performance under variable operating conditions. Specifically, when applying the optimal control gain in coordinated regulation, it also includes: dynamically adjusting the control deadband based on the standard deviation of real-time reactive power fluctuations; and segmentally scaling the optimal control gain based on the amplitude of the real-time reactive power deviation to form the actual control gain applied under different deviations.

[0185] Dynamically adjust the control dead zone to achieve the following: control dead zone ε q It is no longer a fixed empirical value, but is adjusted according to the real-time cleanliness of the system. The adjustment rule can be expressed as: q (t)=ε q 0*(1+β*σ q (t)); where ε q (t) is the actual control dead zone at the current moment; ε q0 is a basic dead zone value (for example, 0.01pu); σ q (t) is the standard deviation of the total reactive power fluctuation in the most recent time window, which reflects the current noise or disturbance level of the system; β is a positive adaptation coefficient. When the system is running smoothly and the noise is small (σq (t) is small), the dead zone is automatically narrowed, which improves the sensitivity and accuracy of control; when there is large noise or high-frequency disturbance in the system (σ q (t) is large), the dead zone is automatically relaxed, avoiding the controller’s excessive response to noise and unnecessary frequent adjustment actions (i.e., control chattering).

[0186] The optimal control gain is scaled piecewise, also known as variable gain or gain scheduling strategy, to match the controller's response strength to the deviation. A preferred piecewise scaling rule is as follows: When the absolute value of the total reactive power deviation |ΔQ| is very small (for example, |ΔQ| < 0.05 pu), it indicates that the system is in the fine adjustment stage. At this time, a smaller gain is used, for example, the actual gain k actual =0.5*k opt This can achieve small deviation, weak control, avoid overshoot, and improve steady-state accuracy. When |ΔQ| is in the normal range (for example, 0.05pu≤|ΔQ|<0.15pu), the optimal gain itself is used, that is, k actual =k opt , in order to give full play to its optimal comprehensive performance. When |ΔQ| is very large (for example, |ΔQ| ≥ 0.15pu), it indicates that the system has suffered a large impact and needs to respond quickly. At this time, a larger gain can be used, such as k actual =1.5*k opt This can achieve large deviations and strong control, bringing the system state back to the target as quickly as possible.

[0187] It can be seen that through this nonlinear adaptive gain strategy, the controller can, like an experienced operator, intelligently adjust its control strength according to the severity of the situation, thereby achieving better control effects 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). It uses the system model to foresee future dynamic responses and pre-determine the optimal control layout, thereby achieving faster and smoother control effects.

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

[0190] In the first step, the total reactive power trajectory is predicted for multiple time steps in the future based on the system model.

[0191] The controller uses a mathematical model that accurately describes the dynamic behavior of the system to deduce the future state of the system.

[0192] The model is preferably a discrete-time state-space model, whose 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 moment; X(k) is the state vector of the system, including Q CAES ,Q thermal ,i f_CAES ,i f_thermal and other key dynamic variables; U(k) is the control input vector, i.e., the regulation instruction applied to the CAES and thermal power unit excitation system; Y(k) is the output of the system, i.e., the total reactive power Q total (k); A d ,B d ,C d 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 through the 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 {Y(k+1|k),Y(k+2|k),...,Y(k+N p |k)}. Here N p Known as the look-ahead horizon, it defines how far into the future the controller can look.

[0194] The second step is to solve the optimal control sequence that minimizes the deviation between the total reactive power trajectory and the total reactive power demand of the system through rolling optimization. After obtaining the foresight of the future, the controller needs to plan a series of optimal future control action sequences {U(k|k),U(k+1|k),...,U(k+N c -1|k)}, where N c It is called the control time domain (N c ≤N p ). The optimal definition is 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 Σ||Y(k+j|k)-Q _target || 2 qy The penalty is that the predicted output trajectory is different from the target value Q in the entire prediction time domain. target The weight matrix Qy reflects the requirements for tracking accuracy at different times. The second part Σ||ΔU(k+j|k)|| 2 Ru The penalty is the size of the change ΔU of the control action 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 the problem that J mpc Minimize the optimal control sequence U opt ={U*(k|k),U*(k+1|k),...}. This solution process is repeated in each control cycle, which is called rolling optimization.

[0195] The third step is to extract the first element of the optimal control sequence to form the feedforward control instruction. opt After that, 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 used as the feedforward control instruction u at the current moment predict In the next control cycle (k+1), the controller remeasures the actual system state and, starting from this new starting point, repeats the entire process of prediction, optimization, and extraction of the first element. This mechanism not only makes the control forward-looking but also continuously utilizes the latest actual measurement information to correct its subsequent planning, making it highly robust to disturbances and model uncertainties.

[0196] To further enhance 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 demagnetization process, the uncertainty in the prediction of the total reactive power trajectory, and the real-time stability margin of the system. This means that the MPC optimization problem itself is not static but rather changes with the circumstances. Specifically, updates are performed based on the different stages of the demagnetization process: in the initial stage of the main reduction phase, speed is more important than accuracy, so the constraints on the predicted output Y(k+j|k) can be appropriately relaxed; in the terminal approach phase at the end, accuracy is the primary goal, so the constraints need to be tightened.

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

[0198] Update according to the real-time stability margin of the system: the stability margin indicator SM calculated in Example 7 can be 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 violent control actions and make the control behavior more conservative and gentle.

[0199] To suppress disturbances and errors not covered by the model, collaborative regulation 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 reactive power collected in real time and the total reactive power demand of the system, a feedback correction instruction is generated; and the feedback correction instruction is synthesized with the feedforward control instruction to form an excitation adjustment instruction sent to the compressed air energy storage motor and the thermal power unit.

[0201] In specific implementation, a conventional feedback controller (e.g., the PI controller described in detail in the eighth embodiment) works in parallel with the MPC controller. The feedback controller is based on the actual measured total reactive power Q total_actual (k) and target value Q target The deviation ΔQ between actual (k), generates a feedback correction instruction u feedback The actual control instruction u is finally sent to the excitation system. final , is the synthesis of feedforward instruction and feedback instruction, that is, 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, achieving a high degree of unity in the speed, stability and robustness of the control system.

[0203] Example 10: This example describes a preferred engineering implementation for improving the overall performance and reliability of the collaborative voltage regulation method of the present invention. This implementation primarily involves a high-fidelity data processing method for the control system input and a smooth control command generation method for the output. 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 ensure accurate and delay-free state information acquisition by the control system, the present invention proposes a high-frequency data acquisition and filtering preprocessing method. Specifically, the control system acquires data at a sampling period higher than that of conventional power system supervisory control and data acquisition (SCADA) systems. For example, the control system acquires the raw reactive power signal Q of the CAES and thermal power units at a sampling period of 5 ms (i.e., a sampling rate of 200 Hz). CAES_raw and Q thermal_raw Since the original signal collected by high-frequency sampling inevitably contains measurement noise and random disturbances, if this signal is directly used for feedback control, it is easy to cause jitter in the control output. Therefore, the present invention preferably uses a Kalman filter to perform online filtering on the original signal.

[0205] The Kalman filter is an optimal linear state estimator that can provide an optimal estimate of the system state in a dynamic system with uncertainty 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. The 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) is the original measurement value at time k; H is the observation matrix; K f is the Kalman gain. Kalman gain K f Dynamic adjustments are made based on the prediction error covariance and measurement noise covariance. Compared to traditional moving average or low-pass filters, the Kalman filter 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, in order to ensure that the regulation command calculated by the control algorithm can be safely and smoothly executed by the excitation system, the present invention proposes a multi-level limiting and smooth output method for the control command. When the cooperative control algorithm (such as the eighth or ninth embodiment) outputs an ideal excitation current regulation rate di f / dt or target excitation current i f_cmdAfter the instruction is sent, it will not be directly sent to the excitation power unit, but will first pass through a command conditioning module. This module performs multi-level limiting operations to protect physical equipment. Change rate limit: The change rate of the instruction must not exceed the maximum dynamic response rate of the excitation system itself, that is, |di f / dt| <k max_dynamic Absolute value limit: The target value that satisfies the instruction will 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 allowed short-time excitation current. Acceleration limit (optional): To further protect the device, the rate of change of the command (i.e. acceleration) can also be limited. 2 i f / dt 2 | max .

[0208] After the clipping process, the instruction sequence may produce non-smooth inflection points at the clipping boundary. To solve this problem, the module further performs smooth output processing. A preferred method is to use cubic spline interpolation. Specifically, a series of discrete instruction points that have been clipped are used as interpolation nodes, and the controller calculates a set of segmented cubic polynomials i f_cmd (t)=a3*t 3 +a2*t 2 +a1*t+a0 to connect these nodes. The characteristics of cubic spline interpolation ensure that not only are the function values ​​continuous at the nodes, but also their first- and second-order derivatives are continuous. This continuous curve, formed by this set of smooth polynomials, serves as the final command signal to the excitation controller. This ensures that the command sent to the power electronics is continuous and has a smooth first-order derivative. This reduces the impact on the power electronics, avoids the generation of high-frequency harmonics, and ensures a more accurate and smoother final excitation current response.

[0209] Example 11: This example describes the online monitoring and anomaly protection mechanisms implemented during the demagnetization process to ensure system safety and stability, as well as the precise timing determination logic designed to achieve seamless disconnection at the end of the process. These mechanisms are key to transitioning the present invention from theory to engineering practice, ensuring robustness and ultimate effectiveness.

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

[0211] Tracking error monitoring: Real-time calculation of actual excitation current i f_CAES_actual ​​With the target trajectory i f *(t) normalized deviation e track =|if_ CAES_actual -i f *(t)| / i f0 If the deviation exceeds the 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 time, the protection mechanism is triggered. A preferred treatment measure is to actively reduce the risk, for example, by reducing the original demagnetization rate k demag (or k fast ) is reduced by 20%, and the subsequent demagnetization trajectory is 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 of transient stability. If dδ / dt exceeds an emergency threshold (e.g. 10° / s), it usually indicates that the system is rapidly sliding towards the loss of step boundary. At this time, the emergency stability control program should be immediately activated, for example, the current demagnetization rate k demag Immediately reduce the power angle to half or even zero, that is, suspend the demagnetization process and prioritize maintaining system synchronization. At the same time, offset the maximum power angle in this event by δ max and power angle recovery time t recover Key indicators such as performance evaluation and model calibration 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 in the actual implementation process to the theoretically safe and feasible demagnetization plan (such as Example 6), thereby achieving the maximum possible maintenance of system stability in the face of unexpected disturbances or system anomalies.

[0214] When the demagnetization process is nearing the end, the present invention starts a zero excitation disconnection timing accurate judgment 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 the CAES is monitored f_CAES_actual The first time it drops to the minimum excitation threshold i f_min (For example, 0.05*i f_rated ) is below, the system enters the off-grid preparation stage. In this stage, the controller continuously monitors the residual reactive power Q of the CAES motor with higher frequency and accuracy. residual .

[0215] In order to achieve accurate judgment, the following three conditions must be met at the same time: Current amplitude condition: i f_CAES_actual f_min , so that the motor is in a deep demagnetization state. Reactive amplitude condition: the absolute value of the residual reactive power |Q residual ​| is 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 off-grid.

[0216] Reactive power stability condition: To prevent misoperation during the oscillation of reactive power crossing zero, a stability criterion is added. Specifically, the residual reactive power Q residual The standard deviation of fluctuation σ Q_residual Require σ Q_residual Less than a preset minimum value (for example, 0.01pu) to confirm that the residual reactive power has stabilized rather than instantaneously crossing zero.

[0217] If and only if the above three conditions are met at a certain moment, the controller will generate the off-grid enable signal ENABLE disconnect In an optional preferred embodiment, when generating ENABLE disconnect After receiving the signal, the controller will not immediately issue a trip command, but will wait for the best electrical opportunity. Specifically, it will monitor the instantaneous value of the three-phase stator current of the motor, and will issue a trip command CMD to the corresponding circuit breaker only when it detects that any phase current naturally crosses the zero point. open This zero-current breaking technology suppresses arcing generated when the circuit breaker contacts open, extending the life of the circuit breaker while further reducing electromagnetic transients generated by the disconnection operation. Through this precise timing and conditional judgment logic, the present invention achieves the final disconnection of the CAES motor from the grid at an electromagnetic static point that is friendly to both the grid and the device itself.

[0218] Example 12. This example provides a specific, non-restrictive numerical calculation example to illustrate in detail the complete calculation process of the multi-constrained 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 the present invention based on the steps shown in this example.

[0219] Assume that in a specific switching scenario, the following initial operating state parameters and preset parameters of the thermal power coupled compressed air energy storage (CAES) power station and the power grid system are obtained through the method of 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 potential E0 calculated by the refined method of Example 3 = 1.1 (pu). Motor and system parameters: d-axis synchronous reactance X d=1.2(pu). Unit inertia 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 50Hz system). Excitation system equivalent gain coefficient K f =1.5(pu / pu). Electromechanical oscillation natural frequency ω n =5.0(rad / s). Control and constraint parameters: Transient stability critical power angle δ critical =75° (i.e. 75*π / 180rad≈1.309rad). The desired minimum system damping ratio ξ=0.707.

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

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

[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 in order to maintain static stability, the rate of decrease of the excitation current should not exceed 0.5818 per second in theory. Calculate the transient stability boundary B transient , unify the angle unit into radians: δ0=π / 6rad≈0.5236rad; δ critical ≈1.309rad. According to formula B transient =(2*H*ω s *(δ critical -δ0)) / (T d 0'*cos(δ0)), substitute the value:

[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 base value of the excitation current and the pu value. For simplicity, pu / s is still used here); it should be noted that the numerical value of the calculated result here is large, indicating that under this specific working condition, the transient stability constraint is relatively loose. This may be due to the large inertia constant H and sufficient power angle margin. In order to maintain the unity and rationality of the units, we re-examine the dimension of this boundary, which should be consistent with the rate of change of the excitation current. Assuming the base value of the excitation current 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] So far, three independent boundary values ​​have been obtained: B static ≈0.5818,B transient ≈2.279,B dynamic ≈2.5923 (units are pu / s).

[0227] In the second step, the three boundaries are fused to determine the preliminary comprehensive boundary (corresponding to the fourth embodiment) and the weighted average boundary B is calculated. weighted , using the exemplary weight given in Example 4: 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 working condition, static stability is the most important limiting factor of the system. Determine the preliminary comprehensive boundary k max_raw :Adopt 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, evaluating historical performance and current margins:

[0229] Retrieve data from the historical database and calculate the comprehensive historical performance index P history =0.9. According to the rule of Example 4, P history ≥0.8 is considered as average or poor historical performance. At the same time, assuming that the current reactive power margin Q obtained from the DCS system of the thermal power unit is margin Low, does not meet the conditions for increasing the safety factor. Dynamically adjust the safety factor α safe :According to P history =0.9, query the piecewise function rule, and choose the basic safety factor as 0.55. Due to insufficient reactive margin, the factor is not adjusted upward. Therefore, the final dynamic safety factor α safe =0.55. Calculate the final comprehensive stability boundary k max :k max =α safe *k max_raw k max =0.55*0.6982≈0.384(pu / s).

[0230] Through the above complete multi-level calculation, the comprehensive stability boundary k under this specific working condition is finally obtained. max ≈0.384 pu / s. This value will be used as the demagnetization rate k of the main drop section when planning the demagnetization trajectory in Example 5. fast The upper bound of the benchmark (e.g., k fast 0.9*k can be taken max ≈0.3456 pu / s). This example demonstrates how the present invention transforms a complex theoretical model and multi-dimensional constraints into a specific control parameter that can guide engineering practice, fully demonstrating the advancement, rigor, and practicality of the present invention.

[0231] Thirteenth Example: This example describes an important, optimal extension of the collaborative voltage regulation method of the present invention, embodying system intelligence and self-learning capabilities. This extension, executed after a complete FC-CAES collaborative switching task, quantitatively evaluates the control performance of the switching process and uses the evaluation results to iteratively optimize control parameters, thereby improving system performance 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 switching is completed includes the following steps: the first step is to collect and calculate the key performance indicators (KPIs) of the switching process. After the CAES motor is successfully disconnected from the grid and the switching task is completed, the control system will retrieve the entire switching process from the historical recorder (from t trigger to t disconnect Based on this data, the system automatically calculates a series of predefined KPIs to comprehensively and quantitatively evaluate the effectiveness of this control.

[0233] Preferred KPIs may include: Voltage stability indicator: Maximum voltage deviation ΔV max :ΔV max =max(|V bus (t)-V nominal |); where V bus (t) is the time series of bus voltage, V nominal is the rated voltage. This indicator directly reflects the suppression effect on the grid voltage shock. Voltage fluctuation duration t fluctuation : Refers to the time when the voltage first deviates from the normal range (for example ±1%*V nominal ) to finally stabilize in this range. This indicator measures how quickly the system returns to stability.

[0234] Reactive power control accuracy index: reactive power total deviation integral E Q_integral :E Q_integral =∫[0,t disconnect ]|Q total (t)-Q target |dt. This indicator measures the tracking accuracy of the total reactive power to the target value during the entire coordinated regulation process. The smaller the value, the more accurate the coordinated control. Coordination smoothness indicator: thermal power reactive power smoothness S smooth :S smooth =max(|dQ thermal / dt|) / Q thermal_rated This indicator measures the maximum rate of change of the response process of the thermal power unit when it undertakes the reactive power compensation task, reflects the smoothness of the coordination process, and avoids excessive impact on the thermal power unit.

[0235] The second step is to make parameter optimization decisions based on the KPI evaluation results. The system compares the calculated KPI with a set of preset performance benchmarks. For example, the performance benchmark is: ΔV max_limit =1%,t fluctuation_limit =500ms. If all KPIs of this switch are better than the performance benchmark, the current control parameters (such as collaborative control gains k1, k2) are considered to be performing well and do not need to be adjusted. If any KPI does not meet the benchmark requirement (for example, this ΔV maxIf the KPI reaches 1.2%, the online parameter optimization program is started. This program fine-tunes 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 , which is a weighted sum of the KPIs exceeding the benchmark. By performing sensitivity analysis on the control system model, J perf Gradient ΔJ of the parameters to be optimized (such as k1, k2) perf (k i ). Further, update the parameters according to the following rules: k i_new =k i_old -η*ΔJ perf (k i ); where k i_new and k i_old are the parameter values ​​before and after the update, respectively; η is a small learning rate, such as 0.01.

[0237] The third step is the maintenance of the historical database and the iterative update of the optimal parameter set. In order to achieve long-term, cross-task learning and evolution, the present invention also includes a historical database DB switching After each switching task is completed, the system will store the feature parameter set of this task in the database. The feature parameter set is a vector that contains the input, process and result of this task, such as {τ demag ,k1,k2,ΔV max ,t fluctuation ,E Q_integral ,S smooth To avoid unlimited data growth, the database preferably adopts a sliding window mechanism, for example, only retaining the latest 100 switching records. Based on the database, the system will periodically or after each update, use the weighted average method to update a set of global optimal control parameters. When calculating the weighted average, a weight can be assigned to each historical record, and the weight is related to the performance of the record (for example, with 1 / J perf This is directly proportional to the performance of the switching case. That is, the better the historical switching case, the greater the weight its parameters will have in calculating the global optimal parameter set. This optimal control parameter set, updated through weighted iteration based on historical data, will serve as the initial default parameters for the controller during the next switching task.

[0238] According to one aspect of the present 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] The 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 time. 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 compared with the target trajectory i stored in the controller memory f *(t) performs real-time comparison 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 recursive least squares method is used to estimate the constant deviation and drift rate of the tracking error model online.

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

[0242] On this basis, the bias compensation and slope compensation are generated according to the constant deviation and drift rate to construct a corrected trajectory under the comprehensive stability boundary constraint to guide subsequent control.

[0243] The generation of compensation follows a hierarchical strategy. f0 When the offset 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 When the slope compensation Δk is generated comp =-0.7·a1. The construction of the correction trajectory adopts 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 To detect the moment that needs correction. After correction, verification is required. f_corrected * / dt|≤0.9·k max If it is not satisfied, 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 show that the trajectory tracking accuracy is improved from ±5% to ±2% after enabling the correction, while maintaining the stability of the system.

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

[0247] This invention solves the technical problem of the inability to make a scientific decision between speed and safety in the switching control strategy by accurately quantifying the comprehensive stability boundary of the switching process online. This effect is achieved by unifying the modeling and calculation of stability constraints at three different time scales: static, transient, and dynamic, and realizing the determined boundary k max_raw The comprehensiveness of the system avoids the potential instability risk caused by a single consideration. history and the current reactive power margin Q of thermal power units margin Adaptive safety factor α safe , so that the determination of the stability boundary is no longer a static, conservative, one-size-fits-all value, but a dynamic intelligent decision that can reflect the system's historical performance and current support capabilities. This method transforms the abstract out-of-sync risk into a specific, operationally guiding rate limit k max This provides a solid theoretical basis for subsequent demagnetization trajectory planning, and can maximize the execution efficiency of the switching process while ensuring safety, thus breaking away from the dilemma of traditional methods that rely on experience and cannot balance speed and safety.

[0248] The present invention solves the problem of bus voltage transient stability caused by instantaneous reactive power surge by constructing a mathematically smooth and dynamic process pre-verified demagnetization trajectory, transforming the instantaneous hard cutoff of the CAES unit into a fully controllable reactive power soft landing process lasting hundreds of milliseconds. On the one hand, it utilizes a segmented trajectory i composed of linear, exponential, and cosine functions. f*(t), each section optimizes the switching speed at the initial stage, the stability at the middle stage, and the accuracy at the final stage. In particular, the terminal approach section achieves zero rate of change of the excitation current when it reaches the target point, and physically achieves a smooth transition. On the other hand, by pre-calculating the power angle and reactive power before trajectory execution, and iteratively optimizing the trajectory that violates the constraints using the gradient projection method, the planned route map i f_optimal *(t) is theoretically safe and feasible. This technology, which actively creates and manages the transition state, eliminates the sudden change in dQ / dt, allowing the grid voltage to remain highly stable throughout the switching process, with the maximum deviation reduced from 3-5% in traditional methods to less than 0.5%.

[0249] The present invention obtains the most realistic dynamic characteristics of the system and designs an optimal and adaptive controller based on this, so that the reactive compensation of the thermal power unit can be accurately, quickly and smoothly taken over during the demagnetization process, thereby ensuring the dynamic balance of the total reactive power of the system. This effect is achieved by injecting PRBS signals into the system and using the RLS algorithm to identify the transfer function model G(s) online, overcoming the performance degradation problem caused by the traditional controller's reliance on offline and solidified model parameters. On this basis, the stable control gain domain is analyzed by the Routh-Hurwitz criterion, and the particle swarm algorithm (PSO) is used to optimize a comprehensive performance indicator J that takes into account accuracy, stability and speed. The obtained control gain (k1 opt ,k2 opt ) theoretically achieves Pareto optimality. By dynamically adjusting the control deadband and the step-wise scaling gain strategy, the controller intelligently adapts to varying disturbance levels and deviation amplitudes. This series of technical features works together to achieve high performance and strong robustness of the collaborative control system.

[0250] This invention integrates the foresight of model predictive control (MPC) with the robustness of traditional feedback control to build a high-performance collaborative control system that can foresee the future and make timely corrections, further improving the response speed and control accuracy of reactive power compensation. The optimal feedforward control instruction u is calculated in advance through rolling optimization. predict This control method, which takes advantage of the deviation before it occurs, shortens the system's response delay. At the same time, the parallel feedback correction loop is responsible for processing unknown disturbances and parameter mismatch errors not covered by the model, generating feedback instructions u feedback By adding u predict and u feedback This dual-loop structure enables the control system to not only actively respond to predictable dynamic processes, but also passively and robustly eliminate unpredictable errors. Its collaborative control performance (such as the integral E of the reactive tracking error) is very good. Q_integral ) is improved compared to a single feedback control strategy.

[0251] By establishing a mathematical relationship based on the total reactive power deviation for bidirectional, real-time, closed-loop regulation of the excitation of the CAES and thermal power units, this invention achieves dynamic balance of the system's total reactive power during the active demagnetization transition state, resolving the problem of secondary voltage fluctuations caused by delayed or uncoordinated compensation. The CAES demagnetization process and the thermal power compensation process are transformed from two independent, time-sequential open-loop actions into an 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. Instead, its real-time response speed and amplitude precisely and inversely match the rate of decrease in CAES reactive power. This dynamic balance, with one decreasing and one increasing, and a constant total amount, makes the two units appear, from the perspective of the power grid, as a single virtual unit with a constant total reactive power output. This maintains high bus voltage stability throughout the switching process and avoids the under- or over-compensation problems that can occur with traditional methods.

[0252] The preferred embodiments of the present invention are described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within 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 scope of protection of the present invention.

Claims

1. A method for collaborative voltage regulation of a thermal power plant coupled with a 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; Execute the target demagnetization trajectory and coordinately adjust 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; When the reactive output of the compressed air energy storage motor meets the preset off-grid conditions, it will be disconnected from the grid.

2. The method according to claim 1, characterized in that Planning and generating the target demagnetization trajectory of the compressed air energy storage motor, including: 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.

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

4. The method according to claim 2, characterized in that The three stability boundaries are fused to determine the comprehensive stability boundary, including: Assign preset credibility weights to the static, transient and dynamic stability boundaries respectively, and perform weighted fusion to obtain the weighted average boundary; The minimum value among the three stable boundaries is selected as the minimum boundary; The comprehensive stability boundary is determined based on the weighted average boundary and the minimum boundary.

5. The method according to claim 1, wherein Planning and generating the target demagnetization trajectory of the compressed air energy storage motor is to construct a piecewise function consisting of the following parts connected in sequence, including: Set the main step-down section to execute the main excitation reduction at a basically constant rate; The transition section used to ease the rate change; and a terminal approach section for making the rate of change of the excitation current approach zero when the excitation current reaches the target value.

6. The method according to claim 5, characterized in that: The main descending segment is composed of the linear function i f (t)=i f0 -k fast *t defines, where i f (t) is the target excitation current corresponding to time t, i f0 is the initial excitation current obtained from the initial operating state data, k fast is the demagnetization rate determined based on the comprehensive stability boundary; The transition section is composed of the exponential function i f (t)=i f1 *exp(-(t-t1) / τ smooth )+i f_target Define, where i f1 and t1 are the excitation current value and time point at the end of the main drop section, τ smooth is the smoothing time constant, i f_target is the transition target excitation 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))) defined, where i f2 t and t2 are the excitation current value and time point at the end of the transition period, respectively. f_min is the preset minimum excitation current, t final is the end time point of the demagnetization process.

7. The method according to claim 1, characterized in that After planning and generating the target demagnetization trajectory of the compressed air energy storage motor, it also includes: Based on the target demagnetization trajectory, the power angle trajectory and reactive power prediction trajectory during the demagnetization process are pre-calculated; Check whether the power angle trajectory and reactive power prediction trajectory meet the preset power angle constraint and reactive power change rate constraint; If any constraint is not satisfied, the parameters of the target demagnetization trajectory are iteratively optimized using the gradient projection method until all constraints are satisfied.

8. The method according to claim 1, characterized in that Co-regulation is based on control gains determined through online system identification, which includes: Pseudo-random binary sequence disturbances are injected into the excitation systems of the compressed air energy storage motor and the thermal power unit respectively. Synchronously record the respective reactive power responses; Based on the disturbance and response, the recursive least squares method is used to identify the transfer function model of the compressed air energy storage motor and the thermal power unit to determine the control gain.

9. The method according to claim 8, characterized in that After the transfer function model is identified, it also includes: 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; The Routh-Hurwitz criterion is applied to analyze the coupled state space model, and the control gain stability domain that ensures the stability of the closed-loop system is analyzed to determine the control gain.

10. The method according to claim 9, characterized in that Determining the control gains involves: Construct a comprehensive performance index consisting of the weighted sum of total reactive power deviation, excitation current variation, and regulation 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 solve the optimal control gain.

Images