Optimization method for controlling a thermo-chemical energy storage device coupled to a photovoltaic plant and application in a power grid

Abstract

The present application relates to the technical field of thermochemical energy storage device control, and discloses a photovoltaic coupled thermochemical energy storage device control optimization method and application in a power grid, the method comprising: after each complete energy storage cycle, an effective activity coefficient is calculated by a single computing node based on the inlet and outlet temperature and flow of the heat exchanger, a unified operation permission domain is generated, and the domain is written into a three-layer controller constraint register as the only constraint source through a write protection module; during the cycle, a protective correction is triggered based on the Mahalanobis distance of the bed temperature field space distribution characteristic vector relative to the baseline, a reduction factor is generated, and the upper limit of the charging power is corrected to tighten the permission domain. The present application eliminates the inconsistency of interlayer constraints and the internal safety blind area of the cycle, and improves the operation safety margin and cycle life in the scenarios of light abandonment and power grid peak regulation.

Description

Technical Field

[0001] This invention relates to the field of thermochemical energy storage equipment control technology, and in particular to a three-layer collaborative control method and system for thermochemical energy storage systems based on a unified operating permit domain of the activity coefficient of energy storage materials. It is particularly applicable to reversible thermochemical energy storage systems containing calcium oxide / calcium hydroxide fixed-bed reactors, specifically a photovoltaic-coupled thermochemical energy storage equipment control optimization method and its application in the power grid. Background Technology

[0002] Thermochemical energy storage utilizes reversible chemical reactions to absorb heat during the charging phase and release heat during the releasing phase, featuring high energy density, long storage period, and suitability for medium- and high-temperature thermal energy management. The calcium oxide / calcium hydroxide system achieves heat storage and release through the dehydration and decomposition of calcium hydroxide and the hydration of calcium oxide. It has been widely studied due to its wide availability of raw materials, high reaction enthalpy, and relatively large theoretical energy density.

[0003] For this type of system, existing technologies typically have the following foundations: First, relatively mature models have been established around the dynamics, heat and mass transfer, and bed pressure drop of fixed-bed reactors; second, model predictive control methods have been developed to adjust variables such as reactor temperature and output power; and third, there are also technical solutions that combine state estimation with controllers for estimating conversion rate within a single cycle.

[0004] On the other hand, calcium oxide / calcium hydroxide energy storage materials are affected by factors such as sintering, pore structure evolution, and reduction of active sites during long-term cycling, leading to an overall decrease in material activity. This decline does not necessarily exhibit a strictly monotonic sequence; the identification value for a single cycle may show non-monotonic fluctuations due to measurement errors, operating condition variations, and local reactivation. However, overall, the effective activity of the material shows a downward trend.

[0005] A common, yet unresolved, problem in existing control architectures is the inconsistency of constraint boundaries between control layers at different timescales. The scheduling and planning layer often arranges energy storage plans based on longer timescales, the power scheduling layer performs power optimization based on rolling forecasts, and the execution layer completes lower-level control based on the real-time status of the reactor. If each layer maintains its capacity and power boundaries based on its own model, empirical parameters, or independent estimates, after material activity decays, situations can easily arise where the upper-level plan is feasible but the lower-level plan is physically inaccessible, leading to plan execution deviations, system performance degradation, and even increased safety risks.

[0006] Especially for calcium oxide / calcium hydroxide fixed-bed reactors, when the material activity decreases, the reaction rate decreases under the same temperature and steam conditions, the reaction endothermic capacity weakens, and the bed temperature rise rate increases under the same heating power. If the scheduling and optimization layers still determine the operating strategy based on the initial or independently estimated capacity boundaries, it is impossible to guarantee consistency with the physical constraints of the execution layer.

[0007] The fundamental flaw in existing technologies, where each control layer independently maintains its liveness estimate or constraint boundaries, lies in the absence of a unified source of constraints with write protection mechanisms. Each layer can modify its own constraint boundaries, resulting in inconsistencies between layers that cannot be eliminated in the architecture.

[0008] In addition to the aforementioned inconsistencies in interlayer constraints, existing control schemes generally suffer from insufficient ability to detect local anomalies within the cycle. During a single charge-discharge cycle in a calcium oxide / calcium hydroxide fixed-bed reactor, rapid deterioration of localized material activity may occur due to factors such as accelerated sintering in localized areas of the bed, shifts in steam distribution, or localized blockages in heat exchange pipes. This deterioration is reflected in abnormal changes in the spatial distribution characteristics of the bed temperature field during the cycle. However, in existing technologies, the identification or estimation of the activity coefficient is typically performed on a complete cycle basis, failing to capture such localized abrupt events midway through the cycle. Within the time window from the occurrence of a localized anomaly within the cycle to the end of the current cycle, the control system still relies on the constraint parameters calibrated at the end of the previous cycle. The actual bed conditions have deviated from the safety envelope corresponding to these parameters, constituting a safety blind zone. The duration of this blind zone is equal to the remaining time of the current cycle, and can reach several hours under typical operating conditions, sufficient to cause localized overheating or irreversible material damage.

[0009] Therefore, it is necessary to propose a new collaborative control method. On the one hand, the material activity decay information is not interpreted separately by each layer, but is uniformly mapped to an operating permission domain shared by the three layers and cannot be modified independently by each layer, so as to eliminate the problem of inconsistent constraints from the system architecture. On the other hand, the spatial distribution characteristics of the bed temperature field are continuously monitored during the cycle. When a local anomaly is detected, the operating permission domain is corrected in real time without waiting for the cycle to end, thus eliminating the safety blind spot within the cycle. Summary of the Invention

[0010] The technical problem to be solved by this invention is that in the existing multi-layer control architecture of thermochemical energy storage systems, each control layer independently maintains the estimated value of material activity or constraint boundary. After the material activity decays, the constraints between layers become inconsistent, and the upper layer is feasible in terms of planning but the lower layer is physically inaccessible. At the same time, the activity identification is updated with a complete cycle as the update cycle, which cannot capture the rapid deterioration of local material activity caused by local sintering, vapor shift, etc. during the cycle, forming a safety blind zone within the cycle.

[0011] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0012] A control optimization method for photovoltaic-coupled thermochemical energy storage equipment, the method being executed based on a system comprising a calcium oxide / calcium hydroxide fixed bed reactor, an electric heater, a heat exchanger, a scheduling plan controller, a power scheduling controller, a reactor bottom layer controller, a single computing node, an industrial real-time communication interface, write-protected constraint registers respectively set in the above three controller layers, a local default memory, and a bed temperature field acquisition device, including the following steps:

[0013] Step S1: Perform laboratory cycle characteristics calibration on the energy storage material and calibrate the system thermal efficiency and safety parameters to obtain benchmark parameters including the lower limit of the asymptotic activity coefficient α_inf, the apparent activation energy E_a, and the decay model parameters, and write them into the read-only configuration area.

[0014] Step S2: After each complete energy storage and release cycle, the single computing node calculates the heat released by the heat exchanger side based on the readings of the inlet and outlet temperature sensors on the heat exchanger side and the readings of the heat exchange medium mass flow sensor, and takes into account the heat loss from the reactor surface to obtain the original value of the effective activity coefficient α_raw(n) for this cycle. After exponential weighted average filtering, α_filt(n) is obtained.

[0015] Steps S3-S5: Steps S3, S4, and S5 are executed by the scheduling plan controller, the power scheduling controller, and the reactor bottom layer controller, respectively; the single computing node centrally generates a unified operating permission domain Ω(n) based on α_filt(n), and synchronously writes it into the constraint registers of the scheduling plan controller, the power scheduling controller, and the reactor bottom layer controller through the industrial real-time communication interface and the write protection module; the three-layer controllers only have read permissions for the constraint registers, and use Ω(n) in the constraint registers as their sole source of constraints;

[0016] Step S6: The scheduling plan controller and the power scheduling controller execute the day-ahead planning and rolling power optimization with Ω(n) as the boundary, respectively. The reactor bottom controller adaptively adjusts the PID gain and feedforward term according to α_filt(n) and implements three-level operating domain management.

[0017] Step S7: Perform three-level hierarchical management of the material operating domain (A / B / C) according to the distance of α_filt(n) relative to α_inf, and generate maintenance decision support output;

[0018] Step S8: During each energy storage cycle, the single computing node constructs a spatial distribution feature vector of the bed temperature field based on the real-time readings of no less than 5 temperature measurement points along the axial direction of the bed temperature field acquisition device, calculates the Mahalanobis distance D_M_bar(t) of the feature vector relative to the sliding baseline covariance, and when D_M_bar(t) exceeds the threshold D_th2 for N_conf consecutive sampling periods, the protective correction of the intra-cycle allowable domain is triggered immediately, generating the intra-cycle reduction coefficient κ(t) and the intra-cycle allowable domain Ω_intra(t) and writing them into the constraint register, which together with the inter-cycle identification in step S2 constitute a dual time-scale activity monitoring architecture;

[0019] Step S9: Implement a three-level response process for communication link failures of the industrial real-time communication interface. When the failure continues to exceed the threshold, the three-layer controller switches to the conservative and safe default value in the local default memory.

[0020] The application of the photovoltaic coupled thermochemical energy storage device control optimization method in the power grid is described. The method is applied to the thermochemical energy storage peak shaving and frequency regulation scenario under photovoltaic output fluctuation conditions. The method takes the photovoltaic available power prediction value and the heat network demand prediction value as input, and the unified operation permission domain as the inviolable boundary. The corresponding optimization and control are executed at the scheduling planning layer, power scheduling layer and bottom control layer respectively.

[0021] Beneficial effects

[0022] By having a single computing node centrally identify the effective activity coefficient and generate a unified operating permission domain after each complete energy storage and release cycle, and then synchronously distribute it to the three-layer controller via the write protection constraint register, the problem of inconsistent constraints caused by the three-layer controllers independently estimating the material state is eliminated, fundamentally ensuring the consistency of inter-layer constraints.

[0023] The real-time correction mechanism for the dynamic permissive domain within the cycle, based on the spatial distribution characteristics of the bed temperature field, introduced in step S8, together with the inter-cycle identification in step S2, forms a dual-timescale activity monitoring architecture. This improves the temporal resolution of activity monitoring from the cycle level to the control sampling level, eliminating the safety blind spot caused by local anomalies within the cycle. Simulation results show that, in scenarios with local anomalies within the cycle, the complete scheme including step S8 reduces the maximum bed temperature shift by 74.3% during anomalies by 74.3% compared to the scheme without step S8, increases the detection rate of local over-temperature events within the cycle from 0% to 100%, and the average delay from anomaly to permissive domain correction taking effect is only 7.3 seconds.

[0024] The three-layer controller passively receives constraint boundaries through constraint registers and does not participate in liveness estimation and anomaly determination calculations. The constraint registers are restricted by write protection mechanisms to allow only a single computing node to write to them, thus eliminating the possibility of multi-source constraint conflicts at the architectural level.

[0025] The three-level response handling for communication link failures and the local conservative safety default value mechanism ensure the system's ability to safely degrade operation under abnormal communication conditions.

[0026] The material operation domain hierarchical management and maintenance decision support function enables operation and maintenance personnel to obtain quantitative replacement suggestions when the material activity is below the threshold, reducing the subjectivity and lag of operation and maintenance decisions. Attached Figure Description

[0027] Figure 1 This is a schematic diagram of the architecture of a three-layer collaborative control system for a thermochemical energy storage system. The labels in the diagram are explained as follows: 1-Calcium oxide / calcium hydroxide fixed-bed reactor; 2-Electric heater; 3-Heat exchanger; 4-Schedule planning controller; 5-Power scheduling controller; 6-Reactor bottom-level controller; 7-Single computing node; 8-Industrial real-time communication interface; 9-Constraint registers (located in 4, 5, and 6 respectively); 10-Local default memory; 11-Write protection module; 12-Bed temperature field acquisition device. Solid arrows in the diagram indicate the direction of heat / material flow, dashed arrows indicate the direction of signal flow, and double arrows indicate the write channel to the constraint register via the write protection module. Detailed Implementation

[0028] I. System Overall Architecture

[0029] like Figure 1 As shown, the three-layer collaborative control system of the thermochemical energy storage system of the present invention includes the following components: a calcium oxide / calcium hydroxide fixed bed reactor 1, an electric heater 2, a heat exchanger 3, a scheduling plan controller 4, a power scheduling controller 5, a reactor bottom layer controller 6, a single computing node 7, an industrial real-time communication interface 8, constraint registers 9 respectively set in the scheduling plan controller 4, the power scheduling controller 5 and the reactor bottom layer controller 6, a local default memory 10, a write protection module 11 and a bed temperature field acquisition device 12.

[0030] Fixed-bed reactor 1 is used for heat storage and release reactions. Electric heater 2 is used to input heat to fixed-bed reactor 1 during the energy charging stage. Heat exchanger 3 is used to output heat to the heating network during the energy release stage. Single computing node 7 is used to identify the effective activity coefficient and centrally generate a unified operating permission domain after each complete energy storage and release cycle. It is also used to perform anomaly detection and protective correction within the unified operating permission domain during the cycle based on the spatial distribution characteristics of the bed temperature field. Industrial real-time communication interface 8 is used to synchronously write the unified operating permission domain into the constraint registers 9 of each controller. The constraint registers 9 are restricted by write protection module 11, allowing only industrial real-time communication interface 8 to write under authorized conditions. Scheduling plan controller 4, power scheduling controller 5, and reactor bottom-level controller 6 only have read permissions.

[0031] The bed temperature field acquisition device 12 has no less than 5 temperature measurement points evenly arranged along the axis of the fixed bed reactor 1, and the acquired data is transmitted to a single computing node 7 in real time.

[0032] Electric heater 2 is connected to fixed-bed reactor 1 to provide charging heat to fixed-bed reactor 1. Fixed-bed reactor 1 is connected to heat exchanger 3 to output released heat. Single computing node 7 is connected to constraint register 9 inside scheduling controller 4, power scheduling controller 5, and reactor bottom controller 6 via industrial real-time communication interface 8. Write protection module 11 is located between industrial real-time communication interface 8 and constraint register 9, or on the write channel of constraint register 9. Local default memory 10 is connected to scheduling controller 4, power scheduling controller 5, and reactor bottom controller 6 for fault degradation.

[0033] II. Complete Steps of Control Method

[0034] A three-layer collaborative control method for thermochemical energy storage systems based on a unified operating permit domain of energy storage material activity coefficients includes the following steps:

[0035] Step S1: System initialization and baseline parameter calibration

[0036] S1.1 Energy Storage Material Cycle Characteristic Calibration

[0037] Unified Notation: α_hat(n) is the raw calculated value of the effective activity coefficient obtained from a single cycle of identification during the laboratory calibration stage (step S1), without filtering; α_raw(n) is the raw calculated value of the effective activity coefficient obtained from a single cycle of identification during the engineering online identification stage (step S2 and thereafter), without filtering. It is the same physical quantity as α_hat(n) but named differently at different stages, and their definitions are consistent (see Equations 104 and 204); α_filt(n) is the effective activity coefficient after exponential weighted average filtering (the value written to the constraint register); α_eq_intra(t) is the equivalent activity coefficient within the cycle, equal to κ(t) × α_filt(n); α_inf is the lower limit of the asymptotic activity coefficient for long-term cycling of the material, determined by laboratory calibration, with a typical value of 0.60~0.70. α_lab(n) is obsolete throughout the text; α_hat(n) is used uniformly in the laboratory stage, and α_raw(n) is used uniformly in the engineering stage.

[0038] Before the system is put into operation, the laboratory cycle characteristics of the calcium oxide material to be filled are calibrated.

[0039] Sampling procedure: Randomly sample at no less than 3 sampling points in the batch to be filled, with no less than 200g at each point. After mixing evenly, take 500g of the mixture for testing.

[0040] Test equipment: a programmed temperature fixed bed test bench, equipped with a thermogravimetric analyzer with an accuracy better than ±0.1 micrograms and a gas mass flow controller with an accuracy of ±0.5% of full scale.

[0041] Standard test protocol:

[0042] Energy charging stage operation: Increase to the target dehydration temperature at a rate not exceeding 10 degrees Celsius per minute, introduce dry nitrogen gas with a purity of not less than 99.99% at this temperature, maintain until the outlet water vapor mole fraction is less than 0.0005, and record the change in dehydration mass.

[0043] Energy release stage operation: Cool to the target hydration temperature at a rate not exceeding 5 degrees Celsius per minute, introduce water vapor, maintain until the bed mass change rate is less than 0.01 mg per minute, and record the hydration mass change.

[0044] The above operations constitute a complete loop, which is executed at least 20 times, preferably 50 times.

[0045] Definition of laboratory calibration of effective activity coefficient:

[0046] α(n)=Δm_h(n) / Δm_h(1) (Equation 101)

[0047] In the formula: α(n) is the effective activity coefficient of the nth cycle, dimensionless; Δm_h(n) is the change in hydration mass of the nth cycle, in grams; Δm_h(1) is the change in hydration mass of the 1st cycle, in grams; n is the cycle number, n=1,2,3,... When n=1, α(1) is exactly equal to 1 by definition.

[0048] Online identification method for engineering: In engineering applications, the change in hydration quality cannot be directly measured online; therefore, measurable quantities on the heat exchanger side are used as a substitute. After the energy release phase of each complete energy storage and release cycle, the heat released on the heat exchanger side for this cycle is calculated by integrating the readings of the heat exchanger inlet and outlet temperature sensors and the mass flow sensor readings of the heat exchange medium throughout the entire energy release phase.

[0049] Q_ex(n)=∫[ts(n),te(n)]ṁ(t)×cp_avg×[T_out(t)-T_in(t)]dt (Equation 102)

[0050] In the formula: Q_ex(n) is the heat released by the heat exchanger side in the nth cycle, in joules; ts(n) is the start time of the energy release stage in the nth cycle, in seconds; te(n) is the end time of the energy release stage in the nth cycle, in seconds; ṁ(t) is the mass flow rate of the heat exchange medium at time t, in kilograms per second; cp_avg is the isobaric specific heat capacity of the heat exchange medium, in joules per kilogram per Kelvin; T_out(t) is the sensor reading of the heat exchanger outlet temperature at time t, in Kelvin; T_in(t) is the sensor reading of the heat exchanger inlet temperature at time t, in Kelvin.

[0051] In the above formula, cp is taken as the average isobaric specific heat capacity of the heat exchange medium within the operating temperature range, used for simplified calibration calculations in the laboratory. In online engineering applications (step S2.1), the temperature-dependent cp value is used for the integral of released heat to improve accuracy, as detailed in Equation 201. The difference between the two values ​​has an impact of no more than 0.5% on the effective activity coefficient identification value and does not affect the determination of the activity decay trend.

[0052] Considering heat loss from the reactor surface, a heat loss correction is introduced to obtain the corrected heat of energy release:

[0053] Q_c(n)=Q_ex(n)+Q_loss(n) (Equation 103)

[0054] In the formula: Q_c(n) is the corrected heat release in the nth cycle, in joules; Q_loss(n) is the heat loss from the reactor surface in the nth cycle, in joules.

[0055] The effective activity coefficient identification value in the laboratory calibration stage is defined by the ratio of corrected heat of release:

[0056] α_hat(n)=Q_c(n) / Q_c(1) (Equation 104)

[0057] In the formula: α_hat(n) is the original calculated value of the effective activity coefficient obtained from the nth laboratory calibration cycle, which is dimensionless; Q_c(n) is the corrected integral value of the exothermic heat from the nth calibration cycle, in joules; Q_c(1) is the corrected integral value of the exothermic heat from the first calibration cycle, which serves as the activity normalization benchmark, in joules. The engineering online identification stage (step S2) adopts the same definition (see Equation 204), and the original calculated value is denoted as α_raw(n).

[0058] Activity decay model fitting:

[0059] α_mod(n)=α_inf+(1-α_inf)×exp(-k_d×n^β) (Equation 105)

[0060] In the formula: α_mod(n) is the effective activity coefficient predicted by the decay model in the nth cycle, dimensionless; α_inf is the residual activity limit value, dimensionless, representing the lower limit of the effective activity coefficient when the number of cycles approaches infinity; k_d is the decay rate constant, dimensionless; β is the decay shape parameter, dimensionless; exp(·) represents an exponential function with the natural constant e as the base.

[0061] The apparent activation energy was determined using the Kissinger method.

[0062] ln(r_h,i / T_p,i²)=−(E_a / R)×(1 / T_p,i)+C_0 (Equation 106)

[0063] In the formula: r_h,i is the heating rate of the i-th calibration experiment, in Kelvin per second; T_p,i is the peak temperature of the differential thermal analysis of the i-th calibration experiment, in Kelvin; E_a is the apparent activation energy, in joules per mole; R is the universal gas constant, with a value of 8.314 joules per mole per Kelvin; C_0 is the fitting intercept constant; ln(·) represents the natural logarithm function.

[0064] Calibration method for the asymptotic activity coefficient lower limit α_inf: After no less than 20 calibration cycles in step S1, the arithmetic mean of the α_hat values ​​from the last 5 cycles is taken as the initial estimate of α_inf; if the material has not yet entered the stable decay period, α_inf is taken as a conservative value based on long-term cycle data provided by the manufacturer or literature data of similar materials, typically ranging from 0.60 to 0.70. Once calibrated, α_inf is written to the read-only configuration area and used as a fixed parameter in Equations 216, 813, etc.

[0065] S1.2 System thermal efficiency and safety parameter calibration

[0066] During the system commissioning period, perform no less than three standard charge-discharge cycles and calibrate the following parameters.

[0067] System thermal efficiency:

[0068] η_sys=∫ṁ(t)×cp((T_out(t)+T_in(t)) / 2)×(T_out(t)-T_in(t))dt / [∫P_heat(t)dt-ΔE_sen] (Equation 107)

[0069] In the formula: η_sys is the system thermal efficiency, dimensionless; P_heat(t) is the heating power of the electric heater at time t, in watts; ΔE_sen is the change in sensible heat storage in the bed, in joules; the meanings of the other symbols are the same as before.

[0070] Changes in bed sensible heat storage:

[0071] ΔE_sen=Σ(j=1 to N_s)[(m_CaO / N_s)×cp_CaO×(T_r,j(t_e)-T_r,j(t_s))](Equation 108)

[0072] In the formula: ΔE_sen is the change in sensible heat storage in the bed, in joules; N_s is the number of axial temperature measurement points, dimensionless; m_CaO is the total mass of calcium oxide, in kilograms; cp_CaO is the specific heat capacity of calcium oxide at constant pressure, in joules per kilogram per Kelvin; T_r,j(t_e) is the temperature of the j-th temperature measurement point at the end of charging at time t_e, in Kelvin; T_r,j(t_s) is the temperature of the j-th temperature measurement point at the beginning of charging at time t_s, in Kelvin.

[0073] In this step, the following parameters are simultaneously calibrated: zone-wide integrated heat transfer coefficient U_j (in watts per square meter per Kelvin), maximum conversion rate η_x (dimensionless), upper limit of safe temperature rise rate dTr_safe (in Kelvin per second), reactor thermal response time constant τ_min (in seconds), non-reactive power ratio γ (dimensionless), and bottom-level controller reference proportional gain K_p0, reference integral gain K_i0, and reference differential gain K_d0.

[0074] S1.3 Calculation of Engineering Capacity Parameters

[0075] Theoretical maximum energy storage capacity:

[0076] E_th=(m_CaO / M_CaO)×ΔH_r (Equation 109)

[0077] In the formula: E_th is the theoretical maximum energy storage capacity, in joules; m_CaO is the total mass of calcium oxide, in kilograms; M_CaO is the molar mass of calcium oxide, with a value of 56.08 grams per mole; ΔH_r is the enthalpy of reaction, with a value of 104 kilojoules per mole.

[0078] Baseline available energy storage capacity:

[0079] E_avail=E_th×η_x×η_sys (Equation 110)

[0080] In the formula: E_avail is the baseline available energy storage capacity, in joules; η_x is the maximum conversion rate, dimensionless; η_sys is the system thermal efficiency, dimensionless.

[0081] Taking a loading capacity of 67,000 kg as an example: E_th is approximately 34,522 kWh, E_avail is approximately 25,481 kWh, and the system configuration capacity E_conf is 10.8 MWh in this embodiment.

[0082] S1.4 Hard Constraint Parameter Acquisition for Heating Network Operation

[0083] Confirm with the heating network operator the hard constraint parameters such as the maximum rate of change of heating power, the maximum change of water supply temperature over 15 minutes, and the minimum continuous heating power stipulated in the contract, and write them into the control system configuration file in a way that cannot be modified online.

[0084] S1.5 Constraint Register Initialization and Conservative Default Value Pre-writing

[0085] Before system commissioning, conservative safety default values ​​are pre-written into the constraint registers of each of the three controllers, including conservative upper limits for charging power, conservative upper limits for state of charge, and conservative temperature reference values. These conservative safety default values ​​are stored independently of the main communication link in the local non-volatile memory of each controller. They are written by the system integration engineer using a dedicated physical key through an offline interface and cannot be modified through the online software interface during operation.

[0086] Step S2: Identification of effective activity coefficients and centralized generation and synchronous distribution of unified operating license domains

[0087] Execution is triggered at the end of each complete energy storage / dissipation cycle.

[0088] S2.1 Exact Integral of Released Heat

[0089] Q_meas(n)=∫[ts(n),te(n)]ṁ(t)×cp((T_out(t)+T_in(t)) / 2)×[T_out(t)-T_in(t)]dt (Equation 201)

[0090] In the formula: Q_meas(n) is the measured heat released in the nth cycle, in joules; the meanings of the other symbols are the same as in formula 102.

[0091] Energy release termination criterion: The temperature difference between the heat exchanger outlet and inlet is below the calibrated threshold ΔT_end for 300 consecutive seconds.

[0092] S2.2 Reactor Surface Heat Loss Correction

[0093] Q_loss(n)=Σ(j=1 to N_z)[U_j×A_j×(T_r_avg,j(n)-T_amb(n))×t_cyc(n)] (Equation 202)

[0094] In the formula: Q_loss(n) is the surface heat loss of the reactor in the nth cycle, in joules; N_z is the number of axial zones of the reactor, dimensionless; U_j is the overall heat transfer coefficient of the jth zone, in watts per square meter per Kelvin; A_j is the external surface area of ​​the jth zone, in square meters; T_r_avg,j(n) is the time-averaged bed temperature of the jth zone in the nth cycle, in Kelvin; T_amb(n) is the average ambient temperature in the nth cycle, in Kelvin; t_cyc(n) is the duration of the nth cycle, in seconds.

[0095] Q_c(n)=Q_meas(n)+Q_loss(n) (Equation 203)

[0096] In the formula: Q_c(n) is the corrected heat release in the nth cycle, in joules.

[0097] S2.3 Original calculated value of effective activity coefficient

[0098] To distinguish the engineering online identification phase from the laboratory calibration phase (α_hat), the raw value of a single cycle identification will be denoted as α_raw(n) in this section and subsequent sections. Its definition is completely consistent with Equation 104, except that the data source is changed from the laboratory test bench to the heat exchanger side sensor at the engineering site.

[0099] α_raw(n)=Q_c(n) / Q_c(1) (Equation 204)

[0100] In the formula: α_raw(n) is the original calculated value of the effective activity coefficient in the nth cycle, which is dimensionless.

[0101] S2.4 Triple Validity Check

[0102] Verification 1 (Physical upper limit verification): If α_raw(n) is greater than 1 plus 3 multiplied by σ_m, then it is truncated to 1.0 and an alarm is triggered.

[0103] α_raw(n)≤1+3×σ_m (Formula 205)

[0104] In the formula: σ_m is the measurement standard deviation of the integral of the released heat, which is dimensionless and determined by calibration experiments.

[0105] Verification 2 (Single Change Amplitude Verification):

[0106] δ_max(n)=max(δ_cal,3×σ_hist(n)) (Equation 206)

[0107] In the formula: δ_max(n) is the upper limit of the single change amplitude in the nth cycle, dimensionless; δ_cal is the maximum physical decay amount determined by the calibration experiment, dimensionless; σ_hist(n) is the standard deviation of the recent historical single change amount, dimensionless.

[0108] σ_hist(n)=sqrt{[1 / (N_h-1)]×Σ(i=n-N_h to n-1)[α_filt(i)-α_filt(i-1)] 2 (Equation 207)

[0109] In the formula: σ_hist(n) is the standard deviation of the historical single change, dimensionless; N_h is the length of the historical window, taken as min(n-1,20); α_filt(i) is the effective activity coefficient after filtering in the i-th cycle, dimensionless; sqrt{·} represents taking the arithmetic square root. When N_h<3, check 2 is skipped, and only check 1 and check 3 are executed.

[0110] Verification 3 (Model Consistency Verification):

[0111] |α_raw(n)-α_mod(n)|≤3×σ_fit (Equation 208)

[0112] In the formula: α_mod(n) is the predicted value of the decay model in the nth cycle, which is dimensionless; σ_fit is the standard deviation of the model fitting residuals, which is dimensionless; |·| represents taking the absolute value. If the condition is not met, the model parameters are updated.

[0113] S2.5 Exponential Weighted Average Filter

[0114] α_filt(n)=λ×α_raw(n)+(1-λ)×α_filt(n-1) (Equation 209)

[0115] In the formula: α_filt(n) is the effective activity coefficient after filtering in the nth cycle, which is dimensionless; α_filt(n-1) is the effective activity coefficient after filtering in the (n-1)th cycle, which is dimensionless; λ is the exponential weighted average smoothing coefficient, which is dimensionless.

[0116] λ=2 / (N_sm+1) (Equation 210)

[0117] In the formula: N_sm is the equivalent smoothing window length, which is dimensionless. Preferably, N_sm = 5, corresponding to λ = 1 / 3.

[0118] S2.6 Unified Operational License Domain Ω(n) Centralized Calculation

[0119] (a) Upper limit of state of charge:

[0120] SOC_max(n)=f_SOC(α_filt(n)) (Equation 211)

[0121] Preferred implementation method:

[0122] SOC_max(n)=α_filt(n) (Equation 212)

[0123] In the formula: SOC_max(n) is the upper limit of the state of charge after the nth cycle, which is dimensionless; f_SOC(·) is the upper limit mapping function of the state of charge.

[0124] In Equation 212, SOC_max(n) takes α_filt(n) as its value, which is based on the following normalization convention: α_filt(n) is defined as the ratio of the effective heat that can be released in the current cycle to the heat that can be released in the rated cycle, and SOC_max(n) is defined as the ratio of the current available energy storage capacity to the nominal rated capacity. The two have the same normalization benchmark in the dimension of heat that can be released, so they can be numerically equivalent, and the consistency of the energy dimension constraint is guaranteed by E_perm(n) = E_conf × SOC_max(n).

[0125] (ii) Maximum energy limit of the permitted domain:

[0126] E_perm(n)=E_conf×SOC_max(n) (Equation 213)

[0127] In the formula: E_perm(n) is the upper limit of the permitted domain energy after the nth cycle, in joules; E_conf is the system configuration capacity, in joules.

[0128] (III) Maximum safe charging power:

[0129] P_ch_max(n)=C_r×dTr_safe+Q_dot_rxn_nom×α_filt(n)+Q_dot_loss_nom (Equation 214)

[0130] In the formula: P_ch_max(n) is the upper limit of safe charging power after the nth cycle, in watts; C_r is the equivalent heat capacity of the reactor bed, in joules per Kelvin; dTr_safe is the upper limit of safe temperature rise rate, in Kelvin per second; Q_dot_rxn_nom is the nominal reaction endothermic rate under the condition that the activity coefficient is equal to 1, in watts; Q_dot_loss_nom is the nominal heat loss rate, in watts.

[0131] Engineering approximation:

[0132] P_ch_max(n)=P_ch_rated×[γ+(1-γ)×α_filt(n)] (Equation 215)

[0133] In the formula: P_ch_rated is the rated charging power, in watts; γ is the non-reactive power ratio, dimensionless.

[0134] In Equation 215, γ is defined as the ratio of non-reactive power (including the sum of bed sensible heat rise power P_sens_nom and reactor surface heat loss power P_loss_nom) to rated charging power under rated operating conditions:

[0135] γ=(P_sens_nom+P_loss_nom) / P_ch_rated (Equation 215a)

[0136] In the formula: P_sens_nom is the sensible heat rise power of the bed under rated operating conditions, which is equal to C_r×dTr_safe and is in watts; P_loss_nom is the heat loss power of the reactor surface under rated operating conditions and is in watts; P_ch_rated is the rated charging power and is in watts. The typical value range of γ is 0.15~0.30. This definition ensures that Equation 215 and Equation 214 are strictly equivalent under rated operating conditions: when α_filt(n)=1, Equation 215 gives P_ch_rated; when α_filt(n)→0, Equation 215 gives γ×P_ch_rated, which corresponds to the sum of pure sensible heat rise and heat loss, and is consistent with the physical meaning of (Q_dot_rxn_nom×0+Q_dot_loss_nom+C_r×dTr_safe) in Equation 214.

[0137] (iv) Reactor temperature reference values ​​and feasible temperature range:

[0138] ΔT_comp(n)=(R×T_nom² / E_a)×ln(1 / α_filt(n)) (Equation 216)

[0139] In the formula: ΔT_comp(n) is the temperature compensation after the nth cycle, in Kelvin; R is the universal gas constant, with a value of 8.314 joules per mole per Kelvin; T_nom is the nominal reaction temperature, in Kelvin; E_a is the apparent activation energy, in joules per mole; ln(·) represents the natural logarithm function.

[0140] Additional explanation: When α_filt(n) approaches α_inf, ΔT_comp(n) becomes large, potentially causing T_ref(n) = T_ref_0 + ΔT_comp(n) to exceed the intrinsically safe temperature limit T_safe_max. In this case, according to Equation 218: T_ref(n) = min(T_ref_0 + ΔT_comp(n), T_safe_max), the temperature reference value is rigidly limited by T_safe_max. Under these circumstances, the system no longer compensates for activity decay by increasing the temperature. Instead, the operating domain degradation mechanism of Equations 307 and 701 takes over, maintaining interlayer constraint consistency by reducing the charging power and the upper limit of the state of charge, ensuring that the reactor operation always remains within the intrinsically safe boundary.

[0141] T_ref(n)=T_nom+ΔT_comp(n) (Equation 217)

[0142] In the formula: T_ref(n) is the temperature reference value after the nth cycle, in Kelvin.

[0143] Ω_T(n)=[T_ref(n),T_ref(n)+ΔT_tol]∩(-∞,T_safe_max] (Equation 218)

[0144] In the formula: Ω_T(n) is the feasible temperature domain after the nth cycle, in Kelvin; ΔT_tol is the temperature tolerance, in Kelvin; T_safe_max is the upper limit of the safe temperature, in Kelvin; ∩ represents the set intersection operation.

[0145] (v) The discrete state D(n) of the operating domain is determined based on α_filt(n) and T_ref(n) according to the following rules:

[0146] When α_filt(n)≥α_A and T_ref(n)≤T_safe_max−T_margin: D(n)=domain A (normal operating domain);

[0147] When α_B≤α_filt(n)<α_A, or T_ref(n) is in the interval (T_safe_max−T_margin,T_safe_max): D(n)=domain B (reduced operation domain), and the domain B judgment condition must be met for 3 consecutive loops to take effect;

[0148] When the unlimited reference temperature T_ref_0(n) + ΔT_comp(n) ≥ T_safe_max before hard limiting according to Equation 218, it is determined to switch to the restricted operating domain C.

[0149] In the formula: α_A is the lower limit threshold of the normal domain, dimensionless, and is taken as 0.8 in this embodiment; α_B is the lower limit threshold of the derating domain, dimensionless, and is taken as 0.5 in this embodiment; T_margin is the safety temperature margin, in Kelvin. The Ω(n) correction rules corresponding to each operating domain are detailed in step S7.2.

[0150] (vi) Unified runtime license domain structure:

[0151] Ω(n)={SOC_max(n),E_perm(n),P_ch_max(n),Ω_T(n),D(n)} (Equation 219)

[0152] S2.7 Rigorous Mathematical Proof of Inter-Layer Constraint Consistency

[0153] Theorem 1 (Scheduling-Execution Feasibility Domain Inclusion Relationship): If all three layers of controllers read constraint parameters from the same Ω(n), then the power plan feasible set F_plan(n) generated by the scheduling planning layer is a subset of the power execution layer physically reachable set R_exec(n), i.e., F_plan(n)⊆R_exec(n).

[0154] prove:

[0155] The arbitrary power trajectory P_c(t) generated by the scheduling planning layer satisfies 0≤P_c(t)≤P_ch_max(n) and 0≤SOC(t)≤SOC_max(n), where both P_ch_max(n) and SOC_max(n) are read from the constraint register.

[0156] The physical reachable set R_exec(n) of the power execution layer is defined as the set of all power trajectories that the reactor can safely execute under the current material physical state. From Equation 214, P_ch_max(n) is precisely the maximum heating power that ensures the temperature rise rate does not exceed the safety limit dTr_safe; from Equation 211, SOC_max(n) = α_filt(n) is precisely the maximum effective state of charge that the current material can achieve. Therefore, any power trajectory satisfying P_ch_max(n) and SOC_max(n) in the constraint register can be safely executed by the reactor, and F_plan(n) ⊆ R_exec(n) holds. Proof complete.

[0157] Complete timing, write protection mechanism and degradation handling of S2.8Ω(n) synchronous distribution

[0158] The implementation methods of constraint register write protection mechanism include hardware write protection circuit (preferred), software access control or controller firmware lock.

[0159] Complete write timing sequence:

[0160] Time T0: The criterion for the end of the nth cycle of energy release is satisfied.

[0161] Time T1: A single computing node completes all calculations from Equation 201 to Equation 218, generating Ω(n). The time from T1 to T0 does not exceed 5 seconds.

[0162] Time T2: A single computing node encodes Ω(n) into a communication frame, appending a version number n and a message authentication code.

[0163] At time T3: The communication frame is sent to the Layer 3 controller via the industrial real-time communication interface 8. The communication interface type is industrial Ethernet or fieldbus, and the communication cycle is no more than 100 milliseconds.

[0164] At time T4: Each controller receives a communication frame, verifies the message authentication code and version number, and after successful verification, the write protection module 11 allows writing and writes each component of Ω(n) into the constraint register 9.

[0165] At time T5: A write completion confirmation signal is returned to the single compute node 7. The time from T5 to T0 shall not exceed 7.1 seconds.

[0166] Three-level response handling for communication link failures:

[0167] Level 1 Fault (Single Frame Loss): A single computing node will automatically retransmit in the next communication cycle, and the constraint register will maintain the previous valid value Ω(n-1).

[0168] Level 2 fault (more than 10 consecutive frames lost): A single compute node triggers an alarm, and the constraint register retains the last valid value Ω(n-1).

[0169] Level 3 fault (communication link interruption for more than 30 seconds): Each controller automatically switches to the conservative safety default value in the local default memory 10, and triggers a system-level alarm.

[0170] Step S3: Scheduling Plan Layer Control Based on Ω(n)

[0171] S3.1 Monthly Thermal Storage Plan Optimization

[0172] Optimization objective: minΣ(t=1 to T){w_1×max(0,P_pv(t)-P_c(t))+w_2×|Q_sup(t)-Q_dem(t)|+w_3×|P_c(t)-P_c(t-1)|} (Equation 301)

[0173] In the formula: T is the total number of time slots in the optimization time domain, dimensionless; w_1 is the curtailment penalty weight coefficient, dimensionless; w_2 is the heating deviation penalty weight coefficient, dimensionless; w_3 is the charging power fluctuation penalty weight coefficient, dimensionless; P_pv(t) is the photovoltaic predicted power in time slot t, in watts; P_c(t) is the charging power decision variable in time slot t, in watts; Q_sup(t) is the heating power decision variable in time slot t, in watts; Q_dem(t) is the heating network demand predicted power in time slot t, in watts; max(·,·) indicates taking the larger of the two; |·| indicates taking the absolute value.

[0174] State transition equation for charged state:

[0175] SOC(t+1)=SOC(t)+[η_c×P_c(t)×Δt-Q_sup(t) / η_d] / E_perm(n) (Equation 302)

[0176] In the formula: SOC(t) is the state of charge in time slot t, dimensionless; η_c is the charging efficiency, dimensionless; η_d is the energy release efficiency, dimensionless; Δt is the time slot duration in seconds; E_perm(n) is the upper limit of the permitted domain energy in joules, read from the constraint register.

[0177] Constraints (all read from the constraint register, not estimated independently at this layer):

[0178] 0≤SOC(t)≤ρ(n)×SOC_max(n) (Equation 303)

[0179] In the formula: ρ(n) is the charging depth correction coefficient, which is dimensionless and defined in Equation 307.

[0180] 0≤P_c(t)≤P_ch_max(n) (Equation 304)

[0181] Q_sup(t)≥Q_con_min×Δt (Equation 305)

[0182] In the formula: Q_con_min is the minimum continuous heating power specified in the contract, and the unit is watts.

[0183] The energy storage and release plan was recently refined and generated within the framework of the monthly thermal energy storage plan.

[0184] S3.2 Dynamic Energy Storage Cost Control Logic

[0185] LCOE(n)={Σ(i=1 to n)[C_mat(i)×(1+r)^(ni)]+C_aux×(n / N_design)} / Σ(i=1 to n)Q_c(i) (Equation 306)

[0186] In the formula: LCOE(n) is the cumulative cost per kilowatt-hour of electricity in the nth cycle, in monetary units per joule; C_mat(i) is the incremental cost of materials in the i-th cycle, in monetary units; r is the discount rate, dimensionless; C_aux is the total cost of auxiliary equipment, in monetary units; N_design is the number of design life cycles, dimensionless.

[0187] The charging depth correction factor ρ(n) is determined in three segments:

[0188] When α_filt(n) ≥ α_A: ρ(n) = 1.0

[0189] When α_B≤α_filt(n)<α_A: ρ(n)=0.85

[0190] When α_filt(n) < α_B: ρ(n) = 0.70 (Equation 307)

[0191] In the formula: ρ(n) is the charging depth reduction coefficient, which is dimensionless and reflects the monotonically non-increasing relationship that the lower the activity, the stricter the constraint; α_A=0.80 is the boundary threshold between the normal operation domain and the derating operation domain; α_B=0.50 is the boundary threshold between the derating operation domain and the restricted operation domain.

[0192] The value of ρ(n) satisfies ρ(n)|α≥α_A>ρ(n)|α_B≤α<α_A>ρ(n)|α<α_B, ensuring that the energy reduction is monotonically aggravated as the material activity decays.

[0193] Step S4: Power optimization layer model predictive control based on Ω(n)

[0194] Optimize the objective function:

[0195] minJ_MPC=Σ(k=0 to H_p-1)[w_1×e_pv 2 (k)+w_2×e_sup 2 (k)+w_3×ΔP_c 2 (k)]+w_4×e_soc_end 2 (Formula 401)

[0196] In the formula: J_MPC is the model predictive control objective function value, dimensionless (after weight normalization); H_p is the prediction time domain length, dimensionless; e_pv(k) is the amount of light wasted at prediction time k, in watts; e_sup(k) is the heating deviation at prediction time k, in watts; ΔP_c(k) is the change in charging power at prediction time k, in watts; e_soc_end is the state of charge deviation at the prediction terminal time, dimensionless; w_1 to w_4 are weighting coefficients, dimensionless.

[0197] Constraints (all read from the constraint register, not estimated independently at this layer):

[0198] 0≤P_c(k)≤P_ch_max(n) (Equation 402)

[0199] 0≤SOC(k)≤SOC_max(n) (Equation 403)

[0200] SOC(k+1)=SOC(k)+[η_c×P_c(k)×Δt-Q_sup(k) / η_d] / E_perm(n) (Equation 404)

[0201] Step S5: Adaptive control of reactor execution layer based on Ω(n)

[0202] S5.1 Feedforward Control Item

[0203] P_ff(t)=C_r×Ṫ_ref+Q_dot_rxn_mod(T_r(t),X(t),α_filt(n))+Q_dot_loss_mod(T_r(t))-Q_dot_sup_plan(t) (Formula 501)

[0204] In the formula: P_ff(t) is the feedforward control power at time t, in watts; C_r is the equivalent heat capacity of the reactor bed, in joules per Kelvin; Ṫ_ref is the time derivative of the temperature reference value, in Kelvin per second; Q_dot_rxn_mod is the reaction endothermic rate calculated based on the physical model, in watts; T_r(t) is the measured bed temperature at time t, in Kelvin; X(t) is the estimated reaction conversion rate at time t, dimensionless; Q_dot_loss_mod is the heat loss rate calculated based on the model, in watts; Q_dot_sup_plan(t) is the planned heating power at time t, in watts.

[0205] S5.2 Adaptive Proportional-Integral-Derivative Feedback Control Term

[0206] e_T(t)=T_ref(n)-T_r(t) (Formula 502)

[0207] In the formula: e_T(t) is the temperature tracking error at time t, in Kelvin.

[0208] P_fb(t)=K_p(n)×e_T(t)+K_i(n)×∫(0 to t)e_T(τ)dτ+K_d(n)×de_T(t) / dt (Equation 503)

[0209] In the formula: P_fb(t) is the feedback control power at time t, in watts; K_p(n) is the adaptive proportional gain, in watts per Kelvin; K_i(n) is the adaptive integral gain, in watts per Kelvin per second; K_d(n) is the adaptive differential gain, in watts per second per Kelvin; τ is the integral variable, in seconds.

[0210] K_p(n)=K_p0 / α_filt(n) (Formula 504)

[0211] K_i(n)=K_i0 / sqrt(α_filt(n)) (Formula 505)

[0212] K_d(n)=K_d0×sqrt(α_filt(n)) (Formula 506)

[0213] In the formula: K_p0, K_i0, and K_d0 are the reference proportional gain, reference integral gain, and reference differential gain tuned in step S1.2, respectively; sqrt(·) represents taking the arithmetic square root.

[0214] S5.3 Composite Control Output

[0215] P_cmd(t)=P_ff(t)+P_fb(t) (Formula 507)

[0216] P_cmd(t) is output to the electric heater after saturation limiting processing:

[0217] P_out(t)=min(max(0,P_cmd(t)),P_ch_max(n)) (Formula 508)

[0218] In the formula: P_cmd(t) is the composite control power command in watts; P_out(t) is the output power after limiting in watts; min(·,·) means taking the smaller of the two values; max(·,·) means taking the larger of the two values.

[0219] Step S6: Smooth switching control of heating network power

[0220] S6.1 Adaptive Calculation of Switching Transition Time

[0221] τ_trans=max(τ_min,|Q_after-Q_before| / (dQ_max×k_sf)) (Formula 601)

[0222] In the formula: τ_trans is the switching transition time in seconds; τ_min is the reactor thermal response time constant in seconds; Q_after is the target heating power after switching in watts; Q_before is the heating power before switching in watts; dQ_max is the maximum allowable rate of change of heating power in the heating network in watts per second; k_sf is the safety margin coefficient, dimensionless.

[0223] In the formula, τ_min is the reactor thermal response time constant, which is calibrated offline in step S1.2: under rated charging conditions, a step power input is applied to the electric heater, and the time taken for the average bed temperature to rise from the initial steady state to the final step value of 63.2% is recorded. This time is taken as τ_min, with a typical value of 60~180 seconds.

[0224] S6.2 S-shaped transition trajectory

[0225] Q_sup(t)=Q_before+(Q_after-Q_before)×f_S(ξ) (Equation 602)

[0226] f_S(ξ)=3ξ 2 -2ξ 3 (Formula 603)

[0227] In the formula: ξ is the normalized transition time, ξ=(t-t_sw) / τ_trans, dimensionless, with a value range of [0,1]; t_sw is the switching start time, in seconds.

[0228] S6.3 Transition Period Adaptive Monitoring

[0229] If the measured rate of change in water supply temperature exceeds the permissible value, the transition period will be extended.

[0230] τ_trans←τ_trans+300 (Formula 604)

[0231] In the formula: 300 is the extension amount, in seconds (i.e., 5 minutes).

[0232] Step S7: Material Operation Domain Hierarchical Management and Maintenance Decisions

[0233] S7.1 Material Operation Domain Three-Level Classification Determination

[0234] Based on the runtime domain determination function D(n) defined in step S2.6 (v), the following engineering implementation details are supplemented in conjunction with this step:

[0235] (a) Hysteresis mechanism for operation domain switching: To avoid frequent operation domain switching caused by the jitter of α_filt(n) near the threshold, a hysteresis bandwidth h=0.02 is introduced: the threshold for switching from A to B is 0.80, and the threshold for switching back from B to A is 0.82; the threshold for switching from B to C is 0.50, and the threshold for switching back from C to B is 0.52.

[0236] (ii) Run domain switching event log: Each time D(n) changes, a single computing node records the cycle number n, the run domain before and after the switch, the value of α_filt(n) and the timestamp in the event log for operation and maintenance auditing.

[0237] (iii) Maintenance decision support output: When D(n)=C and the number of continuous cycles exceeds N_C_max=10, the system generates a material replacement suggestion message. The message content includes the estimated number of remaining expected available cycles, the suggested replacement window period, and the current value of α_filt(n).

[0238] S7.2 Ω(n) Correction for Each Operating Domain

[0239] Domain A: Ω(n) is calculated in full according to step S2.6 without any additional corrections.

[0240] Domain B (Decreased Operation Domain):

[0241] E_perm_B(n)=β_d×E_conf×SOC_max(n) (Equation 702)

[0242] P_ch_max_B(n)=β_d×P_ch_rated×[γ+(1-γ)×α_filt(n)] (Equation 703)

[0243] In the formula: β_d is the derating coefficient, which is dimensionless, and is taken as 0.70 in this embodiment.

[0244] Domain C (Restricted Runtime Domain):

[0245] P_ch_max_C=0 (Equation 704)

[0246] P_dis_max_C(n)=β_C×P_dis_rated×α_filt(n) (Equation 705)

[0247] In the formula: P_dis_rated is the rated energy release power in watts; β_C is the domain C conservative coefficient, which is dimensionless and is taken as 0.50 in this embodiment.

[0248] In S7.2, β_d is applied to the upper limit of the permitted energy and the upper limit of the charging power, and is related in series with ρ(n) in Equation 307 of S3.2. The final effective reduction factor is β_d·ρ(n). In this embodiment, the effective reduction factor of the domain B interval is 0.70×0.85=0.595.

[0249] The thresholds α_A and α_B in Equation 307 and Equation 701 are completely shared, ensuring that the switching of discrete states in the operating domain and the switching of the energy depth reduction coefficient occur synchronously at the same active threshold, thus avoiding the situation where the timing of state switching and constraint tightening is inconsistent.

[0250] S7.3 Material Replacement Decision Support

[0251] Material replacement recommendations are provided when technical conditions T1 to T3 must be met simultaneously:

[0252] T1: α_filt(n) < α_B remains unchanged for more than N_rep iterations;

[0253] T2: α_mod(n) predicts that the value of α will remain lower than α_B for the next N_pred iterations;

[0254] T3: LCOE(n) increases continuously for more than L_rep cycles.

[0255] In the formula: N_rep is the minimum number of consecutive cycle confirmations, which is 10 in this embodiment; N_pred is the number of prediction and evaluation cycles, which is 50 in this embodiment; L_rep is the number of consecutive cost increase confirmations, which is 5 in this embodiment.

[0256] Step S8: Real-time correction of the dynamic allowable domain within the cycle based on the spatial distribution characteristics of the bed temperature field

[0257] The effective activity coefficient identification and unified operating permission domain update described in step S2 are performed after each complete energy storage and release cycle, and the update cycle is equal to the cycle duration. During the cycle, the bed may experience abnormal temperature field spatial distribution due to local accelerated sintering, vapor distribution shift, or local blockage of heat exchange pipelines. This abnormality reflects the rapid deterioration of local material activity before the end of the cycle, but it has not yet been captured by the inter-cycle identification mechanism in step S2.

[0258] In step S8, the spatial distribution characteristics of the bed temperature field are continuously monitored at a fixed sampling period during the loop. When the distribution characteristics are detected to deviate from the normal baseline by more than a preset threshold, the loop-based immediate correction of the unified operating permission domain is triggered without waiting for the loop to end. The correction is then synchronously applied to the three-layer controller via the write protection constraint register.

[0259] Steps S8 and S2 constitute a dual-timescale activity monitoring architecture: Step S2 is responsible for identifying the global effective activity coefficients and updating the entire permitted domain between loops; Step S8 is responsible for detecting local anomalies and protective correction of the permitted domain within the loop. Both use a single computing node 7 as the sole computing and distribution entity. The three-layer controller does not participate in anomaly detection within the loop and does not independently modify constraint boundaries.

[0260] S8.1 Definition of Spatial Distribution Characteristics of Bed Temperature Field

[0261] N_s temperature measuring points are uniformly arranged along the axis of the fixed-bed reactor, with N_s being no less than 5. The temperature reading of the j-th temperature measuring point at time t is denoted as T_r,j(t), where j = 1, 2, ..., N_s.

[0262] Average axial temperature of the bed:

[0263] T_r_bar(t)=[1 / N_s]×Σ(j=1 to N_s)T_r,j(t)(Equation 801)

[0264] In the formula: T_r_bar(t) is the average axial temperature of the bed at time t, in Kelvin; N_s is the number of axial temperature measurement points, dimensionless; T_r,j(t) is the temperature of the j-th temperature measurement point at time t, in Kelvin.

[0265] Standard deviation of bed axial temperature space:

[0266] σ_T(t)=sqrt{[1 / (N_s-1)]×Σ(j=1 to N_s)[T_r,j(t)-T_r_bar(t)] 2 (Equation 802)

[0267] In the formula: σ_T(t) is the spatial standard deviation of the bed axial temperature at time t, in Kelvin; sqrt{·} represents taking the arithmetic square root.

[0268] Maximum axial temperature gradient of the bed:

[0269] G_T_max(t)=max(j=1 to N_s-1){|T_r,j+1(t)-T_r,j(t)| / Δz_j}(Equation 803)

[0270] In the formula: G_T_max(t) is the maximum axial temperature gradient of the bed at time t, in Kelvin per meter; T_r,j+1(t) is the temperature of the (j+1)th temperature measuring point at time t, in Kelvin; Δz_j is the axial distance between the j-th and (j+1)-th temperature measuring points, in meters; |·| represents taking the absolute value; max(·) represents taking the maximum value among all calculated values.

[0271] Extreme temperature difference in the bed:

[0272] ΔT_range(t)=T_r_max(t)-T_r_min(t) (Equation 804)

[0273] In the formula: ΔT_range(t) is the temperature range of the bed at time t, in Kelvin; T_r_max(t) is the maximum temperature value among all temperature measurement points, in Kelvin; T_r_min(t) is the minimum temperature value among all temperature measurement points, in Kelvin.

[0274] Feature vector of spatial distribution of bed temperature field:

[0275] φ(t)=[σ_T(t),G_T_max(t),ΔT_range(t)]' (Equation 805)

[0276] In the formula: φ(t) is the spatial distribution characteristic vector of the bed temperature field at time t, which is a three-dimensional column vector; the superscript (') indicates the vector transpose.

[0277] Online establishment of the normal operating baseline in S8.2

[0278] Within the first N_b sampling periods after the start of the charging or releasing phase of each cycle, the normal operating baseline for this cycle is established using the sequence of spatial distribution characteristic vectors of the temperature field during this period. The value of N_b ranges from 30 to 120.

[0279] Baseline mean vector:

[0280] φ_base_bar=[1 / N_b]×Σ(k=1 to N_b)φ(t_k) (Equation 806)

[0281] In the formula: φ_base_bar is the baseline mean vector, which is a three-dimensional column vector; N_b is the number of baseline establishment samples, which is dimensionless; φ(t_k) is the feature vector at the k-th sampling time; t_k is the k-th sampling time, in seconds.

[0282] Baseline covariance matrix:

[0283] Σ_base=[1 / (N_b-1)]×Σ(k=1 to N_b)[φ(t_k)-φ_base_bar]×[φ(t_k)-φ_base_bar]' (Equation 807)

[0284] In the formula: Σ_base is the baseline covariance matrix, which is a 3x3 symmetric positive definite matrix; the superscript (') indicates the vector transpose.

[0285] When the minimum eigenvalue of Σ_base is less than the regularization lower bound ε_reg, it is replaced by Σ_base plus ε_reg multiplied by I_3:

[0286] Σ_base←Σ_base+ε_reg×I_3 (Equation 808)

[0287] In the formula: ε_reg is the lower bound of regularization, which is dimensionless and takes the value of 10 to the power of negative 4; I_3 is a 3x3 identity matrix.

[0288] S8.3 Intra-loop anomaly detection indicators

[0289] At each control sampling time t after the baseline is established, calculate the Mahalanobis distance of the current temperature field spatial distribution eigenvector relative to the baseline:

[0290] D_M(t)=sqrt{[φ(t)-φ_base_bar]'×Σ_base_inv×[φ(t)-φ_base_bar]} (Equation 809)

[0291] In the formula: D_M(t) is the Mahalanobis distance at time t, which is dimensionless; Σ_base_inv is the inverse matrix of the baseline covariance matrix Σ_base; sqrt{·} denotes taking the arithmetic square root; the superscript (') indicates the vector transpose. The Mahalanobis distance comprehensively considers the correlation between the three distribution characteristics and their respective normal fluctuation amplitudes.

[0292] Apply a sliding window averaging filter to the Mahalanobis distance sequence to suppress transient noise:

[0293] D_M_bar(t)=[1 / N_w]×Σ(i=0 to N_w-1)D_M(ti×T_s)(Equation 810)

[0294] In the formula: D_M_bar(t) is the filtered Mahalanobis distance at time t, which is dimensionless; N_w is the sliding window length, which is dimensionless and ranges from 5 to 15 seconds. In this invention, N_w is preferably 15 seconds; T_s is the control sampling period, which is in seconds.

[0295] S8.4 Intra-loop exception triggering conditions and graded responses

[0296] Define two levels of exception triggering thresholds D_th1 and D_th2, where D_th1 is less than D_th2.

[0297] D_th1=sqrt(χ 2 _3_99) (Equation 811)

[0298] In the formula: D_th1 is the warning threshold, which is dimensionless and has a value of approximately 3.37; χ²_3_99 is the 99th percentile value of the chi-square distribution with 3 degrees of freedom, which has a value of approximately 11.34; sqrt(·) represents taking the arithmetic square root.

[0299] D_th2=sqrt(χ 2 _3_999) (Equation 812)

[0300] In the formula: D_th2 is the protective correction trigger threshold, which is dimensionless and has a value of approximately 4.03; χ 2 _3_999 is the 99.9th percentile value of the chi-square distribution with 3 degrees of freedom, which is approximately 16.27.

[0301] Level 1 Response (Early Warning): When D_M_bar(t) is greater than or equal to D_th1 and less than D_th2, a single computing node generates a temperature field anomaly early warning log, records the current φ(t) and D_M_bar(t), does not correct the unified running permission domain, and does not write to the constraint register.

[0302] Secondary response (protective correction trigger): When D_M_bar(t) is greater than or equal to D_th2 for at least N_conf consecutive sampling periods, a single computing node immediately performs a protective correction of the permitted domain within the loop.

[0303] In the formula: N_conf is the number of consecutive confirmed samples, which is dimensionless and ranges from 3 to 10. In this embodiment, it is 5, corresponding to a 5-second confirmation time.

[0304] S8.5 Calculation method for protective modification of permission domain within the loop

[0305] After the second-level response is triggered, a single computing node calculates the intra-cycle activity reduction factor κ(t) based on the current filtered Mahalanobis distance:

[0306] κ(t)=min(1,max(κ_min,1-[(D_M_bar(t)-D_th2) / (D_sat-D_th2)]×(1-κ_min))) (Equation 813)

[0307] In the formula: κ(t) is the in-cycle activity reduction factor, dimensionless, with a value range of [κ_min,1]; κ_min is the minimum reduction factor, dimensionless, with a value range of 0.3 to 0.5, and is 0.4 in this embodiment; D_sat is the saturated Mahalanobis distance, dimensionless, and is twice D_th2, with a value of 8.06 in this embodiment; D_M_bar(t) is the current filtered Mahalanobis distance, dimensionless; D_th2 is the protective correction trigger threshold, dimensionless; max(·,·) indicates taking the larger of the two values.

[0308] Endpoint verification: When D_M_bar(t)≤D_th2, the expression within the square brackets [(D_M_bar(t)−D_th2) / (D_sat−D_th2)] is not positive, and the value of (1−[·]×(1−κ_min)) is not less than 1. After max(κ_min,·), it is still not less than 1. After truncation by min(1,·), κ(t)=1, meaning that the loop correction is not triggered. When D_M_bar(t)≥D_sat, the expression within the square brackets is not less than 1, and κ(t)=κ_min, reaching the maximum protective reduction.

[0309] The components of the unified runtime license domain after loop correction are calculated as follows:

[0310] Intra-cycle adjustment of charging power limit:

[0311] P_ch_max_intra(t)=κ(t)×P_ch_max(n) (Equation 814)

[0312] In the formula: P_ch_max_intra(t) is the upper limit of the charging power corrected within the loop, in watts; P_ch_max(n) is the upper limit of the current loop charging power calculated in step S2.6 (from the most recent inter-loop update), in watts. When the intra-loop running domain determination rule defined in the latter part of this section (S8.5) determines that D_intra(t) = domain C, the result of Equation 814 is further forced to overwrite P_ch_max_intra(t) = 0, and this overwrite takes precedence over the value calculated by Equation 814.

[0313] Intra-loop correction of the upper limit of the state of charge:

[0314] SOC_max_intra(t)=κ(t)×SOC_max(n) (Equation 815)

[0315] In the formula: SOC_max_intra(t) is the upper limit of the state of charge after correction within the loop, which is dimensionless; SOC_max(n) is the upper limit of the current state of charge calculated in step S2.6, which is dimensionless.

[0316] Upper bound of feasible temperature range within the loop:

[0317] T_upper_intra(t)=T_ref(n)+κ(t)×ΔT_tol (Equation 816)

[0318] In the formula: T_upper_intra(t) is the upper bound of the temperature feasible region after correction within the loop, in Kelvin; T_ref(n) is the temperature reference value calculated in step S2.6, in Kelvin; ΔT_tol is the temperature tolerance, in Kelvin. The lower bound of the temperature feasible region, T_ref(n), remains unchanged.

[0319] Modify the upper limit of the permitted domain energy within the loop:

[0320] E_perm_intra(t)=E_conf×SOC_max_intra(t) (Equation 817)

[0321] In the formula: E_perm_intra(t) is the upper limit of the permitted domain energy after the loop correction, in joules.

[0322] The correction of the discrete state D_intra(t) in the runtime domain within the loop is determined according to the following rules:

[0323] When κ(t)≥0.8, keep the current D(n) unchanged;

[0324] When 0.5≤κ(t)<0.8, a forced switch to domain B (reduced throttling domain) is performed without waiting for confirmation from three consecutive loops;

[0325] When κ(t) < 0.5, a forced switch to domain C (restricted operating domain) is performed, and P_ch_max_intra(t) is corrected to zero.

[0326] The modified unified runtime permission domain structure within the loop:

[0327] Ω_intra(t)={SOC_max_intra(t),E_perm_intra(t),P_ch_max_intra(t),Ω_T_intra(t),D_intra(t)} (Equation 818)

[0328] In the formula: Ω_intra(t) is the unified running permission domain structure after the loop correction; Ω_T_intra(t) is the temperature feasible domain after the loop correction, and Ω_T_intra(t) = [T_ref(n), T_upper_intra(t)] is the intersection of (-∞, T_safe_max].

[0329] Synchronous Distribution and Write Protection with In-Loop Correction in S8.6

[0330] The synchronous distribution mechanism of Ω_intra(t) is completely consistent with the distribution mechanism of Ω(n) described in step S2.8, using the same industrial real-time communication interface, the same communication frame format, the same message authentication code verification, and the same write protection constraint register write channel.

[0331] The communication frame version number field uses a composite encoding (n, m) of the loop count n and the correction sequence number m within the loop, where m increments from 1. After receiving the Ω_intra(t) communication frame, the three-layer controller verifies the frame integrity using the same verification process as when receiving Ω(n). If the verification passes, it writes each component of Ω_intra(t) into the constraint register, overwriting the previous Ω(n) or the last Ω_intra value.

[0332] The three-layer controller does not participate in the calculation of exception determination within the loop, and does not independently decide whether to accept or reject the correction within the loop.

[0333] The communication latency requirement for loop-based correction is that the maximum allowable latency Δt_intra from the second-level response trigger confirmation to the completion of the write to the third-level controller constraint register does not exceed 500 milliseconds. This latency requirement is stricter than the latency requirement for inter-loop updates in step S2.8.

[0334] S8.7 Coordination rules for intra-loop modifications and inter-loop updates

[0335] When the loop correction is active (i.e., κ(t) < 1), if the current loop terminates normally and triggers the inter-loop identification and full update of Ω(n+1) in step S2, then the following four rules apply:

[0336] Rule 1: The inter-loop identification in step S2 is executed normally, outputting α_filt(n+1) and Ω(n+1).

[0337] Rule 2: After Ω(n+1) is written to the constraint register, the correction state within the loop is automatically reset, κ is reset to 1, Ω_intra is invalidated, and the system runs with the new Ω(n+1) as the constraint.

[0338] Rule 3: If the peak Mahalanobis distance D_M_bar_peak after the correction trigger within the loop exceeds 1.5 times D_th2, then a penalty weight is applied to the original calculated value of the effective activity coefficient for this loop in the exponential weighted filtering of step S2.5.

[0339] α_filt(n+1)=λ_p×α_raw(n+1)+(1-λ_p)×α_filt(n) (Equation 819)

[0340] In the formula: α_filt(n+1) is the effective activity coefficient after filtering in the (n+1)th cycle under penalty weight, which is dimensionless; α_raw(n+1) is the original calculated value of the effective activity coefficient in the (n+1)th cycle, which is dimensionless; α_filt(n) is the effective activity coefficient after filtering in the nth cycle, which is dimensionless; λ_p is the weight coefficient after penalty, which is dimensionless.

[0341] λ_p=min(1,λ×D_th2 / D_M_bar_peak) (Equation 820)

[0342] In the formula: λ is the standard weight coefficient defined in formula 209 of step S2.5, which is dimensionless; D_th2 is the protective correction trigger threshold, which is dimensionless; D_M_bar_peak is the peak value of the filtered Mahalanobis distance during the correction period in this loop, which is dimensionless; min(·,·) means taking the smaller of the two.

[0343] The more severe the anomaly within the loop (the larger D_M_bar_peak), the smaller λ_p becomes. This reduces the filter's confidence in the current loop's identification value, retaining more of the previous filter value to prevent excessive disturbances to the overall decay trend estimation caused by the abnormal loop's identification value. Endpoint verification: When D_M_bar_peak equals D_th2, λ_p equals λ, degenerating into standard filtering with no penalty effect; when D_M_bar_peak is much larger than D_th2, λ_p approaches zero, almost completely retaining the previous filter value.

[0344] Rule 4: If the current loop fails to complete the energy release process normally after the correction within the loop is triggered (e.g., the loop is forced to stop because the upper limit of charging power is corrected to zero due to the limited operating domain), then the inter-loop identification in step S2 skips this loop, α_filt remains α_filt(n) and is not updated, the loop counter n is not incremented, and Ω(n) remains unchanged. The next complete loop that completes normally resumes the normal identification process of step S2.

[0345] Instantaneous adaptation of the underlying controller gain during S8.8 in-cycle correction

[0346] When the in-loop modification is active, define the equivalent activity coefficient within the loop:

[0347] α_eq_intra(t)=κ(t)×α_filt(n) (Equation 821)

[0348] In the formula: α_eq_intra(t) is the equivalent activity coefficient within the loop, dimensionless; κ(t) is the activity reduction factor within the loop, dimensionless; α_filt(n) is the filtered effective activity coefficient updated in the most recent inter-loop update, dimensionless.

[0349] The proportional, integral, and derivative gains of the reactor's bottom-level controller are adaptively adjusted in real time based on the equivalent activity coefficient within the cycle.

[0350] K_p_intra(t)=K_p0 / α_eq_intra(t) (Equation 822)

[0351] In the formula: K_p_intra(t) is the adaptive proportional gain within the loop, in watts per Kelvin; K_p0 is the base proportional gain, in watts per Kelvin.

[0352] K_i_intra(t)=K_i0 / sqrt(α_eq_intra(t)) (Equation 823)

[0353] In the formula: K_i_intra(t) is the adaptive integral gain within the loop, in watts per Kelvin per second; K_i0 is the baseline integral gain, in watts per Kelvin per second; sqrt(·) represents taking the arithmetic square root.

[0354] K_d_intra(t)=K_d0×sqrt(α_eq_intra(t)) (Equation 824)

[0355] In the formula: K_d_intra(t) is the adaptive differential gain within the cycle, in watt-seconds per Kelvin; K_d0 is the baseline differential gain, in watt-seconds per Kelvin.

[0356] Equations 822 to 824 are completely consistent with the adjustment directions of Equations 504 to 506 in step S5.2: the proportional gain increases monotonically as the equivalent activity coefficient decreases, the integral gain increases weakly in the same direction, and the differential gain decreases monotonically.

[0357] Gain rate of change limiting: To prevent rapid changes in κ(t) from causing a gain step and resulting in abrupt changes in control commands, rate of change limiting is applied to the three gains:

[0358] |dK_p_intra / dt|≤K_dot_p_max (Equation 825)

[0359] In the formula: dK_p_intra / dt is the time rate of change of the proportional gain, in watts per Kelvin per second; K_dot_p_max is the maximum allowable rate of change of the proportional gain, in watts per Kelvin per second; |·| indicates taking the absolute value. In this embodiment, K_dot_p_max = 0.2 × K_p0 per second.

[0360] |dK_i_intra / dt|≤K_dot_i_max (Equation 826)

[0361] In the formula: dK_i_intra / dt is the time rate of change of the integral gain; K_dot_i_max is the maximum allowable rate of change of the integral gain. In this embodiment, K_dot_i_max = 0.1 × K_i0 per second.

[0362] |dK_d_intra / dt|≤K_dot_d_max (Equation 827)

[0363] In the formula: dK_d_intra / dt is the time rate of change of the differential gain; K_dot_d_max is the maximum allowable rate of change of the differential gain. In this embodiment, K_dot_d_max = 0.1 × K_d0 per second.

[0364] After the intra-loop correction is reset (κ is reset to 1), the gain smoothly recovers to the inter-loop gain value driven by α_filt(n) in step S5.2 according to the above-mentioned rate of change limiting condition, without a step jump.

[0365] Synchronization compensation of feedforward control terms during S8.9 intra-cycle correction

[0366] When the in-cycle correction is active, the reaction endothermic compensation power of the feedforward control term in step S5.1 is replaced by the in-cycle equivalent activity coefficient α_eq_intra(t) instead of α_filt(n):

[0367] Q_dot_rxn_intra(t)=(m_CaO / M_CaO)×ΔH_r×k_0×exp(-E_a / (R×T_r(t)))×α_eq_intra(t)×(1-X(t))^m×g(P_steam,T_r(t)) (Equation 828)

[0368] In the formula: Q_dot_rxn_intra(t) is the corrected feedforward reaction endothermic rate within the cycle, in watts; m_CaO is the total mass of calcium oxide, in kilograms; M_CaO is the molar mass of calcium oxide, in grams per mole; ΔH_r is the enthalpy of reaction, in joules per mole; k_0 is the frequency factor, in seconds; exp(·) represents an exponential function with the natural constant e as the base; E_a is the apparent activation energy, in joules per mole; R is the constant of the flow rate. The gas constant is expressed in joules per mole per Kelvin; T_r(t) is the measured bed temperature at time t, in Kelvin; α_eq_intra(t) is the equivalent activity coefficient within the cycle, dimensionless; X(t) is the estimated reaction conversion rate at time t, dimensionless; m is the geometric exponent of the contraction core model, dimensionless; g(P_steam,T_r(t)) is a function of water vapor partial pressure and temperature, dimensionless; P_steam is the water vapor partial pressure, in Pascals.

[0369] Corrected feedforward control power:

[0370] P_ff_intra(t)=C_r×Ṫ_ref+Q_dot_rxn_intra(t)+Q_dot_loss_mod(T_r(t))-Q_dot_sup_plan(t) (Equation 829)

[0371] In the formula: the meaning of each symbol is the same as in formula 501, and Q_dot_rxn_mod is replaced with Q_dot_rxn_intra(t).

[0372] The revised upper limit of the output saturation limit for the composite control is now synchronously changed to P_ch_max_intra(t):

[0373] P_out_intra(t)=min(max(0,P_ff_intra(t)+P_fb_intra(t)),P_ch_max_intra(t)) (Equation 830)

[0374] In the formula: P_out_intra(t) is the output power after limiting during the correction period within the loop, in watts; P_fb_intra(t) is the feedback control power calculated by K_p_intra(t), K_i_intra(t), and K_d_intra(t), in watts.

[0375] S8.10 Intra-cycle Correction Constraint Transmission to Power Optimization Layer and Scheduling Layer

[0376] Since Ω_intra(t) is synchronously distributed via the write-protected constraint register, the power optimization layer and scheduling layer can automatically receive the corrected constraint boundaries within the loop without additional logic modifications.

[0377] At the next rolling optimization trigger moment, the power scheduling controller 5 reads the currently valid P_ch_max_intra(t) and SOC_max_intra(t) from the constraint register and substitutes them into the constraint conditions of equations 402 and 403 in step S4 to automatically tighten the feasible region of the rolling optimization problem.

[0378] At the next plan revision trigger time, the scheduling plan controller 4 reads the current valid value from the constraint register and substitutes it into the constraint conditions of equations 303 and 304 in step S3. If the revision within the loop causes P_ch_max_intra(t) to be lower than the charging power target value of the current planned time slot, the scheduling plan controller automatically truncates the charging power target value based on the constraint register value, without needing to independently determine the reason for truncation.

[0379] The protective modification of the permitted domain within the loop only tightens (does not relax) the constraint boundary, without destroying the inter-layer constraint inclusion relationship established in step S2.7. The proof is as follows:

[0380] From Equation 814, P_ch_max_intra(t) = κ(t) × P_ch_max(n), where 0 < κ(t) ≤ 1, therefore P_ch_max_intra(t) ≤ P_ch_max(n). From the proof of Theorem 1 in step S2.7, any power trajectory satisfying P_c(t) ≤ P_ch_max(n) belongs to the execution layer's physically reachable set R_exec(n). Since a power trajectory satisfying P_c(t) ≤ P_ch_max_intra(t) ≤ P_ch_max(n) necessarily belongs to R_exec(n), Theorem 1 continues to hold in the loop-corrected scenario.

[0381] Relationship between S8.11 In-cycle correction and bed temperature rise rate safety constraints

[0382] From Equation 214, P_ch_max(n) is calibrated according to the activity coefficient α_filt(n) to make the bed temperature rise rate exactly equal to dTr_safe. When the local activity within the cycle is reduced to κ(t)·α_filt(n), the rate of decrease in the reaction endothermic rate in the abnormal region is (1-κ(t))·Q_dot_rxn_nom. From Equation 814, P_ch_max_intra(t)=κ(t)×P_ch_max(n), and the corresponding charging-side input power also decreases proportionally to κ(t). Substituting the temperature rise balance relationship in Equation 214:

[0383] C_r×(dT_r / dt)_intra_max=P_ch_max_intra(t)-κ(t)·Q_dot_rxn_nom-Q_dot_loss_nom (Equation 831)

[0384] Substituting P_ch_max_intra(t)=κ(t)·(C_r·dTr_safe+Q_dot_rxn_nom+Q_dot_loss_nom) into the equation and rearranging, we get:

[0385] (dT_r / dt)_intra_max=κ(t)·dTr_safe+[(κ(t)-1)·Q_dot_loss_nom / C_r]≤dTr_safe (Equation 832)

[0386] In the formula: (dT_r / dt)_intra_max is the upper bound of the bed temperature rise rate after correction within the cycle, in Kelvin per second; the definitions of other symbols are the same as before.

[0387] It should be noted that the derivation of Equation 832 is based on the following simplifying assumption: the equivalent reduction ratio of the reaction endothermic rate within the anomalous region is consistent with the global reduction coefficient κ(t). Strictly speaking, κ(t) is obtained by mapping the Mahalanobis distance D_M_bar(t) of the global bed temperature field (Equation 813), while the proportion of the impact of local activity degradation on the global reaction endothermic rate depends on the volume ratio ξ_loc of the anomalous region. When the volume ratio of the anomalous region is small (ξ_loc≤20%), κ(t) shows a conservative deviation relative to the true equivalent reduction coefficient, meaning that the upper bound of the temperature rise rate given by Equation 832 is still a safe upper bound; when the volume ratio of the anomalous region is large, this conservatism is enhanced. Therefore, the inequality in Equation 832 always holds in engineering practice, with the only trade-off being a higher degree of conservatism when the anomalous region is extremely large. This trade-off is allowed by the design to achieve model simplicity and real-time computational efficiency.

[0388] Equation 832 shows that the maximum temperature rise rate corresponding to the upper limit of the charging power after the protective correction within the cycle strictly does not exceed the upper limit of the safe temperature rise rate, and the safety constraint continues to be guaranteed in the scenario of correction within the cycle.

[0389] In the cycle-corrected activation state, the reaction endothermic rate in the feedforward compensation term is reduced proportionally to α_eq_intra(t)=κ(t)×α_filt(n). This reduction relationship is based on the engineering equivalence assumption: κ(t) itself is a global dimensionless reduction factor based on the spatial distribution characteristics of the bed temperature field. Using it to reduce the reaction endothermic rate will slightly overestimate the decrease in effective endothermic capacity. The resulting deviation is a conservative deviation (i.e., the compensation amount is too small and the bed temperature rise is too fast, leaving a safety margin), so it does not affect the system safety boundary.

[0390] Based on the above steps, the theoretical upper limit of the total response delay from the occurrence of an exception within a loop to the completion of the constraint register update is:

[0391] τ_total_max=T_s / 2+T_win+N_conf×T_s+T_calc+T_comm+T_jitter (Equation 833)

[0392] In the formula: T_s is the control sampling period, taken as 1 second; T_win is the Mahalanobis distance sliding window alignment delay, taken as 15 seconds; N_conf is the number of consecutive acknowledgment sampling periods, taken as 5; T_calc is the computation delay of a single computing node executing κ(t), taken as 1 second; T_comm is the industrial real-time communication interface write delay, taken as 0.5 seconds; T_jitter is the scheduling jitter margin of communication and computing tasks, taken as 0.5 seconds. Substituting, we get:

[0393] τ_total_max=0.5+15+5+1+0.5+0.5=22.5 seconds.

[0394] Impact of S8.12 Intra-cycle Correction on Power Grid Dispatch Coordination

[0395] When this invention is applied to a photovoltaic coupled grid dispatch scenario, the intra-cycle correction has the following specific impacts on the coordinated operation of the grid side.

[0396] Impact on real-time power dispatch: When the intra-loop correction is triggered, the power dispatch controller 5 reads P_ch_max_intra(t) from the constraint register as the upper limit of the frequency regulation response power. If the charging power frequency regulation command P_AGC(t) issued by the automatic generation control system exceeds P_ch_max_intra(t), the power dispatch controller 5 will truncate the response power to P_ch_max_intra(t) and report the current adjustable capacity update information to the grid energy management system through the reverse communication channel.

[0397] P_report(t) = P_ch_max_intra(t) (Equation 834)

[0398] In the formula: P_ch_max_intra(t) is the upper limit of the charging power after correction within the loop, in watts. P_report(t)_val: the power value field in the communication message, in watts; P_report(t)_ts: the timestamp field in the communication message; P_report(t)_crc: the checksum field in the communication message. The complete message P_report(t)={P_report(t)_val,P_report(t)_ts,P_report(t)_crc} is encapsulated by a single computing node and sent through the industrial real-time communication interface.

[0399] This information enables the power grid dispatch center to promptly learn of temporary reductions in the adjustable capacity of energy storage units due to intracycle anomalies, thus avoiding the continuous issuance of unexecutable frequency regulation commands based on outdated capacity information.

[0400] Impact on frequency response: When the grid frequency is lower than the rated frequency and the frequency deviation exceeds the frequency threshold, the available margin for reducing charging power by the power dispatch controller 5 is constrained by P_ch_max_intra(t). When κ(t) is small, i.e., the anomaly is more severe, the charging power itself is already at a low level or zero, and the amount of grid load that can be released is correspondingly reduced; the power dispatch controller 5 automatically reports the available capacity of the reduced frequency response to the grid.

[0401] Impact on economic optimization operation: When the duration of the correction within the cycle exceeds one scheduling time slot (15 minutes), the scheduling plan controller 4 will solve the problem in the next time slot with the tightened P_ch_max_intra(t) and SOC_max_intra(t) as constraints. The charging energy during the low electricity price period will be automatically reduced to avoid safety risks caused by forced charging under abnormal material conditions.

[0402] Summary of S8.13 In-loop Correction Parameters

[0403] The complete set of parameters and their preferred values ​​for step S8 are summarized below.

[0404] The number of axial temperature measurement points N_s is dimensionless, with a preferred value of 5. Note: This is consistent with step S1.2.

[0405] The baseline establishment sampling number N_b is dimensionless, with a preferred value of 60. This means that the baseline is established in 60 seconds.

[0406] The lower bound for covariance regularization, ε_reg, is dimensionless, with a preferred value of 10 to the power of -4. This is to prevent singularities in the covariance matrix.

[0407] The sliding window length N_w is dimensionless, with a preferred value of 15 seconds. This corresponds to a filtering window with 15 sampling points.

[0408] The warning Mahalanobis distance threshold D_th1 is dimensionless, with a preferred value of 3.37. It corresponds to the square root of the 99th percentile of the chi-square distribution with 3 degrees of freedom.

[0409] The protective correction trigger threshold D_th2 is dimensionless, with a preferred value of 4.03. It corresponds to the square root of the 99.9th percentile of the chi-square distribution with 3 degrees of freedom.

[0410] The number of consecutive confirmed samples N_conf is dimensionless, with a preferred value of 5, which corresponds to a 5-second confirmation time.

[0411] The minimum reduction factor κ_min is dimensionless, with an optimal value of 0.4. This means that the maximum reduction is 60%.

[0412] The saturated Mahalanobis distance D_sat is dimensionless, with a preferred value of 8.06. Note: It is taken as twice D_th2.

[0413] The maximum delay Δt_intra for distribution within the loop is in milliseconds, with a preferred value of 500, indicating a high time urgency requirement.

[0414] The penalty weight triggers the Mahalanobis distance multiple, which is dimensionless and has a preferred value of 1.5. Note: This is the triggering condition for Rule 3 of S8.7.

[0415] The maximum rate of change of proportional gain is K_dot_p_max, measured in watts per Kelvin per second; the maximum rate of change of integral gain is K_dot_i_max, measured in watts per Kelvin per second squared; and the maximum rate of change of differential gain is K_dot_d_max, measured in watts per Kelvin. The typical numerical ranges for these three values ​​are K_dot_p_max ∈ [0.5, 2.0], K_dot_i_max ∈ [0.05, 0.2], and K_dot_d_max ∈ [1.0, 5.0], respectively. These values ​​are tuned according to the specific reactor's thermal time constant during the engineering commissioning phase.

[0416] III. Examples and Simulation Verification

[0417] Unified naming conventions for the scheme

[0418] To facilitate comparison and verification of the technical contributions of each component step of the present invention, the following naming convention is uniformly adopted for all embodiments in this specification, and the entire text follows this convention without being changed due to different embodiments:

[0419] Solution A: The complete solution of this invention, with steps S1 to S8 fully enabled, serves as a performance benchmark.

[0420] Scheme A′: The ablation comparison scheme of the present invention enables steps S1 to S7, disables step S8 (i.e. disables the real-time correction mechanism of the dynamic permission domain based on the spatial distribution characteristics of the bed temperature field within the cycle), and the rest of the configuration is completely consistent with Scheme A, which is used to isolate and evaluate the technical contribution of step S8.

[0421] Option B: Fixed parameter comparison scheme, disables inter-cycle online identification in step S2 and intra-cycle correction in step S8, uses the initially calibrated activity coefficient and permission domain throughout, and the rest of the configuration is the same as Option A, used to evaluate the technical contribution of inter-cycle identification.

[0422] Option C: Independent hierarchical comparison scheme, where the three controllers independently maintain the activity estimation and constraint boundaries, there is no unified operating license domain, and a single computing node does not participate in constraint distribution, used to compare with existing technical architectures.

[0423] Example 1: Simulation Platform and Common Parameters

[0424] The simulation platform was built using MATLAB / Simulink R2022a and Aspen Custom Modeler V12. The fixed-bed reactor employed a one-dimensional unsteady-state distributed parameter model, discretized along the axial direction into 50 finite volume elements with a time step of 1 second. Eight temperature measurement points were uniformly arranged along the reactor's axial direction for the bed temperature field acquisition device. The photovoltaic power output curve was based on the measured power curve of a typical sunny day from a 50MW photovoltaic power plant (sampling interval of 1 minute). The heating network load curve was based on the typical daily heating load curve of an industrial park in northern China.

[0425] Common parameters: reactor loading 67000kg; electric heater rated power 1.2MW; heat exchanger heat exchange area 80m²; initial effective activity coefficient α_filt(0)=1.00; asymptotic activity coefficient lower limit α_inf=0.65; scheduling plan cycle 24 hours; power scheduling rolling cycle 15 minutes; bottom control sampling cycle 1 second; bed temperature field anomaly detection sampling cycle T_s=1 second; Mahalanobis distance sliding window length 15 seconds; operating domain threshold α_A=0.80, α_B=0.50; Mahalanobis distance threshold D_th1=3.37, D_th2=4.03, D_sat=8.06.

[0426] Example 2a: Performance verification under long-term cyclic conditions

[0427] This embodiment verifies the performance of the complete scheme of the present invention compared with the fixed parameter scheme and the independent layered scheme under long-term cyclic conditions (200 consecutive complete energy storage and release cycles without human intervention or local abnormal events within the cycle).

[0428] The schemes being compared are Scheme A, Scheme B, and Scheme C. The naming conventions for these schemes are defined in the "Unified Naming Conventions for Schemes" section at the beginning of this chapter. In this embodiment, because no local exception events are injected within the loop, the exception monitoring module in step S8 of Scheme A continues to run, but the protective correction is not triggered.

[0429] The simulation results are as follows:

[0430] Indicator J1 (interlayer constraint consistency index at the end of 200 cycles): defined as the maximum relative deviation of the P_ch_max values ​​in the constraint registers of the three layers of controllers at the same time. Scheme A is 0.0% (completely consistent across all three layers), Scheme B is 0.0% (consistent despite being fixed), and Scheme C is 7.8% (significantly inconsistent across layers).

[0431] Indicator J2 (relative value of cumulative wasted light over 200 cycles): 1.00 for Scheme A, 1.18 for Scheme B (the schedulable capacity is underestimated because the constraints are not updated with activity decay), and 1.09 for Scheme C.

[0432] Indicator J3 (number of times the event "upper-level plan is feasible but lower-level physical inaccessibility" occurs in 200 iterations): Option A is 0 times, Option B is 0 times (conservative constraints prevent out-of-bounds behavior), and Option C is 23 times (inconsistency between levels directly leads to plan unexecutability).

[0433] Indicator J4 (activity estimation deviation, defined as the absolute deviation between the activity coefficient estimate used by each control layer and the actual effective activity of the material): In scheme A, each layer independently maintains the initial activity parameter α_init=1.00, while the actual effective activity of the material α_actual=0.82, so J4=|α_init−α_actual|=0.18; In the scheme of this invention, α_filt(n)=0.823 is identified and uniformly distributed by a single computing node, so J4≈|0.823−0.82|=0.003; Scheme B is 0.18 (obtained from |α_init−α_actual|=|1.00−0.82|, scheme B does not include online identification, and this deviation accumulates monotonically, reaching more than 12 times that of scheme A by the 87th cycle); Scheme C is 0.087.

[0434] The above results show that Scheme A, through a unified operating permission domain mechanism, simultaneously achieves inter-layer constraint consistency (J1, J3) and full utilization of schedulable capacity (J2) under long-term cyclic operating conditions, while the relatively independent layered architecture (Scheme C) fundamentally eliminates the inter-layer inconsistency problem.

[0435] Example 2b: Performance verification under local abnormal conditions within the loop

[0436] This embodiment verifies the performance difference between the complete solution of the present invention and the ablation solution under the condition of local sintering mutation event injected in the cycle, and focuses on verifying the technical contribution of step S8.

[0437] The schemes used for comparison are Scheme A and Scheme A′. The naming conventions for the schemes are defined in the "Unified Naming Conventions for Schemes" at the beginning of this chapter. The simulation was executed for a total of 50 complete cycles. During the charging phase of the 10th, 20th, 30th, and 40th cycles, a local sintering mutation event was injected. The event lasted for about 5 to 10 minutes, and the event intensity corresponded to D_M_bar_peak between 5.5 and 7.5.

[0438] The evaluation indicators and comparison results are as follows:

[0439] Evaluation index I1 (maximum bed temperature shift during abnormal periods): Scheme A′ is 18.3K, Scheme A is 4.7K, a reduction of 74.3%.

[0440] Evaluation index I2 (integrated mean of spatial non-uniformity of bed temperature field within the circulation system): defined as the time integral mean of the spatial standard deviation of the bed temperature field during the abnormal event. Scheme A is 0.029, and Scheme A′ is 0.033; Scheme A represents a 12.1% reduction compared to Scheme A′.

[0441] Evaluation index I3 (cumulative duration of local over-temperature and over-limit within the cycle): defined as the cumulative number of seconds during which the temperature at any temperature measurement point in the bed exceeds the safe upper limit during an abnormal event. Scheme A is 0.0 seconds, and Scheme A′ is 0.8 seconds; Scheme A completely eliminates local over-temperature and over-limit compared to Scheme A′.

[0442] Evaluation metric I4 (Normalized Cumulative Energy Deviation): Defined as the normalized value of the cumulative energy execution deviation relative to the ideal planned value. Scheme A is 1.00 (baseline), and Scheme A′ is 1.02; Scheme A is 2.0% lower than Scheme A′. The difference in this metric is small because although the ablation scheme lacks in-cycle correction, constraint consistency can still be restored through inter-cycle identification after the abnormal event ends. Therefore, the difference in cumulative energy deviation mainly comes from short-term execution deviations during the abnormal event.

[0443] Evaluation metric I5 (Detection rate of local over-temperature events within the cycle): defined as the proportion of injected local abnormal events that are correctly identified by the control system and trigger a protective response before the end of the cycle, out of the total number of injected events. Scheme A is 100% (4 / 4), and Scheme A' is 0% (0 / 4); Scheme A represents a 100 percentage point improvement over Scheme A'.

[0444] Evaluation metric I6 (average delay from anomaly to permissioned domain correction taking effect): defined as the average delay between the actual start time of the anomaly event and the completion time of the constraint register update. Scheme A has a delay of 7.3 seconds, while Scheme A′ is not applicable (no in-loop correction mechanism). Note: I6 = 7.3 seconds is significantly less than the 22.5-second upper limit of response delay given in section S8.11. This is because the anomalous event injected in this simulation involves a large abrupt change in the bed temperature field, and the Mahalanobis distance D_M(t) exceeds 1.5 times the D_th2 threshold within approximately 3 sampling periods after the anomaly occurs, resulting in a significantly faster trigger response time than the critical case. The theoretical upper limit of 22.5 seconds in section S8.11 corresponds to the critical case where the anomaly amplitude is exactly equal to D_th2, and is a conservative estimate.

[0445] I7 (Maximum bed temperature exceeding the upper limit of safety during abnormal periods): Scheme A′ is 4.8%, Scheme A is 1.3%, and Scheme B is 4.8%, with a relative decrease of approximately 72.9%. This is on the same order of magnitude as the 74.3% decrease in inter-layer constraint violation reflected in I1, verifying the inhibitory effect of inter-layer constraint consistency on plan-execution deviation.

[0446] Based on the above seven indicators, Scheme A (the complete scheme of this invention, including step S8) shows a significant performance advantage over Scheme A′ (ablation scheme, with step S8 disabled) under local abnormal conditions within the cycle:

[0447] (1) Safety aspects: Scheme A reduces the maximum bed temperature deviation in the cycle by 74.3% (I7), reduces the cumulative duration of local over-temperature events from 0.8 seconds to 0 seconds (I3), and increases the detection rate of local over-temperature events in the cycle from 0% to 100% (I5), fundamentally eliminating the safety blind spot caused by local anomalies in the cycle.

[0448] (2) Control performance: Scheme A reduces the mean integral of the spatial non-uniformity of the bed temperature field by 12.1% (I2) and reduces the normalized cumulative energy deviation by 2.0% (I4), indicating that the dynamic correction mechanism within the cycle in step S8 does not significantly sacrifice energy execution accuracy while ensuring safety.

[0449] (3) Response timeliness: The average delay from the abnormality to the effective correction of the permission domain in Scheme A is only 7.3 seconds (I6), which is far lower than the theoretical upper limit of 22.5 seconds given in S8.11, and meets the timeliness requirements of real-time protection of abnormality in the loop.

[0450] The above comparison demonstrates that the real-time correction mechanism for the dynamic permissive domain within the loop, based on the Mahalanobis distance of the characteristic vector of the bed temperature field spatial distribution, introduced in step S8, is a necessary technical means to eliminate the safety blind zone of local anomalies within the loop. Disabling this step would cause the control system to be unable to respond to local sudden events during the loop, significantly increasing the safety risk. The dual-timescale activity monitoring architecture formed by the inter-loop identification in steps S8 and S2 is the key to achieving the dual technical effects of "consistent constraints between loops + no blind zone of safety within the loop" in this invention.

[0451] Example 3: Detailed time sequence of photovoltaic coupling operation in a single day (24 hours)

[0452] This embodiment provides the complete timeline from 0:00 to 0:00 the following day on a typical sunny day (using scheme A), focusing on demonstrating the response process of step S8 under abnormal events within the loop. The loop number n=87 for this day, and at the beginning of the loop α_filt(87)=0.82, the corresponding running domain is located in domain A (α_filt≥α_A=0.80). The unified running permission domain Ω(87) has been generated by step S2 at the end of the 86th loop and written into the three-level controller constraint register.

[0453] Phase 1 (0:00 to 6:00, nighttime energy release phase):

[0454] The energy release phase continues, with the heat network load supplied by the energy storage system. The reactor bottom controller 6 performs bottom control according to P_dis_max(87) and T_lower(87) in Ω(87). The Mahalanobis distance D_M_bar(t) in step S8 has an average value of 0.8 in this phase, which is much lower than D_th1=3.37, and there is no abnormal triggering.

[0455] Phase 2 (6:00 AM to 11:00 AM, photovoltaic output increases and charging starts):

[0456] At 6:00 AM, the photovoltaic output began to increase, and the scheduling controller 4 formulated a charging plan according to SOC_max(87) and P_ch_max(87) in Ω(87). The power scheduling controller 5 read the latest constraint value from the constraint register every 15 minutes at the optimization trigger time. At 8:30 AM, the system switched from the energy release mode to the energy charging mode, and performed an S-shaped transition according to Equations 601 to 604. The entire process was completed within 7.1 seconds (consistent with T5−T0≤7.1 seconds given in S2.8). At 10:00 AM, the photovoltaic output reached the peak of the day, and the electric heater 2 input heat to the reactor at full load with P_ch_max(87).

[0457] Phase 3 (11:00 to 15:00, local sintering mutation event within the cycle):

[0458] 11:15:10: A sudden change in local sintering occurred in the area of ​​the 3rd temperature measurement point, and the temperature rise rate at this point was significantly faster than that of the adjacent measurement points.

[0459] From 11:15:10 to 11:15:25 (Hampel filtering and sliding window alignment delay T_win = 15 seconds): D_M_bar(t) gradually increased from the normal level to 5.8, exceeding D_th2 = 4.03 for the first time.

[0460] From 11:15:25 to 11:15:30 (continuous acknowledgment delay N_conf×T_s=5×1=5 seconds): D_M_bar(t) remains above D_th2 within this window, satisfying the continuous acknowledgment criterion of S8.4, triggering a level 2 response.

[0461] From 11:15:30 to 11:15:31 (κ(t) calculation delay T_calc ≤ 1 second): a single computing node calculates κ(t) according to Equation 813:

[0462] κ(t)=min(1,max(κ_min,1-[(5.8−4.03) / (D_sat−4.03)]×(1−0.4)))

[0463] Taking D_sat = 8.06 (= 2 × D_th2), we calculate κ(t) ≈ 0.7365.

[0464] From 11:15:31 to 11:15:31.5 (communication write delay T_comm≤0.5 seconds): Ω_intra(t) communication frame is sent and written to the three-layer controller constraint register.

[0465] Considering the offset of T_s / 2 = 0.5 seconds introduced by the misalignment of sampling time and the scheduling jitter margin T_jitter = 0.5 seconds, the equivalent total delay from the significant increase of D_M_bar(t) (approximately 11:15:10) to the completion of the constraint register update (approximately 11:15:31.5) is 22.5 seconds, which is exactly equal to the theoretical upper limit τ_total_max = 22.5 seconds given by Equation 833 in S8.11, thus verifying the rationality of the delay design of the loop response link in this invention.

[0466] Since κ(t) = 0.7365, satisfying 0.5 ≤ κ(t) < 0.8, the loop execution domain D_intra(t) is switched to domain B.

[0467] From 11:15:32: the reactor bottom controller 6 reads P_ch_max_intra(t) from the constraint register, and updates the upper limit of the saturation limit to 0.7365 times the original value. The gain is adjusted according to equations 822 to 824 with α_eq_intra(t) = 0.7365 × α_filt (87), and the rate of change is limited by equations 825 to 827. The feedforward control term is calculated according to equation 828, replacing α_filt (87) with α_eq_intra(t).

[0468] From 11:15:32 to approximately 11:18: the local sintering mutation continued to develop. D_M_bar(t) continued to rise after the activation of the secondary response, reaching a peak value of D_M_bar_peak = 6.9 during this anomaly at approximately 11:18, and then gradually decreased. A single computing node continuously updates κ(t) according to the control sampling period T_s = 1 second: when D_M_bar(t) = 6.9, κ(t) ≈ 0.5727 is calculated according to Equation 813, and the intra-loop operating domain D_intra(t) is still located in the domain B interval [0.5, 0.8); the rate of change of κ(t) is constrained by the gain change rate limit specified in S8.8, and smoothly transitions to the new value.

[0469] From 11:18 to 11:20: D_M_bar(t) gradually decreased from its peak, and κ(t) rose accordingly.

[0470] Phase 4 (3 PM to 6 PM, charging ends and switching begins):

[0471] At 3 PM, the photovoltaic output dropped below the demand of the heating network, and the dispatch plan controller 4 triggered a switch to the energy charging and releasing mode. Step S6 performed an S-shaped transition according to equations 601 to 604. During the transition, step S8 continued to run, and no abnormalities were triggered during this stage.

[0472] Phase 5 (6 PM to midnight the following day, nighttime energy release phase):

[0473] The energy release phase begins, and this cycle is completed normally.

[0474] At the end of the loop, step S2 performs inter-loop identification: Since this loop previously triggered an in-loop correction and D_M_bar_peak is greater than 1.5 times D_th2, the filter uses penalized weights.

[0475] λ_p=min(1,λ×D_th2 / D_M_bar_peak)

[0476] =min(1,(1 / 3)×4.03 / 6.9)

[0477] ≈min(1,0.1947)

[0478] =0.1947

[0479] Filtering results:

[0480] α_filt(88)=0.1947×α_raw(88)+0.8053×α_filt(87)

[0481] The identification value for this abnormal loop is assigned a low confidence level. After generating Ω(88) and writing it to the constraint register, the correction state within the loop is completely reset, and the three-layer controller enters the next loop.

[0482] Example 4: Parallel expansion of multiple reactors

[0483] When the system contains M parallel fixed-bed reactors (M≥2), the extension of the present invention is as follows.

[0484] Each reactor is independently equipped with a reactor bottom controller 6 and an independently equipped bed temperature field acquisition device 12. A single computing node 7 is the sole computing and distribution entity for all M reactors.

[0485] Step S2, the inter-cycle identification, is performed independently for each reactor, generating M sets of independent effective activity coefficients α_filt_q(n), q=1,2,...,M.

[0486] In the formula: q is the reactor number, q=1,2,...,M; α_filt_q(n) is the effective activity coefficient of the q-th reactor after the nth cycle after filtering, which is dimensionless.

[0487] The components within the unified operating permission domain are aggregated according to the reactor's physical characteristics. The system-level safe charging power limit is aggregated as a sum of all reactors.

[0488] P_ch_max_sys(n)=Σ(q=1 to M)P_ch_max_q(n) (Equation 835)

[0489] In the formula: P_ch_max_sys(n) is the upper limit of the system-level safe charging power, in watts; P_ch_max_q(n) is the upper limit of the safe charging power of the q-th reactor after the nth cycle, in watts.

[0490] The upper limit of the system-level state of charge is the minimum of the most conservative values ​​among all reactors:

[0491] SOC_max_sys(n) = min(q=1 to M)SOC_max_q(n) (Equation 836)

[0492] In the formula: SOC_max_sys(n) is the upper limit of the system-level state of charge, which is dimensionless; SOC_max_q(n) is the upper limit of the state of charge of the q-th reactor after the nth cycle, which is dimensionless; min(·) represents taking the minimum value among all values.

[0493] The in-loop anomaly detection in step S8 is performed independently for each reactor. When the q-th reactor triggers a secondary response, a single compute node only modifies the upper limit of the charging power P_ch_max_q_intra(t) for that reactor, and the system-level upper limit of charging power is updated accordingly:

[0494] P_ch_max_sys_intra(t)=Σ(q=1 to M)P_ch_max_q_eff(t) (Equation 837)

[0495] In the formula: P_ch_max_sys_intra(t) is the upper limit of the system-level safe charging power after the loop correction, in watts; P_ch_max_q_eff(t) is the upper limit of the current effective charging power of the q-th reactor, in watts. If the q-th reactor is in the loop correction activation state, then P_ch_max_q_eff(t) = P_ch_max_q_intra(t); otherwise, P_ch_max_q_eff(t) = P_ch_max_q(n).

[0496] The system-level unified operating license domain Ω_sys_intra(t) is written to the constraint registers of the scheduling plan controller 4 and the power scheduling controller 5. The constraint registers of each reactor's underlying controller 6 are written to their respective individual machine license domains. This extension method maintains the core architecture of the present invention, which is "single computing node + write protection constraint register + three-layer shared license domain", and only adds parallel aggregation rules based on reactor number at the license domain aggregation layer. Therefore, the two core technical effects of the present invention, namely, inter-layer constraint consistency and no blind zone within the cycle, also hold true in multi-reactor parallel systems.

Claims

1. A control optimization method for a photovoltaic-coupled thermochemical energy storage device, the method being executed based on a system comprising: A calcium oxide / calcium hydroxide fixed-bed reactor (1), an electric heater (2), a heat exchanger (3), a scheduling plan controller (4), a power scheduling controller (5), a reactor bottom-level controller (6), a single computing node (7), an industrial real-time communication interface (8), constraint registers (9) respectively set in the scheduling plan controller, the power scheduling controller and the reactor bottom-level controller, local default memories (10) respectively connected to the three controllers, a write protection module (11) for writing protection of the constraint registers (9), and a bed temperature field acquisition device (12) uniformly arranged along the axial direction of the fixed-bed reactor and having no less than 5 temperature measurement points, characterized in that it includes the following steps: After each complete energy storage and release cycle, the single computing node calculates the corrected heat release for this cycle based on the readings of the heat exchanger inlet and outlet temperature sensors and the readings of the heat exchange medium mass flow sensor, and obtains the original value of the effective activity coefficient α_raw(n) for this cycle. After exponential weighted average filtering, α_filt(n) is obtained. The single computing node centrally generates a unified operating permission domain Ω(n) based on α_filt(n). The Ω(n) includes at least the upper limit of state of charge SOC_max(n), the upper limit of permission domain energy E_perm(n), the upper limit of safe charging power P_ch_max(n), and the temperature feasible domain Ω_T(n) calculated based on α_filt(n), where P_ch_max(n) is a monotonically increasing function of α_filt(n). The single computing node synchronously writes Ω(n) into the constraint registers (9) of the scheduling plan controller, power scheduling controller and reactor bottom layer controller through the write protection module (11) via the industrial real-time communication interface (8); the three-layer controllers only have read permissions for the constraint registers, and use Ω(n) in the constraint registers as their only constraint source to perform scheduling plan optimization, rolling power optimization and bottom layer closed-loop control respectively; During each energy storage and release cycle, the single computing node constructs a spatial distribution feature vector φ(t) of the bed temperature field based on real-time readings from no less than 5 temperature measurement points of the bed temperature field acquisition device, and calculates the Mahalanobis distance D_M_bar(t) between φ(t) and the mean of the sliding baseline and the baseline covariance matrix. When D_M_bar(t) exceeds the protective correction trigger threshold D_th2 for N_conf consecutive sampling periods, the single computing node generates an in-cycle activity reduction factor κ(t) based on D_M_bar(t) according to a preset mapping relationship, where 0 < κ(t) ≤ 1, and uses κ(t) to tighten each component of Ω(n) to obtain the in-cycle correction permission domain Ω_intra(t), which is synchronously written to the constraint register of the three-layer controller through the write protection module, overwriting the previous Ω(n) or the previous Ω_intra. In the in-cycle correction permission domain Ω_intra(t), P_ch_max_intra(t) = κ(t) × P_ch_max(n).

2. The method according to claim 1, characterized in that, The write protection module (11) is implemented as at least one of hardware write protection circuit, software permission control or controller firmware lock, so that the constraint register (9) can only be written by the single computing node (7) through the industrial real-time communication interface (8); the three-layer controller does not participate in the identification of effective activity coefficients, does not participate in the judgment of abnormality in the loop, and cannot independently modify the constraint components in the constraint register.

3. The method according to claim 1, characterized in that, The in-cycle activity reduction factor κ(t) is calculated using the following formula: κ(t)=min(1,max(κ_min,1−[(D_M_bar(t)−D_th2) / (D_sat−D_th2)]×(1−κ_min))) Where κ_min is the minimum reduction factor, with a value ranging from 0.3 to 0.5; D_sat is the saturated Mahalanobis distance, which is twice D_th2; when D_M_bar(t)≤D_th2, κ(t)=1, and the in-loop correction is not triggered; when D_M_bar(t)≥D_sat, κ(t)=κ_min.

4. The method according to claim 1, characterized in that, The coordination of the intra-loop correction and inter-loop identification is performed according to the following rules: (a) When the loop ends normally and inter-loop identification is completed, and Ω(n+1) is written to the constraint register by the write protection module (11), the intra-loop correction state is automatically reset, κ is reset to 1, and Ω_intra is invalidated; (b) If the number of intra-loop corrections triggered within the validity period of Ω(n) exceeds the threshold N_trig, then an additional penalty filter is applied to α_filt(n+1) in the next inter-loop identification; (c) If the current loop fails to end normally, the current inter-loop identification update is skipped, the loop count n is not incremented, Ω(n) remains unchanged, and item (a) is not executed. The original κ(t) and Ω_intra continue to be valid in the new loop until the next normal inter-loop identification is completed.

5. The method according to claim 1, characterized in that, When the in-cycle correction is active, the reactor bottom controller (6) replaces α_filt(n) with the in-cycle equivalent activity coefficient α_eq_intra(t)=κ(t)×α_filt(n) and adjusts the closed-loop control gain in real time according to the following relationship: K_p_intra(t)=K_p0 / α_eq_intra(t); K_i_intra(t)=K_i0 / sqrt(α_eq_intra(t)); K_d_intra(t)=K_d0×sqrt(α_eq_intra(t)); Furthermore, upper limits K_dot_p_max, K_dot_i_max, and K_dot_d_max are applied to the time rates of change of K_p_intra, K_i_intra, and K_d_intra, respectively; the reaction endothermic compensation rate in the feedforward control term is calculated by replacing α_filt(n) with α_eq_intra(t).

6. The method according to claim 1, characterized in that, The safe charging power limit P_ch_max(n) is calculated using the following formula: P_ch_max(n)=P_ch_rated×[γ+(1-γ)×α_filt(n)] Where P_ch_rated is the rated charging power, and γ is the non-reaction power ratio, defined as the ratio of the sum of the bed sensible heat rise power and the reactor surface heat loss power under rated operating conditions to the rated charging power, with γ ranging from 0.15 to 0.30; this formula guarantees that P_ch_rated is given when α_filt(n)=1, and γ×P_ch_rated is given when α_filt(n)→0.

7. The method according to claim 1, characterized in that, The single computing node (7) divides the material operating state into three levels based on α_filt(n): normal operating domain A, derating operating domain B, and restricted operating domain C, and switches with hysteresis according to the following rules: switching from A to B is based on α_filt(n) < α_A, and switching back from B to A must satisfy α_filt(n) ≥ α_A + h; switching from B to C is based on α_filt(n) < α_B, and switching back from C to B must satisfy α_filt(n) ≥ α_B + h; where α_A > α_B, and h is the hysteresis bandwidth; Ω corresponding to different operating domains (n) The correction is applied to P_ch_max(n) and E_perm(n) = E_conf × SOC_max(n) with different reduction coefficients β_d and ρ(n) respectively, and ρ(n) × SOC_max(n) is used as the upper bound of SOC operation. In domain A, β_d = 1 and ρ = 1; in domain B, β_d ∈ [0.7, 0.9] and ρ ∈ [0.85, 0.95]; in domain C, β_d ∈ [0.4, 0.6] and ρ ∈ [0.6, 0.8]; the hysteresis band h = 0.02 is set for inter-domain switching to avoid repeated jumps at the boundary.

8. The method according to claim 1, characterized in that, It also includes a three-level response process for communication link failures of the industrial real-time communication interface (8): when a single frame is lost, the single computing node retransmits it in the next communication cycle, and the constraint register maintains the previous valid Ω value unchanged; When more than a preset number of frames are lost consecutively, an alarm is triggered and the last valid Ω value is maintained; when the communication link is interrupted for more than a preset duration, the three-layer controller automatically switches to the conservative and safe default value in the local default memory (10) connected to it. The conservative and safe default value is pre-written through the offline physical interface and cannot be modified through the online interface during operation.

9. A control and optimization system for a photovoltaic-coupled thermochemical energy storage device, the system being used to execute the control and optimization method for a photovoltaic-coupled thermochemical energy storage device as described in claim 1, characterized in that, include: A calcium oxide / calcium hydroxide fixed-bed reactor (1) is connected to an electric heater (2) and a heat exchanger (3), wherein the electric heater is used to input heat into the fixed-bed reactor during the charging phase and the heat exchanger is used to output the heat of reaction to the heat network during the energy release phase; A bed temperature field acquisition device (12) is arranged with no less than 5 temperature measurement points along the axial direction of the fixed bed reactor (1); The scheduling plan controller (4), power scheduling controller (5) and reactor bottom controller (6) each have built-in constraint registers (9); The single computing node (7) is the sole entity for identifying the effective activity coefficient within the system, generating the unified operating permit domain, and determining the anomaly within the cycle. On the one hand, the single computing node is connected to the inlet and outlet temperature sensors and the heat exchange medium mass flow sensor of the heat exchanger (3) to generate the unified operating permit domain Ω(n) between cycles after each complete energy storage and release cycle. On the other hand, it is connected to the bed temperature field acquisition device (12) to generate the corrected permit domain Ω_intra(t) within the cycle based on the Mahalanobis distance D_M_bar(t) during the cycle. An industrial real-time communication interface (8) is connected between the single computing node (7) and the constraint register (9) of the three-layer controller, and is used to synchronously write Ω(n) and Ω_intra(t) into the constraint register of the three-layer controller. Write protection module (11) is set between the industrial real-time communication interface (8) and the constraint register (9) or on the write channel of the constraint register. It is configured to allow only the single computing node (7) to write to the constraint register through the industrial real-time communication interface (8). The three-layer controller has only read permission for the constraint register and uses the permission field in the constraint register as its only source of constraint. The local default memory (10) is connected to the three-layer controller and stores conservative and safe default values ​​for degraded operation when the communication interruption of the industrial real-time communication interface (8) exceeds the preset time.

10. The application of the method according to any one of claims 1 to 8 in a photovoltaic-coupled thermochemical energy storage grid-connected scenario, characterized in that, With the photovoltaic available power forecast and the heating network heat demand forecast as input, the scheduling plan controller (4), the power scheduling controller (5) and the reactor bottom layer controller (6) respectively use the unified operating permission domain Ω(n) or the cycle correction permission domain Ω_intra(t) in the constraint register (9) as the inviolable boundary to perform the day-ahead energy storage and release plan preparation, day-ahead rolling power optimization and bottom layer closed-loop control; when the photovoltaic output fluctuation or the heating network load change causes the bed temperature field spatial distribution to be abnormal, the single computing node (7) triggers the cycle permission domain protective correction without waiting for the current cycle to end, so that the constraint boundary of the three-layer controller is consistent with the actual reachable state of the material in the second-level time resolution within the cycle. The total response delay τ_total from the occurrence of local anomaly in the bed temperature field to the completion of the three-layer synchronization of Ω_intra(t) is no more than 25 seconds, and much less than a complete charge / release cycle.

Images