Flywheel energy storage grid frequency support method based on model predictive control

Abstract

This invention provides a method for supporting grid frequency for flywheel energy storage based on model predictive control, comprising: S1. Synchronous acquisition of grid operating parameters and flywheel energy storage system state parameters; S2. Preprocessing of acquired parameters; S3. Construction and parameter identification of a dual-model system with the DC bus as the intermediate link; S4. Implementation of a grid-side power prediction and collaborative control algorithm; S5. Implementation of a turbine-side flywheel output power prediction and collaborative control algorithm. This invention achieves accurate and rapid grid frequency support by constructing a multi-model collaborative system, designing a multi-algorithm bidirectional closed-loop interactive collaborative control strategy, and combining full-process data preprocessing and parameter adaptive closed-loop update mechanisms. Simultaneously, it maintains stable DC bus voltage, ensures the safe and reliable operation of the flywheel energy storage system, and improves the overall system's adaptive adjustment capability and disturbance rejection performance under all operating conditions.

Description

Technical Field

[0001] This invention relates to the field of power system energy storage control and grid frequency regulation technology, specifically a flywheel energy storage grid frequency support method based on model predictive control. Background Technology

[0002] As the proportion of new energy power generation in the power system continues to increase, the intermittent, volatile, and low-inertia characteristics of new energy power generation such as wind power and photovoltaics after their connection to the grid can lead to a significant decrease in the rotational inertia and damping characteristics of the power grid. This makes the grid frequency highly susceptible to large and rapid fluctuations, seriously threatening the safe and stable operation of the power system. Frequency support is a key technology for maintaining grid frequency stability in a power system. Flywheel energy storage systems, with their fast response speed, high power density, long charge-discharge cycle life, and pollution-free operation, are ideal energy storage methods for smoothing grid frequency fluctuations and providing frequency support, and have become a key area of ​​research and application in the field of power system frequency control.

[0003] Current control methods for flywheel energy storage grid frequency support mostly employ a single algorithm or a simple series of multiple algorithms, such as standalone virtual synchronous generator (VSG) control, PID control, or model predictive control (MPC), or unidirectional series control of VSG with PID or MPC. However, these methods have several technical drawbacks: First, while VSG control can provide virtual inertia and damping for the grid, it lacks global optimization capabilities and struggles to balance control performance under multiple constraints. Second, while model predictive control offers advantages in global optimization under multiple constraints, it is susceptible to prediction errors caused by internal and external disturbances to the grid, reducing control accuracy. Third, while PID control is simple in structure and fast in response, it has poor adaptability and struggles to adapt to changes in grid and flywheel operating conditions. Fourth, when multiple algorithms are combined, they are mostly unidirectional signal transmissions without bidirectional closed-loop interaction, failing to fully leverage the advantages of each algorithm and easily leading to control objective conflicts. Fifth, grid-side power control and turbine-side flywheel control lack deep coordination, resulting in poor power matching with the DC bus as an intermediate link, which can easily lead to drastic fluctuations in DC bus voltage and affect system stability. Summary of the Invention

[0004] To address the problem that existing flywheel energy storage grid frequency support methods cannot meet the actual accuracy and stability requirements under complex operating conditions, this invention provides a flywheel energy storage grid frequency support method based on model predictive control. This method can achieve accurate and rapid grid frequency support under complex operating conditions, while maintaining DC bus voltage stability, ensuring the safe and reliable operation of the flywheel energy storage system, and improving the overall system's adaptive adjustment capability and anti-disturbance performance under all operating conditions.

[0005] The technical solution of this invention is as follows: a flywheel energy storage grid frequency support method based on model predictive control, characterized by comprising the following steps: S1. Synchronously collect power grid operating parameters and flywheel energy storage system status parameters; The power grid operating parameters include: real-time power grid frequency, power grid voltage amplitude, grid-connected current of energy storage system, power grid frequency change rate, and power grid frequency deviation. The status parameters of the flywheel energy storage system include: real-time flywheel speed, actual flywheel charging and discharging power, flywheel motor stator current, DC bus voltage, flywheel mechanical loss power, and flywheel energy storage capacity. S2. Perform preprocessing on all parameters acquired and calculated in step S1, and output the full set of valid data, which are the filtered values ​​of all parameters after processing. The preprocessing operations include: outlier detection and correction, smoothing and denoising, and validity verification. S3. Based on the aforementioned valid data, construct a dual-model system with the DC bus as the intermediate link, and perform parameter identification; The dual-model system includes: a grid-side power prediction model and a flywheel energy storage device-side dynamic model; The grid-side power prediction model outputs a reference value for grid-connected power prediction to support grid frequency. The dynamic model output of the flywheel energy storage machine side is used to predict the flywheel output power for DC bus voltage support, as well as the charging and discharging power constraint boundary and the current constraint boundary. S4. Based on effective data, grid-side power prediction model and charging / discharging power constraint boundary, integrate model predictive control, disturbance observation and virtual synchronous generator control to construct a grid-connected power prediction cooperative control algorithm and output the optimal grid-connected power prediction value; S5. Based on effective data, the dynamic model of the flywheel energy storage machine, the charging and discharging power constraint boundary, the current constraint boundary, and the optimal grid-connected power prediction value, the flywheel output power prediction collaborative control algorithm is constructed by integrating model predictive control, PID control, and fuzzy control. The algorithm outputs the optimal control quantity of the flywheel motor and drives the flywheel motor to perform the corresponding charging, discharging, or standby actions.

[0006] Its further features are: It also includes step S6, which calculates the actual deviation of the control effect based on the actual operating data after S1 to S5, compares the actual deviation with the target deviation to complete the adjustment effect evaluation; updates the weights of the dual-model system, the grid-connected power prediction collaborative control algorithm and the flywheel output power prediction collaborative control algorithm in a closed loop according to the evaluation results, and synchronizes the updated parameters and weights to the corresponding models and algorithms as inputs for the next acquisition cycle to form closed-loop adjustment, repeating S1 to S6 until the grid frequency stabilizes within the rated range; In step S2, the outlier detection and correction are implemented using the 3σ criterion; the smoothing and denoising are implemented using a moving average filter; the validity check is implemented by pre-setting a corresponding physical reasonable range for each parameter, and invalid data that exceeds the valid range is replaced by the valid filter value of the previous period. In step S3, the grid-side power prediction model integrates the second-order linear model of the grid-connected filter, the second-order disturbance observer, and the second-order model of the virtual synchronous generator to describe the mapping relationship between the control input and the grid-connected output power and to observe the total internal and external disturbances that affect the model's accuracy; the grid-side power prediction model is: X g (k+1)=A gd ·X g (k)+B gd ·U g (k)+d total_g (k); In the formula, X g (k+1) represents the predicted value of the network-side state variable in the (k+1)th sampling period; X g (k) represents the grid-side state variables in the kth sampling period, including the grid-side converter output current, the grid-connected current, and the filter capacitor voltage; A gd B is the discretized network-side state matrix. gd U is the discretized network-side input matrix; g (k) represents the grid-side converter output voltage control quantity during the kth sampling period; d total_g (k) represents the total internal and external disturbances of the network-side model during the k-th sampling period; The dynamic model of the flywheel energy storage device is a second-order linear model used to describe the mapping relationship between control input and flywheel output power. The model parameters are identified and updated online based on effective data from a preset number of sampling periods using the recursive least squares method. The dynamic model of the flywheel energy storage device is as follows: ω f_pred (k+1)=ω f_filt (k)+(T s / J f )·(1.5·p n ·ψ f ·i sq_ref (k)-B ω ·ω f_filt (k)-T load_filt (k)); P f_pred (k+1)=1.5·p n ·ψ f ·i sq_ref (k)·ω f_pred (k+1); U dc_pred(k+1)=U dc_filt (k)+(T s / C dc )·(P dc_g (k)-P dc_f (k)-P dc_loss (k)) / (U dc_filt (k)); In the formula, ω f_pred (k+1) is the predicted flywheel speed in the (k+1)th sampling period; ω f_filt (k) represents the filtered value of the flywheel real-time speed during the kth sampling period; i sq_ref (k) represents the reference value of the flywheel motor stator q-axis current during the kth sampling period; T load_filt (k) represents the flywheel load torque filter value in the kth sampling period; T s For the sampling period, J f p is the moment of inertia of the flywheel rotor. n Let ψ be the number of pole pairs of the flywheel motor. f For the permanent magnet flux linkage of the flywheel motor; P f_pred (k+1) is the predicted flywheel output power value for the (k+1)th sampling period; U dc_pred (k+1) is the predicted value of the DC bus voltage in the (k+1)th sampling period; U dc_filt (k) represents the DC bus voltage filter value for the kth sampling period; C dc P is the capacitance value of the DC bus capacitor; dc_g (k) represents the power flowing from the grid-side converter to the DC bus during the k-th sampling period; P dc_loss (k) represents the DC bus line and switching loss power in the kth sampling period; P dc_f (k) represents the power flowing from the DC bus to the machine-side converter during the kth sampling period; In step S3, the charging and discharging power constraint boundary is determined based on the rated power of the flywheel motor, the current flywheel energy storage capacity, the real-time speed of the flywheel, and the safe range of the DC bus voltage. P f_min (k)≤P f (k)≤P f_max (k); Among them, P f_max (k) represents the maximum allowable charging power in the kth sampling period, which is determined by the rated power of the flywheel motor, the current flywheel energy storage capacity filter value, and the upper limit of the DC bus voltage; P f_min (k) represents the maximum allowable discharge power in the kth sampling period, which is determined by the rated power of the flywheel motor, the current real-time flywheel speed filter value, and the lower limit of the DC bus voltage; P f (k) represents the actual charging and discharging power of the flywheel motor port at the current moment; The current constraint boundary is determined based on the rated stator current of the flywheel motor and the motor demagnetization characteristics; i s_min ≤i sq_ref (k)≤i s_max ; Among them, i s_max i represents the positive upper limit of the q-axis current of the motor stator, corresponding to the maximum charging torque; s_min The negative lower limit of the q-axis current of the motor stator corresponds to the maximum discharge torque, both determined by the motor's rated current and demagnetization characteristics; i sq_ref (k) is the reference value of the q-axis current of the flywheel motor stator in the kth sampling period; In step S4, the grid-connected power prediction and coordinated control algorithm is implemented through a two-way, two-way closed-loop interaction of three types of algorithms: model predictive control, disturbance observation, and virtual synchronous generator control. The specific process is as follows: The virtual synchronous generator control provides a power reference trajectory that meets the grid frequency support requirements for model predictive control. The model predictive control corrects the power reference trajectory in the feasible region based on the charging and discharging power constraint boundary, and updates the corrected reference value in reverse to the input of the virtual synchronous generator control. The formula for calculating the power reference trajectory in the prediction time domain is: ; In the formula, P g_ref (k+i|k) represents the grid-connected power reference value for the k+i sampling period in the k-th sampling period; P0 represents the grid-connected steady-state active power reference value of the flywheel energy storage system; k p The active-frequency droop factor; Δf g_filt (k) represents the power grid frequency deviation filter value for the kth sampling period; k i The integral regulation coefficient for virtual synchronous generator control; Δf g_filt (k+m|k) is the predicted frequency deviation value for the (k+m)th sampling period; The formula for correcting the feasible region is: P g_ref_corr (k+i|k)=sat(P g_ref (k+i|k),P g_min_corr ,P g_max_corr )); In the formula, P g_ref_corr (k+i|k) is the corrected grid-connected power reference value predicted at time k for the (k+i)th sampling period; sat(·) is the saturation limiting function; P g_max_corr P is the corrected maximum allowable grid-connected power. g_min_corr The revised minimum allowable grid-connected power; The corrected formula for updating the virtual synchronous generator control with the reference value in reverse is: P m_ref (k+1)=P g_ref_corr (k+1|k); In the formula, P m_ref (k+1) is the virtual mechanical power reference value of the virtual synchronous generator in the (k+1)th sampling period; The total disturbance observation value output by the disturbance observer is embedded into the prediction model of the model predictive control in real time to complete error correction, and the prediction deviation of the model predictive control is input back into the disturbance observer to update the total disturbance observation value. The formula for calculating the prediction deviation is: e g (k)=Y g_filt (k)-C gd ·X g (k|k-1); In the formula, e g (k) represents the prediction bias of the network-side model in the kth sampling period; Y g_filt (k) is the filtered value of the grid-connected output in the kth sampling period, C gd X is the discretized grid-side output matrix in the grid-side power prediction model. g (k|k-1) represents the predicted state variable value for the k-th period obtained from the model predictive control during the (k-1)-th sampling period; The virtual synchronous generator control adaptively adjusts the virtual inertia and damping coefficient based on the output of the disturbance observer and the grid frequency change rate. The disturbance observer observes and compensates for the internal disturbances of the grid-side prediction model caused by parameter adjustments in real time. The adaptive adjustment formula is: J vsg (k)=J vsg0 +k j ·|d total_g_obs (k)|; D vsg (k)=D vsg0 +k d ·|ROCOF filt (k)|; Among them, J vsg (k) represents the virtual moment of inertia adaptively adjusted in the kth sampling period; D vsg (k) is the virtual damping coefficient adaptively adjusted in the kth sampling period; J vsg0 D is the nominal reference value for the virtual moment of inertia. vsg0 The nominal reference value for the virtual damping coefficient; k j k is the virtual inertia adjustment factor. d ROCOF is the virtual damping adjustment factor. filt (k) represents the filtered value of the power grid frequency change rate in the kth sampling period; dtotal_g_obs (k) represents the total network-side disturbance observation value output by the disturbance observer; In step S4, the cost function of model predictive control is constructed, the constrained quadratic programming optimization problem is solved, and the optimal control quantity and the optimal grid-connected power prediction value are output based on the rolling time-domain optimization principle. The cost function formula on the network side is: ; In the formula, J g Let N be the cost function of the predictive control of the network-side model, and the optimization objective is to minimize the value of this function; p For prediction in the time domain, N c To control the time domain; Q g Y is the weighting coefficient for the output tracking error; g_ref (k+i|k) is the output reference value generated based on the corrected power reference trajectory; Y g (k+i|k) represents the predicted output value for the (k+i)th cycle obtained from the network-side prediction model; R g The weighting coefficient for the rate of change of the control quantity; △U g (k+i|k) represents the change in the control variable; ρ g The penalty weight for the relaxation factor; ε g These are slack variables; The constraints corresponding to the cost function include: prediction model constraints, control quantity constraints, output quantity constraints, and power constraints. The constraint of the prediction model: X g (k+i|k)=A gd ·X g (k+i-1|k)+B gd ·U g (k+i-1|k)+d total_g_obs (k); In the formula, X g (k+i|k) represents the predicted value of the network-side state variable for the (k+i)th sampling period from the kth sampling period; d total_g_obs Output the total disturbance observation value for the second-order disturbance observer; U g (k+i-1|k) is the predicted value of the grid-side converter output voltage control quantity for the k+i-1th period in the kth sampling period; A gd B gd These are the state matrix and input matrix of the network-side discrete prediction model, respectively; the value of i ranges from 1 to N. p ; The control constraint: U g_min ≤U g (k+i|k)≤U g_max ; In the formula, the value of i ranges from 0 to N.c -1; U g_min U g_max These are the upper and lower limits of the grid-side converter output voltage; The output constraint: i g_min ≤i g (k+i|k)≤i g_max ; In the formula, the value of i ranges from 1 to N. p i g_min i g_max These are the upper and lower limits of the grid-connected current; i g (k+i|k) is the instantaneous value of the grid-connected phase current of the flywheel energy storage system in the (k+i)th sampling period predicted at the kth sampling time; The power constraint: P g_min_corr ≤P g (k+i|k)≤P g_max_corr , In the formula, the value of i ranges from 1 to N. p ;P g_max_corr P is the corrected maximum allowable grid-connected power. g_min_corr P is the revised minimum allowable grid-connected power. g (k+i|k) is the predicted active power of the flywheel energy storage system at the k+i sampling period, obtained at the k-th sampling time. In step S5, the flywheel output power prediction and cooperative control algorithm is implemented through a two-way, two-way closed-loop interaction of three types of algorithms: model predictive control, PID control, and fuzzy control. The specific process is as follows: The fuzzy control dynamically outputs the weighting coefficients of the model predictive control and the PID control based on the DC bus voltage deviation and flywheel speed deviation. The DC bus voltage tracking error of the model predictive control and the flywheel speed tracking error of the PID control are input inversely to the fuzzy control, thereby completing the adaptive adjustment of the fuzzy control quantization factor and the universe of discourse; specifically as follows: In the positive interaction, the weighting coefficients of the fuzzy control output are directly used for the weighted fusion of the two control quantities to obtain the final optimal stator current reference value of the flywheel motor. The weighted fusion formula is as follows: i sq_opt (k)=α mpc (k)·i sq_mpc (k|k)+α pid (k)·i sq_pid_sat (k); Among them, i sq_opt (k) represents the optimal stator current reference value finally output by the machine side in the kth sampling period, α mpc (k) and α pid(k) represents the precise weighting coefficients transformed from the fuzzy inference results; i sq_mpc (k|k) represents the optimal control quantity for machine-side model predictive control; i sq_pid_sat (k) is the current reference value after the output of the PID control is limited by the stator current constraint boundary output by S3. Based on the optimal stator current reference value, the corresponding optimal flywheel output power reference value is calculated using the following formula: P f_opt (k)=1.5·p n ·ψ f ·i sq_opt (k)·ω f_filt (k); Among them, P f_opt (k) represents the optimal flywheel output power reference value for the kth sampling period; p n ψ is the number of pole pairs of the flywheel motor; f For the permanent magnet flux linkage of the flywheel motor; ω f_filt (k) is the filtered value of the real-time speed of the flywheel in the kth sampling period; In the reverse interaction, the formula for adaptive adjustment of the quantization factor is: △U dc_max (k)=max(|△U dc (k)|,|△U dc (k-1)|,|△U dc (k-2)|); k u =6 / △U dc_max (k); Among them, △U dc_max (k) represents the maximum absolute value of the DC bus voltage deviation over the past three sampling periods; k u The quantization factor for the DC bus voltage deviation, when ΔU dc_max (k) takes effect when it is not 0, realizing adaptive adjustment of the fuzzy control domain, improving the adjustment accuracy when the deviation is small, and improving the response speed when the deviation is large; The model predictive control generates a power feedforward based on the optimal grid-connected power prediction value, providing feedforward compensation for PID control. The PID control generates a speed safety constraint based on the flywheel speed safety boundary, which in turn corrects the constraint conditions of the model predictive control. The specific details are as follows: The feedforward correction formula is: ; Where, △ω f_ff (k) represents the speed feedforward compensation amount for the kth sampling period; P g_opt (k+1) is the optimal grid-connected power prediction value output by S4; J f(k) represents the latest value of the flywheel moment of inertia obtained from S3 parameter identification; T s Sampling period Based on the feedforward compensation, the input speed deviation of the PID control is corrected. The corrected speed deviation formula is as follows: △ω f_corr (k)=△ω f (k)+△ω f_ff (k); Where, △ω f_corr (k) represents the corrected rotational speed deviation; ω f (k) represents the real-time rotational speed of the flywheel during the kth sampling period; In the reverse interaction, the PID control corrects the constraints of the model predictive control. The constraint correction formula is as follows: ω f_min_safe (k)=max(ω f_min ,ω f_filt (k)+△ω f_pid_min (k)); ω f_max_safe (k)=min(ω f_max ,ω f_filt (k)+△ω f_pid_max (k)); Where, ω f_min_safe (k), ω f_max_safe (k) represents the upper and lower limits of the modified flywheel speed safety constraint; △ω f_pid_min (k), △ω f_pid_max (k) represents the upper and lower safety limits of the rotational speed output by the PID control; ω f_min ω is the minimum permissible operating speed of the flywheel. f_max This refers to the maximum permissible operating speed of the flywheel; Simultaneously, the charging and discharging power constraint boundary calculated in real time on the machine side is fed back into the grid-connected power prediction and collaborative control algorithm to correct the feasible region of the grid-side power reference trajectory; specifically, it includes the following: The corrected formula is: P g_max_corr (k)=min(P g_max ,P f_max (k)·η f ·η g ); P g_min_corr (k)=max(P g_min ,P f_min (k)·η f ·η g ); In the formula, P g_max_corr (k) and P g_min_corr(k) represents the modified maximum and minimum allowable grid-connected power; P g_max P g_min P represents the upper and lower limits of the rated power of the grid-side converter; f_max (k), P f_min (k) represents the upper and lower limits of the flywheel charging and discharging power constraint calculated in real time on the machine side; η g η is the rated efficiency of the grid-side converter. f This refers to the rated efficiency of the machine-side converter; In step S5, the input variables of the fuzzy control are the DC bus voltage deviation and the flywheel speed deviation. The fuzzy subsets of both input variables include: negative large, negative medium, negative small, zero, positive small, positive medium, and positive large. The output variables are the weighting coefficients of the model predictive control and the weighting coefficients of the PID control. The sum of the two weighting coefficients is always 1. Based on a 7×7 combination of two input variables, 49 fuzzy rules were designed. The Mamdani inference method was used to complete the fuzzy inference, and the centroid method was used to complete the defuzzification, so as to obtain the accurate weight coefficients. In step S6, the actual operating data includes the actual filtered values ​​of grid frequency, DC bus voltage and flywheel speed, as well as the actual grid-connected power data output after S4 is executed, and the actual charging and discharging power data output by the flywheel motor after S5 is executed; the actual deviations include: actual deviation of grid frequency, actual deviation of DC bus voltage, actual deviation of flywheel speed, grid-side model prediction error and machine-side model prediction error; The specific rules for updating parameters and weights are as follows: If the actual deviation of the power grid frequency exceeds the preset target deviation and the sampling period is continuously set, the model parameters of the dual-model system are updated using the proportional adjustment method, and the weight domains of the model predictive control and the PID control are adjusted at the same time. The methods for updating network-side model parameters include: 1) Increase the virtual inertia and virtual damping coefficient of the virtual synchronous generator to enhance the system's inertial support capability, suppress rapid frequency fluctuations, and accelerate frequency deviation convergence; 2) Increase the output tracking weight in the model predictive control cost function to improve the tracking accuracy of the power reference trajectory and accelerate the convergence speed of the frequency deviation; 3) Increase the bandwidth of the disturbance observer to improve the response speed of disturbance observation and quickly correct model errors caused by large disturbances on the power grid side; Methods for updating aircraft-side model parameters include: 1) Reduce the forgetting factor of recursive least squares method, improve the tracking speed of parameter identification for model parameter changes, and quickly adapt to model parameter fluctuations caused by changes in working conditions. 2) The proportional adjustment method is used to directly correct the core parameters of the machine-side model and quickly compensate for the model prediction error; The weighted universe of discourse method for machine-side fuzzy control includes: 1) When the power grid frequency deviation is large and it is necessary to prioritize ensuring the power grid frequency support performance; α mpc_max (k+1)=min(1.0,α mpc_max (k)+k α_up ·E(k)); α pid_max (k+1)=max(0.0,α pid_max (k)-k α_down ·E(k)); Where, α mpc_max (k+1) represents the upper bound of the universe of discourse for the updated model prediction control weight coefficients; α pid_max (k+1) represents the upper bound of the universe of discourse for the updated PID control weight coefficients; k α_up and k α_down The weighted universe is used to adjust the coefficients; α is always satisfied. mpc_max (k)+α pid_max The constraint is (k)=1; E(k) is the power grid frequency deviation evaluation index. E(k) = |Δf g_filt (k)|,△f g_filt (k) represents the power grid frequency deviation filter value in the kth sampling period; 2) When the flywheel speed deviation is large and close to the safety boundary, and it is necessary to prioritize equipment safety: α pid_max (k+1)=min(1.0,α pid_max (k)+k α_pid_up ·E ω (k)); α mpc_max (k+1)=max(0.0,α mpc_max (k)-k α_mpc_down ·E ω (k)); Where, k α_pid_up k α_mpc_down E is the weighting universe adjustment coefficient, ranging from 0.1 to 0.3, always satisfying the constraint that the sum of the weighting coefficients is 1; ω (k) represents the flywheel speed tracking error index. E ω (k)=|ω f_ref (k)-ω f_filt (k)|; Where, ω f_ref (k) is the reference value for flywheel speed; ω f_filt(k) is the filtered value of the flywheel real-time speed in the kth sampling period output by S2; If the actual deviation of the power grid frequency is within the preset target deviation range, only the model parameters are fine-tuned or the parameters are kept unchanged; The termination condition for the closed-loop regulation is set as follows: within 100 consecutive acquisition cycles, the actual deviation of the power grid frequency is within ±0.05 Hz and the actual deviation of the power grid frequency does not exceed the target deviation of 0.02 Hz. After the termination condition is met, the adaptive update of parameters and weights is stopped, the current control strategy is maintained, and the closed-loop regulation is restarted until the power grid frequency deviates from the rated range again.

[0007] The flywheel energy storage grid frequency support method based on model predictive control provided in this application has at least the following beneficial effects: Significantly improves the accuracy and speed of grid frequency support: By integrating the global optimization advantages of model predictive control and the inertia support advantages of virtual synchronous generators, combined with the real-time error correction capability of disturbance observation, it effectively suppresses large and rapid fluctuations in grid frequency, enhances the grid's inertia and damping characteristics, and achieves precise and rapid regulation of grid frequency. Achieves global coordinated control between the grid side and the generator side: Constructing a dual-model system with the DC bus as the intermediate link, and through bidirectional power feedback constraints between the grid side and the generator side, it ensures a precise match between grid-connected power demand and flywheel charging and discharging capabilities, effectively suppressing DC bus voltage fluctuations and improving the operational stability of the entire flywheel energy storage system after grid connection. Fully balances grid frequency support performance and flywheel equipment safety: Through fuzzy control, it achieves adaptive weight allocation between model predictive control and PID control, prioritizing frequency support when grid frequency deviation is large, and prioritizing equipment safety protection when flywheel speed approaches the safety boundary, avoiding problems such as flywheel overcharging and over-discharging, and mechanical damage, achieving an optimized balance of multiple control objectives. Significantly enhances the system's anti-disturbance capability and all-condition adaptive performance: Constructs an adaptive closed-loop update mechanism for model parameters and algorithm weights based on control effect evaluation. This mechanism can adjust relevant parameters and weights in real time according to grid disturbances and flywheel operating conditions, ensuring the system always maintains optimal control and adapts to different grid operating conditions and flywheel energy storage states. Improves model prediction accuracy and control command reliability: Through full-process data preprocessing including outlier detection and correction, smoothing and denoising, and validity verification, outliers, random noise, and invalid data in the original collected parameters are effectively eliminated. This provides a reliable data foundation for model construction and control algorithms, improving model prediction accuracy and control command accuracy from the source. Achieves full closed-loop stable control of grid frequency: Designs a clear closed-loop regulation termination and restart mechanism. After the grid frequency stabilizes, parameter updates are stopped and a steady-state control strategy is maintained. When the grid frequency deviates from the rated range again, the adaptive regulation process is quickly restarted, ensuring long-term grid frequency stability while reducing system computational burden and improving control efficiency. Enhance the integration and complementarity of control algorithms: Break away from the traditional control architecture of multiple algorithms in a unidirectional series, design a pairwise bidirectional closed-loop interaction mechanism for all integrated algorithms, so that each algorithm can provide feedback and correction to each other, give full play to the technical advantages of each algorithm, make up for the inherent defects of a single algorithm, and achieve a control effect of 1+1>2. Attached Figure Description

[0008] Figure 1 This is a schematic diagram illustrating the steps of a model predictive control-based flywheel energy storage system's grid frequency support method.

[0009] Figure 2 This is a flowchart illustrating the grid frequency support method for flywheel energy storage systems based on model predictive control.

[0010] Figure 3This is a schematic diagram of the execution flow of the grid-connected power prediction and collaborative control algorithm.

[0011] Figure 4 This is a schematic diagram of the execution flow of the engine-side flywheel output power prediction and collaborative control algorithm. Detailed Implementation

[0012] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0013] like Figures 1-2 As shown, this embodiment provides a grid frequency support method for flywheel energy storage systems based on model predictive control. With model predictive control as the core, it integrates algorithms such as disturbance observation, virtual synchronous generator, PID, and fuzzy control to construct a collaborative control strategy with a grid-side-machine-side dual-model system and multi-algorithm bidirectional closed-loop interaction. It also designs a full-process parameter adaptive closed-loop adjustment mechanism to achieve accurate grid frequency support for the flywheel energy storage system.

[0014] S1. Synchronous acquisition of grid operating parameters and flywheel energy storage system state parameters; This step is the data foundation of the entire frequency support method, and its core function is to provide accurate and synchronous raw operating data for all subsequent control links. This step completes the acquisition of all parameters at fixed time intervals. This time interval is called the fixed sampling period of the system discrete control, denoted by T. s Indicated, the unit is seconds. In this scheme, T s The value is set to 100 microseconds. This setting ensures real-time control and enables rapid response to fluctuations in the power grid frequency, while avoiding the computational burden and hardware cost pressure caused by excessively high sampling frequencies. All raw parameters collected and calculated in this step are used as inputs for the subsequent preprocessing step S2 to ensure the continuity and consistency of the data link.

[0015] First, the power grid operating parameters are collected and calculated. The real-time frequency of the power grid is the actual frequency at the grid connection point at the current moment, denoted by the symbol f. g (k) indicates that k is the index of the sampling period, representing the k-th sampling period, in Hertz. The real-time frequency of the power grid is synchronously acquired through a power grid frequency sensor or a software phase-locked loop. Both acquisition methods can be flexibly selected according to the on-site hardware configuration and both can meet the accuracy requirements. The rated frequency of the power grid is the nominal standard operating frequency of the power grid, denoted by the symbol f. N This indicates that the value is 50 Hz, which is the standard rated frequency of the power grid. No additional data acquisition is required; it is directly preset as a fixed parameter in the system. Based on the acquired real-time power grid frequency, the frequency deviation is calculated. This frequency deviation reflects the degree to which the current power grid frequency deviates from the rated value and is the core basis for triggering frequency support actions, denoted by the symbol Δf. g (k) represents the unit Hertz, and the calculation formula is: △fg (k)=f g (k)-f N ; Among them, f g (k) is the real-time frequency of the power grid in the k-th sampling period, f N This is the rated frequency of the power grid.

[0016] Next, we calculate the grid frequency change rate, which reflects how quickly the grid frequency changes and is an important indicator for judging the severity of grid disturbances. It is denoted by ROCOF(k) and its unit is Hertz per second. It is calculated using a first-order difference, and the formula is as follows: ROCOF(k)=(f g (k)-f g (k-1)) / T s ; In the formula, f g (k-1) is the real-time frequency of the power grid in the (k-1)th sampling period, T s This is the fixed sampling period for the discrete control of the system.

[0017] The grid voltage amplitude is the voltage amplitude at the grid connection point at the current moment, denoted by the symbol U. g (k) represents the voltage in volts, acquired by a voltage sensor. The grid-connected current of the energy storage system is the instantaneous value of the phase current on the grid-connected side of the flywheel energy storage system at the current moment, denoted by the symbol i. g (k) indicates that the unit is ampere, which is obtained by a Hall current sensor.

[0018] Simultaneously, the state parameters of the flywheel energy storage system are collected and calculated. The real-time speed of the flywheel is the mechanical angular velocity of the flywheel rotor at the current moment, denoted by the symbol ω. f (k) represents the unit radians per second, acquired through a photoelectric speed sensor or a sensorless observation algorithm of the motor. The two acquisition methods can be flexibly selected based on the type of flywheel motor and hardware configuration. The rated speed of the flywheel is the nominal standard operating angular velocity of the flywheel rotor, denoted by the symbol ω. fN The value is expressed in radians per second and is determined by the flywheel device's factory parameters, directly preset as a fixed parameter in the system. Based on the collected real-time flywheel speed, a per-unit speed value is calculated. This per-unit speed value visually reflects the ratio of the flywheel speed to the rated speed and is an auxiliary indicator for judging the flywheel's energy storage capacity status, denoted by the symbol λ. ω (k) represents the calculation formula: λ ω (k)=(ω f (k)) / ω fN ; In the formula, ω f (k) represents the real-time rotational speed of the flywheel during the kth sampling period, ωfN This is the rated speed of the flywheel.

[0019] The actual value of flywheel charging and discharging power is the actual charging and discharging power of the flywheel motor port at the current moment, denoted by the symbol P. f (k) represents the power in watts, acquired by a power sensor. This scheme defines charging as positive and discharging as negative; this convention is applied throughout the entire scheme to ensure consistent power flow direction determination. The flywheel motor stator current is the instantaneous value of the phase current of the flywheel motor stator at the current moment, denoted by the symbol i. s (k) indicates that the unit is ampere, which is obtained by a Hall current sensor.

[0020] The mechanical power loss of a flywheel is the sum of the power losses of the flywheel rotor due to mechanical friction and wind resistance at the current moment, denoted by the symbol P. loss (k) indicates that the unit is watts, calculated using a preset loss model without the need for additional hardware data acquisition. The calculation formula for the loss model is: P loss (k)=B ω ·ω f (k) 2 ; In the formula, B ω This is the mechanical damping coefficient of the flywheel, measured in Newton-meter-seconds per radian. This value is determined by the factory parameters of the flywheel equipment and is directly preset as a fixed parameter in the system.

[0021] Flywheel energy storage capacity is the state of charge (SOC) of the flywheel energy storage system at the current moment. f (k) indicates that the unit is per unit, calculated based on the flywheel speed and moment of inertia, and reflects the amount of energy currently stored in the flywheel. The calculation formula is: SOC f (k)=△(k) / E N =(J f ·ω f (k) 2 ·1 / 2) / ( J f ·ω fN 2 ·1 / 2)=((ω f (k)) / ω fN ) 2 =λ ω (k) 2 ; In the formula, ω f (k) represents the real-time rotational speed of the flywheel during the kth sampling period, ω fN λ is the rated speed of the flywheel. ω (k) is the per-unit value of the rotational speed; J fThe moment of inertia of the flywheel rotor, expressed in kilograms per square meter, is determined by the factory specifications of the flywheel equipment and is preset as a fixed parameter in the system. The absolute kinetic energy currently stored in the flywheel: Δ(k) = J f ·ω f (k) 2 · 1 / 2; Maximum rated kinetic energy of the flywheel: E N =J f ·ω fN 2 ·1 / 2.

[0022] The DC bus voltage is the instantaneous voltage across the DC bus capacitor of the flywheel energy storage system at the current moment, denoted by the symbol U. dc (k) represents the voltage in volts, acquired by a voltage sensor. The rated voltage of the DC bus is the nominal standard operating voltage of the DC bus, denoted by the symbol U. dcN This indicates that the unit is volts. This value is determined by the system hardware design and is directly preset in the system as a fixed parameter.

[0023] S2. Preprocessing of Acquired Parameters: The core function of this step is to purify the raw parameters acquired in S1, eliminating outliers, random noise, and invalid data to provide accurate and reliable valid data for subsequent model building and control algorithms. This step sequentially executes three sub-steps: outlier detection and correction, smoothing and denoising, and validity verification. These three sub-steps are interconnected, with the output of the previous sub-step serving as the input for the next. The final output of all valid data serves as the input data source for all subsequent steps S3 to S6, ensuring the accuracy of subsequent model calculations and control outputs.

[0024] The first step is outlier detection and correction, using the 3σ criterion to identify and correct outliers for all parameters. Outliers may appear in the raw data due to sensor malfunctions, electromagnetic interference, etc. If not corrected, they will severely affect the prediction accuracy of subsequent models and the stability of control algorithms. For each acquired raw parameter sequence, the raw data from the first 20 sampling periods are taken as a statistical sample. The sample mean and sample standard deviation are calculated. Using the sample mean as the center and three times the sample standard deviation as the boundary, it is determined whether the raw value of the current sampling period is an outlier. The real-time frequency f of the power grid is used as the reference. g For example, in (i), the sample mean is represented by the symbol μ. f The calculation formula is as follows: ; In the formula, i is the index of the sample data, k is the index of the current sampling period, and f g(i) represents the original sample data used in the calculation. In this embodiment, k is set to 20. Based on the principle of fixed-window moving average, the value of i ranges from k-19 to k, achieving real-time online rolling average to meet the engineering requirements of real-time control systems such as flywheel energy storage. Flywheel energy storage is an online real-time system that requires rapid calculation of the average value within each control cycle (e.g., 1ms / 10ms). The sliding window has a fixed computational load, eliminating the need to accumulate all historical data and perfectly suited to the real-time computing requirements of embedded controllers. During flywheel operation, parameters may drift (e.g., bearing wear, changes in vacuum, increased losses). The sliding window uses only the most recent 20 data points, enabling it to track dynamic system changes and reflect recent operating conditions. In contrast, the full average is affected by the initial old data, which "submerges" the influence of new data, making it unable to track parameter changes. Real-time control only cares about the state near the current moment, not the global state from startup to the present; using the most recent fixed window data is more in line with engineering practice. As the full average increases with k, the weight of new data decreases (e.g., when k=1000, new data only accounts for 1 / 1000 of the weight), eventually causing the average value to hardly change and losing its filtering significance. The sliding window, however, has a fixed weight of 1 / 20 for new data, always reflecting recent changes.

[0025] The standard deviation of a sample is represented by the symbol σ. f The calculation formula is as follows: ; The raw values ​​of the current sampling period are evaluated. If the absolute value of the difference between the current value and the sample mean is greater than three times the sample standard deviation, it is considered an outlier, and the filtered value from the previous period is used to replace the current outlier. If the value does not exceed the threshold, the original value is retained. All collected raw parameters undergo outlier detection and correction according to the above 3σ criterion to obtain the corrected parameter sequence.

[0026] The second step is to perform smoothing and denoising. A moving average filter is used to smooth the corrected parameters, with a filter window length of 5. Even after outlier correction, the original data may still contain random noise. The moving average filter, by averaging data over multiple consecutive sampling periods, effectively suppresses random noise while preserving the data's trend. Taking the real-time grid frequency as an example, the real-time grid frequency value after outlier correction is denoted by the symbol f. g_corr (i) indicates that the filtered real-time frequency value of the power grid obtained after filtering is represented by the symbol f. g_filt (k) represents the calculation formula: ; In the formula, i is the index of the data within the filtering window, ranging from k-4 to k, and in this implementation, k is set to 5. All parameters that have undergone outlier correction are smoothed and denoised using the above moving average filtering to obtain the filtered values ​​of the corresponding parameters.

[0027] The third step is to perform validity verification. For each parameter's filter value, a physically reasonable range is set, and the validity is verified one by one. Invalid data that exceeds the valid range is replaced with the valid filter value from the previous cycle. Finally, all valid data that has been verified is output. The physically reasonable range is based on the safe operation requirements of the power grid and flywheel energy storage system, as well as the rated parameter settings of the equipment, to ensure that the valid data conforms to the actual physical meaning and avoids malfunctions in the control algorithm due to invalid data. The valid ranges of all parameters in this embodiment are as follows: 1. The effective range of the real-time frequency filter value for the power grid is 49.5 Hz to 50.5 Hz; 2. The effective range of the power grid frequency deviation filter value is -0.5 Hz to 0.5 Hz; 3. The effective range of the power grid frequency change rate filter value is -2 Hz to 2 Hz. 4. The effective range of the grid voltage amplitude filter value is 0.9 times the grid rated voltage amplitude to 1.1 times the grid rated voltage amplitude; 5. The effective range of the grid-connected current filtering value of the energy storage system is from -1 times the rated current amplitude of the grid-connected converter to 1 times the rated current amplitude of the grid-connected converter; 6. The effective range of the per-unit value of the flywheel real-time speed filter value is 0.33 to 1.1; 7. The effective range of the actual value of the flywheel charging and discharging power filter is from -1 times the rated charging and discharging power of the flywheel energy storage system to 1 times the rated charging and discharging power of the flywheel energy storage system; 8. The effective range of the flywheel motor stator current filter value is -1 times the rated stator current amplitude of the flywheel motor to 1 times the rated stator current amplitude of the flywheel motor; 9. The effective range of the flywheel energy storage capacity filter value is 0.1 to 1.2; 10. The effective range of the DC bus voltage filter value is 0.8 times the rated DC bus voltage to 1.2 times the rated DC bus voltage.

[0028] The mechanical power loss of the flywheel is a calculated quantity, not a directly collected quantity, and its calculation parameters have been constrained, meaning it is indirectly constrained. Furthermore, mechanical loss is a passive physical quantity that only reflects wind resistance plus friction loss, and will not cause overvoltage, overcurrent, or overspeed of the equipment, so no independent safety threshold is required.

[0029] S3. Construction and Parameter Identification of a Dual-Model System with the DC Bus as the Intermediate Link: The core function of this step is to construct a mathematical model system that supports subsequent control algorithms, providing an accurate prediction basis for model predictive control. This scheme adopts a dual-model system with the DC bus as the intermediate link, including a grid-side power prediction model and a flywheel energy storage machine-side dynamic model. The two models are deeply correlated through the DC bus power balance equation, ensuring the matching between the grid-side power demand and the flywheel charging and discharging power on the machine side. The grid-side prediction model (mathematical model / state-space equation), total disturbance observations, and constraint boundaries output in this step provide inputs for the S4 grid-side collaborative algorithm, while the machine-side dynamic model and power and current constraint boundaries provide inputs for the S5 machine-side collaborative algorithm. The model parameters can be updated through the closed-loop feedback in S6, achieving adaptive optimization of the model.

[0030] First, the core constraint of the dual-model system is clarified, namely the DC bus power balance equation. This equation is the core link between the grid-side and machine-side models, describing the relationship between the power flowing from the grid-side converter to the DC bus, the power flowing from the DC bus to the machine-side converter, the DC bus power loss, and the change in DC bus capacitor energy storage. The calculation formula is as follows: P dc_g (k)=P dc_f (k)+P dc_loss (k)+C dc ·U dc_filt (k)·(dU dc ) / dt; In the formula, P dc_g (k) represents the power flowing from the grid-side converter to the DC bus during the k-th sampling period, in watts; P dc_f (k) represents the power flowing from the DC bus to the machine-side converter during the k-th sampling period, in watts; P dc_loss (k) represents the DC bus line and switch loss power in the kth sampling period, in watts, calculated using a preset loss model, requiring no additional data acquisition; C dc This refers to the capacitance value of the DC bus capacitor, measured in farads. It is determined by the system hardware design and is preset as a fixed parameter in the system. dc_filt (k) is the DC bus voltage filter value for the kth sampling period, obtained from the output of step S2; U dc The DC bus voltage is given by dU, where t is time. dc ) / dt is the rate of change of the DC bus voltage.

[0031] The first part completes the construction of the grid-side power prediction model. The grid-side power prediction model integrates the second-order linear model of the grid-connected filter, the second-order disturbance observer, and the second-order model of the virtual synchronous generator. The input is the effective grid-side data output by S2 (the grid operation parameters collected by S1, and the filtered values ​​obtained after preprocessing by S2). The outputs are the grid-connected power prediction reference trajectory, the grid-side discrete prediction model, the total disturbance observation value, and the grid-connected converter constraint boundary.

[0032] First, a second-order virtual synchronous generator model is constructed. The core of this model is the simulation of the generator's rotor motion equations, providing virtual inertia and damping support for the power grid. This model is the core of the power grid frequency support. In traditional power systems, the rotor inertia of a synchronous generator can suppress rapid changes in the power grid frequency. The virtual synchronous generator simulates this characteristic through algorithms, improving the frequency stability of power systems with a high proportion of renewable energy. The rotor motion equations in the continuous domain are: J vsg ·(dω vsg ) / dt=T m -T e -D vsg ·(ω vsg -ω N ); J vsg ω represents the virtual moment of inertia of the virtual synchronous generator, expressed in kilograms per square meter. This value is an adjustable parameter used to control the magnitude of the virtual inertia. vsg The virtual electric angular velocity of the virtual synchronous generator, in radians per second, (dω) vsg ω / dt is the rate of change of the virtual electric angular velocity; N The rated electric angular velocity of the power grid is calculated using the formula ω. N =2π·f N The unit is radians per second; T m The virtual mechanical torque of the virtual synchronous generator is expressed in Newton-meters (N·m); T e The virtual electromagnetic torque of the virtual synchronous generator is expressed in Newton-meters (N·m); D vsg This is the virtual damping coefficient of the virtual synchronous generator, measured in Newton-meter-seconds per radian. This value is an adjustable parameter used to control the magnitude of the virtual damping.

[0033] Near the rated angular velocity, the torque equation is converted into a power form to obtain the rotor motion equation in the power domain, which is easier to integrate with the power control system. The formula is as follows: J vsg ·ω N ·(d△ω vsg ) / dt=P m_ref -P e -D vsg ·ω N ·△ωvsg ; In the formula, △ω vsg The virtual angular velocity deviation of the virtual synchronous generator is calculated using the formula Δω. vsg =ω vsg -ω N The unit is radians per second; P m_ref This is a virtual mechanical power reference value for a virtual synchronous generator, in watts; P e This refers to the virtual electromagnetic power of the virtual synchronous generator, i.e., the active power output to the grid, measured in watts.

[0034] The virtual mechanical power reference value is generated from the active power-frequency droop characteristic. The input is the grid frequency deviation filter value output in step S2, realizing the function of automatically adjusting the grid-connected power according to the grid frequency deviation. The formula is: P m_ref =P0+k p ·△f g_filt (k); In the formula, P0 is the reference value of the grid-connected steady-state active power of the flywheel energy storage system, in watts. Under steady-state conditions, its value is 0 to ensure that the flywheel energy storage system does not participate in active power regulation when the grid frequency is stable; k p The active power-frequency droop factor, measured in watts per hertz, is determined by the voltage regulation accuracy requirements supported by the grid frequency. A larger value results in greater power regulation for the same frequency deviation. △f g_filt (k) is the power grid frequency deviation filter value for the kth sampling period, which is obtained from the output of S2.

[0035] The virtual synchronous generator model outputs grid-connected power reference values ​​and virtual electric angular velocity, which serve as inputs to the grid-connected filter model and as reference targets for grid-side model predictive control.

[0036] Then, a second-order linear model of the grid-connected filter is constructed. This scheme uses an LCL filter on the grid-connected side, which has better high-frequency harmonic suppression capability and can meet the harmonic requirements of the grid-connected current. A second-order linear state-space model is constructed to describe the mapping relationship between the control input and the grid-connected output power, grid-connected point voltage, and current. The state-space equations in the continuous domain are: (dX g ) / dt=A g ·X g +B g ·U g +D g ·d g ; Y g =C g ·X g ; In the formula, the state variable X gThis is a column vector containing the grid-side converter output current, the grid-connected current, and the LCL filter capacitor voltage; all units are international standard units. A g The network-side state matrix is ​​determined by the parameters of the LCL filter; B g The input matrix on the network side is determined by the parameters of the LCL filter; C g The output matrix on the network side is determined by the physical quantities to be output; D g The network-side perturbation matrix is ​​determined by the parameters of the LCL filter; U g This refers to the output voltage control quantity of the grid-side converter, measured in volts; d g For the grid voltage disturbance term, the value is the filtered value of the grid voltage amplitude output from S2; Y g The output of the model includes grid-connected current, grid-connected active power, and grid-connected reactive power.

[0037] Using the zero-order hold method, with a sampling period T s To discretize the step size, the continuous-domain state-space equation is transformed into a discrete-domain form, resulting in a discretized network-side filter model, which is easier to implement in a digital controller. The formula is as follows: X g (k+1)=A gd ·X g (k)+B gd ·U g (k)+D gd ·d g (k); Y g (k)=C gd ·X g (k); In the formula, A gd B gd C gd D gd These are the discretized state matrix, input matrix, output matrix, and perturbation matrix, respectively; X g (k) represents the network-side state variable during the k-th sampling period; U g (k) represents the output voltage control quantity of the grid-side converter in the kth sampling period; d g (k) represents the grid voltage disturbance term in the kth sampling period; X g (k+1) represents the predicted value of the network-side state variable in the (k+1)th sampling period; Y g (k) represents the model output during the k-th sampling period. This model describes the current control input (the grid-side converter output voltage control quantity U). g The mapping relationship between (k) and the grid-connected state variables and output quantities in the next cycle provides a predictive basis for grid-side model predictive control.

[0038] A second-order disturbance observer is then constructed to observe the total internal and external disturbances of the grid-side model, including all disturbance terms affecting model accuracy such as model parameter mismatch, grid voltage fluctuations, and nonlinear disturbances, thereby improving the accuracy of the prediction model. The accuracy of model predictive control is highly dependent on the accuracy of the prediction model. Various disturbances existing in actual operation will cause deviations between the model prediction values ​​and the actual values. The disturbance observer can observe these disturbances in real time and compensate for them in the prediction model.

[0039] First, the discretized network-side model is rewritten in a form containing the total perturbation, denoted as the network-side discrete prediction model: X g (k+1)=A gd ·X g (k)+B gd ·U g (k)+d total_g (k); In the formula, d total_g (k) represents the total internal and external disturbances of the network-side model during the kth sampling period.

[0040] Construct a second-order linear extended state observer, expanding the total disturbance into new state variables, and define the extended state vector Z. g (k) is a column vector containing the observed values ​​of the state variables z. 1g First-order difference observations of state variables z 2g Total disturbance observation value z 3g The iterative formula for a second-order extended state observer in the discrete domain is: z 1g (k+1)=z 1g (k)+T s ·(z 2g (k)+β 1g ·(Y g_filt (k)-z 1g (k))); z 2g (k+1)=z 2g (k)+T s ·(z 3g (k)+β 2g ·(Y g_filt (k)-z 1g (k))+A gd ·X g (k)+B gd ·U g (k)); z 3g (k+1)=z 3g (k)+T s ·(β 3g ·(Y g_filt (k)-z1g (k))); In the formula, z 1g (k) represents the observed state variable value in the kth sampling period; z 2g (k) represents the first-order difference observation of the state variable in the kth sampling period; z 3g (k) represents the total disturbance observation value in the kth sampling period; β 1g β 2g β 3g The gain coefficient of the observer is designed based on the observer bandwidth, which is 3 to 5 times the network-side control bandwidth to ensure that the observer's response speed is faster than the control loop, while avoiding the introduction of excessive noise; Y g_filt (k) is the filtered value of the grid-connected output in the kth sampling period. The filtered value of the grid-connected output includes the filtered values ​​of the grid-connected current, grid-connected active power and grid-connected reactive power after preprocessing by S2.

[0041] The second-order disturbance observer outputs the total disturbance observation value d. total_g_obs (k)=z 3g (k) is used to correct the network-side prediction model in real time and eliminate prediction errors caused by disturbances.

[0042] The second part completes the construction and parameter identification of the dynamic model of the flywheel energy storage machine side. The dynamic model of the machine side is a second-order linear model, which describes the mapping relationship between the control input (reference value of the q-axis current of the flywheel motor stator) and the flywheel output power, DC bus voltage and flywheel speed. The input is the effective data of the flywheel side output by S2 (preprocessed state parameters of the flywheel energy storage system). The output is the flywheel output power prediction model, the charging and discharging power constraint boundary and the current constraint boundary and the machine side discrete prediction model.

[0043] First, a dynamic model of the flywheel motor is constructed. In this scheme, the flywheel is driven by a permanent magnet synchronous motor (PMSM), which has advantages such as high efficiency and high power density, making it suitable for flywheel energy storage system applications. A second-order linear model is then constructed in the dq rotating coordinate system, with the core components being the electromagnetic torque equation and the mechanical motion equation. The electromagnetic torque equation describes the relationship between the stator current and the electromagnetic torque, and the formula is: T ef (k)=1.5·p n ·ψ f ·i sq_filt (k); In the formula, T ef (k) represents the electromagnetic torque of the flywheel motor in the kth sampling period, in Newton-meters; p n The number of pole pairs of the flywheel motor is determined by the motor's factory parameters and is preset as a fixed parameter in the system; ψ fThe flux linkage of the permanent magnet in the flywheel motor is measured in Weber, determined by the motor's factory specifications and preset as a fixed parameter in the system; i sq_filt (k) is the filtered value of the q-axis component of the flywheel motor stator current in the kth sampling period, in amperes, which is obtained by coordinate transformation from the filtered value of the flywheel motor stator current output by S2.

[0044] The mechanical equations of motion for the flywheel are the core second-order equations of the machine-side model, describing the relationship between electromagnetic torque, load torque, and flywheel speed. Their continuous-domain expression is: J f ·(dω f ) / dt=T ef (k)-B ω ·ω f (k)-T load (k); In the formula, J f T represents the moment of inertia of the flywheel rotor, expressed in kilograms per square meter. This value is a parameter to be identified and will be updated online in subsequent parameter identification stages. load (k) represents the flywheel load torque in the kth sampling period, in Newton-meters, which is determined by mechanical losses under steady-state conditions.

[0045] With sampling period T s Using a step size, the mechanical motion equations are discretized to obtain a discrete-domain machine-side speed prediction model, which is convenient for implementation in a digital controller. The formula is: ω f (k+1)=ω f (k)+(T s / J f )·(T ef (k)-B ω ·ω f (k)-T load (k)); In the formula, ω f (k+1) is the predicted flywheel speed for the (k+1)th sampling period; ω f (k) represents the real-time rotational speed of the flywheel during the kth sampling period; B ω is the flywheel mechanical damping coefficient.

[0046] Based on the DC bus power balance equation, a generator-side power prediction model is constructed. The input power of the generator-side converter comes from the DC bus, and the output power is the electromagnetic power of the flywheel motor. Considering the converter efficiency, the power conversion relationship is obtained: P dc_f (k)=(P ef (k)) / η f +P loss_inv (k); In the formula, Pef (k) represents the electromagnetic power of the flywheel motor in the kth sampling period, calculated using the formula: P ef (k)=T ef (k)·ω f (k), the unit is watts; η f The rated efficiency of the generator-side converter is determined by the equipment's factory parameters and is directly preset as a fixed parameter in the system; P loss_inv (k) represents the switching and conduction losses of the machine-side converter in the kth sampling period, in watts, which are calculated using a preset loss model and do not require additional data acquisition.

[0047] Based on the above equations, a complete prediction model for the generator side is obtained: the flywheel energy storage generator side dynamic model, which describes the mapping relationship between the current control input and the flywheel speed, flywheel output power, and DC bus voltage in the next cycle. ω f_pred (k+1)=ω f_filt (k)+(T s / J f )·(1.5·p n ·ψ f ·i sq_ref (k)-B ω ·ω f_filt (k)-T load_filt (k)); P f_pred (k+1)=1.5·p n ·ψ f ·i sq_ref (k)·ω f_pred (k+1); U dc_pred (k+1)=U dc_filt (k)+(T s / C dc )·(P dc_g (k)-P dc_f (k)-P dc_loss (k)) / U dc_filt (k); In the formula, ω f_pred (k+1) is the predicted flywheel speed for the (k+1)th sampling period; ω f_filt (k) is the filtered value of the flywheel real-time speed in the kth sampling period, obtained from the output of S2; i sq_ref (k) represents the reference value of the flywheel motor stator q-axis current in the kth sampling period, i.e., the machine-side control input, which serves as the target command output by the controller for closed-loop tracking with the actual q-axis current; T load_filt (k) represents the flywheel load torque filter value in the kth sampling period, and the flywheel mechanical loss power filter value P obtained from the output of S2.loss_filt (k) and flywheel real-time speed filter value ω f_filt (k) Calculation yields: T load_filt (k)=P loss_filt (k) / ω f_filt (k); P f_pred (k+1) is the predicted flywheel output power value for the (k+1)th sampling period; U dc_pred (k+1) is the predicted value of the DC bus voltage in the (k+1)th sampling period; U dc_filt (k) is the DC bus voltage filter value for the kth sampling period, obtained from the output of S2; C dc P is the capacitance value of the DC bus capacitor; dc_g (k) represents the power flowing from the grid-side converter to the DC bus during the k-th sampling period; P dc_loss (k) represents the DC bus line and switch loss power in the kth sampling period.

[0048] Then, the parameters of the dynamic model on the machine side are identified. The parameters to be identified include: flywheel moment of inertia, mechanical damping coefficient, and load torque coefficient. The recursive least squares method is used for online parameter identification. The input is the effective data output by S2 in the first 20 sampling periods. The objective function is to minimize the sum of squared errors between the actual output power of the flywheel and the model prediction. The objective function expression is: ; In the formula, J ident The objective function identified for the parameters; P f_filt (i) is the filtered value of the actual flywheel charging and discharging power in the i-th sampling period, obtained from the output of S2; P f_pred (i) is the predicted output power of the flywheel in the i-th sampling period; in this embodiment, k is 20; i is the index of the sample data, ranging from k-19 to k. The power data of the most recent 20 sampling periods are used for model matching to achieve real-time online rolling calculation, which is suitable for the engineering requirements of real-time control systems such as flywheel energy storage.

[0049] First, the discretized mechanical motion equations are rewritten in linear regression form to facilitate the application of the recursive least squares method: y(k)= T (k)·θ(k); Observed output y(k)=ω f_filt (k)-ω f_filt (k-1); Regression vector (k) is a column vector containing electromagnetic torque, rotational speed, and load torque; the parameter vector to be identified, θ(k), is a column vector containing moment of inertia, damping coefficient, and load torque coefficient; ωf_filt (k) is the filtered value of the real-time speed of the flywheel in the kth sampling period.

[0050] The iterative formula for the recursive least squares method is: K(k)=(P(k-1)· (k)) / (λ+ T (k)·P(k-1)· (k)); P(k)=(1 / λ)· (IK(k)· T (k))·P(k-1); θ(k) = θ(k-1) + K(k)·(y(k)- T (k)·θ(k-1)); In the formula, K(k) is the gain matrix of the kth sampling period; P(k) is the covariance matrix of the kth sampling period; λ is the forgetting factor, with a value of 0.98, which is used to weaken the influence of historical data and improve the ability to track time-varying parameters; I is the identity matrix; y(k) is the observation output of the kth sampling period. θ(k) is the regression vector for the kth sampling period; θ(k) is the vector of parameters to be identified for the kth sampling period.

[0051] The parameter identification process is iteratively updated every 20 sampling periods. The updated parameters are directly substituted into the machine-side dynamic model to ensure the accuracy of the model's predictions. During actual system operation, online parameter identification updates the parameter estimates time-by-time without storing all historical data or recalculating. The recursive least squares method is used to update the parameter estimates with new data each time, avoiding repeated processing of all historical data. This method is computationally efficient and more suitable for real-time control systems.

[0052] Finally, based on the effective flywheel data output by S2, combined with the equipment's rated parameters and safe operating range, the charging and discharging power constraint boundaries and current constraint boundaries are calculated and output. These constraint boundaries are used for subsequent optimization of the control algorithm to ensure that control commands remain within the equipment's safe operating range and to prevent equipment damage. The constraint boundaries are as follows: 1. Flywheel charging and discharging power constraint boundary: P f_min (k)≤P f (k)≤P f_max (k), where P f_max (k) represents the maximum allowable charging power in the kth sampling period, which is determined by the rated power of the flywheel motor, the current flywheel energy storage capacity filter value, and the upper limit of the DC bus voltage; P f_min(k) represents the maximum allowable discharge power in the kth sampling period, which is determined by the rated power of the flywheel motor, the current real-time speed filter value of the flywheel, and the lower limit of the DC bus voltage. 2. Flywheel motor stator current constraint boundary: i s_min ≤i sq_ref (k)≤i s_max i s_max i represents the positive upper limit of the q-axis current of the motor stator, corresponding to the maximum charging torque; s_min The negative lower limit of the q-axis current of the motor stator corresponds to the maximum discharge torque, both of which are determined by the motor's rated current and demagnetization characteristics. 3. DC bus voltage constraint boundary: U dc_min ≤U dc (k)≤U dc _ max U dc_min U is the minimum permissible operating voltage for the DC bus. dc_max This is the highest permissible operating voltage of the DC bus, determined by the rated parameters of the DC bus capacitors and the withstand voltage rating of the converter components.

[0053] S4. Implementation of Grid-Side Power Prediction and Coordinated Control Algorithm; The core function of this step is to generate the optimal control quantity for the grid-side converter, accurately respond to grid frequency deviation requirements, and achieve active support for grid frequency. (Reference) Figure 3 This step deeply integrates three types of algorithms: model predictive control, disturbance observation, and virtual synchronous generator control. It breaks away from the traditional unidirectional series control architecture, achieving two-way, two-way closed-loop interaction between the three algorithms. This fully leverages the inertial support advantages of the virtual synchronous generator, the multi-constraint global optimization advantages of model predictive control, and the real-time error correction advantages of the disturbance observer. The inputs to this step include the full set of valid grid-side data output from S2, the grid-side discrete prediction model and total disturbance observations constructed in S3, and the grid-connected converter constraint boundaries output from S3. Simultaneously, it receives the flywheel charging and discharging power constraint boundaries from the S5 generator-side algorithm in real time. The final outputs are the optimal control quantity for the grid-side converter and the optimal grid-connected power prediction value. The optimal grid-connected power prediction value serves as the core input to the S5 generator-side collaborative algorithm, achieving global coordination between grid-side frequency support requirements and generator-side power execution capabilities.

[0054] First, clarify the core logic of the two-way interaction of the three types of algorithms: 1. Two-way interaction between virtual synchronous generator control and model predictive control: Virtual synchronous generator control provides model predictive control with a power reference trajectory and state reference target that meet the grid frequency support requirements; Model predictive control, based on grid-side equipment constraints and generator-side feedback power capability constraints, corrects the feasible region of the reference trajectory and updates the corrected reference value in reverse to the input of virtual synchronous generator control, avoiding the generation of power commands that exceed the equipment's execution capability, and achieving closed-loop matching between support requirements and execution capabilities.

[0055] 2. Two-way interaction between disturbance observer and model predictive control: The total disturbance observation value output by the disturbance observer is embedded into the prediction model of the model predictive control in real time to correct the prediction bias and improve the prediction accuracy; the deviation between the prediction value of the previous cycle output by the model predictive control and the actual filtered value output by S2 is input back to the disturbance observer as an error driving term to update the total disturbance observation value in real time, thereby realizing closed-loop optimization of disturbance observation and predictive control.

[0056] 3. Two-way interaction between virtual synchronous generator control and disturbance observer: The total disturbance observation value output by the disturbance observer and the grid frequency change rate output by S2 jointly drive the virtual inertia and damping coefficient of the virtual synchronous generator to adaptively adjust, thereby improving the frequency support capability under large disturbances; the internal disturbance of the grid-side prediction model caused by the adaptive adjustment of virtual synchronous generator parameters is captured in real time by the disturbance observer and compensated through the total disturbance observation value, ensuring that the model prediction accuracy is not affected by the dynamic adjustment of parameters.

[0057] The first step is to complete the state prediction in the prediction time domain. Based on the grid-side discrete prediction model built on S3, and combined with the total disturbance observation value output by the disturbance observer, a corrected multi-step prediction model is constructed to achieve accurate prediction of the grid connection state in the future time domain. In this scheme, the prediction time domain is represented by the symbol N. p This indicates that the value is 3, and the time domain is controlled by the symbol N. c This indicates that a value of 1 indicates that this parameter combination ensures both the foresight of the prediction and keeps the computational load within an acceptable range for real-time control, meeting the requirements of microsecond-level control cycles. The corrected multi-step prediction model is as follows: X g (k+i|k)=A gd ·X g (k+i-1|k)+B gd ·U g (k+i-1|k)+d total_g_obs (k); In the formula, X g (k+i|k) represents the predicted value of the network-side state variable for the (k+i)th sampling period in the kth sampling period; i is the step size index in the prediction time domain, ranging from 1 to N. p A gdB gd U represents the discretized network-side state matrix and input matrix in S3; g (k+i-1|k) represents the predicted value of the network-side control quantity for the k+i-1th sampling period from the kth sampling period; d total_g_obs (k) represents the total network-side disturbance observation value output by the second-order disturbance observer in S3.

[0058] The second step completes the bidirectional interaction between virtual synchronous generator control and model predictive control, achieving a closed-loop match between frequency support requirements and equipment execution capabilities. First, the virtual synchronous generator control generates a grid-connected power reference trajectory in the predictive time domain. This trajectory is generated based on the active power-frequency droop characteristics and rotor motion equations, taking into account both rapid adjustment of frequency deviation and grid inertia support. The formula for calculating the power reference trajectory in the predictive time domain is: ; In the formula, P g_ref (k+i|k) represents the grid-connected power reference value for the k+i sampling period in the k-th sampling period; P0 is the grid-connected steady-state active power reference value of the flywheel energy storage system, which is 0 under steady-state conditions; k p The active-frequency droop factor; k i The integral regulation coefficient for virtual synchronous generator control, whose specific value is determined using conventional VSG parameter tuning methods in this field, is used to eliminate frequency steady-state deviation; Δf g_filt (k) is the power grid frequency deviation filter value output by S2 in the kth sampling period; △f g_filt (k+m|k) is the predicted frequency deviation value for the (k+m)th sampling period, obtained using a known linear prediction method in the short prediction time domain based on the rotor motion equation of the virtual synchronous generator in S3. The expression is: △f g_filt (k+m|k)=△f g_filt (k)+ROCOF filt (k)·m·T s Where m is the prediction step size, taking values ​​of 1, 2, or 3; T s This is the system sampling period.

[0059] Model predictive control, based on grid-connected power constraint boundaries, performs feasible region limiting correction on the generated reference trajectory. The constraint boundaries simultaneously include the grid-side converter rated constraints output by S3 and the flywheel power capability constraints fed back by the machine-side algorithm in S5, ensuring that the corrected reference trajectory remains within the equipment's executable range. The correction formula is: P g_ref_corr (k+i|k)=sat(P g_ref (k+i|k),P g_min_corr ,P g_max_corr ); In the formula, P g_ref_corr (k+i|k) represents the corrected grid-connected power reference value predicted at time k for the (k+i)th sampling period; sat(·) is the saturation limiting function, which outputs the corresponding boundary value when the input value exceeds the upper and lower limits; P g_max_corr The revised maximum allowable grid-connected power is determined jointly by the rated power of the grid-side converter and the maximum charging power of the flywheel; P g_min_corr The revised minimum allowable grid-connected power is determined by the rated power of the grid-side converter and the maximum discharge power of the flywheel, both of which are updated in real time by the S5 machine-side algorithm.

[0060] The corrected power reference trajectory updates the input parameters of the virtual synchronous generator control in reverse, achieving closed-loop interaction between the two. The update formula is as follows: P m_ref (k+1)=P g_ref_corr (k+1|k); In the formula, P m_ref (k+1) represents the virtual mechanical power reference value of the virtual synchronous generator in the (k+1)th sampling period. Through this reverse update, the feasible power reference value corrected by model predictive control is directly used as the input for the virtual synchronous generator control in the next period, avoiding the virtual synchronous generator from generating power commands exceeding the flywheel's execution capability due to large grid disturbances, while ensuring optimal execution of frequency support requirements. The corrected power reference trajectory is also converted into reference trajectories for grid-connected current and grid-connected point voltage, serving as the tracking target for the model predictive control cost function.

[0061] The third step completes the bidirectional interaction between the disturbance observer and model predictive control, achieving real-time correction of prediction errors and closed-loop optimization of disturbance observations. In the forward interaction, the total disturbance observation value output by the disturbance observer is directly embedded into each iteration of the prediction model, as shown in the multi-step prediction model in the first step, correcting the model prediction value in real time and eliminating prediction biases caused by internal and external disturbances. In the reverse interaction, the deviation between the predicted value of the model predictive control in the previous cycle and the actual filtered value output by S2 serves as the error driving term of the disturbance observer, which is input back into the second-order extended state observer constructed by S3 to update the total disturbance observation value d in the next cycle. total_g_obs (k), the formula for calculating prediction bias is: e g (k)=Y g_filt (k)-C gd ·X g (k|k-1); In the formula, e g (k) represents the prediction bias of the network-side model in the kth sampling period; Y g_filt(k) represents the filtered output value of the grid-connected quantity in the kth sampling period of S2, specifically including: the filtered value of the grid-connected current of the energy storage system (grid-connected current of the energy storage system), the filtered value of the grid-connected active power, and the filtered value of the grid-connected reactive power; C gd X is the discretized network-side output matrix in S3; g (k|k-1) represents the predicted state variable value for the k-th period obtained from the model predictive control during the (k-1)-th sampling period. The deviation e g (k) Directly substitute the second-order extended state observer iteration formula of S3 to update the extended state variables and the total disturbance observation value, so as to continuously optimize the disturbance observation accuracy with the control process. Specifically, the original system state variables are: converter side current, grid current, and filter capacitor voltage; the extended state variables are: based on the original state variables, the total system disturbance is added as a new "virtual state" to construct the extended state column vector; the total disturbance observation value is updated by the linear extended state observer (LESO) based on the closed-loop recursion of the observation error, that is, the difference between the actual grid output filter value and the observed output value is used to correct the predicted value of the extended state variables, so as to obtain the real-time updated total disturbance observation value, which is used for disturbance compensation in subsequent model predictive control.

[0062] The fourth step completes the two-way interaction between the virtual synchronous generator control and the disturbance observer, achieving adaptive inertia damping and full compensation for internal disturbances within the model. In the forward interaction, the total disturbance observation value output by the disturbance observer and the grid frequency change rate output by S2 jointly drive the adaptive adjustment of the virtual inertia and damping coefficient of the virtual synchronous generator. This enhances the inertial support capability during large grid disturbances and reduces parameters to avoid system oscillations in steady state. The adaptive adjustment formula is: J vsg (k)=J vsg0 +k j ·|d total_g_obs (k)|; D vsg (k)=D vsg0 +k d ·|ROCOF filt (k)|; Among them, J vsg (k) represents the virtual moment of inertia adaptively adjusted in the kth sampling period; D vsg (k) is the virtual damping coefficient adaptively adjusted in the kth sampling period; J vsg0 D is the nominal reference value for the virtual moment of inertia. vsg0 The nominal reference value for the virtual damping coefficient; k j k is the virtual inertia adjustment factor. d ROCOF is the virtual damping adjustment factor. filt (k) is the filtered value of the power grid frequency change rate output by S2 in the kth sampling period.

[0063] In the reverse interaction, the adaptive adjustment of the virtual synchronous generator parameters will cause changes in the internal parameters of the grid-side prediction model, generating internal disturbances in the model. These internal disturbances will be captured in real time by the second-order disturbance observer and included in the total disturbance observation value in the output. Full compensation will be completed through the disturbance correction stage of the prediction model, ensuring that the model prediction accuracy is not affected by the dynamic adjustment of the virtual inertia and damping coefficient, thus achieving decoupling between adaptive adjustment and prediction accuracy.

[0064] The fifth step involves constructing the cost function for model predictive control, completing the constrained optimization solution, and outputting the optimal control quantity and the optimal grid-connected power prediction value. The cost function comprehensively considers the tracking accuracy of frequency support, the smoothness of the control quantity, and the penalty for constraint violations, ensuring that the optimization results meet both the grid frequency support requirements and the equipment safety operation requirements. The cost function formula is as follows: ; In the formula, J g Let Q be the cost function of the predictive control of the network-side model, and let Q be the optimization objective. g The weighting coefficients for the output tracking error are a 3rd order diagonal matrix. The diagonal elements correspond to the tracking weights of active power, reactive power, and grid connection point voltage, respectively. The larger the weight, the higher the tracking priority of the corresponding physical quantity. The weights are determined by conventional MPC parameter tuning methods in this field. Y g_ref (k+i|k) is the output reference value generated based on the corrected power reference trajectory; Y g (k+i|k) represents the predicted output value for the (k+i)th cycle obtained from the network-side prediction model; R g The weighting coefficient for the rate of change of the control quantity is used to suppress drastic fluctuations in the control quantity, avoid frequent operation of converter devices, and ensure system stability. It is determined by conventional model prediction control parameter tuning methods in this field (common range: 10~100); △U g (k+i|k) represents the change in the control quantity, and the calculation formula is: △U g (k+i|k)=U g (k+i|k)-U g (k+i-1|k); ρ g This is the penalty weight for the relaxation factor, used to amplify the cost of constraint violation. Its value is a much larger positive number than the tracking weight and the control smoothing weight; in this implementation, it is set to 1000. ε gThese are slack variables, which are non-negative optimization variables in the model predictive control optimization process. They are automatically obtained by the solver. For example, when the constraint is slightly violated, the solver automatically gives a very small positive number (such as 0.01 or 0.1). When the constraint is not violated, the value is 0. They are used to soften hard constraints and avoid the optimization problem being unsolvable under extreme conditions.

[0065] The constraints corresponding to the cost function are as follows, ensuring that all optimization results meet physical constraints and equipment safety requirements: 1. Prediction model constraints: X g (k+i|k)=A gd ·X g (k+i-1|k)+B gd ·U g (k+i-1|k)+d total_g_obs (k), In the formula, X g (k+i|k) represents the predicted value of the network-side state variable for the (k+i)th sampling period from the kth sampling period; d total_g_obs Output the total disturbance observation value for the second-order disturbance observer; U g (k+i-1|k) is the predicted value of the grid-side converter output voltage control quantity for the k+i-1th period in the kth sampling period; A gd B gd These are the state matrix and input matrix of the network-side discrete prediction model, respectively. The value of i ranges from 1 to N. p This ensures that the prediction process strictly follows the physical characteristics of the network-side system; 2. Control quantity constraint: U g_min ≤U g (k+i|k)≤U g_max , In the formula, the value of i ranges from 0 to N. c -1, when i is 0, indicates that only the control quantity at the current moment is constrained; U g_min U g_max The upper and lower limits of the grid-side converter output voltage are determined by the DC bus voltage and the converter modulation strategy. 3. Output constraint: i g_min ≤i g (k+i|k)≤i g_max , In the formula, the value of i ranges from 1 to N. p i g_min i g_max The upper and lower limits of the grid-connected current are determined by the rated current of the grid-connected converter; i g (k+i|k) is the instantaneous value of the grid-connected phase current of the flywheel energy storage system predicted at the k-th sampling time and in the k+i-th sampling period; 4. Power constraint: P g_min_corr ≤P g (k+i|k)≤P g_max_corr , In the formula, the value of i ranges from 1 to N. p This ensures that the optimized grid-connected power is within the feasible range on the generator side. P g_max_corr P is the corrected maximum allowable grid-connected power. g_min_corr The revised minimum allowable grid-connected power is determined jointly by the rated power of the grid-side converter and the maximum discharge power of the flywheel; P g (k+i|k) is the predicted active power of the flywheel energy storage system at the k-th sampling time and the predicted value of the flywheel energy storage system at the k+i-th sampling period.

[0066] Within each sampling period, based on the cost function formula and four constraints, the cost function is minimized. Under the constraints of the four constraints, a constrained quadratic programming optimization problem is solved to obtain the optimal control quantity sequence in the control time domain. {U g (k|k), U g (k+1|k)...,U g (k+N c -1|k)}; Based on the principle of rolling time-domain optimization, only the first control variable U in the sequence is considered. g (k|k) serves as the actual output control quantity of the grid-side converter, driving the grid-side converter to perform the corresponding switching action, i.e., N c =1, only the optimal control quantity U at the current moment. g (k|k) is the actual output control quantity; U g (k|k) is the reference vector for the output voltage of the grid-side converter in the synchronous rotating dq coordinate system. It naturally contains frequency and phase information synchronized with the grid. After SVPWM modulation, it generates three-phase PWM drive pulses to control the grid-side converter output to achieve three-phase grid-connected AC power that is completely synchronized with the grid frequency, voltage, and phase. Simultaneously, based on the optimal control quantity and the grid-side power prediction model after total disturbance correction, the optimal grid-connected power prediction value P is calculated. g_opt (k+1) is the predicted value of the grid-connected active power in the (k+1)th sampling period under the action of the optimal control quantity. This value is directly used as the core input of the S5 machine-side collaborative control algorithm to realize the accurate transmission of grid-side frequency support requirements to the machine-side power execution link.

[0067] S5. Implementation of the flywheel output power prediction and collaborative control algorithm on the machine side; the core function of this step is to accurately execute the grid-connected power demand from the grid side while maintaining DC bus voltage stability, ensuring the safe and reliable operation of the flywheel equipment. (Reference) Figure 4 This step deeply integrates three types of algorithms: model predictive control, PID control, and fuzzy control. It breaks away from the traditional unidirectional architecture of master-slave control, achieving bidirectional closed-loop interaction between the three algorithms. This fully leverages the advantages of model predictive control in DC bus voltage regulation and power point tracking, the flywheel speed safety protection advantage of PID control, and the multi-objective adaptive weight allocation advantage of fuzzy control. The inputs to this step include the full set of valid flywheel-side data output from S2, the optimal grid-connected power prediction value output from S4, and the machine-side dynamic model and constraint boundaries constructed in S3. Simultaneously, the flywheel charging and discharging power constraint boundaries calculated in real-time on the machine side are fed back to the grid-side algorithm in S4, achieving global coordination between machine-side execution capabilities and grid-side support requirements. Finally, the optimal control quantity for the flywheel motor is output, driving the flywheel motor to precisely execute charging, discharging, or standby actions.

[0068] First, clarify the core logic of the two-way interaction of the three types of algorithms: 1. Two-way interaction between fuzzy control and model predictive control and PID control: Fuzzy control is based on the DC bus voltage deviation and flywheel speed deviation, and dynamically outputs the weight coefficients of model predictive control and PID control. The sum of the two weight coefficients is always 1, so as to achieve an adaptive balance of three objectives: grid frequency support, DC bus voltage stabilization and flywheel safety protection. The tracking error of model predictive control and PID control is input back into fuzzy control, dynamically adjusting the quantization factor and domain range of fuzzy control, and optimizing the adaptive adjustment capability of fuzzy control.

[0069] 2. Two-way interaction between model predictive control and PID control: Model predictive control takes DC bus voltage stability and grid-connected power tracking as its core objectives. Its output power feedforward provides feedforward compensation for PID control, improving the dynamic response speed of PID control. PID control takes flywheel speed safety as its core objective. Its output speed safety constraint corrects the constraint boundary of model predictive control in reverse, ensuring that the optimal solution of model predictive control is always within the safe operating range of the flywheel and avoiding overcharging and over-discharging of the equipment.

[0070] 3. Two-way interaction between machine-side and grid-side algorithms: In the forward interaction, the optimal grid-connected power prediction value output by S4 serves as the core input for machine-side model predictive control, ensuring that the machine-side power output and the grid-side grid-connected power demand are accurately matched; in the reverse interaction, the flywheel charging and discharging power constraint boundary calculated in real time by the machine-side algorithm is fed back to the S4 grid-side algorithm to correct the feasible region of the grid-side power reference trajectory and prevent the grid side from generating power commands that exceed the flywheel's capacity.

[0071] First, two core control deviations are defined as the basic inputs for the three types of algorithms. Both deviations are calculated based on the valid data output by S2. The first is the DC bus voltage deviation, which reflects the degree to which the DC bus voltage deviates from its rated value. It is the core control objective for maintaining DC bus stability and is denoted by ΔU. dc (k) represents the calculation formula: △U dc (k)=U dc_ref -U dc_filt (k); In the formula, U dc_ref The rated reference voltage of the DC bus is taken as U, which is the rated voltage of the DC bus. dcN ;U dc_filt (k) is the DC bus voltage filter value of the kth sampling period output by S2.

[0072] The second is the flywheel speed deviation, which reflects the degree to which the actual speed of the flywheel deviates from the safe reference value. It is the core control target for ensuring the safety of flywheel equipment and is represented by the symbol △ω. f (k) represents the calculation formula: △ω f (k)=ω f_ref (k)-ω f_filt (k); In the formula, ω f_ref (k) is the reference value of the flywheel speed in the kth sampling period, which is determined by the safety target of the flywheel energy storage capacity and the rated parameters of the equipment; ω f_filt (k) is the filtered value of the flywheel real-time speed in the kth sampling period output by S2.

[0073] The first step is to construct the generator-side model predictive control module. The core objective is to maintain a stable DC bus voltage and accurately track the optimal grid-connected power demand (predicted optimal grid-connected power value) output by S4. This ensures real-time matching between the grid-side power demand and the generator-side flywheel charging and discharging power, preventing converter malfunctions caused by drastic fluctuations in DC bus voltage. The prediction time domain of the generator-side model predictive control is denoted by N. pf This indicates that the value is 3, and the time domain is controlled by the symbol N. cf This indicates that a value of 1 is used to balance predictive foresight and real-time computational efficiency. Based on the discrete dynamic model of the machine side constructed using S3, and combined with the DC bus power balance equation, a multi-step prediction model of the machine side is constructed. The model parameters adopt the latest results obtained from S3 parameter identification to ensure prediction accuracy. The multi-step prediction model of the machine side is as follows: ω f (k+i|k)=ω f (k+i-1|k)+(T s / J f (k))·(T ef (k+i-1|k)-Bω (k)·ω f (k+i-1|k)-T load_filt (k)); ; ; In the formula, ω f (k+i|k) represents the predicted flywheel speed for the (k+i)th sampling period from the kth sampling period; i is the step size index in the prediction time domain, ranging from 1 to N. pf J f (k), B ω (k) represents the latest flywheel moment of inertia and mechanical damping coefficient obtained from the S3 parameter identification in the kth sampling period; T ef (k+i-1|k) is the predicted value of the electromagnetic torque in the (k+i-1)th sampling period; T load_filt (k) is the flywheel load torque filter value of the kth sampling period output by S2; U dc (k+i|k) represents the predicted DC bus voltage for the (k+i)th sampling period from the kth sampling period; P g_opt (k+i|k) is the sequence of optimal grid-connected power prediction values ​​in the prediction time domain output by S4; P dc_f (k+i-1|k) represents the predicted DC bus-side generator power value for the (k+i-1)th sampling period; i sq_ref (k+i-1|k) is the machine-side control input, i.e., the reference value of the flywheel motor stator q-axis current; η f P represents the rated efficiency of the generator-side converter. loss_inv (k) represents the power loss of the machine-side converter.

[0074] Based on the machine-side multi-step prediction model, a cost function for machine-side model predictive control is constructed. The core optimization objective is accurate tracking of the DC bus voltage, while also considering the smoothness of the control quantity and constraint satisfaction. The cost function formula is as follows: ; In the formula, J f_mpc Let Q be the cost function of the machine-side model predictive control, and let Q be the optimization objective. dc This is the weighting coefficient for DC bus voltage tracking error, and it is set to a large positive number (1000 in this embodiment) to ensure the priority of DC bus voltage stability; R f The weighting coefficient for the rate of change of the control quantity is used to suppress drastic fluctuations in the stator current reference value and protect the flywheel motor and the machine-side converter; △i sq_ref (k+i|k) represents the change in the control quantity, and the calculation formula is: △i sq_ref (k+i|k)=i sq_ref(k+i|k)-i sq_ref (k+i-1|k); ρ f The penalty weight for the relaxation factor is a large positive number; ε f These are slack variables used to soften hard constraints and prevent the optimization problem from becoming unsolvable.

[0075] The constraints corresponding to the cost function are as follows, ensuring that the optimization results meet the safe operation requirements of the flywheel equipment and the DC system: 1. Machine-side prediction model constraints: In the above flywheel speed prediction equation, DC bus voltage prediction equation, and DC bus-side machine-side power prediction equation, the value of i ranges from 1 to N. pf This ensures that the prediction process strictly follows the physical characteristics of the machine-side system; 2. Stator current constraint: i s_min ≤i sq_ref (k+i|k)≤i s_max The value of i ranges from 0 to N. cf -1, the constraint boundary is determined by the motor stator current constraint boundary output by S3; 3. Flywheel speed constraint: ω f_min ≤ω f (k+i|k)≤ω f_max The value of i ranges from 1 to N. pf ω f_min ω represents the minimum permissible operating speed of the flywheel, corresponding to a lower limit of 0.1 for the flywheel's energy storage capacity. f_max This represents the maximum permissible operating speed of the flywheel, corresponding to a maximum flywheel energy storage capacity of 1.2. 4. Charge / discharge power constraint: P f_min (k)≤P f (k+i|k)≤P f_max (k), where i ranges from 1 to N. pf The constraint boundary is determined by the flywheel charging and discharging power constraint boundary output by S3.

[0076] Within each sampling period, based on the cost function J f_mpc Given the formula and four constraints, the cost function is minimized. Under these constraints, a constrained quadratic programming optimization problem is solved to obtain the optimal control quantity i for the machine-side model predictive control. sq_mpc (k|k), and correspondingly obtain the flywheel power control quantity P of the model predictive control output. f_mpc (k).

[0077] The second step involves constructing a machine-side PID control module. The core objective is to ensure the safe and stable operation of the flywheel speed, preventing overcharging, over-discharging, and mechanical damage caused by the flywheel speed exceeding the safe range. The input to the PID control is the flywheel speed deviation. A position-based discrete PID algorithm is used, and the output is a reference value for the flywheel motor stator current. The algorithm formula is as follows: ; In the formula, i sq_pid (k) represents the stator current reference value output by the PID control in the kth sampling period; K p K is the proportional coefficient for PID control, used for rapid response to speed deviation; i K represents the integral coefficient of the PID control, used to eliminate steady-state speed deviation; d These are the derivative coefficients for PID control, used to suppress speed fluctuations. The three coefficients are tuned based on the control bandwidth and stability requirements of the flywheel speed loop.

[0078] The output of the PID control is limited by the stator current constraint boundary output by S3 to obtain the limited current reference value i. sq_pid_sat (k), and correspondingly obtain the flywheel power control quantity P of the PID control output. f_pid_sat (k) ensures that the output control commands meet the requirements for safe operation of the motor.

[0079] The third step involves constructing a fuzzy control module. Its core function is to dynamically allocate the weight coefficients of model predictive control and PID control based on DC bus voltage deviation and flywheel speed deviation, achieving an adaptive balance between grid frequency support performance and flywheel equipment safety. When the grid frequency deviation is large and DC bus voltage fluctuations are significant, the weight of model predictive control is increased to prioritize grid-connected power tracking and DC bus stability, meeting grid frequency support requirements. When the flywheel speed approaches the safety boundary and the speed deviation is large, the weight of PID control is increased to prioritize flywheel equipment safety and prevent equipment damage.

[0080] First, clarify the input and output variables of fuzzy control, and complete the quantization and fuzzification of the input variables: 1. Input variable 1 is the DC bus voltage deviation ΔU dc (k), the quantized universe of discourse is [-6,6], the fuzzy subset contains 7 subsets: negative large, negative medium, negative small, zero, positive small, positive medium, and positive large, and the membership function adopts the triangular membership function; 2. Input variable 2 is the flywheel speed deviation Δω f (k), the quantized universe of discourse is [-6,6], the fuzzy subset contains 7 subsets: negative large, negative medium, negative small, zero, positive small, positive medium, and positive large, and the membership function adopts the triangular membership function; 3. Output variable 1 is the weighting coefficient α of the model predictive control.mpc (k), with a universe of discourse of [0,1], and a fuzzy subset containing five elements: minimal, small, medium, large, and maximal. The membership function is a Gaussian membership function. 4. Output variable 2 is the weighting coefficient α of the PID control. pid (k), satisfying α mpc (k)+α pid (k)=1, ensuring that the weighted sum of the two control variables is full control, with no control gaps or control superposition redundancy.

[0081] Based on the 7×7 input fuzzy subset combination, a complete set of 49 fuzzy rules is designed to cover all input combination scenarios and ensure adaptive weight allocation under all working conditions. The fuzzy rules are shown in Table 1.

[0082] Table 1. Fuzzy Control Rule Table

[0083] Fuzzy inference employs the Mamdani inference method to map input fuzzy values ​​to output fuzzy values; defuzzification uses the centroid method to convert the membership function of the fuzzy output into precise weight coefficients α. mpc (k) and α pid (k), the centroid method for fuzzification is: ; In the formula, x j For discrete points within the output universe; μ j represents the membership degree value of the corresponding discrete point; n is the total number of discrete points in the universe of discourse. The weight coefficients calculated by this formula strictly satisfy α. mpc (k)+α pid The constraint (k)=1.

[0084] The fourth step involves completing the bidirectional interaction between the three types of algorithms, achieving deep integration and global collaboration of the algorithms, rather than simple weighted summation. At the same time, it completes the bidirectional interaction between the machine side and the network side, realizing the collaborative control of the entire system.

[0085] The first interaction is a two-way interaction between fuzzy control and model predictive control / PID control. In the forward interaction, the weighting coefficients of the fuzzy control output are directly used for the weighted fusion of the two control quantities to obtain the final optimal stator current reference value for the flywheel motor. The weighted fusion formula is as follows: i sq_opt (k)=α mpc (k)·i sq_mpc (k|k)+α pid (k)·i sq_pid_sat (k); Among them, i sq_opt(k) is the optimal stator current reference value of the final output of the machine side in the kth sampling period. This value is directly used as the control input of the machine side converter. It is converted into a converter switching signal through the vector control algorithm to drive the flywheel motor to perform the corresponding charging and discharging actions.

[0086] Based on the optimal stator current reference value, the corresponding optimal flywheel output power reference value is calculated using the following formula: P f_opt (k)=1.5·p n ·ψ f ·i sq_opt (k)·ω f_filt (k); Among them, P f_opt (k) is the optimal flywheel output power reference value for the kth sampling period, which is used for the subsequent closed-loop effect evaluation of S6.

[0087] In the reverse interaction, the DC bus voltage tracking error of model predictive control and the flywheel speed tracking error of PID control are used as input variables for fuzzy control in the next sampling period. Simultaneously, based on changes in the tracking error, the quantization factor of the fuzzy control input variables is dynamically adjusted to optimize the sensitivity and adjustment accuracy of the fuzzy control. The adaptive adjustment formula for the quantization factor is: △U dc_max (k)=max(|△U dc (k)|,|△U dc (k-1)|,|△U dc (k-2)|); k u =6 / △U dc_max (k); Among them, △U dc_max (k) represents the maximum absolute value of the DC bus voltage deviation over the past three sampling periods; k u The quantization factor for the DC bus voltage deviation, when ΔU dc_max (k) takes effect when it is not 0, enabling adaptive adjustment of the fuzzy control domain. This improves adjustment accuracy when the deviation is small and response speed when the deviation is large. Using the same method, the quantization factor of the speed deviation is adaptively adjusted based on the flywheel speed tracking error to achieve optimal performance of fuzzy control under all operating conditions.

[0088] The second interaction involves a two-way interaction between model predictive control and PID control, achieving deep synergy between power point tracking performance and equipment safety protection. In the forward interaction, model predictive control generates a power feedforward based on the grid-connected power demand output by S4, correcting the speed deviation input of PID control, compensating in advance for speed fluctuations caused by changes in grid-connected power, and improving the dynamic response capability of PID control. The feedforward correction formula is: ; Where, △ω f_ff (k) represents the speed feedforward compensation amount for the kth sampling period; P g_opt (k+1) represents the optimal grid-connected power prediction value output by S4; J f (k) represents the latest value of the flywheel moment of inertia obtained from the S3 parameter identification.

[0089] Based on the feedforward compensation, the input speed deviation of the PID control is corrected. The corrected speed deviation formula is as follows: △ω f_corr (k)=△ω f (k)+△ω f_ff (k); In the formula, △ω f_corr (k) represents the corrected speed deviation; The corrected speed deviation is input into the PID control algorithm, which effectively reduces flywheel speed fluctuations when grid-connected power changes rapidly, and improves the synchronization of power point tracking and system stability.

[0090] In the reverse interaction, PID control generates speed safety constraints based on the flywheel speed safety boundary, and then reversely corrects the constraints of model predictive control to ensure that the optimal solution of model predictive control is always within the safe operating range of the flywheel. This avoids model predictive control generating control commands that exceed the flywheel's safety capabilities in order to meet grid-connected power requirements. The constraint correction formula is: ω f_min_safe (k)=max(ω f_min ,ω f_filt (k)+△ω f_pid_min (k)); ω f_max_safe (k)=min(ω f_max ,ω f_filt (k)+△ω f_pid_max (k)); Where, ω f_min_safe (k), ω f_max_safe (k) represents the upper and lower limits of the modified flywheel speed safety constraint; △ω f_pid_min (k), △ω f_pid_max (k) represents the upper and lower safety limits of the PID control output speed, calculated based on the flywheel speed deviation and safety boundary. The corrected speed constraints are updated in the constraints of the machine-side model predictive control to achieve reverse constraints of power control for equipment safety protection, ensuring the long-term reliable operation of the flywheel equipment.

[0091] The third interaction is a bidirectional interaction between the machine-side algorithm and the grid-side algorithm, achieving a global match between grid-side frequency support requirements and machine-side execution capabilities. In the forward interaction, the optimal grid-connected power prediction value output by S4 is directly used as the core input of the machine-side model predictive control, embedded in the DC bus voltage prediction equation of the machine-side multi-step prediction model, ensuring that the machine-side power output accurately matches the grid-side grid-connected power requirements and maintaining DC bus voltage stability. In the reverse interaction, the machine-side algorithm, based on the flywheel effective data output by S2, calculates the flywheel charging and discharging power constraint boundary in real time, and feeds this constraint boundary back to the S4 grid-side algorithm to correct the feasible region of the grid-side power reference trajectory and the power constraint of the model predictive control. The correction formula is as follows: P g_max_corr (k)=min(P g_max ,P f_max (k)·η f ·η g ); P g_min_corr (k)=max(P g_min ,P f_min (k)·η f ·η g ); In the formula, P g_max P g_min P represents the upper and lower limits of the rated power of the grid-side converter; f_max (k), P f_min (k) represents the upper and lower limits of the flywheel charging and discharging power constraint calculated in real time on the machine side; η g This represents the rated efficiency of the grid-side converter. The corrected power constraint is directly used for limiting the virtual synchronous generator reference trajectory in S4 and for constraining the model predictive control, ensuring that the grid-connected power command generated by the grid side is always within the executable range of the flywheel energy storage system, thus avoiding control failure problems caused by demand-capacity mismatch.

[0092] S6. Closed-loop performance evaluation and adaptive updating of model parameters and algorithm weights: The core function of this step is to quantitatively evaluate the frequency support control effect based on the actual operating data after S1 to S5, and to complete the adaptive closed-loop update of model parameters, algorithm weights, and constraint boundaries, thereby achieving continuous optimization of control performance. The inputs to this step are the full amount of valid data output from S2, the control quantities and predicted values ​​output from S4 and S5, and the model and parameters constructed in S3. All updated parameters and weights are synchronized to the corresponding modules in S3, S4, and S5 as inputs for the next sampling period, forming the closed-loop control loop of the entire system. The system will repeat the entire process from S1 to S6 until the grid frequency stabilizes within the rated range, achieving full closed-loop optimized control of grid frequency support.

[0093] First, quantitative evaluation indicators for control effectiveness are defined, including one core evaluation indicator and four auxiliary evaluation indicators. All indicators are calculated based on the effective data output by S2 and the actual control results, providing a quantitative basis for adaptive parameter updates.

[0094] The core evaluation index is the power grid frequency deviation evaluation index, which is used to quantify the core control effect of power grid frequency support. It is represented by the symbol E(k), and the calculation formula is: E(k) = |Δf g_filt (k)|; Where E(k) is the frequency deviation evaluation index for the k-th sampling period; △f g_filt (k) is the power grid frequency deviation filter value output by S2 in the kth sampling period. In this scheme, the target control value of the frequency deviation is 0.02 Hz, that is, when E(k) ≤ 0.02 Hz, the frequency control effect is judged to be satisfactory.

[0095] The auxiliary evaluation indicators are used to quantify the control performance and model accuracy of each component of the system, providing a basis for updating the parameters of each component. They are defined as follows: 1. DC bus voltage deviation index, used to evaluate the stability control effect of the DC bus, denoted by E. dc (k) represents the calculation formula: E dc (k)=|U dc_ref -U dc_filt (k)|; Among them, U dc_ref U is the rated reference voltage of the DC bus. dc_filt (k) is the DC bus voltage filter value of the kth sampling period output by S2.

[0096] 2. Flywheel speed tracking error index, used to evaluate the control effect and safety protection performance of the flywheel speed, denoted by the symbol E. ω (k) represents the calculation formula: E ω (k)=|ω f_ref (k)-ω f_filt (k)|; Where, ω f_ref (k) is the reference value for flywheel speed; ω f_filt (k) is the filtered value of the flywheel real-time speed in the kth sampling period output by S2.

[0097] 3. The prediction error index of the network-side model is used to evaluate the accuracy of the network-side prediction model and to provide a basis for updating the parameters of the network-side model and observers. It is denoted by E. model_g (k) represents the calculation formula: E model_g (k)=|Y g_filt(k)-Y g (k|k-1)|; Among them, Y g_filt (k) is the filtered value of the grid-connected output during the kth sampling period of S2; Y g (k|k-1) is the predicted output value of the k-th cycle obtained by the network-side model predictive control in the (k-1)-th sampling cycle.

[0098] 4. Aircraft-side model prediction error index, used to evaluate the accuracy of the aircraft-side dynamic model and provide a basis for updating aircraft-side model parameters, denoted by E. model_f (k) represents the calculation formula: E model_f (k)=|ω f_filt (k)-ω f (k|k-1)|; Where, ω f_filt (k) is the filtered value of the flywheel real-time speed in the kth sampling period output by S2; ω f (k|k-1) is the predicted flywheel speed value for the kth cycle obtained by the machine-side model predictive control in the (k-1)th sampling cycle.

[0099] Based on the above evaluation indicators, the adaptive update of model parameters and algorithm weights is completed in three parts. All update rules have clearly defined trigger conditions and update magnitudes to ensure smooth and stable parameter adjustment and avoid system oscillation caused by frequent parameter fluctuations.

[0100] The first part completes the adaptive update of the grid-side model and algorithm parameters. Based on the frequency deviation evaluation index and the grid-side model prediction error index, the virtual synchronous generator parameters, model prediction control weights, and disturbance observer parameters are updated to optimize the grid-side frequency support performance and model prediction accuracy.

[0101] When E(k) > 0.02 Hz and this state persists for more than 3 sampling periods, the power grid frequency is determined to be in a state of large disturbance, and frequency support capability needs to be strengthened. The following parameters are updated: 1. Increase the virtual inertia and virtual damping coefficient of the virtual synchronous generator to enhance the system's inertial support capability, suppress rapid frequency fluctuations, and accelerate frequency deviation convergence. The updated formula is as follows: J vsg (k+1)=J vsg (k)·(1+k j_up ·E(k)); D vsg (k+1)=D vsg (k)·(1+k d_up ·E(k)); Among them, J vsg (k+1),D vsg(k+1) represents the updated virtual moment of inertia and virtual damping coefficient; k j_up k d_up This is a positive adjustment coefficient, with a value ranging from 0.1 to 0.5, to ensure smooth parameter adjustment and avoid system oscillation; at the same time, an upper limit is set for the parameter to prevent excessive parameter values ​​from causing system overshoot and oscillation.

[0102] 2. Increase the output tracking weight in the model predictive control cost function to improve the tracking accuracy of the power reference trajectory and accelerate the convergence speed of the frequency deviation. The updated formula is as follows: Q g (k+1)=Q g (k)·(1+k q_up ·E(k)); Among them, Q g (k+1) represents the updated output tracking weight coefficients; k q_up This is the weight adjustment coefficient, with a value range of 0.05 to 0.2. At the same time, upper and lower limits for the weight are set to avoid drastic fluctuations in the control quantity due to weight imbalance.

[0103] 3. Increase the bandwidth of the disturbance observer to improve the response speed of disturbance observation and quickly correct model errors caused by large disturbances on the grid side. The update formula is as follows: ω o_g (k+1)=min(ω o_g_max ,ω o_g (k)·(1+k o_up ·E model_g (k))); Where, ω o_g (k+1) is the updated observer bandwidth; ω o_g_max The upper limit of the observer bandwidth is determined by the system's noise suppression capability; k o_up This is the observer bandwidth adjustment factor, with a value ranging from 0.1 to 0.3.

[0104] When E(k) ≤ 0.02 Hz, the system is in steady-state condition, and only when the network-side model prediction error E model_g (k) When the set threshold is exceeded, the observer gain and model weight are fine-tuned while the other parameters remain unchanged to avoid frequent parameter adjustments that could cause steady-state fluctuations in the system and ensure the stability of steady-state operation.

[0105] The second part completes the adaptive update of the aircraft-side model parameters. Based on the aircraft-side model prediction error index and frequency deviation evaluation index, the aircraft-side dynamic model parameters and parameter identification forgetting factor are updated to improve the prediction accuracy and parameter tracking capability of the aircraft-side model.

[0106] When E(k) > 0.02 Hz, and E model_f (k)>0.01·ω fNIf the prediction error of the machine-side model is large, it is necessary to speed up parameter identification, correct the model parameters, and perform the following updates: 1. Reduce the forgetting factor of recursive least squares method to improve the speed of parameter identification in tracking changes in model parameters, and quickly adapt to model parameter fluctuations caused by changes in operating conditions. The update formula is as follows: λ(k+1)=max(λ min ,λ(k)-k λ_down ·E model_f (k)); Where λ(k+1) is the updated forgetting factor; λ min The lower limit of the forgetting factor is set at 0.90 to avoid excessively small forgetting factors that could lead to divergence in parameter identification; k λ_down This is the forgetting factor adjustment coefficient, with a value ranging from 0.01 to 0.05.

[0107] 2. The proportional adjustment method is used to directly correct the core parameters of the machine-side model, quickly compensate for model prediction errors, and improve prediction accuracy. The update formula is as follows: ; ; Among them, J f (k+1), B ω (k+1) represents the updated flywheel moment of inertia and mechanical damping coefficient; k j_f k b_f This is the adjustment coefficient for the machine-side parameters, ranging from 0.05 to 0.2, to ensure the stability of parameter adjustment and avoid sudden parameter changes. When E(k) ≤ 0.02 Hz, and E... model_f (k)≤0.01·ω fN At the same time, keep the machine-side model parameters unchanged, and only use the recursive least squares method in S3 to complete the regular parameter iteration, maintain the steady-state accuracy of the model, and avoid frequent parameter adjustments from affecting system stability.

[0108] The third part completes the adaptive update of the domain of the machine-side fuzzy control weights. Based on the frequency deviation evaluation index, DC bus voltage deviation index, and flywheel speed deviation index, the domain range of the weight coefficients is adjusted to optimize the dynamic allocation capability of the control weights and realize the adaptive adjustment of the priority of multiple control objectives.

[0109] When E(k) > 0.02 Hz, the grid frequency deviation is large, and priority needs to be given to ensuring the grid frequency support performance. In this case, the following universe of discourse update is performed to increase the upper limit of the weight of the model predictive control and enhance the DC bus voltage regulation and grid-connected power tracking capabilities: α mpc_max (k+1)=min(1.0,α mpc_max (k)+k α_up ·E(k)); α pid_max (k+1)=max(0.0,α pid_max (k)-k α_down ·E(k)); Where, α mpc_max (k+1) represents the upper bound of the universe of discourse for the updated model prediction control weight coefficients; α pid_max (k+1) represents the upper bound of the universe of discourse for the updated PID control weight coefficients; k α_up k α_down The weighting universe adjustment coefficient ranges from 0.1 to 0.3, always satisfying α. mpc_max (k)+α pid_max The constraint (k)=1.

[0110] When E ω (k)>0.05·ω fN When the flywheel speed deviation is large and close to the safety boundary, and it is necessary to prioritize equipment safety, the following domain update is performed to increase the upper limit of the PID control weight and strengthen the flywheel speed safety protection capability: α pid_max (k+1)=min(1.0,α pid_max (k)+k α_pid_up ·E ω (k)); α mpc_max (k+1)=max(0.0,α mpc_max (k)-k α_mpc_down ·E ω (k)); Where, k α_pid_up k α_mpc_down The weighting universe adjustment coefficients range from 0.1 to 0.3, always satisfying the constraint that the sum of the weighting coefficients is 1. This applies when E(k) ≤ 0.02 Hz and E ω (k)≤0.05·ω fN When the system is in steady state, the domain of fuzzy control and the quantization factor remain unchanged to ensure the smoothness of steady-state control and avoid system fluctuations caused by frequent weight switching.

[0111] Finally, the termination condition and restart mechanism of the closed-loop regulation are clarified to achieve steady-state maintenance after the grid frequency stabilizes and rapid response after disturbances. The termination condition of the closed-loop regulation is as follows: when the grid frequency deviation filter value is within ±0.05 Hz for 100 consecutive sampling periods, and the frequency deviation evaluation index E(k) is less than or equal to 0.02 Hz for 100 consecutive sampling periods, the grid frequency is determined to be stable within the rated range. At this time, the adaptive update of model parameters and algorithm weights is stopped, the current control strategy and parameters are maintained, and the system enters the steady-state maintenance mode.

[0112] When the grid frequency deviates from the rated range again, i.e., the absolute value of the grid frequency deviation filter value exceeds 0.05 Hz, the steady-state holding mode is immediately exited, the parameter adaptive update and closed-loop adjustment process is restarted, and the entire process from S1 to S6 is repeated to quickly respond to grid frequency disturbances and provide continuous frequency support.

Claims

1. A flywheel energy storage grid frequency support method based on model predictive control, characterized in that, It includes the following steps: S1. Synchronously collect power grid operating parameters and flywheel energy storage system status parameters; The power grid operating parameters include: real-time power grid frequency, power grid voltage amplitude, grid-connected current of energy storage system, power grid frequency change rate, and power grid frequency deviation. The status parameters of the flywheel energy storage system include: real-time flywheel speed, actual flywheel charging and discharging power, flywheel motor stator current, DC bus voltage, flywheel mechanical loss power, and flywheel energy storage capacity. S2. Perform preprocessing on all parameters acquired and calculated in step S1, and output the full set of valid data, which are the filtered values ​​of all parameters after processing. The preprocessing operations include: outlier detection and correction, smoothing and denoising, and validity verification. S3. Based on the aforementioned valid data, construct a dual-model system with the DC bus as the intermediate link, and perform parameter identification; The dual-model system includes: a grid-side power prediction model and a flywheel energy storage device-side dynamic model; The grid-side power prediction model outputs a reference value for grid-connected power prediction to support grid frequency. The dynamic model output of the flywheel energy storage machine side is used to predict the flywheel output power for DC bus voltage support, as well as the charging and discharging power constraint boundary and the current constraint boundary. S4. Based on effective data, grid-side power prediction model and charging / discharging power constraint boundary, integrate model predictive control, disturbance observation and virtual synchronous generator control to construct a grid-connected power prediction cooperative control algorithm and output the optimal grid-connected power prediction value; S5. Based on effective data, the dynamic model of the flywheel energy storage machine, the charging and discharging power constraint boundary, the current constraint boundary, and the optimal grid-connected power prediction value, the flywheel output power prediction collaborative control algorithm is constructed by integrating model predictive control, PID control, and fuzzy control. The algorithm outputs the optimal control quantity of the flywheel motor and drives the flywheel motor to perform the corresponding charging, discharging, or standby actions.

2. The flywheel energy storage grid frequency support method based on model predictive control according to claim 1, characterized in that: It also includes step S6, which calculates the actual deviation of the control effect based on the actual operating data after S1 to S5, compares the actual deviation with the target deviation to complete the adjustment effect evaluation, and updates the weights of the dual-model system, the grid-connected power prediction collaborative control algorithm and the flywheel output power prediction collaborative control algorithm in a closed loop according to the evaluation results. The updated parameters and weights are synchronized to the corresponding models and algorithms as inputs for the next acquisition cycle to form closed-loop adjustment. S1 to S6 are repeated until the grid frequency stabilizes within the rated range.

3. The flywheel energy storage grid frequency support method based on model predictive control according to claim 1, characterized in that: In step S2, the outlier detection and correction are implemented using the 3σ criterion; the smoothing and denoising are implemented using a moving average filter; and the validity check is implemented by pre-setting a corresponding physical reasonable range for each parameter, replacing invalid data that exceeds the valid range with the valid filter value of the previous period.

4. The flywheel energy storage grid frequency support method based on model predictive control according to claim 1, characterized in that: In step S3, the grid-side power prediction model integrates the second-order linear model of the grid-connected filter, the second-order disturbance observer, and the second-order model of the virtual synchronous generator to describe the mapping relationship between the control input and the grid-connected output power and to observe the total internal and external disturbances that affect the model's accuracy; the grid-side power prediction model is: X g (k+1)=A gd ·X g (k)+B gd ·U g (k)+d total_g (k); In the formula, X g (k+1) represents the predicted value of the network-side state variable in the (k+1)th sampling period; X g (k) represents the grid-side state variables in the kth sampling period, including the grid-side converter output current, the grid-connected current, and the filter capacitor voltage; A gd B is the discretized network-side state matrix. gd U is the discretized network-side input matrix; g (k) represents the grid-side converter output voltage control quantity during the kth sampling period; d total_g (k) represents the total internal and external disturbances of the network-side model during the k-th sampling period; The dynamic model of the flywheel energy storage device is a second-order linear model used to describe the mapping relationship between control input and flywheel output power. The model parameters are identified and updated online using a recursive least squares method based on valid data from a preset number of sampling periods. The dynamic model of the flywheel energy storage device is as follows: oh f_pred (k+1)=ω f_filt (k)+(T s / J f )·(1.5·p n ·ψ f ·i sq_ref (k)-B ω ·oh f_filt (k)-T load_filt (k)); P f_pred (k+1)=1.5·p n ·ψ f ·i sq_ref (k)·ω f_pred (k+1); U dc_pred (k+1)=U dc_filt (k)+(T s / C dc )·(P dc_g (k)-P dc_f (k)-P dc_loss (k)) / (U dc_filt (k)); In the formula, ω f_pred (k+1) is the predicted flywheel speed for the (k+1)th sampling period; ω f_filt (k) represents the filtered value of the flywheel real-time speed during the kth sampling period; i sq_ref (k) represents the reference value of the flywheel motor stator q-axis current during the kth sampling period; T load_filt (k) represents the flywheel load torque filter value in the kth sampling period; T s For the sampling period, J f p is the moment of inertia of the flywheel rotor. n Let ψ be the number of pole pairs of the flywheel motor. f For the permanent magnet flux linkage of the flywheel motor; P f_pred (k+1) is the predicted flywheel output power value for the (k+1)th sampling period; U dc_pred (k+1) is the predicted value of the DC bus voltage in the (k+1)th sampling period; U dc_filt (k) represents the DC bus voltage filter value for the kth sampling period; C dc P is the capacitance value of the DC bus capacitor; dc_g (k) represents the power flowing from the grid-side converter to the DC bus during the k-th sampling period; P dc_loss (k) represents the DC bus line and switching loss power in the kth sampling period; P dc_f (k) represents the power flowing from the DC bus to the machine-side converter during the kth sampling period.

5. The flywheel energy storage grid frequency support method based on model predictive control according to claim 1, characterized in that: In step S3, the charging and discharging power constraint boundary is determined based on the rated power of the flywheel motor, the current flywheel energy storage capacity, the real-time speed of the flywheel, and the safe range of the DC bus voltage. P f_min (k)≤P f (k)≤P f_max (k); Among them, P f_max (k) represents the maximum allowable charging power in the kth sampling period, which is determined by the rated power of the flywheel motor, the current flywheel energy storage capacity filter value, and the upper limit of the DC bus voltage; P f_min (k) represents the maximum allowable discharge power in the kth sampling period, which is determined by the rated power of the flywheel motor, the current real-time flywheel speed filter value, and the lower limit of the DC bus voltage; P f (k) represents the actual charging and discharging power of the flywheel motor port at the current moment; The current constraint boundary is determined based on the rated stator current of the flywheel motor and the motor demagnetization characteristics; i s_min ≤i sq_ref (k)≤i s_max ; Among them, i s_max i represents the positive upper limit of the q-axis current of the motor stator, corresponding to the maximum charging torque; s_min The negative lower limit of the q-axis current of the motor stator corresponds to the maximum discharge torque, both determined by the motor's rated current and demagnetization characteristics; i sq_ref (k) is the reference value of the q-axis current of the flywheel motor stator in the kth sampling period.

6. The flywheel energy storage grid frequency support method based on model predictive control according to claim 4, characterized in that: In step S4, the grid-connected power prediction and coordinated control algorithm is implemented through a two-way, two-way closed-loop interaction of three types of algorithms: model predictive control, disturbance observation, and virtual synchronous generator control. The specific process is as follows: The virtual synchronous generator control provides a power reference trajectory that meets the grid frequency support requirements for model predictive control. The model predictive control corrects the power reference trajectory in the feasible region based on the charging and discharging power constraint boundary, and updates the corrected reference value in reverse to the input of the virtual synchronous generator control. The formula for calculating the power reference trajectory in the prediction time domain is: ； In the formula, P g_ref (k+i|k) represents the grid-connected power reference value for the k+i sampling period in the k-th sampling period; P0 represents the grid-connected steady-state active power reference value of the flywheel energy storage system; k p The active-frequency droop factor; Δf g_filt (k) represents the power grid frequency deviation filter value for the kth sampling period; k i The integral regulation coefficient for virtual synchronous generator control; Δf g_filt (k+m|k) is the predicted frequency deviation value for the (k+m)th sampling period; The formula for correcting the feasible region is: P g_ref_corr (k+i|k)=sat(P g_ref (k+i|k),P g_min_corr ,P g_max_corr )); In the formula, P g_ref_corr (k+i|k) is the corrected grid-connected power reference value predicted at time k for the (k+i)th sampling period; sat(·) is the saturation limiting function; P g_max_corr P is the corrected maximum allowable grid-connected power. g_min_corr The revised minimum allowable grid-connected power; The corrected formula for updating the virtual synchronous generator control with the reference value in reverse is: P m_ref (k+1)=P g_ref_corr (k+1|k); In the formula, P m_ref (k+1) is the virtual mechanical power reference value of the virtual synchronous generator in the (k+1)th sampling period; The total disturbance observation value output by the disturbance observer is embedded into the prediction model of the model predictive control in real time to complete error correction, and the prediction deviation of the model predictive control is input back into the disturbance observer to update the total disturbance observation value. The formula for calculating the prediction deviation is: e g (k)=Y g_filt (k)-C gd ·X g (k|k-1); In the formula, e g (k) represents the prediction bias of the network-side model in the kth sampling period; Y g_filt (k) is the filtered value of the grid-connected output in the kth sampling period, C gd X is the discretized grid-side output matrix in the grid-side power prediction model. g (k|k-1) represents the predicted state variable value for the k-th period obtained from the model predictive control during the (k-1)-th sampling period; The virtual synchronous generator control adaptively adjusts the virtual inertia and damping coefficient based on the output of the disturbance observer and the grid frequency change rate. The disturbance observer observes and compensates for the internal disturbances of the grid-side prediction model caused by parameter adjustments in real time. The adaptive adjustment formula is: J vsg (k)=J vsg0 +k j ·|d total_g_obs (k)|; D vsg (k)=D vsg0 +k d ·|ROCOF filt (k)|; Among them, J vsg (k) represents the virtual moment of inertia adaptively adjusted in the kth sampling period; D vsg (k) is the virtual damping coefficient adaptively adjusted in the kth sampling period; J vsg0 D is the nominal reference value for the virtual moment of inertia. vsg0 The nominal reference value for the virtual damping coefficient; k j k is the virtual inertia adjustment factor. d ROCOF is the virtual damping adjustment factor. filt (k) represents the filtered value of the power grid frequency change rate in the kth sampling period; d total_g_obs (k) represents the total network-side disturbance observation value output by the disturbance observer.

7. The flywheel energy storage grid frequency support method based on model predictive control according to claim 6, characterized in that: In step S4, the cost function of model predictive control is constructed, the constrained quadratic programming optimization problem is solved, and the optimal control quantity and the optimal grid-connected power prediction value are output based on the rolling time-domain optimization principle. The cost function formula on the network side is: ； In the formula, J g Let N be the cost function of the predictive control of the network-side model, and the optimization objective is to minimize the value of this function; p For predicting the time domain, N c To control the time domain; Q g Y is the weighting coefficient for the output tracking error; g_ref (k+i|k) is the output reference value generated based on the corrected power reference trajectory; Y g (k+i|k) represents the predicted output value for the (k+i)th cycle obtained from the network-side prediction model; R g The weighting coefficient for the rate of change of the control quantity; △U g (k+i|k) represents the change in the control variable; ρ g The penalty weight for the relaxation factor; ε g These are slack variables; The constraints corresponding to the cost function include: prediction model constraints, control quantity constraints, output quantity constraints, and power constraints. The constraint of the prediction model: X g (k+i|k)=A gd ·X g (k+i-1|k)+B gd ·U g (k+i-1|k)+d total_g_obs (k); In the formula, X g (k+i|k) represents the predicted value of the network-side state variable for the (k+i)th sampling period from the kth sampling period; d total_g_obs Output the total disturbance observation value for the second-order disturbance observer; U g (k+i-1|k) is the predicted value of the grid-side converter output voltage control quantity for the k+i-1th period in the kth sampling period; A gd B gd These are the state matrix and input matrix of the network-side discrete prediction model, respectively; the value of i ranges from 1 to N. p ; The control constraint: U g_min ≤U g (k+i|k)≤U g_max ; In the formula, the value of i ranges from 0 to N. c -1; U g_min U g_max These are the upper and lower limits of the grid-side converter output voltage; The output constraint: i g_min ≤i g (k+i|k)≤i g_max ; In the formula, the value of i ranges from 1 to N. p i g_min i g_max These are the upper and lower limits of the grid-connected current; i g (k+i|k) is the instantaneous value of the grid-connected phase current of the flywheel energy storage system in the (k+i)th sampling period predicted at the kth sampling time; The power constraint: P g_min_corr ≤P g (k+i|k)≤P g_max_corr , In the formula, the value of i ranges from 1 to N. p ;P g_max_corr P is the corrected maximum allowable grid-connected power. g_min_corr P is the revised minimum allowable grid-connected power. g (k+i|k) is the predicted active power of the flywheel energy storage system at the k+i sampling period, obtained at the k-th sampling time.

8. The flywheel energy storage grid frequency support method based on model predictive control according to claim 7, characterized in that: In step S5, the flywheel output power prediction and cooperative control algorithm is implemented through a two-way, two-way closed-loop interaction of three types of algorithms: model predictive control, PID control, and fuzzy control. The specific process is as follows: The fuzzy control dynamically outputs the weighting coefficients of the model predictive control and the PID control based on the DC bus voltage deviation and flywheel speed deviation. The DC bus voltage tracking error of the model predictive control and the flywheel speed tracking error of the PID control are input inversely to the fuzzy control, thereby completing the adaptive adjustment of the fuzzy control quantization factor and universe of discourse; specifically as follows: In the positive interaction, the weighting coefficients of the fuzzy control output are directly used for the weighted fusion of the two control quantities to obtain the final optimal stator current reference value of the flywheel motor. The weighted fusion formula is as follows: i sq_opt (k)=α mpc (k)·i sq_mpc (k|k)+α pid (k)·i sq_pid_sat (k); Among them, i sq_opt (k) represents the optimal stator current reference value finally output by the machine side in the kth sampling period, α mpc (k) and α pid (k) represents the precise weighting coefficients transformed from the fuzzy inference results; i sq_mpc (k|k) represents the optimal control quantity for machine-side model predictive control; i sq_pid_sat (k) is the current reference value after the output of the PID control is limited by the stator current constraint boundary output by S3. Based on the optimal stator current reference value, the corresponding optimal flywheel output power reference value is calculated using the following formula: P f_opt (k)=1.5·p n ·ψ f ·i sq_opt (k)·ω f_filt (k); Among them, P f_opt (k) represents the optimal flywheel output power reference value for the kth sampling period; p n ψ is the number of pole pairs of the flywheel motor; f For the permanent magnet flux linkage of the flywheel motor; ω f_filt (k) is the filtered value of the real-time speed of the flywheel in the kth sampling period; In the reverse interaction, the formula for adaptive adjustment of the quantization factor is: △U dc_max (k)=max(|△U dc (k)|,|△U dc (k-1)|,|△U dc (k-2)|); k u =6 / △U dc_max (k); Among them, △U dc_max (k) represents the maximum absolute value of the DC bus voltage deviation over the past three sampling periods; k u The quantization factor for the DC bus voltage deviation, when ΔU dc_max (k) takes effect when it is not 0, realizing adaptive adjustment of the fuzzy control domain, improving the adjustment accuracy when the deviation is small, and improving the response speed when the deviation is large; The model predictive control generates a power feedforward based on the optimal grid-connected power prediction value, providing feedforward compensation for PID control. The PID control generates a speed safety constraint based on the flywheel speed safety boundary, which in turn corrects the constraint conditions of the model predictive control. The specific details are as follows: The feedforward correction formula is: ； Where, △ω f_ff (k) represents the speed feedforward compensation amount for the kth sampling period; P g_opt (k+1) represents the optimal grid-connected power prediction value output by S4; J f (k) represents the latest value of the flywheel moment of inertia obtained from S3 parameter identification; T s Sampling period Based on the feedforward compensation, the input speed deviation of the PID control is corrected. The corrected speed deviation formula is as follows: △ω f_corr (k)=△ω f (k)+△ω f_ff (k); Where, △ω f_corr (k) represents the corrected rotational speed deviation; ω f (k) represents the real-time rotational speed of the flywheel during the kth sampling period; In the reverse interaction, the PID control corrects the constraints of the model predictive control. The constraint correction formula is as follows: oh f_min_safe (k)=max(ω f_min ,oh f_filt (k)+△ω f_pid_min (k)); oh f_max_safe (k)=min(ω f_max ,oh f_filt (k)+△ω f_pid_max (k)); Where, ω f_min_safe (k), ω f_max_safe (k) represents the upper and lower limits of the modified flywheel speed safety constraint; △ω f_pid_min (k), △ω f_pid_max (k) represents the upper and lower safety limits of the rotational speed output by the PID control; ω f_min ω is the minimum permissible operating speed of the flywheel. f_max This refers to the maximum permissible operating speed of the flywheel; Simultaneously, the charging and discharging power constraint boundary calculated in real time on the machine side is fed back into the grid-connected power prediction and collaborative control algorithm to correct the feasible region of the grid-side power reference trajectory; specifically, it includes the following: The corrected formula is: P g_max_corr (k)=min(P g_max ,P f_max (k)·h f ·or g ); P g_min_corr (k)=max(P g_min ,P f_min (k)·h f ·or g ); In the formula, P g_max_corr (k) and P g_min_corr (k) represents the modified maximum and minimum allowable grid-connected power; P g_max P g_min P represents the upper and lower limits of the rated power of the grid-side converter; f_max (k), P f_min (k) represents the upper and lower limits of the flywheel charging and discharging power constraint calculated in real time on the machine side; η g η is the rated efficiency of the grid-side converter. f This refers to the rated efficiency of the machine-side converter.

9. The flywheel energy storage grid frequency support method based on model predictive control according to claim 1, characterized in that: In step S5, the input variables of the fuzzy control are the DC bus voltage deviation and the flywheel speed deviation. The fuzzy subsets of both input variables include: negative large, negative medium, negative small, zero, positive small, positive medium, and positive large. The output variables are the weighting coefficients of the model predictive control and the weighting coefficients of the PID control. The sum of the two weighting coefficients is always 1. Based on a 7×7 combination design of two input variables, 49 fuzzy rules were designed. The Mamdani inference method was used to complete the fuzzy inference, and the centroid method was used to complete the defuzzification, so as to obtain the accurate weight coefficients.

10. The flywheel energy storage grid frequency support method based on model predictive control according to claim 2, characterized in that: In step S6, the actual operating data includes the actual filtered values ​​of grid frequency, DC bus voltage and flywheel speed, as well as the actual grid-connected power data output after S4 is executed, and the actual charging and discharging power data output by the flywheel motor after S5 is executed; the actual deviations include: actual deviation of grid frequency, actual deviation of DC bus voltage, actual deviation of flywheel speed, grid-side model prediction error and machine-side model prediction error; The specific rules for updating parameters and weights are as follows: If the actual deviation of the power grid frequency exceeds the preset target deviation and the sampling period is continuously set, the model parameters of the dual-model system are updated using the proportional adjustment method, and the weight domains of the model predictive control and the PID control are adjusted at the same time. The methods for updating network-side model parameters include: 1) Increase the virtual inertia and virtual damping coefficient of the virtual synchronous generator to enhance the system's inertial support capability, suppress rapid frequency fluctuations, and accelerate frequency deviation convergence; 2) Increase the output tracking weight in the model predictive control cost function to improve the tracking accuracy of the power reference trajectory and accelerate the convergence speed of the frequency deviation; 3) Increase the bandwidth of the disturbance observer to improve the response speed of disturbance observation and quickly correct model errors caused by large disturbances on the power grid side; Methods for updating aircraft-side model parameters include: 1) Reduce the forgetting factor of recursive least squares method, improve the tracking speed of parameter identification for model parameter changes, and quickly adapt to model parameter fluctuations caused by changes in working conditions. 2) The proportional adjustment method is used to directly correct the core parameters of the machine-side model and quickly compensate for the model prediction error; The weighted universe of discourse method for machine-side fuzzy control includes: 1) When the power grid frequency deviation is large and it is necessary to prioritize ensuring the power grid frequency support performance; a mpc_max (k+1)=min(1.0,a mpc_max (k)+k α_up ·E(k); α pid_max (k+1)=max(0.0,α pid_max (k)-k α_down ·E(k)); Where, α mpc_max (k+1) represents the upper bound of the universe of discourse for the updated model prediction control weight coefficients; α pid_max (k+1) represents the upper bound of the universe of discourse for the updated PID control weight coefficients; k α_up and k α_down The weighted universe is used to adjust the coefficients; α is always satisfied. mpc_max (k)+α pid_max The constraint is (k)=1; E(k) is the power grid frequency deviation evaluation index. E(k) = |Δf g_filt (k)|,△f g_filt (k) represents the power grid frequency deviation filter value in the kth sampling period; 2) When the flywheel speed deviation is large and close to the safety boundary, and it is necessary to prioritize equipment safety: a pid_max (k+1)=min(1.0,a pid_max (k)+k α_pid_up ·E ω (k)); α mpc_max (k+1)=max(0.0,α mpc_max (k)-k α_mpc_down ·E ω (k)); Where, k α_pid_up k α_mpc_down E is the weighting universe adjustment coefficient, ranging from 0.1 to 0.3, always satisfying the constraint that the sum of the weighting coefficients is 1; ω (k) represents the flywheel speed tracking error index. E ω (k)=|ω f_ref (k)-ω f_filt (k)|; Where, ω f_ref (k) is the reference value for flywheel speed; ω f_filt (k) is the filtered value of the flywheel real-time speed in the kth sampling period output by S2; If the actual deviation of the power grid frequency is within the preset target deviation range, only the model parameters are fine-tuned or the parameters are kept unchanged; The termination condition for the closed-loop regulation is set as follows: within 100 consecutive acquisition cycles, the actual deviation of the power grid frequency is within ±0.05 Hz and the actual deviation of the power grid frequency does not exceed the target deviation of 0.02 Hz. After the termination condition is met, the adaptive update of parameters and weights is stopped, the current control strategy is maintained, and the closed-loop regulation is restarted until the power grid frequency deviates from the rated range again.

Images