Energy storage power station operation mode switching method based on super capacitor and intelligent algorithm

Abstract

The invention provides an energy storage power station operation mode switching method based on a super capacitor and an intelligent algorithm, which relates to the technical field of energy storage power stations, and comprises the steps of obtaining real-time state parameters and load demand information of the super capacitor, constructing a coupling state space of electric energy, heat energy and aging dimensions, calculating a state evolution track, and predicting a state transition result. A dynamic planning framework is combined to optimize a switching strategy, optimal matching of a super capacitor state track and a power grid demand track is realized, a coupling weight coefficient is dynamically adjusted through a feedback mechanism, the operation efficiency of an energy storage power station is improved, the service life of the super capacitor is prolonged, and the stability of a power grid is enhanced.

Description

TECHNICAL FIELD

[0001] The present application relates to energy storage power station technology, and particularly to an energy storage power station operation mode switching method based on supercapacitors and intelligent algorithms. BACKGROUND

[0002] With the widespread application of renewable energy and the rapid development of smart grids, energy storage technology plays an increasingly important role in power systems. As a high-power-density, long-cycle-life energy storage device, supercapacitors exhibit unique advantages in grid peak shaving and frequency regulation, new energy grid connection, and microgrid applications due to their fast charging and discharging characteristics and high efficiency. Energy storage power stations can effectively improve the stability and reliability of the grid, smooth load fluctuations, and promote the consumption of renewable energy by properly scheduling and managing supercapacitors and other energy storage devices.

[0003] Currently, the operation modes of energy storage power stations typically include charging mode, discharging mode, standby mode, and other states. According to the grid load demand and the state of energy storage devices, energy storage power stations need to switch between different modes to meet the dynamic operation requirements of the system. Traditional energy storage power station operation mode switching methods are mainly based on simple rule logic or fixed schedules, which cannot fully adapt to complex and changing grid demands and the dynamic characteristics of energy storage devices.

[0004] Traditional methods often only focus on the electrical energy state of energy storage devices, ignoring the impact of thermal energy state and aging process on device performance, and cannot fully reflect the actual operating state of supercapacitors, resulting in one-sidedness in mode switching decisions. Existing technologies lack accurate prediction capabilities for the evolution of energy storage device states, making it difficult to effectively estimate device state changes over a period of time, resulting in a lack of foresight and scientificity in mode switching decisions. Existing switching methods typically use static decision strategies, making it difficult to adaptively adjust based on actual operating conditions, and cannot maximize the service life of supercapacitors and other energy storage devices while meeting grid demands, making it difficult to achieve economic and efficient operation of energy storage power stations. SUMMARY

[0005] The embodiments of the present application provide an energy storage power station operation mode switching method based on supercapacitors and intelligent algorithms, which can solve the problems in the prior art.

[0006] In a first aspect, the embodiments of the present application provide an energy storage power station operation mode switching method based on supercapacitors and intelligent algorithms, comprising:

[0007] Acquire real-time state parameters of a super capacitor in an energy storage power station and load demand information of a power grid side, based on the real-time state parameters, calculate state evolution trajectories of the super capacitor along a coupling state space of the super capacitor in an electric energy dimension, a thermal energy dimension and an aging dimension under different operating modes by constructing the coupling state space, and obtain a state transition prediction result;

[0008] According to the state transition prediction result and the load demand information of the power grid side, determine a candidate operating mode sequence of the energy storage power station within a preset time window; according to the candidate operating mode sequence, search within the preset time window through a dynamic programming framework, so that a mode switching time combination with a minimum deviation degree of a super capacitor state trajectory and a power grid demand trajectory is obtained, introduce a long-term performance degradation penalty term calculated based on the trajectory feature vector into a cost function of the dynamic programming framework, and generate a switching strategy;

[0009] Control the energy storage power station to perform a mode switching action based on the switching strategy, collect actual state trajectory data after the switching as feedback information, use a deviation between the feedback information and the state transition prediction result to adjust coupling weight coefficients between the electric energy dimension, the thermal energy dimension and the aging dimension in the calculation of the state evolution trajectory through a gradient iteration algorithm, and return the adjusted coupling weight coefficients to the state analysis step for updating subsequent state evolution trajectory calculation.

[0010] Based on the real-time state parameters, calculate state evolution trajectories of the super capacitor along a coupling state space of the super capacitor in an electric energy dimension, a thermal energy dimension and an aging dimension under different operating modes by constructing the coupling state space, and obtain a state transition prediction result.

[0011] Map the real-time state parameters to the coupling state space to obtain a current state point containing an electric energy dimension state component, a thermal energy dimension state component and an aging dimension state component, and construct a local linearization region in the coupling state space based on the current state point, and obtain coupling weight coefficients by calculating a partial derivative matrix of each dimension state component with respect to other dimension state components in the local linearization region;

[0012] Based on the coupling weight coefficients, establish a state increment calculation rule, the state increment calculation rule converts the electric energy dimension state component at a current time into a thermal energy dimension state increment through a first coupling weight coefficient, and converts the thermal energy dimension state increment into an aging dimension state increment through a second coupling weight coefficient, and uses the state increment calculation rule to iteratively deduce the evolution process of each dimension state component from the current state point in a time direction to form a state evolution trajectory containing time sequence correlation.

[0013] constructing a local linearization region in the coupling state space based on the current state point, and obtaining coupling weight coefficients by calculating a partial derivative matrix of each dimension state component to other dimension state components within the local linearization region comprises:

[0014] constructing an adaptive local linearization region based on the current state point, calculating a rate of change of the current state point in each dimension, determining a linearization radius of each dimension according to the rate of change, and constructing an ellipsoid-shaped local linearization region taking the linearization radius as a principal axis; within the ellipsoid-shaped local linearization region, a hierarchical importance sampling method is used to generate a plurality of sampling point groups;

[0015] performing recursive principal component analysis on the plurality of sampling point groups, constructing nonlinear coupling characteristics between state components by constructing a high-order matrix decomposition, and using singular value decomposition to reduce and reconstruct the high-order matrix; based on the reconstructed characteristics, dynamic coupling weight coefficients are obtained by calculating conditional partial derivatives of each dimension state component to other dimension state components under different time scales.

[0016] According to the candidate operation mode sequence, a dynamic programming framework is used to search within the preset time window, so that the combination of mode switching time points with the minimum deviation degree of the super capacitor state trajectory and the grid demand trajectory comprises:

[0017] Based on the candidate operation mode sequence, a state transition graph is constructed, a hierarchical beam search method is used to prune the state transition graph, and a preset number of paths with the minimum deviation at each layer are retained;

[0018] The pruned state transition graph is applied to dynamic programming search, the preset time window is divided into a plurality of time segments, and the cumulative deviation of all feasible state transitions within each time segment is calculated; based on the cumulative deviation, an optimal state transition matrix is constructed, and the combination of operation mode switching time points that minimizes the deviation degree of the super capacitor state trajectory and the grid demand trajectory is determined by backtracking the optimal state transition matrix.

[0019] Based on the candidate operation mode sequence, a state transition graph is constructed, a hierarchical beam search method is used to prune the state transition graph, and a preset number of paths with the minimum deviation at each layer are retained;

[0020] The state transition graph is divided into a plurality of search layers in time sequence, the deviation value of each layer node to the target state is calculated, and a priority queue is established for each layer node based on the deviation value;

[0021] Perform hierarchical beam search on the nodes in the priority queue, in each layer of search process, construct a local state transition subgraph based on the state transition gradient between each node and its adjacent nodes, calculate the path connectivity between nodes according to the local state transition subgraph, and take the weighted combination of the path connectivity and the node bias value as the node score, and retain the preset number of nodes with the optimal score as the retained nodes of the current layer; the path of the retained node in the state transition graph is taken as a candidate path, the cumulative bias of the candidate path is calculated, and the candidate path whose cumulative bias exceeds the dynamic threshold is removed.

[0022] Introduce a long-term performance degradation penalty term calculated based on the trajectory feature vector into the cost function of the dynamic programming framework, and the generation of the switching strategy includes:

[0023] Calculate the performance degradation feature within a sliding time window based on the trajectory feature vector, extract the change rate of the trajectory feature vector within adjacent time windows, and calculate the degradation speed of each performance indicator according to the change rate; segmentally linearly fit the main direction component under multiple time scales, take the projection value of each performance indicator on the main direction component as the corresponding weight coefficient, obtain the comprehensive performance degradation rate, and construct a long-term performance degradation constraint based on the comprehensive performance degradation rate;

[0024] Combine the long-term performance degradation constraint and the state transition cost to obtain a stage cost calculation method; for each decision stage, determine the cumulative cost of different switching schemes based on the stage cost calculation method; solve the optimal switching sequence through a recursive step so that the sum of the cumulative costs is minimized, and generate a switching strategy according to the optimal switching sequence.

[0025] Adjust the coupling weight coefficients between the electrical energy dimension, the thermal energy dimension and the aging dimension in the state evolution trajectory calculation by using the deviation between the feedback information and the state transition prediction result through a gradient iteration algorithm, including:

[0026] Calculate the deviation between the feedback information and the state transition prediction result, and decompose the deviation according to the electrical energy dimension, the thermal energy dimension and the aging dimension;

[0027] Construct a coupling sensitivity matrix between the state dimensions for the deviation decomposition result, and perform segmented calculation on the coupling sensitivity matrix according to the correlation of the state changes in different time scales to obtain the gradient change direction between the dimensions;

[0028] Iteratively adjust the coupling weight coefficients based on the gradient change direction, and determine the iteration step size by calculating the change rate of the state response between adjacent dimensions in each iteration process, and complete the adjustment of the weight coefficients when the state deviation is less than a preset deviation threshold.

[0029] In a second aspect, the present application provides an electronic device, comprising:

[0030] a processor;

[0031] a memory for storing processor-executable instructions;

[0032] wherein the processor is configured to invoke the instructions stored in the memory to perform the method described above.

[0033] In a third aspect, the present application provides a computer-readable storage medium having stored thereon computer program instructions, which when executed by a processor, implement the method described above.

[0034] The present application has the following advantages:

[0035] The present application can accurately predict the state evolution trajectory under different operating modes by constructing a coupling state space of supercapacitors in the dimensions of electrical energy, thermal energy and aging, and comprehensively considering the multi-dimensional operating state of supercapacitors, thereby ensuring the normal operation of the energy storage power station while improving the service life and efficiency of the supercapacitors.

[0036] The present application adopts a dynamic programming framework to search for the optimal mode switching time combination, and introduces a long-term performance degradation penalty term calculated based on the trajectory feature vector into the cost function, so that the generated switching strategy can not only meet the load demand of the power grid side, but also effectively reduce the performance degradation of the supercapacitors, thereby achieving optimal control of the operation of the energy storage power station.

[0037] The present application uses actual state trajectory data as feedback information, dynamically adjusts the coupling weight coefficients between the dimensions of electrical energy, thermal energy and aging by using the deviation between the feedback information and the state transition prediction result, forms a closed-loop self-learning mechanism, and continuously optimizes the accuracy of state evolution trajectory calculation, so that the mode switching decision is more intelligent and adaptive. BRIEF DESCRIPTION OF DRAWINGS

[0038] Figure 1 Fig. 1 is a flowchart of an energy storage power station operating mode switching method based on supercapacitors and intelligent algorithms according to an embodiment of the present application;

[0039] Figure 2 Fig. 2 is a hierarchical beam search state transition graph pruning flowchart based on a candidate operating mode sequence according to an embodiment of the present application. DETAILED DESCRIPTION

[0040] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the protection scope of the present application.

[0041] The technical solutions of the present application will be described in detail below with specific embodiments. The following specific embodiments can be combined with each other, and some embodiments may not be described again for the same or similar concepts or processes.

[0042] Figure 1 The flowchart of the method for switching the operation mode of the energy storage power station based on the super capacitor and the intelligent algorithm in the embodiments of the present application is shown in FIG. 1, which comprises the following steps. Figure 1

[0043] The real-time state parameters of the super capacitor in the energy storage power station and the load demand information of the power grid side are acquired, and based on the real-time state parameters, the state evolution trajectory of the super capacitor along the coupling state space in the energy dimension, the thermal energy dimension and the aging dimension under different operation modes is calculated by constructing the coupling state space of the super capacitor in the energy dimension, the thermal energy dimension and the aging dimension, to obtain a state transition prediction result.

[0044] According to the state transition prediction result and the load demand information of the power grid side, a candidate operation mode sequence of the energy storage power station in a preset time window is determined; according to the candidate operation mode sequence, a mode switching time combination in which the deviation degree of the super capacitor state trajectory and the power grid demand trajectory is minimum is searched in the preset time window through a dynamic programming framework, a long-term performance degradation penalty term calculated based on the trajectory feature vector is introduced into a cost function of the dynamic programming framework, and a switching strategy is generated.

[0045] Based on the switching strategy, the energy storage power station is controlled to perform a mode switching action, and the actual state trajectory data after the switching is collected as feedback information. By using the deviation between the feedback information and the state transition prediction result, the coupling weight coefficients between the energy dimension, the thermal energy dimension and the aging dimension in the state evolution trajectory calculation are adjusted through a gradient iteration algorithm, and the adjusted coupling weight coefficients are fed back to the state analysis step for updating the subsequent state evolution trajectory calculation.

[0046] In an optional embodiment, based on the real-time state parameters, the state evolution trajectory of the super capacitor along the coupling state space in the energy dimension, the thermal energy dimension and the aging dimension under different operation modes is calculated by constructing the coupling state space of the super capacitor in the energy dimension, the thermal energy dimension and the aging dimension, which comprises the following steps.

[0047] ​mapping the real-time state parameters into a coupling state space to obtain a current state point comprising an electrical energy dimension state component, a thermal energy dimension state component and an aging dimension state component, and constructing a local linearization region in the coupling state space based on the current state point, and obtaining coupling weight coefficients by calculating a partial derivative matrix of each dimension state component with respect to other dimension state components within the local linearization region;

[0048] establishing a state increment calculation rule based on the coupling weight coefficients, the state increment calculation rule converting the electrical energy dimension state component at a current time into a thermal energy dimension state increment through a first coupling weight coefficient, and converting the thermal energy dimension state increment into an aging dimension state increment through a second coupling weight coefficient, and using the state increment calculation rule to iteratively deduce an evolution process of each dimension state component from the current state point in a time direction to form a state evolution trajectory comprising time sequence correlation.

[0049] The mapping process of the real-time state parameters establishes a three-dimensional coordinate system. The electrical energy dimension state component is calculated by multiplying the voltage square by the capacitance value and then dividing by two to obtain the instantaneous stored electrical energy, and then dividing by the rated capacity to normalize to the range of zero to one. The thermal energy dimension state component uses a distributed temperature sensor array, with platinum resistance sensors arranged on the positive electrode, negative electrode and shell, respectively, with a weight distribution of 40%, 40% and 20% for weighted average, and the final temperature is normalized by subtracting the ambient temperature and dividing by the maximum allowed temperature rise. The aging dimension state component combines internal resistance measurement and cycle count, the internal resistance is measured by a 1 kHz alternating current impedance method, the ratio of the current internal resistance to the initial internal resistance is used as an aging indicator, and the ratio of the cycle count to the rated life is used as another indicator, and the two indicators are linearly weighted to obtain the final aging state component.

[0050] The current state point determines three normalized state components as coordinate values to locate the position in the coupling state space, the coordinate accuracy is kept to four decimal places, and the update frequency is synchronized with the sensor sampling frequency to 1 kHz. The effectiveness check limits the electrical energy dimension to not more than 1.05, the thermal energy dimension to not more than 1.2, and the aging dimension to not less than 0.5. The local linearization region is constructed by expanding the neighborhood radius centered on the current state point, the electrical energy dimension radius is adaptively set to 3% to 8% according to the energy storage state, the thermal energy dimension radius is set to 2% to 6% according to the temperature state, and the aging dimension radius is set to 1% to 4% according to the health state.

[0051] The partial derivative matrix calculation adopts a five-point central difference method, applies a perturbation of one-tenth of the neighborhood radius to each state component, and observes the response changes of other state components. The partial derivative of electrical energy to thermal energy is calculated by multiplying the current square by the internal resistance and dividing by the thermal capacity. The partial derivative of thermal energy to aging is based on the Arrhenius equation exponential relationship. The partial derivative of aging to electrical energy adopts a linear relationship of capacity attenuation. The partial derivative of electrical energy to aging considers the influence of deep discharge and adopts a piecewise linear model. The partial derivative of thermal energy to electrical energy adopts a quadratic polynomial fitting. The partial derivative of aging to thermal energy is calculated by the internal resistance change rate. The matrix update period is set to one hundred milliseconds, and the element precision is maintained to six decimal places.

[0052] The coupling weight coefficient obtains the singular value decomposition of the partial derivative matrix. The first coupling weight coefficient is equal to the absolute value of the matrix first row second column element multiplied by a normalization coefficient, with a numerical range of zero point one to zero point nine and a default value of zero point five. The second coupling weight coefficient is equal to the matrix second row third column element after hyperbolic tangent transformation, with a numerical range of zero point zero one to zero point five and a default value of zero point one. The weight coefficient adopts a double buffer storage mechanism, with the main and standby buffers being used alternately, and the effectiveness check ensures that the value is within a reasonable range.

[0053] The state increment calculation rule establishes the causal relationship chain between dimensions. The state decrement of the electrical energy dimension is equal to the power loss multiplied by the time step and divided by the total capacity. The power loss is the current square multiplied by the internal resistance. The state increment of the thermal energy dimension includes the electrical energy conversion part and the environmental heat dissipation part. The conversion efficiency is determined by the first coupling weight coefficient, and the heat dissipation is calculated using Newton's cooling law. The state increment of the aging dimension is equal to the state increment of the thermal energy multiplied by the second coupling weight coefficient and then multiplied by the temperature index acceleration factor. The time step is set to one millisecond, and the calculation order strictly follows the dependence relationship of electrical energy, thermal energy, and aging.

[0054] The iterative deduction process starts from the current state point to establish a time axis forward recursive loop. Each iteration step includes state increment calculation, state update, boundary check, and result storage. The deduction depth is set to three hundred sixty thousand calculation steps corresponding to one hour, and multi-thread parallelization is adopted to share the calculation load. The state space is divided into sub-regions to allocate threads. The calculation process maintains a ring queue history cache with a length of ten thousand history points, supporting two retrieval methods of time stamp and serial number.

[0055] The state evolution trajectory forms a continuous curve by connecting the time sequence state points. The trajectory data structure includes time and three state dimensions, and each trajectory point stores four-tuple data. The trajectory point connection adopts linear interpolation, the density is adjusted adaptively according to the curvature, and the smoothness is guaranteed by cubic B-spline fitting. The trajectory feature extraction includes geometric features such as length, curvature, and twist, and dynamic features such as change rate, acceleration, and energy dissipation rate. The feature vector dimension is twenty and is normalized. The storage format adopts a hierarchical compression strategy, with incremental encoding in the time dimension and quantization encoding in the state dimension, with a compression rate of seventy-five percent.

[0056] Data case set initial state energy dimension zero point eight, thermal energy dimension zero point three, aging dimension zero point nine five, corresponding energy storage percentage eighty, temperature thirty-five degrees Celsius, health status ninety-five percent. The first coupling weight coefficient is zero point four, the second coupling weight coefficient is zero point zero eight, and the time step is one millisecond. After one hour of evolution, the energy dimension decreases to zero point seven two, the thermal energy dimension increases to zero point four five, and the aging dimension decreases to zero point nine four eight. The total length of the trajectory is zero point two three units, the maximum curvature three point five occurs at the twenty-second minute, and the energy dissipation rate decreases from zero point zero five to zero point zero three. The boundary condition processing includes starting the energy saving mode when the energy is lower than zero point zero five, triggering the thermal protection when the thermal energy is more than zero point eight five, and sending a maintenance reminder when the aging is lower than zero point eight.

[0057] In an optional implementation, a local linearization region is constructed based on the current state point in the coupling state space, and the coupling weight coefficients are obtained by calculating the partial derivative matrix of each dimensional state component to other dimensional state components within the local linearization region, comprising:

[0058] An adaptive local linearization region is constructed based on the current state point, the change rate of the current state point in each dimension is calculated, the linearization radius of each dimension is determined according to the change rate, and the linearization radius is used as the principal axis to construct an ellipsoidal local linearization region; within the ellipsoidal local linearization region, a plurality of sampling point groups are generated by using a stratified importance sampling method;

[0059] Recursive principal component analysis is performed on the plurality of sampling point groups, the nonlinear coupling characteristics between state components are constructed by constructing high-order matrix decomposition, and the high-order matrix is reconstructed by dimension reduction using singular value decomposition; based on the reconstructed characteristics, dynamic coupling weight coefficients are obtained by calculating the conditional partial derivative of each dimensional state component to other dimensional state components under different time scales.

[0060] An adaptive local linearization region is constructed based on the current state point, which is realized by real-time monitoring of the instantaneous state change rate of the energy dimension, the thermal energy dimension and the aging dimension in the three-dimensional coupling state space of the battery energy storage system. The energy dimension change rate is measured by a power sensor, and the normalized change amount of the energy state per unit time is obtained by dividing the total capacity of the battery, and the change rate value range is zero point zero zero one to zero point five per second.

[0061] The thermal energy dimension rate of change monitors the temperature gradient through a distributed temperature sensor array, calculates the derivative of temperature with respect to time, and divides by the maximum allowed temperature rise to obtain a normalized thermal energy rate of change with a range of zero point zero zero five to zero point two per second. The aging dimension rate of change combines the internal resistance growth rate and the cycle count increment, with the internal resistance rate of change obtained through continuous impedance measurements and the cycle count rate of change calculated from the charge and discharge detection algorithm, and is normalized to a range of zero point zero zero zero one to zero point zero one per second after comprehensive weighting.

[0062] The linearization radius of each dimension is determined using a rate adaptive algorithm, with the electrical energy dimension linearization radius equal to the electrical energy rate of change multiplied by the response time constant and then multiplied by the safety factor, the response time constant set to ten seconds, and the safety factor in the range of one point five to three point zero, resulting in an electrical energy linearization radius in the final range of zero point zero one five to one point five. The thermal energy dimension linearization radius takes into account the thermal inertia effect, and is equal to the thermal energy rate of change multiplied by the thermal time constant and then multiplied by the expansion coefficient, with the thermal time constant set to sixty seconds and the expansion coefficient in the range of two point zero to five point zero, resulting in a thermal energy linearization radius in the final range of zero point zero three to six point zero.

[0063] The aging dimension linearization radius is based on the irreversible nature of the aging process, and is equal to the aging rate of change multiplied by the aging time constant and then multiplied by the prediction coefficient, with the aging time constant set to three thousand six hundred seconds and the prediction coefficient in the range of ten point zero to fifty point zero, resulting in an aging linearization radius in the final range of zero point zero three six to one point eight.

[0064] The ellipsoidal local linearization region is constructed by taking the linearization radii of the three dimensions as the major axis radii of the ellipsoid and the current state point as the center of the ellipsoid to establish a three-dimensional ellipsoid region. The ellipsoid equation uses the standard form, with the electrical energy dimension coordinate minus the square of the current electrical energy state divided by the square of the electrical energy linearization radius, plus the thermal energy dimension coordinate minus the square of the current thermal energy state divided by the square of the thermal energy linearization radius, and then plus the aging dimension coordinate minus the square of the current aging state divided by the square of the aging linearization radius, with the sum less than or equal to one. The ellipsoid volume calculation uses the three-dimensional ellipsoid volume formula, which is equal to four times the circumference ratio divided by three multiplied by the product of the three radii, with a volume range of zero point zero zero one to fifty cubic units. The ellipsoid boundary check uses the Mahalanobis distance from the point to the center of the ellipsoid, with a distance less than one indicating that the point is inside the ellipsoid.

[0065] The stratified importance sampling method divides the ellipsoidal region into three concentric ellipsoidal layers according to the Mahalanobis distance from the center. The inner layer has a Mahalanobis distance of zero to zero point three three, the middle layer has a Mahalanobis distance of zero point three three to zero point six seven, and the outer layer has a Mahalanobis distance of zero point six seven to one point zero. The sampling density is highest in the inner layer, and the number of sampling points accounts for fifty percent of the total number. The number of sampling points in the middle layer accounts for thirty percent, and the number of sampling points in the outer layer accounts for twenty percent. The importance weight distribution adopts an exponential decreasing mode, the inner layer weight coefficient is two point zero, the middle layer weight coefficient is one point zero, and the outer layer weight coefficient is zero point five. The sampling point generation adopts Latin hypercube sampling to ensure uniform distribution, and the sampling points generated in each ellipsoidal layer are mapped to the three-dimensional state space through ellipsoidal coordinate transformation.

[0066] The multiple sampling point groups are formed by grouping the sampling points in the ellipsoidal region according to state similarity clustering. The K-means clustering algorithm is used to set the number of clusters to eight to sixteen. The similarity measure uses weighted Euclidean distance, with a weight of zero point four for the electrical energy dimension, a weight of zero point four for the thermal energy dimension, and a weight of zero point two for the aging dimension. The K-means++ algorithm is used for cluster center initialization, and the iterative convergence condition is set to a cluster center movement distance less than zero point zero zero one or an iteration number greater than one hundred. Each sampling point group contains eight to thirty-two sampling points, and the standard deviation of the distance between the sampling points in the group is less than zero point zero five, and the distance between the groups is greater than zero point two.

[0067] Recursive principal component analysis is performed on each sampling point group, and the covariance matrix of the sampling point group is calculated. The covariance matrix has a dimension of three by three corresponding to the three state dimensions. Principal component extraction uses eigenvalue decomposition, and the principal components with a cumulative contribution rate exceeding ninety-five percent are retained, usually the first two to three principal components are retained. Recursive update uses a sliding window mechanism, and the window length is set to one thousand time steps. When a new sampling point is added, the oldest sampling point is removed to keep the window length constant. The principal component loading matrix is updated using an incremental learning algorithm, and the learning rate is set to zero point zero one. The convergence criterion is that the norm of the change in the load vector is less than zero point zero zero one.

[0068] High-order matrix construction extends the original three-dimensional state vector using quadratic and cubic cross terms of state components. The quadratic term contains nine elements corresponding to the product of each two state components, and the cubic term contains twenty-seven elements corresponding to the product of each three state components. The extended high-order feature vector has a dimension of thirty-nine, including three linear terms, nine quadratic terms, and twenty-seven cubic terms. Nonlinear coupling feature extraction identifies important coupling relationships by calculating the variance contribution of high-order terms. High-order terms with a variance contribution greater than five percent are retained as key coupling features. Feature selection uses a recursive feature elimination algorithm, which removes the lowest ten percent of features in each iteration until the number of remaining features stabilizes at ten to fifteen.

[0069] The singular value decomposition performs dimension reduction reconstruction on the high-order matrix with dimension of sampling points multiplied by high-order feature dimension, usually two hundred to five hundred rows multiplied by thirty-nine columns. The singular value decomposition generates a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix, with singular values arranged in descending order. The dimension reduction reconstruction retains singular vectors corresponding to the top eight to twelve largest singular values, with reconstruction error controlled within five percent of the original data variance. The reconstructed feature matrix dimension is reduced to the number of sampling points multiplied by the number of retained singular values, achieving effective compression of the feature space.

[0070] Different time scales include a short time scale of one to ten seconds, a medium time scale of ten to sixty seconds, and a long time scale of sixty to six hundred seconds. The short time scale mainly reflects the transient response of electrochemical reactions, the medium time scale corresponds to the heat conduction process, and the long time scale involves the cumulative effect of aging. Each time scale uses a corresponding filter window, with one second sliding average for the short time scale, ten second sliding average for the medium time scale, and sixty second sliding average for the long time scale. The time scale weight distribution is dynamically adjusted according to the current operating conditions, with the short time scale weight increasing to zero point six under high power operating conditions and the long time scale weight increasing to zero point five under steady state operating conditions.

[0071] The conditional partial derivative calculation uses the finite difference method to observe the response changes of other state components by applying a small perturbation to each state component in the reconstructed feature space. The perturbation amplitude is set to one thousandth of the current state value, and the perturbation direction includes positive and negative directions to improve numerical accuracy. The conditional partial derivatives of the electrical energy state component with respect to the thermal energy and aging state components are calculated by fixing the aging state and observing the thermal energy response, and by fixing the thermal energy state and observing the aging response. The conditional partial derivatives of the thermal energy state component with respect to the electrical energy and aging state components use a similar method, fixing the third state component to calculate the partial derivative relationship between the other two state components. The conditional partial derivatives of the aging state component with respect to the electrical energy and thermal energy state components consider the one-way nature of the aging process, and use historical data fitting to obtain approximate values of the partial derivatives.

[0072] The dynamic coupling weight coefficients are obtained by weighted fusion of conditional partial derivatives at different time scales. The first coupling weight coefficient corresponds to the coupling strength between the electrical energy and thermal energy dimensions, and is equal to the absolute value of the conditional partial derivative of electrical energy with respect to thermal energy multiplied by the short-time scale weight, plus the absolute value of the conditional partial derivative of thermal energy with respect to electrical energy multiplied by the medium-time scale weight, with a value ranging from 0.1 to 0.9. The second coupling weight coefficient corresponds to the coupling strength between the thermal energy and aging dimensions, and is equal to the absolute value of the conditional partial derivative of thermal energy with respect to aging multiplied by the medium-time scale weight, plus the absolute value of the conditional partial derivative of aging with respect to thermal energy multiplied by the long-time scale weight, with a value ranging from 0.05 to 0.5. The third coupling weight coefficient corresponds to the coupling strength between the electrical energy and aging dimensions, and is calculated by integrating the cross-correlation function of electrical energy and aging state components at different time scales, with a value ranging from 0.01 to 0.3.

[0073] The data case uses operational data from a lithium battery energy storage system. The initial state is set with an electrical energy dimension of 0.75, a thermal energy dimension of 0.4, and an aging dimension of 0.92. The rate of change of electrical energy is 0.08 per second, the rate of change of thermal energy is 0.03 per second, and the rate of change of aging is 0.002 per second. The corresponding linearization radii are 0.24 for electrical energy, 0.36 for thermal energy, and 0.216 for aging. Two hundred and fifty-six sampling points are generated within an ellipsoidal region, divided into twelve sampling point groups. Recursive principal component analysis retains the first three principal components, with a cumulative contribution rate of 97%. High-order matrix singular value decomposition retains the first ten singular values, with a reconstruction error of 3.5%. The calculated conditional partial derivatives at different time scales are: short-term electrical energy partial derivative with respect to thermal energy of 0.63, medium-term thermal energy partial derivative with respect to aging of 0.27, and long-term aging partial derivative with respect to electrical energy of 0.15. The final dynamic coupling weight coefficients are 0.58 for the first coupling weight coefficient, 0.32 for the second coupling weight coefficient, and 0.18 for the third coupling weight coefficient. The weight coefficient update cycle is 50 milliseconds, and the calculation precision is maintained to four decimal places.

[0074] In one optional implementation, based on the candidate operating mode sequence, a dynamic programming framework is used to search within the preset time window for the combination of mode switching times that minimizes the deviation between the supercapacitor state trajectory and the grid demand trajectory. This combination includes:

[0075] A state transition graph is constructed based on the candidate operating mode sequence. The state transition graph is pruned using a hierarchical bundle search method to retain a preset number of paths with the smallest deviation at each level.

[0076] The state transition graph after pruning is applied to a dynamic programming search, the preset time window is divided into multiple time segments, and the cumulative deviation of all feasible state transitions in each time segment is calculated; based on the cumulative deviation, an optimal state transition matrix is constructed, and the combination of operation mode switching time points that minimizes the deviation degree of the super capacitor state trajectory and the grid demand trajectory is determined by backtracking the optimal state transition matrix.

[0077] According to the system operating environment and constraint conditions, a candidate operation mode sequence is obtained, which includes different working states such as charging mode, discharging mode, standby mode, etc. Each operation mode corresponds to different super capacitor charging and discharging strategies and power curves. For example, in the charging mode, the super capacitor obtains energy from the grid at a preset power; in the discharging mode, the super capacitor provides energy to the load at a specific power; in the standby mode, the super capacitor keeps the current energy state unchanged.

[0078] Based on the obtained candidate operation mode sequence, a state transition graph is constructed, which represents the operation modes adopted by the system at different time points and their mutual conversion relationship within a preset time window. The nodes in the state transition graph represent the operation mode at a specific time point, and the edges represent the transition process from one mode to another. Each edge is also associated with a deviation value, which represents the deviation degree between the super capacitor state trajectory and the grid demand trajectory caused by the state transition.

[0079] After the state transition graph is constructed, a hierarchical beam search method is used to prune the state transition graph to improve the search efficiency. The core idea of hierarchical beam search is that in the search process, only the preset number of paths with the smallest deviation are retained at each layer, and other paths with lower possibility are discarded. In specific implementation, for each time layer in the preset time window, the system calculates the deviation values of all state transitions, and sorts them in ascending order of deviation values, and only retains the top K state transitions with the smallest deviation, where K is a preset beam width parameter. For example, if the preset beam width K = 5, only 5 paths with the smallest deviation are retained in each time layer to enter the next layer of search.

[0080] In practical applications, the beam width parameter K can be adjusted according to the calculation resources and accuracy requirements. A larger K value can improve the search accuracy but increase the calculation complexity, and a smaller K value can speed up the search but may miss the global optimal solution. In this embodiment, through experimental comparison, K = 10 is selected as the beam width parameter, which achieves a good balance between calculation efficiency and search accuracy.

[0081] After pruning the state transition graph, a dynamic programming search algorithm is applied to the pruned state transition graph. The core of the dynamic programming search is to decompose a complex problem into a series of sub-problems, and gradually build a global optimal solution by solving the sub-problems. In this embodiment, the preset time window is divided into a plurality of time segments, for example, a 24-hour time window is divided into 144 time segments of 10 minutes. For each time segment, the system calculates the cumulative deviation of all feasible state transitions.

[0082] The instantaneous deviation caused by the state transition of the current time segment, and the historical cumulative deviation from the starting time to the current time segment. The instantaneous deviation represents the difference between the state of the super capacitor in a certain operating mode and the grid demand within the current time segment. The historical cumulative deviation represents the total deviation that the system has generated before the current time segment. Adding these two parts of deviation together gives the cumulative deviation of the current state transition.

[0083] In a specific implementation, it is assumed that the preset time window is 24 hours, divided into 144 time segments, each with a length of 10 minutes. The candidate operating modes include three: charging mode, discharging mode and standby mode. For time segment t and operating mode m, the system calculates the instantaneous deviation abs(st,m-dt) between the state st,m of the super capacitor in this mode and the grid demand dt. Then, the instantaneous deviation is added to the historical cumulative deviation to obtain the cumulative deviation value.

[0084] Based on the calculated cumulative deviation value, the system constructs an optimal state transition matrix, which records the minimum cumulative deviation value and its corresponding predecessor state from the starting time to each time segment in each operating mode. For example, for time segment t and operating mode m, the optimal state transition matrix records dp[t][m] (indicating the minimum cumulative deviation value) and prev[t][m] (indicating the operating mode of the previous time segment to reach this minimum value).

[0085] After completing the construction of the optimal state transition matrix, the optimal operating mode switching time combination is determined by backtracking the matrix. The backtracking process starts from the end of the preset time window, and gradually backtracks to the starting time according to the predecessor state recorded in the optimal state transition matrix, to obtain the complete optimal operating mode sequence and its switching time.

[0086] In practical applications, taking an electric vehicle charging station as an example, the preset time window is 24 hours, the power grid demand curve is a typical daily load curve, including morning and evening peak and midnight valley characteristics. The initial energy state of the super capacitor is 50% of the rated capacity. By applying the above method, the system obtains the optimal operation mode switching time combination: 0:00-6:00 uses the charging mode (power grid load valley period), 6:00-9:00 uses the discharging mode (morning peak period), 9:00-17:00 uses the standby mode (daytime stable period), 17:00-22:00 uses the discharging mode (evening peak period), and 22:00-24:00 uses the charging mode (night low valley period starts). This mode switching combination reduces the total deviation degree of the super capacitor state trajectory and the power grid demand trajectory by 37.8%, significantly improving the system operation efficiency compared with the traditional fixed time period switching strategy.

[0087] Through the above method, the system can efficiently find the operation mode switching time combination that minimizes the deviation degree of the super capacitor state trajectory and the power grid demand trajectory within the preset time window, improving the operation efficiency of the energy management system and the stability of the power grid.

[0088] In an optional implementation, a state transition graph is constructed based on the candidate operation mode sequence, a hierarchical beam search method is used to prune the state transition graph, and a preset number of paths with the smallest deviation in each layer are retained, including:

[0089] The state transition graph is divided into multiple search layers in chronological order, the deviation value of each layer node from the target state is calculated, and a priority queue is established for each layer node based on the deviation value;

[0090] The hierarchical beam search is performed on the nodes in the priority queue, in each layer search process, a local state transition subgraph is constructed based on the state transition gradient between each node and its adjacent node, the path connectivity between nodes is calculated according to the local state transition subgraph, and the weighted combination of the path connectivity and the node deviation value is taken as the node score, and a preset number of nodes with the optimal score are retained as the retained nodes of the current layer; the path of the retained node in the state transition graph is taken as a candidate path, the cumulative deviation of the candidate path is calculated, and the candidate path with a cumulative deviation exceeding a dynamic threshold is eliminated.

[0091] As shown in Figure 2 , the method comprises:

[0092] A candidate operation mode sequence is obtained, which contains multiple operation state points. Taking an air conditioning system as an example, the candidate operation mode sequence contains parameter combinations such as temperature set point, wind speed level, and operation mode. For example, a candidate sequence can contain state points such as {(temperature 25℃, wind speed 3 level, cooling mode), (temperature 26℃, wind speed 2 level, cooling mode), (temperature 24℃, wind speed 4 level, cooling mode)}.

[0093] Based on the above candidate sequence construction state transition diagram, the nodes in the graph represent different running states, and the edges represent the transition relationship between states. The system divides the state transition diagram into multiple search layers according to time sequence, and each layer corresponds to a time point in the time sequence. Taking a 10-hour running cycle as an example, if the sampling interval is 1 hour, the state transition diagram is divided into 10 layers, and each layer contains all the running state nodes at this time point.

[0094] For each node in each layer, calculate its deviation value from the target state. The deviation value calculation considers multiple dimensions, including energy consumption deviation, comfort deviation, etc. Taking energy consumption as an example, if the target energy consumption is 2 kWh, and the expected energy consumption of a certain node is 2.5 kWh, then the energy consumption deviation is 0.5 kWh. The system normalizes and weights the sum of the dimension deviations to obtain the comprehensive deviation value. For example, the energy consumption deviation weight is 0.6, and the comfort deviation weight is 0.4, then the comprehensive deviation value is 0.6 x energy consumption deviation + 0.4 x comfort deviation.

[0095] Based on the calculated deviation value, a priority queue is established for each layer node. The priority queue is arranged in ascending order of deviation value, and the node with smaller deviation value has higher priority. For example, the first layer has 5 nodes with deviation values of 0.3, 0.5, 0.2, 0.7, and 0.4, respectively. The priority queue is sorted as: 0.2, 0.3, 0.4, 0.5, 0.7.

[0096] Perform hierarchical beam search on the nodes in the priority queue. In each layer search process, for each node in the priority queue, the system constructs a local state transition subgraph based on the state transition gradient between the node and its adjacent nodes. The state transition gradient represents the difficulty or cost of transitioning from one state to another. For example, the gradient value of the air conditioner adjusting from 25°C to 26°C is 0.1, while the gradient value of adjusting from 25°C to 22°C is 0.3, indicating that the latter state change is more drastic.

[0097] According to the local state transition subgraph, calculate the path connectivity between nodes. Path connectivity reflects the reachability and smoothness from the current node to the next layer node. Connectivity calculation considers the number of reachable nodes and the size of the transition gradient. For example, if the current node can connect to 3 nodes in the next layer with an average transition gradient of 0.2, then its path connectivity can be represented as 3 x (1-0.2) = 2.4.

[0098] The path connectivity and the node bias value are combined by weighting to obtain a node score, and the score calculation formula is: node score = a x (1-normalized bias value) + b x normalized connectivity, wherein a and b are weight coefficients, and a+b=1. In actual application, a=0.7 and b=0.3 can be set. For example, the normalized bias value of a certain node is 0.4, and the normalized connectivity is 0.8, then the node score = 0.7 x (1-0.4) + 0.3 x 0.8 = 0.66.

[0099] According to the node score, a preset number of nodes with the optimal score are reserved as the reserved nodes of the current layer in each layer, and the preset number can be set according to the computing resource and the accuracy requirement, for example, 5 nodes with the highest score are reserved. If the node scores of the second layer are 0.75, 0.66, 0.82, 0.58, 0.93, 0.79 and 0.62, and the preset reserved number is 5, then the system reserves the nodes with the scores of 0.93, 0.82, 0.79, 0.75 and 0.66.

[0100] The path of the reserved node in the state transition graph is taken as a candidate path, and for each candidate path, the system calculates the cumulative bias thereof. The cumulative bias is the weighted sum of the biases of all the nodes on the path, and the weight of a recent node is higher and the weight of a long-term node is lower. For example, a path contains 3 nodes, and the biases are 0.2, 0.3 and 0.4, and the weights are 0.5, 0.3 and 0.2, then the cumulative bias = 0.5 x 0.2 + 0.3 x 0.3 + 0.2 x 0.4 = 0.28.

[0101] A dynamic threshold is set to eliminate the candidate paths with too large cumulative bias, and the dynamic threshold is determined according to the cumulative bias distribution of the current reserved path, and can be set as the average cumulative bias plus 1.5 times of the standard deviation. For example, the average value of the cumulative bias of the current reserved path is 0.35, and the standard deviation is 0.1, then the dynamic threshold = 0.35 + 1.5 x 0.1 = 0.5. The system eliminates the candidate paths with the cumulative bias exceeding 0.5.

[0102] Through the above hierarchical beam search process, the state transition path is constructed and optimized layer by layer, and finally the operation mode sequence meeting the optimization target and the smooth state transition is obtained. For example, in a 10-layer state transition graph, the optimal path finally reserved is represented as: (temperature 25℃, wind speed 3 levels, refrigeration mode) → (temperature 25℃, wind speed 3 levels, refrigeration mode) → (temperature 26℃, wind speed 2 levels, refrigeration mode) →...→ (temperature 24℃, wind speed 2 levels, refrigeration mode). The path not only meets the energy consumption and comfort requirements, but also ensures the smooth and feasible state transition.

[0103] In an optional embodiment, a long-term performance degradation penalty term calculated based on the trajectory feature vector is introduced into a cost function of the dynamic programming framework, and generating the switching strategy comprises:

[0104] The performance degradation features are calculated based on the trajectory feature vectors within a sliding time window, the change rates of the trajectory feature vectors within adjacent time windows are extracted, and the degradation speeds of each performance index are calculated according to the change rates; the degradation speeds are piecewise linearly fitted at multiple time scales to obtain a principal direction component, the projection values of each performance index on the principal direction component are taken as corresponding weight coefficients, a comprehensive performance degradation rate is obtained, and a long-term performance degradation constraint is constructed based on the comprehensive performance degradation rate;

[0105] The long-term performance degradation constraint is combined with the state transition cost to obtain a stage cost calculation mode; for each decision stage, the cumulative cost of different switching schemes is determined based on the stage cost calculation mode; the optimal switching sequence is solved through a recursive step, so that the sum of the cumulative costs is minimized, and a switching strategy is generated according to the optimal switching sequence.

[0106] The performance degradation feature calculation is realized based on statistical analysis of the trajectory feature vectors within a sliding time window, the length of the sliding time window is set to sixty sampling points, and the window sliding step is fifteen sampling points. The trajectory feature vector includes three dimensions of battery capacity attenuation rate, internal resistance growth rate, and thermal cycle loss rate. All trajectory feature vector samples in the current time window are obtained, the weighted average value of each dimension is calculated, the weight coefficient adopts an exponential decay function, and the decay coefficient is 0.95. The weighted average result is taken as the performance degradation feature of the current time window.

[0107] The change rate of the trajectory feature vector within the adjacent time window is extracted by using the forward difference method, the trajectory feature vectors of the current time window and the previous time window are obtained. The difference between the feature vector of the subsequent window and the feature vector of the current window is calculated. The change rate is obtained by dividing the difference by the time interval, and the time interval is equal to the actual time corresponding to the window sliding step. Anomaly detection is performed on each dimension of the change rate vector, and the three-sigma rule is adopted. The values exceeding the range are replaced by the median of the adjacent window change rate.

[0108] The three-dimensional values of the change rate vector are obtained, the square of each dimension value is calculated, and the square root of the sum of the three square values is taken to obtain the module length of the change rate vector. The module length result is normalized by dividing by the historical maximum module length to obtain a degradation speed in the range of zero to one. The degradation speed is smoothed by moving average, and the arithmetic mean of the latest five degradation speed samples is taken.

[0109] The multi-time scale piecewise linear fitting processes the degradation speed time series in short-term, medium-term and long-term three levels. The short-term level takes the latest five minutes of data, the medium-term level takes the latest one hour of data, and the long-term level takes the latest twenty-four hours of data. The piecewise linear fitting is performed on the degradation speed data of each time level. The sliding window variance analysis is used to detect the segmentation point, and when the variance change exceeds the threshold, it is determined as the segmentation boundary. The least square method is used to calculate the linear fitting parameters for each segment to obtain the slope and intercept. The root mean square value of the fitting error is calculated, and it is compared with ten percent of the standard deviation of the original data to verify the fitting quality.

[0110] The principal direction component extraction is realized by principal component analysis. A fitting result matrix is constructed, with rows corresponding to time points and columns corresponding to the slope values of the three time scales. The covariance matrix of the fitting result matrix is calculated. The eigenvalues and eigenvectors of the covariance matrix are solved. The eigenvectors are sorted by eigenvalue size, and the first two to three eigenvectors are selected as the principal direction components. The variance contribution rate of each principal direction component is calculated to ensure that the cumulative contribution rate reaches 85%. The principal direction components are normalized to unit vectors.

[0111] The current performance indicator vector is obtained, which contains three dimensions of degradation speed values. Vector inner product operation is performed on each principal direction component. The performance indicator vector and the principal direction component vector are multiplied element by element and then summed. The projection value is normalized by dividing the length of the principal direction component to ensure that the projection value range is between -1 and 1. The absolute value of the projection value is taken as the weight coefficient. All weight coefficients are normalized to ensure that the sum of the weight coefficients is equal to one.

[0112] The degradation speed of each performance indicator and the corresponding weight coefficient are obtained. The degradation speed of each performance indicator is multiplied by the corresponding weight coefficient. The sum of all weighted results is obtained to get the initial comprehensive degradation rate. The regularization term is added, and the initial comprehensive degradation rate is added to the regularization coefficient multiplied by the square term of the degradation rate. The exponential moving average smoothing is performed on the comprehensive degradation rate, and the current value is multiplied by the smoothing coefficient and added to the historical value multiplied by one minus the smoothing coefficient.

[0113] A preset threshold of the comprehensive performance degradation rate is set, and the difference between the current comprehensive degradation rate and the threshold is calculated. When the difference is greater than zero, the difference is multiplied by a proportion coefficient to obtain a penalty cost. A time window for constraint verification is set, and the constraint range is updated in a rolling time domain. Based on the current degradation trend, the degradation level in the future time window is predicted, and when the predicted value exceeds the threshold, the constraint violation is triggered.

[0114] The state transition cost is calculated, including the weighted sum of three components of device start-stop loss, switching delay, and control loss. The long-term performance degradation constraint cost is calculated, based on the product of the penalty cost and the constraint weight coefficient. The state transition cost is multiplied by one minus the constraint weight coefficient, and the constraint cost is multiplied by the constraint weight coefficient. The two parts of the cost are added together to obtain the stage cost. The stage cost is normalized by dividing by the historical maximum stage cost and multiplying by ten.

[0115] The cumulative cost calculation performs a forward recursion for each decision stage, sets the decision stage time interval and the planning horizon length. All feasible switching schemes for the current stage are enumerated. The stage cost is calculated for each switching scheme. The current stage cost is added to the minimum cumulative cost of the subsequent stage, and the future cost is discounted using a time discount factor. The cumulative cost value corresponding to each switching scheme is stored.

[0116] The recursive step solution uses a backward recursion algorithm, initialized from the last stage of the planning horizon, with the cumulative cost of the terminal stage set to zero. The recursion is performed from the second last stage to the first. For each state of each stage, all switching decisions are enumerated. The stage cost corresponding to each decision is calculated and added to the cumulative cost of the next stage. The decision that minimizes the total cost is selected as the optimal decision. The optimal decision and the corresponding cumulative cost are stored in the optimal strategy table.

[0117] The cumulative cost sum minimization is achieved through global optimization, setting the optimization objective function as the sum of all stage costs within the planning horizon. The constraint conditions include device operation limits and performance degradation limits. The dynamic programming algorithm is used to solve the global optimal solution. For each state and stage combination, a feasibility check is performed to eliminate switching sequences that do not meet the constraints. The optimal switching sequence is obtained through backward recursion, and the optimal switching decision of each stage is recorded.

[0118] The switching strategy generation is based on the optimal switching sequence, extracting the optimal switching decision of each stage from the optimal strategy table. The execution time of each switching is determined, based on the decision stage time and the switching advance amount calculation. The switching object is determined, including the device identifier and the switching type that need to be switched. The switching target state is determined, based on the optimal decision corresponding device operation parameters. The switching execution parameters are set, including the switching rate, transition time, and safety check points. The switching strategy is encoded into a state machine form, supporting both conditional triggering and time-triggered execution modes.

[0119] The data case adopts a four-hour energy storage operation cycle verification, and the initial values of the three dimensions of the trajectory feature vector are set to zero point zero five, zero point zero eight and zero point one two. The sliding window length is sixty sampling points, and the sliding step is fifteen sampling points. The calculated change rate vector is zero point zero zero two, zero point zero zero three and zero point zero zero one. The degradation speed calculation result is zero point zero zero four. The multi-time scale fitting obtains a short-term slope of zero point zero zero one, a medium-term slope of zero point zero zero zero eight and a long-term slope of zero point zero zero zero three. The principal component analysis extracts two main direction components, and the characteristic values are zero point sixty two and zero point twenty eight. The projection value calculation obtains the weight coefficients of zero point four, zero point three five and zero point two five. The comprehensive performance degradation rate calculation result is zero point zero two three. The long-term constraint threshold is set to zero point zero five, and the penalty cost is zero. The dynamic programming sets twenty-four decision stages, and each stage has eight switching schemes. The backward recursion solution obtains the optimal switching sequence, and the cumulative cost is twelve point seven, including sixteen device switching operations. The switching strategy generation includes a complete description of the switching time, device object, target state and execution parameters. The strategy execution verification shows a success rate of ninety-eight point five percent, and the average execution delay is two point three seconds.

[0120] In an optional implementation, the coupling weight coefficients between the electric energy dimension, the thermal energy dimension and the aging dimension in the state evolution trajectory calculation are adjusted by a gradient iteration algorithm based on the deviation between the feedback information and the state transition prediction result, including:

[0121] The deviation between the feedback information and the state transition prediction result is calculated, and the deviation is decomposed according to the electric energy dimension, the thermal energy dimension and the aging dimension;

[0122] A coupling sensitivity matrix between the state dimensions is constructed according to the deviation decomposition result, the coupling sensitivity matrix is calculated in sections according to the correlation of the state changes in different time scales, and the gradient change direction between the dimensions is obtained;

[0123] The coupling weight coefficients are iteratively adjusted based on the gradient change direction, and the iteration step is determined by calculating the change rate of the state response between adjacent dimensions in each iteration process, and the adjustment of the weight coefficients is completed when the state deviation is less than a preset deviation threshold.

[0124] The deviation between the feedback information and the state transition prediction result is calculated by vector difference operation. The feedback information is represented by a three-dimensional vector, including the actual state value of the electrical energy dimension, the actual state value of the thermal energy dimension, and the actual state value of the aging dimension. The state transition prediction result is also represented by a three-dimensional vector, including the predicted state value of the corresponding three dimensions. The deviation vector is equal to the feedback information vector minus the prediction result vector, resulting in a three-dimensional deviation vector. The deviation calculation precision is maintained to six decimal places, with a numerical range of negative five to positive five. The deviation vector storage uses a floating-point array structure, and a single deviation vector occupies twenty-four bytes of memory. The deviation calculation frequency is synchronized with the feedback information update frequency, usually one to ten times per second.

[0125] The deviation is decomposed into the electrical energy dimension, the thermal energy dimension, and the aging dimension by vector component extraction. The electrical energy dimension deviation is equal to the first component of the deviation vector, reflecting the prediction error of the battery power state. The thermal energy dimension deviation is equal to the second component of the deviation vector, reflecting the prediction error of the battery temperature state. The aging dimension deviation is equal to the third component of the deviation vector, reflecting the prediction error of the battery health state. Each dimension deviation value is independently stored and processed, supporting individual analysis and adjustment. The deviation decomposition result is stored in a structure, including four fields: dimension identifier, deviation value, timestamp, and validity identifier. Abnormal deviation detection uses a statistical threshold method, and deviation values exceeding three times the standard deviation from the historical mean are marked as abnormal and trigger the abnormal handling process.

[0126] The coupling sensitivity matrix between state dimensions is constructed based on the correlation analysis of the deviation decomposition result. The coupling sensitivity matrix is a three-by-three symmetric matrix, and the matrix elements represent the correlation strength between different dimension deviations. The diagonal elements of the matrix are set to one, representing complete correlation of each dimension with itself. The non-diagonal elements are calculated by the Pearson correlation coefficient, reflecting the linear correlation degree of different dimension deviations. The correlation coefficient calculation uses a sliding window method, with a window length of one hundred to five hundred deviation samples. The matrix calculation precision is maintained to four decimal places, and the update frequency is once every fifty deviation samples processed. The matrix storage uses a two-dimensional array structure, occupying seventy-two bytes of memory.

[0127] The correlation calculation of state changes in each dimension at different time scales is achieved through multi-level correlation analysis. The time scale is divided into short-term, medium-term, and long-term, corresponding to time windows of one minute, ten minutes, and one hour. The covariance matrix of state changes between dimensions is calculated within each time scale. The covariance matrix reflects the synchronization and strength of state changes between different dimensions. The short-term correlation weight is set to zero point five, the medium-term correlation weight is set to zero point three, and the long-term correlation weight is set to zero point two. The comprehensive correlation is calculated by weighted average, with a weight sum of one to ensure normalization. The correlation calculation uses incremental update method to avoid repeated calculation of historical data and improve calculation efficiency.

[0128] The segmented calculation of the coupling sensitivity matrix is realized based on the weighted combination of the time scale correlation. The segmentation strategy determines the segmentation boundary according to the significant degree of correlation change. The time point when the correlation change exceeds zero point one is taken as the segmentation point. The coupling sensitivity matrix in each time period is calculated by linear combination of the corresponding time scale correlation matrix. The segmented matrix storage adopts a time index structure, which supports fast searching of the sensitivity matrix of a specific time period. The segmented calculation result is verified by the cross-validation method to ensure the numerical stability and prediction accuracy of the segmented matrix. The segmented matrix update adopts an inertial calculation strategy, which only performs calculation operations when needed.

[0129] The gradient change direction between each dimension is obtained by eigenvalue decomposition of the coupling sensitivity matrix. The eigenvalue decomposition decomposes the sensitivity matrix into eigenvalues and eigenvectors, and the eigenvectors represent the main direction of gradient change. The eigenvector corresponding to the largest eigenvalue is selected as the main gradient direction, and the eigenvector corresponding to the second largest eigenvalue is selected as the secondary gradient direction. The gradient direction is normalized to a unit vector to ensure the numerical stability of the direction calculation. The calculation precision of the gradient direction is kept to five decimal places, and the calculation frequency is synchronized with the sensitivity matrix update frequency. The gradient direction storage adopts a vector array structure, which supports parallel processing of multiple gradient directions.

[0130] The iterative adjustment of the coupling weight coefficient is realized based on the gradient descent algorithm. The initial coupling weight coefficient is set to equal weight, with the electrical energy dimension weight being zero point three three three, the thermal energy dimension weight being zero point three three three, and the aging dimension weight being zero point three four. The weight adjustment direction is equal to the negative gradient direction, ensuring adjustment in the direction of decreasing deviation. The weight adjustment amplitude is controlled by the iteration step, which is calculated based on the product of the gradient modulus and the learning rate. The learning rate is set to be in the range of zero point zero one to zero point one, and is dynamically adjusted according to the convergence speed and stability requirements. The weight coefficient constraint ensures that the sum of all weights is equal to one, and the range of a single weight value is zero to one.

[0131] The change rate of state response between adjacent dimensions is calculated to determine the iteration step. The state response is defined as the sensitivity of the state value of each dimension to the change of the weight coefficient, which is calculated by numerical differentiation method. The change rate is equal to the state response increment divided by the weight coefficient increment, and the increment size is set to one thousandth of the current weight value. The change rate between the electrical energy dimension and the thermal energy dimension, the change rate between the thermal energy dimension and the aging dimension, and the change rate between the aging dimension and the electrical energy dimension are calculated respectively. The central difference method is used for change rate calculation to improve numerical accuracy, and the average value of the forward difference and backward difference results is taken as the final change rate. The change rate anomaly detection is based on the historical statistical distribution, and the change rate exceeding the reasonable range is replaced by the historical median.

[0132] The iteration step length determination is based on the adaptive algorithm of the state response change rate between adjacent dimensions. The step length calculation formula is the learning rate divided by the absolute value of the change rate multiplied by the adjustment factor. The adjustment factor ranges from 0.5 to 2.0. The upper limit of the step length is set to 10% of the current weight value, and the lower limit is set to 1 / 1000 of the current weight value. The step length calculation includes a momentum term, and the momentum coefficient is set to 0.9 to accelerate convergence using historical gradient information. The step length is adaptively adjusted based on convergence speed monitoring. When the deviation decreases by less than a threshold value for five consecutive iterations, the step length is increased. When the deviation increases for three consecutive iterations, the step length is decreased. The step length storage uses a historical buffer to save the step length values of the last 20 iterations for trend analysis.

[0133] The iteration process is implemented through a loop structure, and each iteration includes four steps: gradient calculation, step length determination, weight update, and deviation evaluation. The gradient calculation is based on the current deviation decomposition result and the sensitivity matrix to perform matrix vector multiplication operation. The weight update uses vector addition, and the current weight vector is added to the step length multiplied by the negative gradient vector. After the weight update, normalization processing is performed to ensure that the weight sum is equal to one. The normalization method is to divide each weight value by the weight sum. The deviation evaluation is achieved by recalculating the deviation between the state prediction result after updating the weight and the feedback information. The iteration termination condition check is based on the comparison of the state deviation and the preset deviation threshold. When the deviation is less than the threshold, the iteration ends.

[0134] The state deviation is compared with the preset deviation threshold using the Euclidean distance metric. The state deviation is equal to the square root of the sum of the squares of the components of the three-dimensional deviation vector, reflecting the overall amplitude of the deviation. The preset deviation threshold is determined according to the system control accuracy requirement, and is usually set to the range of 0.01 to 0.1. The deviation comparison accuracy is maintained to four decimal places, and the comparison frequency is synchronized with the iteration frequency. The deviation history record uses a ring buffer to store the deviation values of the last 1000 iterations for convergence analysis and parameter tuning. The deviation threshold is adaptively adjusted based on the convergence behavior analysis. When oscillatory convergence occurs, the threshold is appropriately relaxed. When monotonic convergence occurs, the threshold is appropriately tightened.

[0135] The determination of the weight coefficient adjustment is based on the achievement of the deviation threshold and the verification of the convergence stability. The achievement of the deviation threshold requires that the state deviation of five consecutive iterations is less than the preset threshold. The verification of the convergence stability requires that the change amplitude of the weight coefficient in the last ten iterations is less than 1 / 1000. The weight coefficient adjustment result storage includes four information: the optimal weight value, the iteration number, the convergence deviation, and the adjustment time consumption. The adjustment process monitoring record includes the weight value, the deviation value, the step length, and the gradient modulus of each iteration, supporting subsequent analysis and debugging. The abnormal termination processing includes the maximum iteration number limit and the divergence detection. When the iteration number exceeds 1000 or the deviation increases continuously for 10 times, the process is forcibly terminated and rolled back to the optimal historical state.

[0136] In a second aspect of the embodiment of the present application, an electronic device is provided, comprising:

[0137] a processor;

[0138] a memory for storing processor-executable instructions;

[0139] wherein the processor is configured to invoke the instructions stored by the memory to perform the method as described above.

[0140] In a third aspect, the present application provides a computer readable storage medium having stored thereon computer program instructions, which when executed by a processor, implement the method as described above.

[0141] The present application can be a method, apparatus, system, and / or computer program product. Computer program products can include computer-readable storage media having computer-readable program instructions loaded thereon for performing various aspects of the present application.

[0142] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for switching operation mode of energy storage power station based on super capacitor and intelligent algorithm, characterized in that, The method comprises: acquiring real-time state parameters of a super capacitor in an energy storage power station and load demand information of a power grid side, based on the real-time state parameters, constructing a coupled state space of the super capacitor in an electric energy dimension, a thermal energy dimension and an aging dimension, calculating state evolution trajectories of the super capacitor along the coupled state space in different operating modes, and obtaining a state transition prediction result; determining a candidate operating mode sequence of the energy storage power station within a preset time window according to the state transition prediction result and the load demand information of the power grid side, searching within the preset time window through a dynamic programming framework according to the candidate operating mode sequence, combining mode switching time points at which a super capacitor state trajectory and a power grid demand trajectory have a minimum deviation, introducing a long-term performance degradation penalty term calculated based on a trajectory feature vector into a cost function of the dynamic programming framework, and generating a switching strategy; controlling the energy storage power station to perform mode switching actions based on the switching strategy, collecting actual state trajectory data after switching as feedback information, adjusting coupling weight coefficients between the electric energy dimension, the thermal energy dimension and the aging dimension in state evolution trajectory calculation through a gradient iteration algorithm based on a deviation between the feedback information and the state transition prediction result, and returning the adjusted coupling weight coefficients to a state analysis step for updating subsequent state evolution trajectory calculation.

2. The method of claim 1, wherein, The method comprises: mapping the real-time state parameters to the coupled state space to obtain a current state point comprising an electric energy dimension state component, a thermal energy dimension state component and an aging dimension state component, and constructing a local linearization region in the coupled state space based on the current state point, obtaining coupling weight coefficients by calculating a partial derivative matrix of each dimension state component with respect to other dimension state components in the local linearization region; establishing a state increment calculation rule based on the coupling weight coefficients, the state increment calculation rule converting the electric energy dimension state component at a current time into a thermal energy dimension state increment through a first coupling weight coefficient, converting the thermal energy dimension state increment into an aging dimension state increment through a second coupling weight coefficient, and iteratively deducing evolution processes of each dimension state component in a time direction from the current state point to form a state evolution trajectory comprising time sequence correlation by using the state increment calculation rule.

3. The method of claim 2, wherein, The method comprises: mapping the real-time state parameters to the coupled state space to obtain a current state point comprising an electric energy dimension state component, a thermal energy dimension state component and an aging dimension state component, and constructing a local linearization region in the coupled state space based on the current state point, obtaining coupling weight coefficients by calculating a partial derivative matrix of each dimension state component with respect to other dimension state components in the local linearization region; constructing an adaptive local linearization region based on the current state point, calculating a rate of change of the current state point in each dimension, determining a linearization radius of each dimension according to the rate of change, and constructing an ellipsoid-shaped local linearization region with the linearization radius as a principal axis; within the ellipsoid-shaped local linearization region, generating a plurality of sampling point groups by using a stratified importance sampling method; performing recursive principal component analysis on the plurality of sampling point groups, constructing a nonlinear coupling feature between state components by constructing a high-order matrix decomposition, and performing dimensionality reduction reconstruction on the high-order matrix by using singular value decomposition; based on the reconstructed feature, obtaining dynamic coupling weight coefficients by calculating conditional partial derivatives of state components in each dimension to state components in other dimensions under different time scales.

4. The method of claim 1, wherein, According to the candidate operation mode sequence, a dynamic programming framework is used to search within the preset time window, so that the combination of mode switching time points with the minimum deviation degree of the super capacitor state trajectory and the power grid demand trajectory includes: Based on the candidate operation mode sequence, a state transition graph is constructed, a stratified beam search method is used to prune the state transition graph, and a preset number of paths with the minimum deviation in each layer are retained; The dynamic programming search is applied to the pruned state transition graph, the preset time window is divided into a plurality of time segments, and the cumulative deviation of all feasible state transitions in each time segment is calculated; based on the cumulative deviation, an optimal state transition matrix is constructed, and the combination of operation mode switching time points that minimizes the deviation degree of the super capacitor state trajectory and the power grid demand trajectory is determined by backtracking the optimal state transition matrix.

5. The method of claim 4, wherein, Based on the candidate operation mode sequence, a state transition graph is constructed, a stratified beam search method is used to prune the state transition graph, and a preset number of paths with the minimum deviation in each layer are retained: The state transition graph is divided into a plurality of search layers in chronological order, the deviation value of each node from the target state is calculated, and a priority queue is established for each node based on the deviation value; The stratified beam search is performed on the nodes in the priority queue, in each layer of search process, a local state transition subgraph is constructed based on the state transition gradient between each node and its adjacent node, the path connectivity between nodes is calculated according to the local state transition subgraph, and the weighted combination of the path connectivity and the node deviation value is taken as the node score, a preset number of nodes with the optimal score are retained as the retained nodes of the current layer; the path of the retained node in the state transition graph is taken as a candidate path, the cumulative deviation of the candidate path is calculated, and the candidate path with a cumulative deviation exceeding a dynamic threshold is eliminated.

6. The method of claim 1, wherein, In the cost function of the dynamic programming framework, a long-term performance degradation penalty term calculated based on the trajectory feature vector is introduced, and the switching strategy is generated. The performance degradation features are calculated in a sliding time window based on a trajectory feature vector, a change rate of the trajectory feature vectors in adjacent time windows is extracted, and a degradation speed of each performance index is calculated according to the change rate; the degradation speed is piecewise linearly fitted in multiple time scales to obtain a main direction component, and a projection value of each performance index on the main direction component is taken as a corresponding weight coefficient to obtain a comprehensive performance degradation rate, and a long-term performance degradation constraint is constructed based on the comprehensive performance degradation rate; The long-term performance degradation constraint and the state transition cost are combined to obtain a stage cost calculation mode; for each decision stage, the cumulative cost of different switching schemes is determined based on the stage cost calculation mode; the optimal switching sequence is solved through a recursive step, so that the sum of the cumulative costs is minimized, and a switching strategy is generated according to the optimal switching sequence.

7. The method of claim 1, wherein, The coupling weight coefficients between the electrical energy dimension, the thermal energy dimension and the aging dimension in the state evolution trajectory calculation are adjusted by a gradient iteration algorithm through the deviation between the feedback information and the state transition prediction result, including: The deviation between the feedback information and the state transition prediction result is calculated, and the deviation is decomposed according to the electrical energy dimension, the thermal energy dimension and the aging dimension; A coupling sensitivity matrix between the state dimensions is constructed for the deviation decomposition result, the coupling sensitivity matrix is calculated in sections according to the correlation of the state changes in each dimension in different time scales to obtain the gradient change direction between each dimension; The coupling weight coefficients are iteratively adjusted based on the gradient change direction, and the iteration step is determined by calculating the change rate of the state response between adjacent dimensions in each iteration process, and the adjustment of the weight coefficients is completed when the state deviation is less than a preset deviation threshold.

8. An energy storage plant operating mode switching system based on supercapacitor and intelligent algorithm for implementing the method of any one of the preceding claims 1-7, characterized in that, It includes: The first unit is configured to obtain real-time state parameters of a super capacitor in an energy storage power station and load demand information of a grid side, calculate state evolution trajectories of the super capacitor along a coupling state space of the super capacitor in the electrical energy dimension, the thermal energy dimension and the aging dimension in different operating modes based on the real-time state parameters by constructing the coupling state space, and obtain a state transition prediction result; The second unit is configured to determine a candidate operating mode sequence of the energy storage power station in a preset time window according to the state transition prediction result and the load demand information of the grid side, search in the preset time window through a dynamic programming framework according to the candidate operating mode sequence, combine mode switching time points of a mode in which a super capacitor state trajectory and a grid demand trajectory have a minimum deviation, and introduce a long-term performance degradation penalty term calculated based on the trajectory feature vector into a cost function of the dynamic programming framework to generate a switching strategy; The third unit is configured to control the energy storage power station to perform a mode switching action based on the switching strategy, collect actual state trajectory data after the switching as feedback information, adjust coupling weight coefficients between the electrical energy dimension, the thermal energy dimension and the aging dimension in the state evolution trajectory calculation through a gradient iteration algorithm through the deviation between the feedback information and the state transition prediction result, and return the adjusted coupling weight coefficients to a state analysis step for updating subsequent state evolution trajectory calculation.

9. An electronic device, comprising: It includes: A processor; a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored by the memory to perform the method of any one of claims 1 to 7.

10. A computer-readable storage medium having stored thereon computer program instructions, wherein, The computer program instructions, when executed by a processor, implement the method of any one of claims 1 to 7.

Images