Optimal scheduling method for wind-solar-thermal storage multi-energy system based on adaptive dynamic programming

Abstract

The present application relates to the technical field of energy optimization scheduling, and specifically proposes a wind-solar-thermal-storage multi-energy system optimization scheduling method based on adaptive dynamic programming. The core of the method is to reversely utilize the coupling effect of rotor thermal stress "heat conduction lag-morphology difference": first, the morphological deviation of the node virtual temperature and the reference sequence is used to complete the heat risk trend quantification; then, the risk adjustment coefficient, the battery polarization impedance and the state of charge constraint term are extracted under the power regulation condition, and the battery inhibition weight is calculated by means of normalization and feature fusion. The system automatically distributes the scheduling instructions to the battery and the solar thermal power station by embedding the dynamic weight into the quadratic programming objective function. The present application takes "radial temperature field deduction-thermal expansion stress prediction-morphological deviation quantification" as the basis for judgment, eliminates the blind area of heat conduction lag caused by single-point threshold, guarantees the accuracy and response speed under the second-level rolling optimization, and realizes the collaborative optimization of rotor thermal impact inhibition, battery aging delay and unit climbing cost.

Description

Technical Field

[0001] This invention relates to the field of energy optimization scheduling technology, specifically to an optimization scheduling method for a multi-energy system based on adaptive dynamic programming, encompassing wind, solar, thermal, and energy storage. Background Technology

[0002] In new power systems dominated by new energy sources, concentrated solar power (CSP) plants, with their built-in thermal storage and steam turbine generators, combine the dispatchability and rotational inertia of thermal power plants, becoming the foundation for wind and solar energy integration. Electrochemical energy storage, with its millisecond-level response speed, works alongside CSPs to mitigate the random fluctuations in new energy sources. Consequently, optimized scheduling methods for multi-energy systems involving wind, solar, thermal, and energy storage have gradually become a research hotspot.

[0003] Existing scheduling technologies only set a fixed ramp rate to constrain turbine output. However, in actual operation, the peak value of rotor thermal stress lags significantly later than power changes due to the lag in radial heat conduction. This spatiotemporal misalignment not only stems from the thermal inertia of the large mass of metal components but is also affected by multiple factors such as the transient change in main steam temperature, ambient temperature fluctuations, and the rate characteristics of electrochemical energy storage. Furthermore, existing wind, solar, thermal, and energy storage systems generally rely on the current wall temperature or SOC single-point threshold to trigger power allocation. Since there are differences in the shape of battery polarization and thermal stress response curves, and the optimization weights depend only on whether the limit is exceeded rather than the future stress trend, the fixed battery suppression weights cannot accurately reflect the local evolution of thermal risks, resulting in a lag in response time during multi-energy optimization scheduling of wind, solar, thermal, and energy storage. Summary of the Invention

[0004] To address the technical problem that the peak value of rotor thermal stress lags significantly later than power change due to radial heat conduction lag, leading to reduced optimization weight accuracy and consequently, a lag in response time during multi-energy system optimization scheduling involving wind, solar, thermal, and energy storage, this invention aims to provide an adaptive dynamic programming-based optimization scheduling method for multi-energy systems involving wind, solar, thermal, and energy storage. The specific technical solution adopted is as follows: This invention proposes an optimal scheduling method for a multi-energy system (wind, solar, thermal, and energy storage) based on adaptive dynamic programming. The method includes: In the optimized scheduling of energy, the rotor node temperature, steam temperature, battery current, and remaining battery power are obtained; The initial temperature of a node is obtained by performing a shutdown-startup attenuation analysis on the node temperature; the corrected temperature of the node is obtained by correcting the initial temperature of the node using the ambient temperature; and the virtual temperature of the node is obtained by performing a recursive fusion analysis on the corrected temperature using the steam temperature. By analyzing the differences in the distribution of virtual temperature at nodes, the predicted value of expansion stress is obtained; based on the difference between the predicted value of expansion stress and the preset benchmark value, the risk adjustment coefficient is obtained. By analyzing battery polarization loss, the polarization voltage difference is obtained; by the difference between the polarization voltage difference and the battery current, the battery impedance benchmark term is obtained; and based on the decay penalty of the remaining battery capacity, the charge constraint term is obtained. Calculate the characteristic fusion value of the risk adjustment coefficient, battery impedance benchmark term, and normalized value of charge constraint term to obtain the battery suppression weight; construct an operation loss assessment function, which includes a battery power term and a solar thermal power plant power change term; use the battery suppression weight as the weight of the battery power term in the operation loss assessment function; optimize the scheduling of the energy system based on the output of the operation loss assessment function.

[0005] Furthermore, the method for obtaining the initial temperature includes: The current shutdown / start status is determined by a preset threshold. When the shutdown duration is greater than or equal to the preset threshold, it is determined to be a cold start, and the initial temperature of the node is the ambient temperature. When the shutdown duration is less than the preset threshold, it is determined to be a hot start. The node's initial temperature is obtained by performing a decay analysis on the node based on the shutdown duration.

[0006] Furthermore, the method for obtaining the attenuation analysis includes: The difference between the node's temperature at the moment of shutdown and the ambient temperature is calculated, and the natural cooling temperature is obtained by negatively weighting the temperature by the shutdown duration. Based on the feature fusion of the natural cooling temperature and the ambient temperature, the attenuation analysis results of the node temperature cooling by the shutdown duration are obtained.

[0007] Furthermore, the method for obtaining the node correction temperature includes: The difference between the current ambient temperature and the predicted temperature on the rotor surface is calculated and weighted using a preset temperature weight to obtain the temperature correction difference at the current moment; the initial temperature is corrected based on the temperature correction difference to obtain the corrected temperature at the current moment.

[0008] Furthermore, the method for obtaining the node's virtual temperature includes: The corrected temperature of each node at the previous time step is weighted and fused with the steam temperature at the current time step according to a preset weighting coefficient, and then recursively analyzed to obtain the virtual node temperature of each node at different times.

[0009] Furthermore, the method for obtaining the predicted expansion stress value includes: Based on the difference between the average virtual temperature of all nodes at the current moment and the virtual temperature of the node at the rotor boundary node, and corrected by the material adjustment coefficient, the predicted value of the rotor's expansion stress is obtained.

[0010] Furthermore, the method for obtaining the risk adjustment coefficient includes: Based on the difference between the predicted value of rotor expansion stress and the preset benchmark value at different times, the coupling matrix is ​​obtained; through dynamic path finding of the coupling matrix, all path trajectories are obtained; based on the maximum value in each path trajectory, the minimum value among all maximum values ​​is obtained and normalized to obtain the trend deviation value of the control cycle; through weighted negative correlation analysis of the trend deviation value, the risk adjustment coefficient at the current time is obtained.

[0011] Furthermore, the method for obtaining the polarization voltage difference includes: The polarization loss of lithium ions in the electrode material is analyzed by considering the differences between the energy storage battery voltage, open circuit voltage, and internal resistance voltage, and the polarization voltage difference of the battery is obtained.

[0012] Furthermore, the method for correcting the degradation of the remaining battery capacity includes: Based on the difference between the remaining battery capacity and the preset threshold, an exponential decay penalty for overcharging and over-discharging is applied to obtain the charge constraint term.

[0013] Furthermore, the operating loss evaluation function includes: in, This represents the result of the current control cycle operating loss assessment function; This indicates the battery suppression weight for the current control cycle; This represents the battery power calculated in the current control cycle; This represents the preset smoothing coefficient, which is set to 100 in this embodiment; This indicates the rate of change in the power output of the solar thermal power plant during the current control period.

[0014] The present invention has the following beneficial effects: This invention provides a data foundation for multi-energy optimized scheduling by acquiring rotor and battery operating status parameters in real time. By performing shutdown-and-restart attenuation analysis on node temperatures and correcting for ambient temperature, the corrected radial temperature of the rotor can be reconstructed from two dimensions: cooling time and residual heat. Combined with a recursive analysis of the influence of steam temperature on rotor nodes, the virtual node temperature is calculated, accurately representing the future temperature state of the rotor nodes after startup. By analyzing the radial changes in the virtual node temperatures of different nodes at the same time, the predicted expansion stress value is calculated, effectively representing the degree of thermal expansion stress change caused by different material expansion changes at different temperatures on the rotor surface. Based on the difference between the predicted expansion stress value and the preset benchmark value, a risk adjustment coefficient is calculated. This step accurately analyzes the impact of ambient temperature on rotor temperature changes, representing the potential risk of rotor damage due to uneven metal expansion. By calculating the battery impedance benchmark term through the difference between the battery electrode voltage loss and battery current, the degree of polarization resistance growth can be quantified, accurately reflecting the risk of lithium-ion diffusion obstruction and accelerated aging. Combining the risk adjustment coefficient, battery impedance benchmark term, and charge constraint term, the battery suppression weight is calculated, accurately analyzing the degree of thermal and electrical coupling degradation. By applying a weighted penalty to the battery power term in the operating loss assessment function through battery suppression weights, and dynamically adjusting the operating loss assessment function, the optimal control command for the control cycle can be obtained through dynamic optimization of weights when changes occur in both thermal and electrical dimensions. This ensures that the battery yields in advance before the peak of thermal stress, achieving optimal lifespan and tracking accuracy in second-level rolling, thereby improving the response speed and reducing the response time when optimizing the scheduling of multiple energy sources such as wind, solar, thermal, and energy storage. Attached Figure Description

[0015] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 The flowchart illustrates an optimized scheduling method for a multi-energy system based on adaptive dynamic programming, encompassing wind, solar, thermal, and energy storage, as provided in one embodiment of the present invention. Detailed Implementation

[0017] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the following, in conjunction with the accompanying drawings and preferred embodiments, details the specific implementation, structure, features, and effects of an adaptive dynamic programming-based multi-energy system scheduling method for wind, solar, thermal, and energy storage systems proposed according to the present invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0018] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0019] The following description, in conjunction with the accompanying drawings, details the specific scheme of the wind-solar-thermal-storage multi-energy system optimization scheduling method provided by the present invention.

[0020] Please see Figure 1 The diagram shows a flowchart of an optimal scheduling method for a multi-energy system based on adaptive dynamic programming, which is provided by an embodiment of the present invention.

[0021] Step S101: In the energy optimization scheduling, obtain the rotor node temperature, steam temperature, battery current, and remaining battery power.

[0022] In wind-solar-thermal-storage systems, the thermal processes of a solar thermal power plant exhibit large inertia on the order of seconds, while the electrical response processes of an energy storage battery management system have rapid characteristics on the order of milliseconds. These two systems differ significantly in sampling frequency and response scale. To achieve coordinated control within a unified time frame, data from different frequencies needs to be sampled and aligned. In this embodiment, the control cycle... For data acquisition, firstly, the main steam temperature and turbine metal wall temperature (hereinafter referred to as ambient temperature) of the control cycle are directly read through the distributed control system (DCS) interface of the solar thermal power plant, in degrees Celsius. Secondly, the battery voltage and current are collected in each cycle via the data bus of the Battery Management System (BMS). The collected current and voltage data are arranged in the order of acquisition to obtain voltage and current sequences. Next, the arithmetic mean of the voltage and current sequences is calculated using arithmetic mean filtering, which serves as the battery terminal voltage and battery current for that control cycle. Then, the remaining battery capacity is read via the BMS data bus. Finally, the grid-connected power for the current control cycle is obtained from the grid dispatch terminal, and the real-time output power for the current control cycle is obtained from the wind and solar inverters. Based on the power balance principle, the net power gap value that the system needs to adjust is obtained by subtracting the real-time output power from the grid-connected power for the current control cycle. This net power gap value can be used to solve the constraints of the operating loss assessment function in subsequent modules.

[0023] In the initial state, the node temperature of each node of the rotor is the ambient temperature.

[0024] Step S102: Obtain the initial temperature of the node by performing a shutdown-startup attenuation analysis on the node temperature; correct the initial temperature of the node by using the ambient temperature to obtain the corrected temperature of the node; and obtain the virtual temperature of the node by performing a recursive fusion analysis on the corrected temperature using the steam temperature.

[0025] During the operation of a solar thermal power plant, the temperature of the turbine rotor inside the turbine typically differs significantly from the ambient temperature. After the turbine stops operating, it undergoes natural cooling. Due to the closed nature of the turbine, heat exchange occurs between the internal components and the metal walls, while the metal walls exchange heat with the outside air; the efficiencies of these two heat exchange methods differ. During startup, the turbine interior is not fully cooled, resulting in a turbine temperature higher than the ambient temperature. Furthermore, the turbine's cooling rate decays exponentially with the shutdown duration (i.e., cooling time). Therefore, by analyzing the decay of the rotor node temperature at the shutdown moment with the shutdown-startup duration, the initial temperature of the current cycle node can be obtained. A higher initial temperature indicates a shorter turbine cooling time (shutdown time).

[0026] Due to parameter errors or external disturbances, the calculated internal temperature distribution of the rotor may gradually deviate from the actual physical state. To ensure the accuracy of the rotor node temperatures, it is necessary to correct the initial temperature of the rotor. Since a steam turbine is a relatively enclosed space, there is a strong correlation between the ambient temperature and the rotor temperature. Therefore, by using the difference between the ambient temperature and the rotor temperature, the initial temperature of each rotor node is corrected to obtain the corrected node temperature. The smaller the difference between this corrected value and the initial temperature, the less affected the rotor node temperature is by external disturbances.

[0027] Traditional methods assess risk solely based on current temperature or stress values, neglecting the inherent hysteresis of thermal stress as a product of large-inertia thermodynamic processes. Predicting future rotor node temperatures reflects the true dynamics of thermal inertia, enabling tracking of the rotor's internal thermal state. Furthermore, since steam temperature directly convects with the rotor surface, it determines the temperature of each rotor node. Therefore, a recursive fusion analysis is performed using the current steam temperature and the correction temperature from the previous control cycle to obtain the virtual node temperature for the next N control cycles. The predicted virtual node temperatures provide a complete temperature field for subsequent thermal stress prediction, improving the accuracy of thermal stress analysis. In this embodiment, N represents the number of predicted control cycles, which can be 60.

[0028] Step S103: Obtain the predicted value of expansion stress by the distribution difference of virtual temperature at the nodes; obtain the risk adjustment coefficient based on the difference between the predicted value of expansion stress and the preset benchmark value.

[0029] Due to internal temperature fluctuations, instantaneous temperature differences arise at each radial node due to varying thermal inertia. The rotor boundary region is thinner than the center, causing the rotor boundary to heat up faster under steam impact, while the center lags behind, creating a temperature gradient that is higher on the outside and lower on the inside. This gradient causes uneven expansion and contraction of the metal fibers; the surface desires to expand but is constrained internally, thus generating thermal expansion stress. Therefore, by analyzing the distribution of virtual temperatures at the nodes, the predicted value of the rotor's expansion stress is obtained. A larger value indicates a more severe temperature deviation across different regions of the rotor surface, a higher risk of thermal shock, and a higher probability of fatigue crack initiation.

[0030] Next, the predicted values ​​of expansion stress for the next N control cycles will be calculated and arranged in chronological order to obtain the expansion stress prediction sequence, which will be used to characterize the future changes in rotor thermal expansion stress.

[0031] Furthermore, traditional methods for analyzing thermal expansion stress only focus on whether the stress exceeds a fixed numerical threshold, making it difficult to distinguish between instantaneous spikes and continuous deterioration—two risk modes. Therefore, a baseline sequence of length N is established, and the risk adjustment coefficient is obtained by comparing the predicted expansion stress value with the preset baseline value. A larger value indicates a smaller deviation between the predicted and ideal states, and a lower risk of damage caused by the rotor's thermal expansion stress. In this embodiment, all elements of the baseline sequence are set to 0, representing the ideal, safest state where the rotor has no thermal expansion stress.

[0032] Step S104: Analyze the battery polarization loss to obtain the polarization voltage difference; obtain the battery impedance reference term based on the difference between the polarization voltage difference and the battery current; obtain the charge constraint term based on the decay penalty of the remaining battery capacity.

[0033] For energy storage batteries, the battery terminal voltage comprises three parts: open-circuit voltage, ohmic voltage drop, and polarization voltage difference. Only the polarization voltage difference can reflect the obstructed diffusion of lithium ions at the electrode / electrolyte interface and the risk of lithium plating, thus characterizing the battery's polarization loss. Therefore, the polarization voltage difference can be obtained by comparing the battery voltage with the open-circuit voltage and ohmic voltage drop. A larger value indicates a larger concentration gradient at the electrode surface, slower ion transport, a smaller usable power margin, and a faster battery aging rate.

[0034] Furthermore, since the battery's internal resistance affects its discharge process, a higher degree of internal polarization indicates higher resistance and a greater impact on discharge. Therefore, a battery impedance benchmark is obtained by analyzing the difference between the polarization voltage difference and the battery current. A larger benchmark value indicates more severe internal polarization, smaller usable power margin, and faster aging. In such cases, the depth of charge and discharge should be actively limited during scheduling to extend cycle life.

[0035] Next, because the active materials undergo drastic phase transitions, internal resistance increases sharply, and aging accelerates nonlinearly in the battery's extremely high or low SOC range, the battery's discharge intensity can be evaluated based on its remaining state of charge. Therefore, a charge constraint term is obtained based on the degradation penalty of the remaining battery capacity. A larger value indicates a smaller charge / discharge margin for the battery, and a higher safety risk and aging rate.

[0036] Step S105: Calculate the characteristic fusion value of the risk adjustment coefficient, battery impedance benchmark term, and normalized value of charge constraint term to obtain the battery suppression weight; construct the operation loss assessment function, which includes the battery power term and the power change term of the solar thermal power plant; use the battery suppression weight as the weight of the battery power term in the operation loss assessment function; optimize the scheduling of the energy system based on the output of the operation loss assessment function.

[0037] In the optimization of the energy system, multiple constraints need to be considered. In the quadratic programming objective function, the losses of the battery and solar thermal components differ significantly in magnitude. Directly using the battery impedance benchmark and charge constraint terms may be ignored or locked by the optimizer. Therefore, for N control cycles, the risk adjustment coefficient, battery impedance benchmark, and charge constraint are normalized respectively. The normalized risk adjustment coefficient, battery impedance benchmark, and charge constraint are then fused to obtain the battery suppression weight. A larger value indicates a higher current battery operating cost, and the system should actively reduce battery output, prioritizing solar thermal regulation to extend electrochemical life and avoid overcharging and over-discharging. In this embodiment, maximum value normalization can be used.

[0038] Furthermore, regarding the construction of the operation loss assessment function, the operation loss assessment function includes a battery power term and a solar thermal power plant power change term. By using the battery suppression weight as the weight of the battery power term in the operation loss assessment function, the operation loss assessment function is adjusted for different control cycles.

[0039] After renewable energy is connected to the grid, real-time source-load imbalance can easily lead to frequency overruns. Dispatch commands must be rigidly executed, and any economic or lifetime optimization cannot come at the expense of power balance; the net power deficit requirement of the grid must be met. Therefore, in a specific implementation of this invention, the power balance constraint is expressed by the formula: In the formula, This represents the battery power calculated in the current control cycle; This indicates the output power of the solar thermal power plant during the current control cycle; This represents the slack variable in the power balance equation; This represents the net power deficit value for the current control cycle. , This represents the minimum and maximum output power, which in this embodiment can be 0 and the battery's rated output power, respectively.

[0040] In one specific implementation of this invention, the objective function is calculated using quadratic programming (QP). The quadratic term matrix and the linear term vector, and the objective function The constraints are input into the controller's built-in QP solver (including but not limited to the effective set method or interior point method solver), and the solver will automatically find a set of optimal solutions. and Next, the calculated optimal battery power command will be... The command is sent to the power storage converter (PCS) via the communication interface to control the charging and discharging of the battery; the optimal solar thermal power command is then transmitted. The digital electro-hydraulic control system (DEH) sent to the solar thermal power plant controls the opening of the turbine regulating valve.

[0041] Finally, the initial temperature of the rotor radial node calculated in the current control cycle is saved to non-volatile memory. This memory is then called when calculating the initial temperature in the next control cycle, ensuring continuous and accurate calculation of the unmeasurable temperature field inside the rotor.

[0042] This invention provides a data foundation for multi-energy optimized scheduling by acquiring rotor and battery operating status parameters in real time. By performing shutdown-and-restart attenuation analysis on node temperatures and correcting for ambient temperature, the corrected radial temperature of the rotor can be reconstructed from two dimensions: cooling time and residual heat. Combined with a recursive analysis of the influence of steam temperature on rotor nodes, the virtual node temperature is calculated, accurately representing the future temperature state of the rotor nodes after startup. By analyzing the radial changes in the virtual node temperatures of different nodes at the same time, the predicted expansion stress value is calculated, effectively representing the degree of thermal expansion stress change caused by different material expansion changes at different temperatures on the rotor surface. Based on the difference between the predicted expansion stress value and the preset benchmark value, a risk adjustment coefficient is calculated. This step accurately analyzes the impact of ambient temperature on rotor temperature changes, representing the potential risk of rotor damage due to uneven metal expansion. By calculating the battery impedance benchmark term through the difference between the battery electrode voltage loss and battery current, the degree of polarization resistance growth can be quantified, accurately reflecting the risk of lithium-ion diffusion obstruction and accelerated aging. Combining the risk adjustment coefficient, battery impedance benchmark term, and charge constraint term, the battery suppression weight is calculated, accurately analyzing the degree of thermal and electrical coupling degradation. By applying a weighted penalty to the battery power term in the operating loss assessment function through battery suppression weights, and dynamically adjusting the operating loss assessment function, the optimal control command for the control cycle can be obtained through dynamic optimization of weights when changes occur in both thermal and electrical dimensions. This ensures that the battery yields in advance before the peak of thermal stress, achieving optimal lifespan and tracking accuracy in second-level rolling, thereby improving the response speed and reducing the response time when optimizing the scheduling of multiple energy sources such as wind, solar, thermal, and energy storage.

[0043] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the initial temperature includes: For the turbine rotor, the boundary to the center is equally divided into M-1 segments along the rotor radius, resulting in M ​​nodes. Since turbine startup is divided into cold start and hot start, the shutdown duration between the current startup time and the last shutdown time is calculated. Here, M represents the number of nodes, which can be 20 in this embodiment.

[0044] Furthermore, when the shutdown duration is greater than or equal to a preset threshold, it indicates that the current startup is a cold start. The ambient temperature is then read as the initial temperature of each node on the rotor. The initial temperatures of the nodes are then distributed from the rotor boundary to the center, constructing an initial temperature sequence of length M. In this embodiment, the preset threshold can be 72 hours, and the ambient temperature is the temperature of the outer metal wall of the turbine.

[0045] When the downtime is less than a preset threshold, it indicates that the current startup is a hot start. At this time, residual temperature remains inside the rotor. The initial temperature of the node is obtained by performing a temperature decay analysis on the node based on the downtime. The larger the value, the higher the residual temperature inside the rotor.

[0046] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the attenuation analysis includes: In estimating the initial temperature of a turbine rotor during hot start-up, the node temperature is susceptible to the combined effects of shutdown duration and natural convection attenuation. Temperature differences deviate from capturing metal heat storage, and shutdown duration deviates from quantifying cooling depth; both jointly map the synchronous decay of shape and value under the coupling of "heat storage-heat dissipation." Therefore, this invention uses the difference between the node shutdown temperature and the ambient temperature, which characterizes the degree of influence of ambient temperature on rotor temperature, and the shutdown duration, which characterizes the superimposed duration of rotor cooling and heat dissipation and shows a negative correlation with node temperature. Thus, the two results are weighted and fused to obtain the natural cooling temperature; subsequently, the natural cooling temperature is fused with the characteristics of the ambient temperature. The smaller this value, the more exponentially the node heat storage has dissipated with the shutdown duration, the closer the rotor radial temperature field is to the environment, and the safer the initial temperature setting for hot start-up, effectively avoiding thermal shock calculation errors.

[0047] In one specific implementation of this invention, the attenuation analysis is expressed by the formula: In the formula, This represents the initial temperature of the i-th node; Indicates ambient temperature; This represents the temperature of the i-th node at the time of the previous shutdown; The rotor's natural cooling coefficient is measured experimentally. Indicates the downtime; Represents an exponential function with the natural constant as its base; This indicates the natural cooling temperature.

[0048] In the above formula, the difference between the temperature at the node shutdown time and the ambient temperature is represented by the difference value. The negative correlation is mapped using an exponential function. The final natural cooling temperature is used as the adjustment amount of the ambient temperature. Feature fusion is achieved through addition.

[0049] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the node correction temperature includes: In the online simulation of thermal stress in a high-inertia steam turbine, the internal temperature of the rotor is affected by both transient shocks from the main steam and contamination from model parameter drift and unknown heat dissipation disturbances. Therefore, the difference between the ambient temperature of the current control cycle and the predicted temperature of the rotor boundary nodes of the previous control cycle is calculated. This difference characterizes the deviation between the predicted boundary temperature and the actual external heat dissipation. This difference is then weighted and amplified using preset temperature weights to obtain the temperature correction difference for the current control cycle. This temperature correction difference characterizes the correction strength of the node temperature. The temperature correction difference is then superimposed onto the initial temperature of each rotor node to obtain the corrected temperature at the current moment. The larger the value, the more significant the deviation of the model from the actual boundary, requiring a stronger pull-back force to ensure the accuracy of thermal stress prediction.

[0050] It should be noted that the rotor boundary predicted temperature is calculated as follows: the collected rotor boundary temperatures are arranged in the order of the control cycles, and the arranged data is used as the output of the EMA (Exponential Moving Average) algorithm, which outputs the predicted temperature of the rotor boundary nodes. The calculation of the EMA algorithm is a well-known technique, and the specific calculation method will not be elaborated here.

[0051] In one specific implementation of this invention, the correction temperature is expressed by the formula: In the formula, This represents the corrected temperature of the i-th node; This represents the initial temperature of the i-th node; This represents the correction equity coefficient, which can be 0.3 in this embodiment; This represents the difference between the ambient temperature of the current control cycle and the predicted temperature of the rotor boundary node in the previous control cycle. This indicates the temperature correction difference for the current control cycle.

[0052] In the above formula, the two features are positively fused by addition.

[0053] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the virtual temperature of a node includes: In the rolling prediction of the temperature field of a high-inertia steam turbine, the virtual node temperature is affected by both the residual metal heat storage from the previous cycle's correction temperature and the impact of new heat flow brought about by the current steam boundary transient. Therefore, the correction temperatures of each node at the previous moment are fused with the steam temperature at the current moment according to a preset weighting coefficient, and recursive analysis is performed to obtain the virtual node temperature of each node at different moments. The larger the value, the more significant the heat accumulation inside the rotor, the higher the risk of subsequent thermal stress peaks, and the stronger the warning level for thermal safety assessment.

[0054] In one specific implementation of this invention, the node virtual temperature is expressed by the formula: In the formula, , These represent the node virtual temperature sequences predicted for the k-th and (k-1)-th control cycles, respectively. Indicates the steam temperature; , These represent the coefficient matrix and coefficient sequence of the rotor material, respectively.

[0055] In the above formula, the two features are positively fused using a product form.

[0056] It should be noted that, , This represents the coefficient matrix and coefficient sequence of each parameter in the discretized one-dimensional unsteady heat conduction differential equation; in this embodiment, the calculation method can be as follows: the stator and rotor radii are... Discretized into If there are 1 node, then the space step size is 1. The virtual simulation time step dt, in this embodiment, uses a control cycle of 1 second, and the thermal diffusivity of the rotor material is a. The Fourier number is defined. Then the coefficient matrix is ​​a diagonal matrix, where the values ​​of each element are: elements on the main diagonal. ; lower diagonal element Top diagonal elements For the coefficient sequence, only the last element is non-zero, and the rest are 0. The value of the last element is... (where h is the surface convection heat transfer coefficient, (where is the thermal conductivity of the material, determined experimentally). The specific discretization principle of the one-dimensional unsteady-state thermal conductivity differential equation is a well-known technique in the field and will not be elaborated here.

[0057] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the predicted value of expansion stress includes: In the online assessment of thermal stress in high-order steam turbines, rotor expansion stress is influenced by both the macroscopic heat capacity determined by the global average temperature and the radial constraint impact generated by the temperature difference between boundary nodes and the center. Therefore, based on the difference between the current virtual average temperature of all nodes and the virtual temperature of the rotor boundary nodes, this difference can characterize the risk of thermal stress concentration caused by uneven thermal expansion inside and outside the rotor. After correction and amplification by a material adjustment coefficient, the predicted value of rotor expansion stress is obtained; the larger the value, the more severe the radial temperature difference, the higher the degree of thermal expansion restriction, and the higher the risk assessment of rotor thermal fatigue damage. In a specific implementation of this invention, the predicted expansion stress value is expressed by the formula: In the formula, This represents the predicted value of the expansion stress for the k-th control cycle. Indicates the elastic modulus of a material; Indicates the coefficient of linear expansion; Indicates Poisson's ratio; This represents the mean of all elements in the virtual temperature sequence of the k-th control cycle node; This represents the predicted virtual node temperature of the rotor boundary node in the k-th control cycle, using the first element in the virtual node temperature sequence. This indicates the material adjustment factor.

[0058] In the above formula, the temperature difference is corrected and amplified by using a material adjustment coefficient in a product form.

[0059] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the risk adjustment coefficient includes: In high-penetration new energy scenarios, rotor thermal stress is affected by both the amplitude impact of transient main steam changes and the morphological drift disturbance caused by the hysteresis of metal thermal conduction. Ideally, the internal temperature of the turbine should gradually converge and stabilize after the control process, and the thermal stress should gradually decrease and tend to calm down over time. The greater the deviation between the predicted state and the ideal state, the higher the risk of thermal fatigue. Therefore, a coupling matrix is ​​constructed by the difference between the predicted rotor expansion stress value and the preset benchmark value at different times. All path trajectories are extracted through dynamic path finding, and normalized using the minimum value among the maximum values ​​of each path to obtain the trend deviation value of the control cycle. This trend deviation value can characterize the degree to which the predicted thermal stress curve deviates from the safe decay pattern. The more significant the deviation of the predicted trajectory from the safe decay pattern, the higher the risk of thermal fatigue. Then, a weighted negative correlation mapping is applied to this deviation value to obtain the risk adjustment coefficient. The larger the value, the lower the unit price of battery output, indicating that batteries can be used to share the pressure of the turbine.

[0060] As an example, in one specific implementation of this invention, the risk adjustment coefficient is expressed by the formula: In the formula, This represents the risk adjustment coefficient for the current control cycle; This represents the preset thermoelectric coupling gain coefficient, which can be 0.5 in this embodiment; This represents the trend deviation value of the current control cycle. In the above formula, a reciprocal form is used for negative correlation mapping. To avoid calculation errors and ensure the subsequent algorithm runs reasonably, a parameter of 1 is added to the denominator to ensure the formula has calculation meaning.

[0061] In the above formula, the two features are negatively fused using the reciprocal form.

[0062] As an example, in a specific implementation of this invention, the expansion stress prediction sequence and a preset benchmark sequence are used as inputs to the Discrete Fréchet distance algorithm, and the algorithm outputs the trend deviation value of the control period. In this embodiment, the Discrete Fréchet distance algorithm is calculated as follows: First, a coupling matrix is ​​constructed using the Discrete Fréchet distance between the expansion stress prediction sequence and the preset benchmark sequence; then, a dynamic programming (DP) algorithm is used to obtain all path trajectories between the lower left corner (1,1) and the upper right corner (N,N) of the coupling matrix; next, the maximum value among each path trajectories is obtained, and the minimum value among all maximum values ​​is obtained. The minimum value is then normalized by comparing it with the maximum value in the expansion stress prediction sequence to obtain the trend deviation value of the control period. The calculation of the Discrete Fréchet distance algorithm and the dynamic programming algorithm are well-known techniques, and the specific calculations will not be described in detail here.

[0063] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the polarization voltage difference includes: Under high-rate charge / discharge conditions in energy storage systems, the battery terminal voltage is affected by both the instantaneous voltage drop due to the ohmic internal resistance and the drift of the open-circuit potential corresponding to the SOC. The superposition of these two factors can easily mask the true diffusion resistance of lithium ions at the electrode / electrolyte interface. Therefore, by measuring the terminal voltage, the open-circuit voltage and internal resistance voltage are separated. That is, the polarization loss of lithium ions in the electrode material is analyzed based on the difference between the energy storage battery voltage, open-circuit voltage, and internal resistance voltage, obtaining the battery's polarization voltage difference. The larger the value, the more significant the polarization loss, the lower the usable power margin, and the higher the risk of battery aging. Therefore, the open-circuit voltage of the battery is obtained by referring to the battery open-circuit voltage characteristic curve table preset in the controller using the battery's SOC value, and the preset battery DC internal resistance is read. In a specific implementation of this invention, the polarization voltage difference is expressed by the formula: In the formula, This indicates the polarization voltage difference of the battery; , These represent the battery terminal voltage and battery current, respectively. Indicates the open-circuit voltage of the battery; This indicates the battery's DC internal resistance. Indicates the absolute value sign. This represents the internal resistance voltage.

[0064] In the above formula, the three features are negatively fused by subtraction.

[0065] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the battery impedance reference term includes: In energy optimization and regulation, the state of the battery electrodes determines the battery's maximum output power, directly determining its maximum usable output power. Electrode aging or increased polarization will limit instantaneous discharge capacity and reduce system scheduling flexibility. Therefore, the battery impedance reference term is obtained by measuring the difference between the polarization voltage difference and the battery current. The larger the value, the more severe the internal polarization of the battery, making it unsuitable for high-power charging and discharging.

[0066] In one specific implementation of this invention, the battery impedance reference term is expressed by the formula: In the formula, Indicates the battery impedance reference term; This represents the normalization scaling factor; This represents a small constant used to prevent the denominator from being zero.

[0067] It should be noted that the normalization scaling factor in this embodiment can be the rated resistance of the battery, and the small constant can be 0.001.

[0068] In the above formula, the two features are positively fused by multiplication, and a small constant is introduced to ensure the meaning of the feature values, which will not be elaborated further.

[0069] Preferably, in some possible implementations of the embodiments of the present invention, the method for correcting the degradation of the remaining battery capacity includes: In energy storage operation scenarios, the charge constraint term is affected by both the cumulative impact of the remaining charge deviating from the 50% range and the aging disturbance amplified nonlinearly by the cycle depth. Therefore, the charge constraint term is obtained by calculating the exponential decay penalty based on the difference between the remaining battery charge and a preset threshold value; the larger the value, the closer the charge is to the upper and lower limits, the higher the assessment of battery life risk, and the stronger the constraint on subsequent output suppression. In a specific implementation of this invention, the charge constraint term is expressed by the formula: In the formula, This represents the charge constraint term for the current control cycle; This represents the penalty coefficient, which can be 10 in this embodiment; This indicates the remaining battery charge in the current control cycle; that is, in the above formula, 0.5 represents the threshold value, used to characterize the battery at 50% capacity, the cubic purpose is to enhance the penalty effect for large deviations, the penalty coefficient is used to adjust the cost weight of constraint violation, and then achieves exponential growth mapping through the exponential function, significantly improving the optimization sensitivity under high-risk conditions.

[0070] Preferably, in some possible implementations of the embodiments of the present invention, the method for obtaining the battery suppression weight includes: During battery power regulation, the battery power is affected by heat, voltage, and remaining charge. Therefore, battery power needs to be suppressed in multiple ways. In a specific implementation of this invention, the battery suppression weight is expressed by the formula: In the formula, This indicates the battery suppression weight for the current control cycle; , , These represent the normalized risk adjustment coefficient, battery impedance benchmark term, and charge constraint term, respectively. , , These represent constraint weights, which can take values ​​of [value] in this implementation. , , .

[0071] It should be noted that the normalization adopts maximum value normalization, and the normalized values ​​of the risk adjustment coefficient, battery impedance reference term, and charge constraint term in the first control cycle are all 1.

[0072] Preferably, in some possible implementations of the embodiments of the present invention, the running loss evaluation function method includes: In wind, solar, thermal, and energy storage systems integrating new energy sources, the operating loss assessment function is affected by both the mechanical lifespan impact caused by the accumulation of thermal stress in the turbine and the polarization and aging disturbances caused by high-rate charging and discharging of the battery. In a specific implementation of this invention, the operating loss assessment function is expressed by the formula: In the formula, This represents the result of the current control cycle operating loss assessment function; This indicates the battery suppression weight for the current control cycle; This represents the battery power calculated in the current control cycle; This represents the preset smoothing coefficient, which is set to 100 in this embodiment; This indicates the rate of change in the power output of the solar thermal power plant during the current control period; Indicates battery power; This represents the power change term of the concentrated solar power (CSP) plant; where the power change rate of the CSP plant in this implementation example is adopted. , This indicates the power output of the solar thermal power plant during the current control cycle. This represents the power output of the solar thermal power plant in the previous control cycle. During initial operation or cold start, the power output of the solar thermal power plant in the previous control cycle is set to 0, and the battery suppression weight can be set to a constant of 1.

[0073] In another implementation of this invention, the operating loss evaluation function is expressed by the formula: In the formula, This represents a very large penalty coefficient, which in this embodiment can take the value of [value missing]. ; This represents the slack variable in the power balance equation; This represents a penalty term that allows for small tracking errors under extreme conditions, ensuring that the solver always has a solution.

[0074] It should be noted that the order of the above embodiments of the present invention is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0075] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

Claims

1. A method for optimal scheduling of a multi-energy system (wind, solar, thermal, and energy storage) based on adaptive dynamic programming, characterized in that, The method includes: In the optimized scheduling of energy, the rotor node temperature, steam temperature, battery current, and remaining battery power are obtained; The initial temperature of a node is obtained by performing a shutdown-startup attenuation analysis on the node temperature; the corrected temperature of the node is obtained by correcting the initial temperature of the node using the ambient temperature; and the virtual temperature of the node is obtained by performing a recursive fusion analysis on the corrected temperature using the steam temperature. By analyzing the differences in the distribution of virtual temperature at nodes, the predicted value of expansion stress is obtained; based on the difference between the predicted value of expansion stress and the preset benchmark value, the risk adjustment coefficient is obtained. By analyzing battery polarization loss, the polarization voltage difference is obtained; by the difference between the polarization voltage difference and the battery current, the battery impedance benchmark term is obtained; and based on the decay penalty of the remaining battery capacity, the charge constraint term is obtained. Calculate the characteristic fusion value of the risk adjustment coefficient, battery impedance benchmark term, and normalized value of charge constraint term to obtain battery suppression weight; construct an operation loss assessment function, which includes battery power term and solar thermal power plant power change term; use the battery suppression weight as the weight of battery power term in the operation loss assessment function; optimize the scheduling of the energy system based on the output of the operation loss assessment function. The method for obtaining the initial temperature includes: The current shutdown and startup status is determined by a preset threshold. When the shutdown duration is greater than or equal to the preset threshold, it is determined to be a cold start, and the initial temperature of the node is the ambient temperature. When the shutdown duration is less than the preset threshold, it is determined to be a hot start. The node's initial temperature is obtained by performing a decay analysis on the node based on the shutdown duration. The method for obtaining the attenuation analysis includes: The difference between the node's temperature at the moment of shutdown and the ambient temperature is calculated, and the natural cooling temperature is obtained by negatively weighting the temperature by the shutdown duration. Based on the feature fusion of the natural cooling temperature and the ambient temperature, the attenuation analysis results of the node temperature cooling by the shutdown duration are obtained. The method for obtaining the node correction temperature includes: Calculate the difference between the current ambient temperature and the predicted rotor surface temperature, and use a preset temperature weight to weight the difference to obtain the current temperature correction difference; correct the initial temperature based on the temperature correction difference to obtain the current corrected temperature. The method for obtaining the virtual temperature of the node includes: The corrected temperature of each node at the previous time step and the steam temperature at the current time step are weighted and fused according to a preset weighting coefficient, and then recursive analysis is performed to obtain the virtual node temperature of each node at different times. The method for obtaining the predicted expansion stress value includes: Based on the difference between the average virtual temperature of all nodes at the current moment and the virtual temperature of the node at the rotor boundary node, and corrected by the material adjustment coefficient, the predicted value of the rotor's expansion stress is obtained. The method for obtaining the risk adjustment coefficient includes: Based on the difference between the predicted value of rotor expansion stress and the preset benchmark value at different times, the coupling matrix is ​​obtained; through dynamic path finding of the coupling matrix, all path trajectories are obtained; based on the maximum value in each path trajectory, the minimum value among all maximum values ​​is obtained and normalized to obtain the trend deviation value of the control cycle; through weighted negative correlation analysis of the trend deviation value, the risk adjustment coefficient at the current time is obtained. The method for obtaining the polarization voltage difference includes: The polarization loss of lithium ions in the electrode material is analyzed based on the differences between the energy storage battery voltage, open circuit voltage, and internal resistance voltage, and the polarization voltage difference of the battery is obtained. The method for correcting the degradation of the remaining battery capacity includes: Based on the difference between the remaining battery capacity and the preset threshold, an exponential decay penalty for overcharging and over-discharging of the battery is applied to obtain the charge constraint term. The operating loss assessment function includes: in, This represents the result of the current control cycle operating loss assessment function; This indicates the battery suppression weight for the current control cycle; This represents the battery power calculated in the current control cycle; This represents the preset smoothing coefficient, with a value of 100. This indicates the rate of change in the power output of the solar thermal power plant during the current control period.

Images