Energy storage spot market coordinated clearing method based on time delay characteristics

Abstract

The application discloses a kind of energy storage spot market collaborative clearing methods based on time delay characteristics, it is related to power market clearing technical field, including the following steps: obtaining the physical operation data of energy storage subject, extracting energy storage subject time delay response characteristics and based on battery temperature evolution and thermal accumulation characteristics, define the climbing limited interval;Multi-dimensional state space including time delay response characteristics is constructed, and based on the climbing limited interval, the safe reachable domain of energy storage subject in multi-dimensional state space is calculated;Dimensionality reduction decomposition and time evolution integration are carried out to the safe reachable domain to obtain non-convex evolution boundary, the evolution boundary is divided into local convex region set and is carried out logic reconstruction, and the segmented time delay climbing constraint is obtained;Based on segmented time delay climbing constraint, pre-established collaborative clearing objective function is solved, and the clearing capacity of the energy storage subject is obtained.The application is used to solve the problem that the prior art does not consider physical time lag and dynamic recession, leading to instruction response mismatch and adjustment resource realization limit.

Description

Technical Field

[0001] This invention relates to the field of electricity market clearing technology, and more specifically, to a collaborative clearing method for the energy storage spot market based on time delay characteristics. Background Technology

[0002] Existing energy storage spot market collaborative clearing schemes mostly focus on multi-market revenue games or macro capacity allocation to meet the system regulation needs under high proportion of new energy access and the economic goals of energy storage entities.

[0003] For example, the invention patent with announcement number CN121073112A discloses a method for energy storage to participate in the spot electricity energy-flexibility market clearing. It uses principal component analysis to construct a multi-dimensional joint fluctuation domain to quantify flexibility demand and establishes a two-level optimization decision model to achieve a balance between energy storage arbitrage profits and system dispatch costs. Another example is the invention patent with announcement number CN121395380A, which discloses a ramp-frequency regulation coordination optimization method considering the uncertainty of energy storage capacity. This method uses chance constraints to characterize the probability distribution boundary of energy storage capacity and linearly reconstructs the ramp-up capability of the units, thereby achieving resource allocation for ramp-up and frequency regulation services under a unified power constraint.

[0004] However, the aforementioned existing technologies still have the following technical problems in practical applications: 1. Energy storage devices often lag in responding to commands, resulting in a disconnect between actual output and dispatch signals during high-frequency dispatch, making it difficult to meet the grid's real-time adjustment needs; 2. In scenarios involving continuous high power or frequent frequency adjustments, the actual ramp-up capability of energy storage devices often declines nonlinearly with the operation process, leading to the risk of unexpected interruption of output during the planned execution, which in turn makes it difficult to fully realize the awarded capacity. 3. If complex physical constraints are directly introduced into the existing model to ensure the feasibility of the plan, it is easy to cause a sharp increase in the clearing calculation time or local deadlock, which makes it difficult to meet the timeliness requirements of high-frequency rolling in the spot market. Summary of the Invention

[0005] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a collaborative clearing method for the energy storage spot market based on time delay characteristics. By extracting the time delay and heat accumulation characteristics of energy storage and reconstructing the non-convex reachable domain into a segmented ramp constraint for clearing optimization solution, the present invention addresses the problems of command response mismatch and limited realization of regulation resources caused by the failure to consider physical time delay and dynamic decay in the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: A collaborative clearing method for the energy storage spot market based on time delay characteristics includes the following steps: The physical operation data of the energy storage unit is obtained, the time delay response characteristics of the energy storage unit are extracted, and the ramp-limited range is defined based on the cell temperature evolution and heat accumulation characteristics. A multidimensional state space including the aforementioned time delay response characteristics is constructed, and the safe reachability domain of the energy storage entity within the multidimensional state space is calculated based on the aforementioned ramp-limited interval. The safe reachable domain is decomposed by dimensionality reduction and time evolution integral to obtain a non-convex evolution boundary. The evolution boundary is divided into a set of locally convex regions and logically reconstructed to obtain the segmented time delay ramping constraint of the energy storage subject. The pre-established collaborative clearing objective function is solved based on the segmented time delay ramping constraint to obtain the clearing capacity of the energy storage entity.

[0007] In a preferred embodiment, the step of extracting the time delay response characteristics of the energy storage entity and defining the ramp-up restricted area based on the cell temperature evolution and heat accumulation characteristics includes: extracting time delay response characteristics based on the time response deviation between the historical scheduling commands of the energy storage entity and the corresponding historical actual output data; acquiring environmental thermal disturbance data and physical operation data and combining them with the time delay response characteristics to deduce the cell temperature evolution trajectory corresponding to the power adjustment of the energy storage entity; and extracting the power adjustment range that makes the evolution trajectory meet the preset heat accumulation constraint as the ramp-up restricted area.

[0008] In a preferred embodiment, the construction of the multidimensional state space including the delay response features includes: constructing a multidimensional state space including the delay response features and the physical operation data, and setting the target boundary set of the multidimensional state space at the deadline of the preset scheduling period according to the preset heat accumulation constraint.

[0009] In a preferred embodiment, calculating the safe reachability domain of the energy storage entity within the multidimensional state space includes: within a preset scheduling period, based on the ramp-limited interval and the environmental thermal disturbance data, performing extreme value calculations on the evolution paths from each state point in the multidimensional state space to the target boundary set to obtain the optimal evolution cost corresponding to each state point; and taking the set of state points whose optimal evolution cost satisfies a preset safety threshold as the safe reachability domain.

[0010] In a preferred embodiment, the dimensionality reduction decomposition of the secure reachable domain includes: performing grid discretization on the multidimensional state space to obtain a numerical array of the secure reachable domain at each grid node; performing tensor decomposition on the numerical array based on a preset rank evaluation truncation criterion to extract local tensor cores characterizing the evolution features of the secure reachable domain and the corresponding mapping basis.

[0011] In a preferred embodiment, obtaining the non-convex evolution boundary includes: taking the scheduling cycle deadline as the starting point and the current time as the ending point, performing reverse reachability propagation on the local tensor core to obtain tensor evolution features characterizing boundary changes; extracting isosurface feature points from the tensor evolution features that satisfy preset boundary conditions, and using the mapping basis to reconstruct and restore the feature points to the multidimensional state space to obtain the non-convex evolution boundary.

[0012] In a preferred embodiment, obtaining the segmented time delay ramping constraint of the energy storage entity includes: transforming the evolution boundary into a set of local convex regions composed of multicellular structures through geometric segmentation; and constructing a set of inequalities characterizing the switching logic between convex regions based on the boundary geometric parameters of the set of local convex regions, as the segmented time delay ramping constraint.

[0013] In a preferred embodiment, obtaining the clearing capacity of the energy storage entity includes: constructing a thermal accumulation sensitivity matrix and extracting features based on the topological relationship and local thermal disturbance data of each energy storage unit within the energy storage entity to obtain the thermal difference adjustment parameters of the energy storage entity; using the thermal difference adjustment parameters to correct the piecewise time delay ramping constraint, and optimizing the cooperative clearing objective function based on the corrected ramping constraint to obtain the clearing capacity.

[0014] In a preferred embodiment, the collaborative clearing objective function includes at least an energy storage lifetime degradation index, operating cost, and dispatch deviation penalty index; the clearing capacity includes the energy market clearing power of the energy storage entity and the ancillary service market clearing reserve capacity.

[0015] In a preferred embodiment, after obtaining the cleared capacity of the energy storage entity, the method further includes: converting the cleared capacity into an active power signal sequence, and performing feedforward correction on the active power signal sequence based on the signal transmission delay in the delay response characteristics, to generate a duty cycle command to drive the energy storage entity to perform power regulation.

[0016] The technical effects and advantages of the collaborative clearing method for the energy storage spot market based on time delay characteristics proposed in this invention are as follows: 1. This invention extracts the time delay response characteristics of the energy storage entity and the ramp-up limited range considering the thermal accumulation characteristics of the battery cells, and constructs a multi-dimensional state space to calculate the safe reachable domain. This effectively avoids the mid-term output disconnection caused by equipment response lag or thermal protection triggering in actual operation, improves the matching degree between control commands and physical responses, and ensures the fulfillment rate and safety of the awarded capacity in actual grid dispatch execution from the source.

[0017] 2. This invention extracts non-convex evolutionary boundaries by performing dimensionality reduction decomposition and time evolution integral on the safe reachable domain, and divides it into a set of locally convex regions for logical reconstruction to generate piecewise time delay ramping constraints. This realizes the transformation of the complex nonlinear dynamic decay boundary of energy storage equipment into a set of convex constraints that are easy for the solver to process, reducing the computational complexity and deadlock risk of joint solution, and meeting the timeliness requirements of high-frequency rolling clearing engine calculation in the spot market environment. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of a collaborative clearing method for the energy storage spot market based on time delay characteristics, provided as an embodiment of the present invention.

[0019] Figure 2 A three-dimensional simulation diagram illustrating the time-delayed contraction of the safe reachable domain boundary, as provided in an embodiment of the present invention.

[0020] Figure 3 This is a schematic diagram showing the simulation comparison of the closed-loop control response waveform provided in an embodiment of the present invention. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0022] Example 1, Figure 1 This invention presents a collaborative clearing method for the energy storage spot market based on time delay characteristics, comprising the following steps: S1. Obtain the physical operation data of the energy storage unit, extract the time delay response characteristics of the energy storage unit, and define the ramp-limited range based on the cell temperature evolution and heat accumulation characteristics, including: S101. Based on the time response deviation between the historical dispatch commands of the energy storage entity and the corresponding historical actual output data, extract the time delay response features, as follows: The system uses a remote monitoring and control terminal of the energy storage converter to collect active power command signals sent from the upper layer and actual output signals fed back by the battery system in real time with a high-frequency sampling period (e.g., 10ms). During the acquisition process, the system uses a high-precision clock synchronization protocol (e.g., IEEE 1588 PTP) to add a uniform absolute timestamp to each frame of acquired signal. Then, in order to convert real-time data into feature samples that can be analyzed by algorithms, a first-in-first-out online sliding data window is built inside the local controller. The recommended length of this window is set to the past 5 minutes, and the active power command signals and actual output signals with timestamps are periodically (e.g., once per minute) continuously pushed into the online sliding data window. The system performs time-series alignment and pairing of the active power command signals and actual output signals according to the timestamps to obtain the "command-output" data column. An online optimization algorithm based on the least squares criterion of the forgetting factor is used to extract the delay response features of the current scheduling cycle from the data series. The online optimization algorithm performs asynchronous shifting of historical scheduling command data on a discrete data sequence within a preset physical allowable delay range (recommended to be 0ms to 1000ms, and equivalently mapped to 0 to 100 discrete sampling steps according to the sampling period) to mathematically simulate a lag process. Then, the weighted sum of squared errors between the shifted command waveform and the actual output waveform is calculated. The optimal shift step number that minimizes the sum of squared errors is multiplied by the sampling period to obtain the current total delay, i.e., the delay response features. The formula for extracting the delay response features by the online optimization algorithm is as follows: (1) (2) in, To find the optimal number of time shift steps identified; The total number of data samples within the online sliding data window is calculated from the time length of the sliding window and the sampling period. For example, if the window length is 5 minutes and the sampling period is 10ms, then the total number of data samples is 30,000. These are discrete sampling time points; The preset forgetting factor is recommended to be between 0.95 and 0.99. For the first Historical actual power output data at each sampling time, in MW; To shift historical scheduling instruction data backward on the timeline Data after each time step, in MW; This represents the current total latency, i.e., the latency response characteristic, in milliseconds (ms). The high-frequency sampling period is (e.g., 10ms).

[0023] S102. Acquire environmental thermal disturbance data and physical operation data, and combine them with the time delay response characteristics to deduce the cell temperature evolution trajectory corresponding to the energy storage main body's adjustment power, as follows: The system collects real-time physical operating data of the main energy storage battery through the Battery Management System (BMS). This physical operating data specifically includes the current cell initial temperature, current state of charge (SOC), and current battery pack terminal voltage. Simultaneously, a temperature sensing sensor network deployed inside and outside the energy storage container collects real-time external ambient temperature as environmental thermal disturbance data. Then, a discrete iterative approximation method is used to extract the dynamic internal resistance parameters of the cells. Specifically, at the time of the calculation... When the thermodynamic state is known, the previous integration iteration step (i.e., time step 1) is used. The cell temperature has been deduced and cached. The system uses the current State of Charge (SOC) as the addressing input parameter to query a pre-calibrated two-dimensional MAP table stored in the controller's memory. This MAP table is composed of a two-dimensional data index grid consisting of discretized temperature calibration points (e.g., the range from -20℃ to 60℃, with 5℃ intervals) and SOC calibration points (e.g., the range from 0% to 100%, with 5% intervals). The matrix elements at the grid intersections store reference internal resistance values ​​obtained by fitting the hybrid pulse power characteristics of a typical battery cell (e.g., a 280Ah lithium iron phosphate cell) through offline experimental measurements. The system addresses and locates four neighboring calibration grid points surrounding the current input parameter in the MAP table. Then, it performs two one-dimensional linear interpolations along the SOC data dimension to obtain two temperature transition internal resistance values. Finally, it performs one linear interpolation along the temperature data dimension on the two temperature transition internal resistance values ​​to obtain the dynamic internal resistance parameter. This dynamic internal resistance parameter is used to calculate the Joule heating power term in the subsequent differential equation. Combining the initial cell temperature, current state of charge, battery pack terminal voltage, external ambient temperature, and the time delay response characteristics, a discrete-time numerical integration algorithm (e.g., fourth-order Runge-Kutta method, with a recommended integration iteration step size of 1 second) is used to solve the pre-constructed thermodynamic state differential equation. By discretely integrating the temperature change rate term in the thermodynamic state differential equation according to the recommended integration iteration step size, the discrete points of the cell evolution temperature are obtained and connected in time sequence to obtain the cell temperature evolution trajectory. The thermodynamic state differential equation is as follows: (3) (4) in, The equivalent heat capacity of the battery cell is pre-calibrated by the battery's factory thermal characteristic test, and the unit is J / ℃; For a moment The cell temperature, i.e., discrete points on the cell temperature evolution trajectory, and Equal to the initial temperature of the battery cell, in °C; The actual cell response current, derived from the adjusted power and taking into account the time delay response characteristics, is expressed in amperes (A). The dynamic internal resistance parameter; This refers to environmental thermal disturbance data, i.e., external ambient temperature, in °C. The equivalent convection and conduction thermal resistance constant from the cell to the external environment is pre-calibrated. In conventional liquid-cooled energy storage systems, this constant is recommended to be calibrated to 0.05℃ / W. For the regulation power of the energy storage unit, The time delay response characteristic, To simulate a constant reference voltage within the time window, since feedforward simulation cannot obtain future measured voltages, the current battery pack terminal voltage from the physical operation data is assigned to... .

[0024] S103. Extract the adjustable power range that makes the evolution trajectory satisfy the preset heat accumulation constraint, as the climbing-limited interval, as follows: Within a preset scheduling period (e.g., the next 15 minutes), a boundary search algorithm based on binary search (with a convergence tolerance of 0.01MW and a maximum number of iterations of 50) is used to solve for the power boundary extreme values ​​of the bidirectional charging and discharging scheduling scenario. The preset heat accumulation constraint includes that the cell evolution temperature does not exceed the preset maximum safe temperature threshold (e.g., 55℃) and the temperature rise rate does not exceed the preset maximum temperature rise rate threshold (e.g., 3℃ / min). Taking the solution of the maximum effective charging power boundary as an example, the system sets the known power value that satisfies the constraint (such as 0MW) as the lower bound of the search and the rated charging power of the energy storage body as the upper bound of the search. In each iteration, the system takes the median of the current upper and lower bounds as the candidate adjustable power input to the thermodynamic state differential equation in step S102 for trajectory feedforward deduction. If the deduction finds that the cell temperature evolution trajectory of the candidate adjustable power in the deduction time window does not satisfy the preset heat accumulation constraint, the candidate adjustable power is determined to be infeasible, and the upper bound of the search is updated to the median. Conversely, if the deduced cell temperature evolution trajectory satisfies the preset heat accumulation constraint, the candidate adjustable power is determined to be feasible, and the lower bound of the search is updated to the median. The system continuously reduces the search interval until the difference between the current upper and lower bounds is less than the convergence tolerance or the maximum number of iterations is reached. At this time, the output lower bound of the search is the maximum effective charging power boundary. For solving the maximum effective discharge power boundary, the same algorithm logic as that used for solving the maximum effective charging power boundary is employed, with the only difference being the mathematical setting of the initial search interval: for the discharge scenario, the known 0MW that satisfies the constraints is set as the upper search bound, and the rated discharge power of the energy storage entity (which is algebraically negative, for example, -10MW) is set as the lower search bound; in each iteration, the median is taken for testing. If it is feasible, the upper search bound is updated to the median; if it is not feasible, the lower search bound is updated to the median. The final converged upper search bound (i.e., the feasible negative power with the largest absolute value) is the maximum effective discharge power boundary. The feasible region closed set enclosed by the maximum effective charge and discharge power boundary is the adjustable power range that makes the evolution trajectory satisfy the preset heat accumulation constraint; the adjustable power range is taken as the ramp-limited interval, and the mathematical set expression of the ramp-limited interval is as follows: (5) in, The set of charge and discharge regulation powers that are allowed to be safely executed, i.e., the ramp-limited range, is measured in MW; The candidate adjustable power is expressed in MW. The time derived from the candidate adjustment power Discrete points on the temperature evolution trajectory of the battery cell, in °C; This refers to the highest safe temperature threshold in the heat accumulation constraint; The maximum temperature rise rate threshold is defined in the heat accumulation constraint.

[0025] This step achieves numerical representation of the dynamic physical constraints of the energy storage entity by online identification of time delay response characteristics and feedforward extrapolation of the thermodynamic evolution trajectory of the battery cell. It eliminates the time series and thermodynamic calculation deviations caused by traditional static models, provides basic state parameters for subsequent reconfiguration of the reachable domain, and thus ensures the physical operation safety of the energy storage entity and the execution reliability of the underlying control commands during continuous scheduling.

[0026] S2. Construct a multi-dimensional state space including the aforementioned time delay response characteristics, and calculate the safe reachability domain of the energy storage entity within the multi-dimensional state space based on the aforementioned ramp-limited interval, including: S201. Construct a multi-dimensional state space including the time delay response characteristics and the physical operation data, and based on the preset heat accumulation constraint, set the target boundary set of the multi-dimensional state space at the deadline of the preset scheduling period, as follows: First, using the initial cell temperature and current state of charge from the physical operation data as the initial coordinates for time-differential evolution, a state vector characterizing the dynamic evolution of the energy storage entity is constructed. The multidimensional state space is a three-dimensional Euclidean space with the state vector as the basic coordinate axis, used to describe the dynamic transfer trajectory of the energy storage entity after the control command is issued. The state vector is defined as containing the dynamic state of charge evolving with time. Dynamic cell temperature And the three-dimensional column vector representing the time delay response features as augmented state dimensions, i.e. Among them, dynamic state of charge The differential evolution follows the Anton-Jackson integral theorem, and its differential equation is: ,in The actual response current of the battery cell is specified (the charging current is taken as a positive value, and the discharging current is taken as a negative value). Rated ampere-hour capacity of the battery pack, in Ah; dynamic cell temperature. The differential evolution follows the thermodynamic state differential equation described in step S102; it should be noted that the time delay response characteristics are regarded as constant system parameters within a single scheduling cycle, and their evolution derivative is always 0; Next, the highest safe temperature threshold in the heat accumulation constraint is extracted, and combined with the battery's inherent physical charge / discharge cutoff range (e.g., a preset state of charge limit range of 10% to 90%), an implicit objective function is established for statically and geometrically partitioning the multidimensional state space; the calculation formula of the implicit objective function is as follows: (6) in, The value of the implicit objective function; This refers to the highest safe temperature threshold. and These are the lower and upper limits of the preset state of charge limit range, respectively. It should be noted that since the maximum temperature rise rate threshold in the heat accumulation constraint is essentially an input coupling constraint on the regulating power, and has been transformed into a range restriction on the regulating power in step S103, the implicit objective function only needs to include independent absolute state threshold constraints. Furthermore, at the deadline of the preset scheduling period (e.g., the next 15 minutes), any coordinate point in the multidimensional state space that causes the dynamic cell temperature to exceed the maximum safe temperature threshold or the dynamic state of charge to exceed the physical charge-discharge cutoff interval will have a negative implicit objective function value. The system measures these regions with negative calculated values ​​as infeasible state subspaces; conversely, it measures regions with positive calculated values ​​as feasible state subspaces. The system further extracts the zero level set of the implicit objective function (i.e., the isosurface with a function value of zero) as the boundary line enclosing the infeasible state subspace, and directly sets the zero level set as the target boundary set.

[0027] S202. Within a preset scheduling period, based on the climbing-limited interval and the environmental thermal disturbance data, the extreme value calculation is performed on the evolution path from each state point in the multidimensional state space to the target boundary set to obtain the optimal evolution cost corresponding to each state point. The set of state points whose optimal evolution cost satisfies a preset safety threshold is taken as the safe reachable domain, as follows: By introducing Hamilton-Jacobi reachability analysis theory, and solving the partial differential equations describing the evolution of the system's safety boundary, the resilience of each state point in the multidimensional state space under extreme disturbances, i.e., the optimal evolution cost, is quantified. The Hamilton-Jacobi partial differential equations and the formulas for calculating the system state space evolution rate function in the equations are as follows: (7) (8) in, To reverse the evolution to time point At that time, the state vector in the multidimensional state space The optimal evolutionary cost function; For reverse evolution time variables (defined) ,in The preset scheduling period end time, (This is positive physical time), measured in seconds; The candidate adjustable power of the energy storage unit is expressed in MW. The set of adjustable powers that can be executed at the current integration time is the climbing-limited interval; The environmental thermal disturbance data variable is the specific adverse temperature selected during the game process, and the unit is °C. This is the set of extreme values ​​of environmental temperature disturbances identified from historical meteorological data (e.g., deviations from historical maximum temperatures). (range in °C), unit is °C; Let be the gradient vector of the three-dimensional partial derivative of the optimal evolutionary cost function with respect to the state vector; is the system state-space evolution rate function, and its first to third rows correspond to the dynamic state of charge, dynamic cell temperature, and the differential rate of change of time-delay response characteristics, respectively. Mathematically, this is represented as the inner product (dot product) of the spatial partial derivative gradient vector and the state change rate vector. This inner product, combined with the minimax optimization operation on the candidate regulation power and the environmental thermal disturbance data, constitutes the numerical Hamiltonian operator of the Hamilton-Jacobi partial differential equation.

[0028] The Hamilton-Jacobi partial differential equation is solved using the level set method, with the following steps: (1) In the multidimensional state space, the dynamic charge state axis (e.g., 10% to 90%) and the dynamic cell temperature axis (e.g., 20°C to 60°C) are set at preset steps (e.g. SOC and 0.2℃) generate a discrete set of state grid points; wherein, since the time delay response feature remains constant in the current scheduling cycle, the system does not need to perform grid division on this time delay dimension, but only needs to perform grid processing on a fixed two-dimensional slice plane intercepted along the time delay response feature dimension. (2) Set a fixed reverse time integration step size (e.g.) ), and at the starting point of reverse evolution (i.e. At point ), the optimal evolution cost function value for each state grid point is... Directly initialize to the implicit objective function described in step S201. The calculated value; through this initialization operation, that is, mathematically using the target boundary set as the terminal boundary constraint condition for calculating the evolution cost of each state point; (3) At each inverse time integration step, perform an algebraic extremum search based on first physics principles for each grid point in the multidimensional state space. The algebraic extremum search specifically includes: finding the optimal evolution strategy by using a nested exhaustive traversal: for each candidate discrete power point in the climbing-limited interval, traverse the discrete temperature points in the set of environmental temperature disturbance extremities, calculate the inner product of the spatial gradient vector and the state change rate vector, and find the environmental disturbance point that maximizes the value of the inner product, so as to simulate the worst strategy (i.e., mathematically) that the environment attempts to maximize the risk of the system going out of bounds. Operator); then, among the maximum inner product values ​​corresponding to all candidate discrete power points, find the minimum minimax inner product value and its corresponding optimal power control point, in order to simulate the optimal regulation strategy of the controller to minimize the risk of exceeding limits (i.e., mathematically). Operators); It should be noted that, in order to capture the non-smooth geometric creases generated during the evolution of the safety boundary and avoid numerical oscillations, this embodiment adopts the fifth-order weighted essential non-oscillatory (WENO-JS) format based on local Lax-Friedrichs numerical Hamiltonian to extract the spatial gradient vector of the optimal evolution cost function at the current grid point. The specific characteristic flux splitting and nonlinear adaptive weight calculation process of the WENO-JS format follows the standard Jiang-Shu discrete operator form. The system uses this discrete operator to automatically assign the dominant weight of the smooth sub-template when crossing the non-differentiable crease region in the multidimensional state space, thereby stably extracting the one-sided spatial partial derivatives and combining them into the spatial gradient vector. (4) The found minimax inner product value is used as the time partial derivative of the partial differential equation at the grid point. The first-order forward Euler method is used to multiply the minimax inner product value with the inverse time integration step size. The product is then added to the optimal evolution cost function value of the current integration step size to update the optimal evolution cost function value of the grid point in the next inverse time integration step size. (5) Repeat the extreme value search and integral update steps until... The calculation terminates when the process proceeds layer by layer to the total duration of the preset scheduling period (e.g., 900 seconds), and the numerical solution of the partial differential equation is output, which is the set of optimal evolution costs corresponding to all discrete state grid points in the multidimensional state space.

[0029] Discrete state grid points in the multidimensional state space whose optimal evolution cost is greater than or equal to a preset safety threshold are extracted, and the set formed by them is defined as the safe reachable region. The preset safety threshold is recommended to be 0 to meet the conservative and high reliability requirements in the actual operation of the power system. Specifically, if the optimal evolution cost obtained by substituting a certain discrete state grid point into the partial differential equation is positive, its positive value attribute physically represents that, within a given preset scheduling period, even if the energy storage system continuously faces the worst disturbance sequence within the set of extreme values ​​of the ambient temperature disturbance, the controller can still optimize to obtain a feasible adjustment power trajectory within the ramp-limited interval, so as to ensure that the evolution trajectory of the dynamic cell temperature and dynamic state of charge of the energy storage body is always enveloped within the safe space defined by the target boundary set. This state point is then measured as a safe reachable state. To verify the ability of the Hamilton-Jacobi partial differential equation to characterize physical boundaries, Figure 2Simulation results based on multidimensional spatial modeling using simulation platforms such as Matlab / Simulink are presented. The figure shows a three-dimensional schematic diagram of the shrinkage of the safe reachable domain boundary of the state space with evolutionary delay. It can be clearly observed in the figure that as the delay response characteristic represented by the vertical axis increases, the geometric boundary of the safe operating domain composed of the state of charge and cell temperature exhibits significant nonlinear collapse and inward contraction. The simulation results intuitively reveal the squeezing effect of communication and execution lag on the physical control margin of energy storage, proving the scientific nature and necessity of the proposed method of broadening the delay characteristics into the state space dimension.

[0030] This step, based on Hamilton-Jacobi reachability analysis and differential game theory, mathematically unifies and decouples the nonlinear thermodynamic mechanism, time delay response characteristics, and external thermal disturbances at the underlying level of the energy storage system. This provides a numerical basis for limiting the risk of the energy storage entity going out of bounds under complex operating conditions and the time delay coupling effect of the underlying hardware.

[0031] S3. Perform dimensionality reduction decomposition and time evolution integration on the safe reachable domain to obtain a non-convex evolution boundary. Divide the evolution boundary into a set of locally convex regions and perform logical reconstruction to obtain the segmented time delay ramping constraint of the energy storage subject.

[0032] In this embodiment, the dimensionality reduction decomposition of the secure reachable domain includes: S301. The multidimensional state space is discretized into a grid to obtain the numerical array of the safe reachable domain at each grid node, as follows: In this embodiment, the multidimensional state space is discretized into a grid. The coordinate axis division range and discrete grid step size are consistent with those in step S202 when solving the Hamilton-Jacobi partial differential equation. That is, the grid is divided along the dynamic charge state axis, the dynamic cell temperature axis, and the time delay response characteristic dimension according to a preset step size. The difference is that the discrete grid in S202 is used to support the spatial partial derivative calculation of the windward difference scheme, while this step aims to combine the optimal evolution cost values ​​corresponding to each grid point in the multidimensional state space calculated in step S2 based on the one-to-one mapping relationship of discrete coordinate points, and according to the preset physical dimension order ( The dynamic state of charge, dynamic cell temperature, and time delay response characteristics are sequentially and fixedly mapped to the first, second, and third modal dimensions of a higher-order tensor, reshaping them into a higher-order numerical tensor structure. As an example parameter configuration, if the dynamic state of charge axis interval is set to 10% to 90% with a preset step size of 0.5%, then 161 discrete nodes are correspondingly divided; if the dynamic cell temperature axis interval is set to 20℃ to 60℃ with a preset step size of 0.2℃, then 201 discrete nodes are correspondingly divided. Combined with a fixed parameter slice set for the time delay response characteristic axis, the numerical array is specifically represented algebraically as a structure with a specification of... A third-order numerical tensor, in which each tensor element stores the optimal evolutionary cost of the corresponding physical state coordinate point.

[0033] S302. Based on a preset rank evaluation truncation criterion, the numerical array is decomposed into tensors to extract local tensor cores and corresponding mapping basis that characterize the evolution of the safe reachable domain, as follows: Using a tensor column decomposition algorithm, matrix singular value decomposition is performed continuously a preset number of times (the number of modal orders of the higher-order numerical tensor minus 1) along the first to third modal order of the higher-order numerical tensor (i.e., corresponding to the dynamic state of charge, dynamic cell temperature, and time delay response characteristics respectively). During each singular value decomposition, the tensor data modulus of the current dimension is expanded into a two-dimensional matrix (denoted as...). And construct the first autocorrelation matrix (i.e. ) and the second autocorrelation matrix (i.e. The eigenvalues ​​and orthogonal eigenvectors of the first and second autocorrelation matrices are calculated respectively. The eigenvalues ​​are then processed by non-negative square root extraction and arranged in descending order of numerical value to construct a singular value diagonal matrix. The orthogonal eigenvectors of the first autocorrelation matrix are combined to obtain the left singular matrix, and the orthogonal eigenvectors of the second autocorrelation matrix are combined to obtain the right singular matrix. The first element of the singular value diagonal matrix is ​​taken as the largest singular value, and the small singular value feature dimensions are truncated according to a preset rank evaluation truncation criterion. The rank evaluation truncation criterion specifically involves setting a desired truncation error threshold (e.g., ...). Using the product of the truncation error threshold and the maximum singular value as the truncation criterion, all non-dominant singular values ​​in the singular value diagonal matrix whose singular values ​​are less than the truncation criterion are identified and removed. Simultaneously, the column vectors and row vectors corresponding to the removed singular values ​​in the left singular matrix and the transpose of the right singular matrix are removed, thereby achieving rank reduction and dimensionality compression in algebraic space. Through the truncation criterion, the original numerical array is decoupled into a set of local tensor cores. These local tensor cores algebraically represent the low-rank approximate subspace matrix of the evolution law of the original high-dimensional safe boundary. As an example, a scale of... After the tensor sequence decomposition and truncation, the third-order numerical tensor yields three low-dimensional local tensor cores, whose algebraic scales are reduced to [values ​​to be filled in]. , as well as (in and This refers to the number of effective singular values ​​that are adaptively retained in the first and second singular value decompositions (i.e., those greater than the truncation benchmark), thereby achieving dimensionality reduction and compression of the original high-dimensional features. The mapping basis is the orthogonal feature vector matrix generated by the singular value decomposition process, which is embedded in the local tensor core and used to restore and map the decomposed low-dimensional subspace features back to the original physical state space. The local tensor core and its contained mapping basis, through the multiplication operation of tensor shrinking, together constitute the low-rank compressed numerical expression of the securely reachable domain. The discretization calculation formula of the tensor column decomposition is as follows: (9) in, The original numerical array obtained by discretization; The total modal order of the higher-order numerical tensor (for example, the total modal order of a third-order numerical tensor is 3). Traverse the dimension index indicator variable (values ​​range from 1 to...) ); For the first The index number of the discrete grid node in each physical dimension; The first one extracted by tensor sequence decomposition The local tensor core, in tensor algebra, has a modular expansion form that orthogonally includes the basis features that map the low-dimensional subspace back to the original physical dimension, i.e., the mapping basis. The index for the internal shrinking summation connecting adjacent local tensor cores; For the truncation error threshold The first one retained after evaluation and screening The rank of the core features of each tensor (as specified by boundary conditions) ); This step uses tensor dimensionality reduction to effectively compress the exponential computational complexity of the high-dimensional state space to the polynomial level, enabling the underlying embedded microprocessor to meet the computing power requirements for online real-time solving.

[0034] It should be noted that in a specific implementation scenario where the time delay response feature is set as a single-parameter slice, the tensor column decomposition is equivalent to the singular value decomposition of a two-dimensional matrix. This embodiment adopts a high-order tensor column decomposition architecture, which aims to build an algorithmic foundation compatible with dimensional expansion: when the actual working condition assessment needs to introduce augmented dimensions or multi-parameter slices of time delay response features, the high-dimensional state tensor after dimensional expansion can be processed without replacing the core computing framework, thereby ensuring the engineering applicability and architectural flexibility of the decomposition logic in complex scenarios.

[0035] In this embodiment, obtaining the non-convex evolution boundary includes: S303. Starting from the end time of the scheduling cycle and ending at the current time, reverse reachability propagation is performed on the local tensor core to obtain the tensor evolution characteristics, as follows: Based on the evolution mechanism of Hamiltonian-Jacobi physics constructed in step S202, the dynamic trajectory of the local tensor core is deduced in the reduced-dimensional algebraic space, thereby capturing and extracting the evolution characteristics of the safety boundary over time. The deduction process is as follows: Using the dynamic Galerkin projection method, the local tensor core is taken as the orthonormal basis operator of the tangent space, for the Hamilton-Jacobi partial differential equation containing... For the strongly nonlinear Hamiltonian part of the operator, a discrete empirical interpolation method and its built-in greedy algorithm are introduced. When iteratively approximating the column space of this strongly nonlinear Hamiltonian, the convergence condition is set as the local approximation residual norm being less than a preset precision threshold (e.g., ...). The iteration count does not exceed the set maximum iteration count (e.g., 50 times). In the iteration loop, the algorithm calculates the interpolation residual vector between the current projected subspace approximation matrix and the real nonlinear column vector one by one. After locating the component index with the largest absolute value in the residual vector, the algorithm extracts the discrete space node with the largest local approximation error and linear independence corresponding to the component index in the tensor grid as the skeleton point. The number of skeleton points extracted is not higher than the preset constant of the dimension of the low-rank tensor quantum space (the preset constant ranges from 15 to 30) to avoid full space traversal calculation during projection integration. Finally, by calculating the minimum residual of the Hamilton-Jacobi partial differential equation in the orthogonal basis direction, the original high-dimensional partial differential equation is transformed into a set of ordinary differential equations with the time partial derivatives of each local tensor core as the dependent variable at the algebraic level. Using the fourth-order Runge-Kutta method, a fixed reverse time integration step size (e.g., 1 second) is set, the same as in step S202. On the reverse time axis, the time representing the end of the scheduling cycle is... Taking the current time as the starting point for integration, and using the current time (at this moment) The total duration of the scheduling cycle is taken as the integration endpoint. Within each reverse time integration step, four evolution slope characteristic values ​​within the same integration step interval are obtained by solving the ordinary differential equation system. Specifically: the tensor state at the current integration starting point is substituted into the ordinary differential equation system to obtain the first initial slope; the system state prediction is shifted to half the integration step using the first initial slope to obtain the second predicted slope; the second predicted slope is used to recalculate the third predicted slope at half the integration step; the system state is then extrapolated to the current integration step endpoint using the third slope to obtain the fourth predicted slope; subsequently, the first to fourth slopes are assigned according to a preset weight ratio (e.g., respectively). , , and The four slopes are weighted and summed to obtain the final equivalent slope. Then, the equivalent slope is multiplied by the reverse time integration step size, and the product is added to the values ​​of the elements inside the local tensor core at the previous integration time, completing one reverse reachability propagation update. When the integration progresses to the current time (i.e., ... After the total duration of the scheduling period is equal to the total duration of the scheduling period, the latest set of local tensor core parameter matrices obtained is the tensor evolution feature characterizing the change of the boundary of the safe reachable domain.

[0036] S304. Extract isosurface feature points that satisfy preset boundary conditions from the tensor evolution features, and reconstruct the feature points back to the multidimensional state space using the mapping basis to obtain the non-convex evolution boundary, as follows: First, a multidimensional matrix multiplication is performed on the local tensor cores carrying the tensor evolution characteristics (i.e., updated by integration in step S303) and their corresponding mapping basis to reconstruct and calculate the optimal evolutionary cost value of each grid node at the current moment. The multidimensional matrix multiplication and reconstruction do not change the value of the evolutionary cost value itself, but rather perform chain multiplication and inner product elimination along the modal order of the higher-order tensors, multiplying adjacent local tensor cores with embedded mapping basis along their shared internal rank index, thereby decoding and restoring the implicit algebraic expression in the dimension-reduced compressed form into an explicit cost value distribution on the full-dimensional grid space. Subsequently, feature extraction is performed on the full-space grid numerical distribution obtained by reconstruction. The optimal evolutionary cost value is equal to 0 (i.e., the physical critical condition that distinguishes between safe and out-of-bounds states) is used as the preset boundary condition. The tensor grid node index items (i.e., the algebraic subscript number of the multidimensional array) that make the optimal evolutionary cost value satisfy the preset boundary condition are extracted and used as isosurface feature points. Finally, based on the one-to-one mapping relationship of discrete coordinate points during the grid discretization process described in step S301, the system reverse maps and restores the algebraic grid index of the extracted isosurface feature points to the actual physical coordinate values ​​in the original multidimensional physical state space (for example, mapping and restoring the grid index to the specific state of charge percentage and cell temperature in degrees Celsius). The restored set of discrete physical coordinate points envelops a closed geometric hypersurface with a nonlinear convex-concave shape in the multidimensional Euclidean geometric space, and this hypersurface is the non-convex evolution boundary.

[0037] In this embodiment, obtaining the segmented time delay ramping constraint of the energy storage body includes: S305. The evolutionary boundary is transformed into a set of local convex regions composed of multicellular shapes through geometric segmentation, as follows: A hyperplane-based Delaunay triangulation algorithm is used to geometrically segment the closed geometric space enclosed by the non-convex evolutionary boundary. The specific steps of this geometric segmentation include: taking the set of discrete physical coordinate points on the non-convex evolutionary boundary as input, performing Delaunay triangulation within the feasible domain space enclosed by the non-convex evolutionary boundary; and generating several basic three-dimensional simplexes based on the empty circumsphere core criterion (i.e., the multidimensional circumsphere of any generated simplex does not contain any other coordinate points from the input point set). The upper limit of the number of basic three-dimensional simplexes is the square-order asymptotic convergence of the total number of input point sets (i.e., theoretically increasing in magnitude). ,in (This refers to the total number of discrete physical coordinate points). The dihedral angle is calculated based on the angle between the outward normal vectors on both sides of the common plane of adjacent simplexes. If the dihedral angle is not greater than 180 degrees, it is determined to satisfy the convexity criterion, and spatial cluster merging is performed on adjacent simplexes that satisfy the convexity criterion. After a single merging generates a new spatial cluster, the dihedral angle between this new spatial cluster and its external adjacent simplexes or the common plane between adjacent spatial clusters is recalculated, and adjacent regions that satisfy the convexity criterion are continuously merged until the dihedral angle between any adjacent independent spatial clusters or simplexes within the internal closed geometric space is greater than 180 degrees. At this point, the internal closed geometric space is divided into... A set of non-overlapping multidimensional convex polycells, which together constitute the set of local convex regions in multidimensional geometric space.

[0038] S306. Based on the boundary geometric parameters of the set of local convex regions, a set of inequalities characterizing the switching logic between convex regions is constructed as the piecewise time delay ramping constraint, as follows: First, extract the unshared external triangular facets from all simplexes belonging to the spatial cluster of the multidimensional convex polytope within the set of local convex regions, those exclusively occupied by a single simplex. Based on the vertex coordinate set of these triangular facets, define the geometrically closed region of the multidimensional convex polytope in the multidimensional state space. Further, based on this geometrically closed region, construct the hyperplane algebraic equations representing the outer boundaries of each convex polytope. (in (where the state vector is in the three-dimensional state space), for each extracted unshared external triangular facet, normal vector mapping and intercept calculation are performed using its vertex coordinates to analytically obtain the geometric normal vector matrix in the hyperplane algebraic equation. Intercept offset column vector This set of parameters directly constitutes the subspace defining the interior of a single convex multicell (i.e., independent geometric constraints). The bottom boundary geometry parameters of ). Subsequently, the construction dimension is State evolution transfer matrix The elements in the matrix are state vectors. The first-order partial derivatives of the dynamic state of charge, dynamic cell temperature, and time delay response characteristics with respect to the regulating power of the energy storage entity to be solved are given. Since the time delay response characteristics remain unchanged within a single cycle, the corresponding matrix element values ​​are 0. Based on the state evolution transfer matrix, the state vector is... Equivalent expression as initial state vector (i.e., the current state of charge, initial cell temperature, and time delay response characteristics) and the regulation power to be solved The linear evolutionary relationship, i.e. Substituting the linear evolution relationship into the independent geometric constraints and performing an algebraic transformation, we obtain the independent linear power constraints. (wherein, the equivalent power normal vector matrix) Equivalent offset column vector ); Next, in order to integrate the independent linear boundary constraints distributed in different local convex regions into a globally solvable mathematical model, the Big M method is used in conjunction with binary discrete switching variables representing the activation states of each convex region to logically reconstruct the independent linear power constraints of each local convex region: by setting a large positive number (recommended value) As a relaxation term, the independent linear boundary constraints are transformed into a system of mixed-integer linear inequalities controlled by 0-1 switching variables, representing the logic of the optimization solver searching and switching between multiple convex regions; the system of mixed-integer linear inequalities is the piecewise time-delay ramping constraint, and the formulas of the inequalities of the piecewise time-delay ramping constraint are as follows: (10) in, The main energy storage regulating power control variable to be solved in the current scheduling cycle is expressed in MW. and The first The equivalent power normal vector matrix and the equivalent offset column vector of a local convex region; The total number of the multidimensional convex polycells; To characterize the first The 0-1 discrete state switching variable of the activation state of a local convex region, when The optimization trajectory of the time representation system is strictly controlled by the first... A convex region, which makes the corresponding linear power rigidly constrained (i.e., ) takes effect; The large positive number set is used in This ensures that the corresponding inequality always holds, thereby achieving the mathematical relaxation transformation of the constraints in the inactive region.

[0039] This step utilizes tensor dimensionality reduction evolution and Big M method logic reconstruction to transform the high-dimensional non-convex dynamic physical security boundary into a mixed-integer linear constraint format that can be directly solved by the optimization solver. This effectively reduces the computational complexity caused by the high-dimensional state space and realizes the efficient embedding of energy storage time-delay coupling characteristics in the grid coordinated dispatch clearing process.

[0040] S4. Solve the pre-established collaborative clearing objective function based on the piecewise time delay ramping constraint to obtain the clearing capacity of the energy storage entity, including: S401. Based on the topological relationship and local thermal disturbance data of each energy storage unit within the energy storage entity, a thermal accumulation sensitivity matrix is ​​constructed and features are extracted to obtain thermal difference adjustment parameters and correct the piecewise time delay ramping constraint, as follows: By reading the device node mapping configuration file pre-stored in the energy storage management system, the series and parallel topology relationships of a large number of battery clusters within the energy storage body are obtained. Addressing the thermal imbalance caused by uneven cooling airflow and physical spatial distribution—for example, battery clusters located at the end of the cooling airflow are more sensitive to heat accumulation caused by the same charging and discharging current compared to those at the air inlet—the system collects the node temperature and charging / discharging power of each energy storage unit in real time as local thermal disturbance data through the CAN bus or Modbus / TCP industrial protocol communication port of the underlying battery management system (BMS). A thermal accumulation sensitivity matrix is ​​then constructed to characterize the impact of unit power disturbance on the temperature evolution of each spatial node. The formula for the thermal accumulation sensitivity matrix is ​​as follows: (11) in, For dimension The thermal accumulation sensitivity matrix, For the first The local cell temperature of an individual battery cluster energy storage unit, in °C. For the injection of the first The active power of each battery cluster node, in MW. The total number of battery cluster energy storage units that serve as independent heat dissipation nodes within the energy storage entity; Using the same algorithmic logic as in step S302 for matrix singular value decomposition, the heat accumulation sensitivity matrix is ​​decomposed into the product of a left singular matrix, a singular value diagonal matrix, and a right singular matrix, with the first element of the singular value diagonal matrix taken as the maximum singular value; a monotonically decreasing negative exponential mapping operator (e.g.) is constructed. ,in To decay the penalty weight scalar, For the maximum singular value, The maximum singular value is converted into a thermal difference adjustment parameter with a value range of (0, 1], based on a preset recommended value range of 0.3 to 0.7 (a thermal safety reduction factor). The thermal difference adjustment parameter is used to perform scalar multiplication on the equivalent offset column vector in the piecewise time delay ramp constraint inequality to complete the correction of the piecewise time delay ramp constraint. The corrected piecewise time delay ramp constraint causes the equivalent offset column vector representing the boundary intercept of the local convex region to shrink proportionally, thereby mapping and coupling the thermal physical constraints of the underlying energy storage unit to the macroscopic operational optimization feasible domain, avoiding the risk of local thermal accumulation caused by the obtained clearing command exceeding the limit.

[0041] S402. Construct a collaborative clearing objective function that includes energy storage lifetime degradation indicators, operating costs, and scheduling deviation penalty indicators, as follows: The energy storage life degradation index The formula is as follows: (12) in, To calculate the unit capacity replacement cost of the energy storage entity by extracting the financial asset ledger and equipment depreciation model during the construction period of the energy storage power station project; Index for discrete scheduling periods; The pre-exponential frequency factor is an empirical constant characterizing the rate of chemical reactions inside the battery cell. It is obtained by conducting accelerated aging cycle charge-discharge experiments on the same batch of energy storage cells in advance and fitting the experimental data using the least squares method. To characterize the activation energy of the chemical reaction barrier inside the battery cell, based on the different material characteristics of lithium iron phosphate battery cells, the recommended typical engineering value range is 40 kJ / mol to 60 kJ / mol. The ideal gas constant is 8.314 J / (mol·K). For time period The estimated absolute temperature within the vector is derived directly from the state vector. The dynamic cell temperature component is calculated by converting Celsius to Kelvin absolute temperature scale. and The time periods to be solved are respectively The cleared power in the domestic electricity market and the cleared reserve capacity in the ancillary services market, in MW; The operating costs The formula is as follows: (13) in, The unit comprehensive operation and maintenance cost coefficient, which takes into account the AC-DC conversion loss of the inverter and the daily inspection cost, is calculated and extracted based on the rated efficiency on the factory nameplate of the energy storage inverter and the historical operation and maintenance details. The time step for discrete scheduling periods (e.g., converting a 15-minute scheduling cycle to an hourly unit, i.e., 0.25 hours). The scheduling deviation penalty index The formula is as follows: (14) in, The deviation penalty price is obtained by reading the spot market settlement rules or grid connection dispatch agreements published by the power trading center, and the unit is yuan / MWh; For the time period The power difference generated within is calculated by subtracting the maximum auxiliary service regulation power that the energy storage entity can actually provide after being constrained by the segmented time delay ramping constraint from the theoretical AGC frequency regulation demand command curve issued by the power grid a day earlier (i.e., ) was calculated; With the goal of minimizing the energy storage lifetime degradation index, operating cost, and dispatch deviation penalty index, a collaborative clearing objective function is constructed, and the calculation formula of the collaborative clearing objective function is as follows: (15) in, To achieve a coordinated clearing objective function, The total number of time periods in the scheduling cycle (for example, if the global scheduling time span is set to 24 hours and the discrete scheduling time period step is 0.25 hours, then the total number of time periods is 96). , , For the preset multi-objective optimization weight coefficients, recommended values ​​are 0.5, 0.3, and 0.2.

[0042] S403. Based on the modified ramping constraint, the collaborative clearing objective function is optimized and solved to obtain the clearing capacity, which includes the clearing power of the electric energy market and the clearing reserve capacity of the ancillary service market for each scheduling period of the energy storage entity, as follows: (1) Establish the energy and reserve coupling balance constraint, as shown in the following formula: (16) in, For the first Energy storage main body regulation power for each scheduling cycle; (2) Establish upper and lower boundary constraints for dynamic charge state, as shown in the following formula: (17) (18) (19) in, and These are the charging power variable and the discharging power variable obtained by decoupling the power clearing from the electricity market, respectively. and All are not less than 0 and satisfy ; For the first Percentage of dynamic state of charge at the end of each time period; This represents the state of charge at the end of the previous time period; and These are the charging efficiency coefficient and the discharging efficiency coefficient, respectively. Their specific values ​​are determined based on the battery's factory cycle test report and the unidirectional conversion loss of the energy storage inverter. The rated total capacity of the energy storage system is expressed in MWh. The time step for discrete scheduling; and These are the preset lower and upper limits of the safe state of charge (e.g., 10% and 90%). (3) Establish the bidirectional maximum effective power constraint for the energy storage inverter (PCS), as shown in the following formula: (20) (twenty one) in, The maximum bidirectional continuous active power allowed by the energy storage inverter equipment under a specific power factor is obtained by reading the factory nameplate calibration value of the PCS hardware on site and the operating log parameters entered by the energy management system (EMS). (4) Market demand capacity limit constraint, the formula is as follows: (twenty two) (twenty three) in, and The first Within each dispatch cycle, the power grid dispatch center issues the lower and upper limits of active power injection to the busbars of the energy storage entities' grid connection points, in MW. The maximum frequency regulation reserve capacity allocated to this energy storage entity by the regional ancillary services market, in MW, is obtained by pulling the day-ahead / intraday ancillary services market demand forecast curve from the power trading center in real time. Based on all constraints, including the modified piecewise time delay ramping constraint, the power clearing power in the electric energy market is... Clearing standby capacity in the ancillary services market Configured as independent decision variables; employing commercial mathematical programming solvers such as CPLEX or Gurobi, calling their mixed-integer nonlinear programming or mixed-integer quadratic constraint solution kernels, the collaborative clearing objective function is solved by branch and bound optimization. The output is the optimal solution set of decision variables that minimizes the global objective function, i.e., the power clearing power of the electricity market during energy transfer in each scheduling cycle. And the clearing-out reserve capacity of the ancillary services market, which is mandatory to suppress grid frequency fluctuations. .

[0043] S404. After obtaining the cleared capacity of the energy storage entity, the cleared capacity is converted into an active power signal sequence, and based on the signal transmission delay in the delay response characteristics, the active power signal sequence is feedforward corrected to generate a duty cycle command to drive the energy storage entity to perform power regulation, as follows: Using a zero-order hold and setting a preset step size (e.g., 5ms) that matches the sampling rate of the underlying controller, the clearing capacity values ​​obtained from the optimization solution for discrete time periods are resampled along the time axis at intervals of the preset step size to obtain a continuously stepped high-frequency discrete active power signal sequence. By reading the time stamp difference between the communication messages and the underlying controller of the energy storage unit, the communication network signal transmission delay in the delay response characteristics is extracted. Subsequently, a feedforward lead compensation controller is adopted, which advances the control command issuance pointer relative to the theoretical timestamp by pre-reading the data to be issued in the command buffer. The step size is adjusted to compensate for the signal transmission delay in advance, resulting in a corrected target power sequence. A low-level proportional-integral closed-loop control algorithm is used to generate a duty cycle command to drive the energy storage unit to regulate power. Specific steps include: using the target power sequence... The set value is the power collected in real time by the sensor at the grid connection point. The instantaneous deviation between the setpoint and the feedback value is used as the algorithm input. The deviation is linearly amplified by a proportional term, and the historical cumulative values ​​of the deviation are summed by an integral term to calculate the dynamic adjustment amount of the output duty cycle. This adjustment amount is then added to a preset duty cycle offset reference term. Above, output the final duty cycle instruction. The calculation formula for generating the duty cycle command is as follows: (twenty four) in, This is a duty cycle instruction, expressed as a percentage (%). The active power measured and fed back at the grid connection point; and These are the proportional gain constant and integral time constant of the underlying closed-loop PI controller, respectively. Their values ​​are preset based on the hardware transfer function of the energy storage inverter output filter. In lithium iron phosphate energy storage systems, their recommended engineering ranges are as follows: , ; The step size is the variable for integration. and These are the historical command power value and the historical feedback power value within the integration domain, respectively. The feedforward duty cycle bias reference term is obtained by pre-querying an offline mapping lookup table of "active power - duty cycle". This offline mapping lookup table is constructed based on the circuit topology of a three-level neutral-point clamped inverter, and its underlying mapping function relationship satisfies the space vector pulse width modulation rule. The specific expression is as follows: ,in The AC side reference phase voltage amplitude under the corresponding target clearing power; The DC bus voltage measurement value is sampled in real time by a voltage transformer configured on the DC side of the inverter to match the preset high-frequency sampling rate (e.g., 10kHz) of the underlying controller's main frequency. The linear modulation region reference coefficient is pre-calibrated based on the inherent losses of the inverter hardware dead time and switching frequency, and the recommended value range is 0.80 to 0.866. Finally, the system sends the duty cycle command to the underlying drive circuit through the high-speed peripheral interface inside the microprocessor (such as hardware PWM register mapping or SPI bus); the underlying drive circuit maps the duty cycle signal into a physical level signal that controls the on and off of the insulated gate bipolar transistor inside the energy storage inverter through an optical fiber isolated transmission channel, thereby realizing the controlled regulation of the active power output of the energy storage body. To demonstrate the closed-loop control performance of this embodiment, Figure 3 A comparison of simulation waveforms of the control response of the energy storage inverter based on the Matlab / Simulink platform is presented. In this simulation environment, communication time delay jitter with an average value of 342ms and typical high-frequency measurement white noise are set. As can be seen from the figure, the conventional closed-loop control (red solid line) lacks time delay feedforward compensation, resulting in phase lag and power overshoot caused by integral saturation when the scheduling command changes abruptly. In contrast, the scheme of this embodiment (green solid line) compensates for the time delay dead zone through feedforward processing, realizing effective following of the scheduling command. Its output waveform rising edge conforms to the physical constraint of the ramp rate of 12.5MW / s of the underlying power device, and the steady-state power ripple is maintained at a low level.

[0044] This step corrects the ramp constraint by using the thermal accumulation sensitivity matrix, optimizes the allocation of electrical energy and reserve capacity using a multi-objective function, and generates duty cycle instructions by combining communication delay feedforward compensation. This improves the accuracy of scheduling instructions from top-level optimization to bottom-level mapping of power devices and effectively reduces response deviations caused by signal lag.

[0045] The above formulas are all dimensionless calculations. The formulas are derived from software simulations using a large amount of collected data, and are the closest to the real situation. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0046] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0047] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0048] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0049] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0050] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A collaborative clearing method for the energy storage spot market based on time delay characteristics, characterized in that, Includes the following steps: The physical operation data of the energy storage unit is obtained, the time delay response characteristics of the energy storage unit are extracted, and the ramp-limited range is defined based on the cell temperature evolution and heat accumulation characteristics. A multidimensional state space including the aforementioned time delay response characteristics is constructed, and the safe reachability domain of the energy storage entity within the multidimensional state space is calculated based on the aforementioned ramp-limited interval. The safe reachable domain is decomposed by dimensionality reduction and time evolution integral to obtain a non-convex evolution boundary. The evolution boundary is divided into a set of locally convex regions and logically reconstructed to obtain the segmented time delay ramping constraint of the energy storage subject. The pre-established collaborative clearing objective function is solved based on the segmented time delay ramping constraint to obtain the clearing capacity of the energy storage entity.

2. The method according to claim 1, characterized in that, The extraction of the time delay response characteristics of the energy storage entity and the definition of the ramp-limited range based on cell temperature evolution and thermal accumulation characteristics include: Based on the time response deviation between the historical dispatch instructions of the energy storage entity and the corresponding historical actual output data, time delay response features are extracted. By acquiring environmental thermal disturbance data and physical operation data and combining them with the time delay response characteristics, the cell temperature evolution trajectory corresponding to the energy storage main body's adjustment power is deduced; The range of adjustable power that makes the evolution trajectory satisfy the preset heat accumulation constraint is extracted as the climbing-limited interval.

3. The method according to claim 2, characterized in that, The construction includes a multi-dimensional state space comprising the time delay response features, including: A multidimensional state space is constructed, including the time delay response characteristics and the physical operation data, and a target boundary set of the multidimensional state space is set at the deadline of the preset scheduling period based on the preset heat accumulation constraint.

4. The method according to claim 3, characterized in that, The calculation of the safe reachability domain of the energy storage entity within the multidimensional state space includes: Within a preset scheduling period, based on the climbing-limited interval and the environmental thermal disturbance data, the extreme value calculation is performed on the evolution path from each state point in the multidimensional state space to the target boundary set, so as to obtain the optimal evolution cost corresponding to each state point. The set of state points whose optimal evolution cost satisfies the preset safety threshold is taken as the safe reachable domain.

5. The method according to claim 4, characterized in that, The dimensionality reduction decomposition of the secure reachable domain includes: The multidimensional state space is discretized into a grid to obtain the numerical array of the safe reachable domain at each grid node; The numerical array is decomposed into tensors based on a preset rank evaluation truncation criterion to extract local tensor cores and corresponding mapping basis that characterize the evolution of the safe reachable domain.

6. The method according to claim 5, characterized in that, The process of obtaining the non-convex evolution boundary includes: Starting from the end time of the scheduling cycle and ending at the current time, reverse reachability propagation is performed on the local tensor core to obtain tensor evolution features that characterize boundary changes. Extract isosurface feature points that satisfy preset boundary conditions from the tensor evolution features, and use the mapping basis to reconstruct and restore the feature points to the multidimensional state space to obtain the non-convex evolution boundary.

7. The method according to claim 6, characterized in that, The step of obtaining the segmented time delay ramping constraint of the energy storage entity includes: The evolutionary boundary is transformed into a set of local convex regions composed of multicellular shapes through geometric segmentation; Based on the boundary geometric parameters of the set of local convex regions, a set of inequalities characterizing the switching logic between convex regions is constructed as the segmented time delay ramping constraint.

8. The method according to claim 7, characterized in that, The process of obtaining the cleared capacity of the energy storage entity includes: Based on the topological relationship of each energy storage unit within the energy storage entity and local thermal disturbance data, a thermal accumulation sensitivity matrix is ​​constructed and features are extracted to obtain the thermal difference adjustment parameters of the energy storage entity. The segmented time delay ramping constraint is corrected using the thermal difference adjustment parameter, and the cooperative clearing objective function is optimized based on the corrected ramping constraint to obtain the clearing capacity.

9. The method according to claim 1 or 8, characterized in that, The collaborative clearing objective function includes at least the energy storage lifetime degradation index, operating cost, and dispatch deviation penalty index; the clearing capacity includes the energy storage entity's clearing power in the electricity market and the clearing reserve capacity in the ancillary services market.

10. The method according to claim 9, characterized in that, After obtaining the clearing capacity of the energy storage entity, the method further includes: The clearing capacity is converted into an active power signal sequence, and based on the signal transmission delay in the delay response characteristics, the active power signal sequence is feedforward corrected to generate a duty cycle command that drives the energy storage entity to perform power regulation.

Images