Virtual power plant demand side response and load management method and system

Abstract

The invention provides a virtual power plant demand side response and load management method and system, and relates to the technical field of power system scheduling, and the method comprises the steps: obtaining a demand side response request and a historical response track of each distributed energy source; performing time-frequency domain decomposition on the historical response trajectory, extracting a transient component and a steady component, constructing a dynamic response model containing the current output power and the load deviation cumulant, and calculating load deviation distribution through Kalman filtering; a Markov decision process is constructed based on load deviation distribution, an optimal power distribution strategy is solved by adopting a value iteration algorithm, and a response execution scheme is generated; collecting a power change sequence after execution of each distributed energy source, converting the power change sequence into a load recovery demand vector, constructing a recovery capacity constraint cone based on the current available adjustment margin of the power grid, and solving an optimal projection point; and performing hierarchical division according to the optimal projection point to obtain a staged recovery sequence. According to the method, the accuracy of virtual power plant demand response and the coordination of load recovery are improved, and the power grid regulation pressure is reduced.

Description

Technical Field

[0001] This invention relates to the field of power system dispatching technology, and in particular to a virtual power plant demand-side response and load management method and system. Background Technology

[0002] With the rapid development of the energy internet and smart grids, virtual power plants (VPS), as a coordinated control technology that aggregates multiple distributed energy sources, are playing an increasingly important role in demand-side response of power systems. VPS integrates various types of distributed energy sources, such as distributed photovoltaics, energy storage systems, and adjustable loads, to form a dispatchable resource pool, participating in grid peak shaving, frequency regulation, and ancillary services. This is significant for improving grid flexibility and renewable energy absorption capacity. In actual operation, VPS needs to respond to grid regulation demands by coordinating and controlling the output power of each distributed energy source to achieve rapid tracking of target regulation power. Simultaneously, after demand response ends, each distributed energy source needs to return to normal operation. This process requires consideration of grid carrying capacity and the recovery priority of each load to avoid causing new grid shocks. Therefore, how to achieve precise demand response control and orderly load recovery management in VPS is a pressing technical problem that needs to be solved. Summary of the Invention

[0003] The embodiments of the present invention provide a method and system for demand-side response and load management of virtual power plants, which can solve the problems in the prior art.

[0004] A first aspect of the present invention provides a virtual power plant demand-side response and load management method, comprising: Obtain the demand-side response requests and historical response trajectories of each distributed energy source, wherein the demand-side response requests include the target adjustment power and the response duration; The historical response trajectory is decomposed in the time-frequency domain to extract transient and steady-state components, and a dynamic response model containing the current output power and the cumulative load deviation is constructed. Based on the dynamic response model, the target adjustment power is used as the control input, and the load deviation distribution at the end of the response duration is calculated by Kalman filtering. A Markov decision process is constructed based on the load deviation distribution. The cumulative load deviation is defined as the decision state set, and the adjustable capacity constraint is defined as the decision action set. The reward function is calculated, and the optimal power allocation strategy is solved using a value iteration algorithm. A response execution plan is generated and issued. Collect the power change sequence of each distributed energy source after executing the response execution scheme, convert the power change sequence into a load recovery demand vector, construct a recovery capacity constraint cone with vertices based on the current available adjustment margin of the power grid, and solve for the optimal projection point of each load recovery demand vector in the constraint cone; The system is divided into layers based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, and then into multiple recovery stages executed sequentially to obtain a phased recovery sequence.

[0005] The historical response trajectory is decomposed in the time-frequency domain to extract transient and steady-state components, and a dynamic response model including the current output power and the cumulative load deviation is constructed, including: Perform wavelet transform decomposition on the historical response trajectory to obtain a multi-scale time-frequency coefficient matrix, calculate the energy distribution of the coefficients at each scale, determine the transient interval, reconstruct the time-frequency coefficients by inverse wavelet transform to obtain the transient component, and reconstruct the time-frequency coefficients outside the transient interval by inverse wavelet transform to obtain the steady-state component. Calculate the transient response time constant of the transient component and establish the transient evolution equation of the current output power; calculate the power deviation sequence of the steady-state component and perform time integration to obtain the evolution data of the load deviation accumulation, fit the accumulation rate function and calculate the steady-state accumulation coefficient, and establish the steady-state evolution equation of the load deviation accumulation. The current output power and the cumulative load deviation are constructed into a state vector. The transient response time constant and the steady-state cumulative coefficient are embedded into the state transition matrix. The transient evolution equation and the steady-state evolution equation are discretized and combined into a dynamic response model in state space form.

[0006] Based on the dynamic response model, using the target regulating power as the control input, the load deviation distribution at the end of the response duration is calculated using Kalman filtering, including: The target regulation power is discretized according to the response duration to generate a control input sequence, which is substituted into the dynamic response model as the control vector. A process noise covariance matrix with a block diagonal structure is constructed, and the time update step of Kalman filtering is performed to calculate the prior state estimate and prior error covariance. Collect real-time power measurements of each distributed energy source, perform measurement update steps and calculate Kalman gain, correct the prior state estimate, and obtain the posterior state estimate and posterior error covariance. Based on the posterior error covariance, calculate the covariance submatrix corresponding to the cumulative load deviation, perform eigenvalue decomposition on the covariance submatrix, identify the eigenvector corresponding to the direction of the cumulative load deviation, construct a multidimensional confidence ellipsoid along the direction of the eigenvector, use its projection boundary on the dimension of the cumulative load deviation as the upper and lower bounds of the distribution, and combine it with the mean of the cumulative load deviation in the posterior state estimation to construct the load deviation distribution.

[0007] A Markov decision process is constructed based on the load deviation distribution. The cumulative load deviation is defined as the decision state set, and the adjustable capacity constraint is defined as the decision action set. The reward function is calculated, and a value iteration algorithm is used to solve for the optimal power allocation strategy. A response execution plan is generated and issued, including: Based on the load deviation distribution, the expected value and variance of each distributed energy source are extracted, and a discretized state set of the cumulative load deviation is generated according to a preset step size. The Cartesian product of these states is then combined to construct the decision state set of the Markov decision process. Based on the adjustable capacity constraints of each distributed energy source, the target adjustable power is decomposed into power allocation increments, a feasible region that satisfies the total power balance and capacity constraints is constructed, and a set of decision actions is generated through grid sampling. Based on the superposition relationship between the cumulative load deviation and the power allocation increment under each decision state, the state transition probability matrix is ​​calculated using the variance of the deviation distribution, and a reward function is constructed that includes the square term of the adjustment power deviation, the weighted sum of squares of the load deviation, and the penalty term of the power allocation change rate. Based on the state transition probability and reward function, a value iteration algorithm is used to update the value function until convergence. According to the optimal value function, the power allocation increment scheme with the largest expected cumulative reward is selected and superimposed on the output power of each distributed energy source at the start of the response, which is then used as the response execution scheme and issued.

[0008] The power change sequence is converted into a load recovery demand vector. A recovery capacity constraint cone with vertices is constructed based on the current available adjustment margin of the power grid. The optimal projection point of each load recovery demand vector within the constraint cone is then solved, including: The peak deviation of each distributed energy source is calculated based on the power change sequence, and the average recovery power demand is calculated in combination with the preset recovery time to generate a load recovery demand vector. The electrical distance between nodes is calculated based on the voltage measurement values ​​of the power grid bus nodes. Multiple voltage control zones are divided, and the reactive and active reserve capacity of each zone is obtained. The reactive reserve capacity is converted into the maximum load recovery power to form a zone reserve capacity vector. The remaining transmission margin is calculated based on the current transmission power and thermal stability limit of the tie line. A cross-zone power transmission constraint matrix is ​​established, and the negative values ​​of the zone reserve capacity vectors are used as vertices and the row vectors of the transmission constraint matrix are used as the cone surface normal vectors to generate the recovery capacity constraint cone. Calculate the line vector connecting each load recovery demand vector to the vertex of the recovery capacity constraint cone. Determine whether the load recovery demand vector is located inside the constraint cone by judging the sign of the inner product of the connecting line vector and the cone surface normal vector. For a load recovery demand vector located outside the constraint cone, an optimization problem is constructed with the square of the distance from the projection point to the original vector as the objective function. The optimal projection point is calculated by solving the Lagrange dual problem.

[0009] Based on the current transmission power and thermal stability limit of the tie line, the remaining transmission margin is calculated, and a cross-regional power transmission constraint matrix is ​​established. Using the negative values ​​of the regional reserve capacity vectors as vertices and the row vectors of the transmission constraint matrix as the cone normal vectors, a recovery capacity constraint cone is generated, including: Calculate the remaining transmission margin of each tie line based on the real-time transmission power and corresponding thermal stability limit of the tie lines between each voltage control zone. Based on the remaining transmission margin and the topological connection relationship of each voltage control zone, a cross-zone power transmission constraint matrix is ​​constructed to characterize the directional constraints of power flow between zones. Each row of the cross-zone power transmission constraint matrix corresponds to the power transmission boundary condition of a tie line, and the absolute value of the matrix element represents the power transmission capacity coefficient of the corresponding zone on the tie line. The negative values ​​of each component of the partition spare capacity vector are used to obtain the constraint cone vertex, which represents the recovery capability benchmark of each partition when there is no cross-partition power support. Extract each row vector of the cross-regional power transmission constraint matrix as a cone normal vector. Each cone normal vector is perpendicular to the hyperplane passing through the vertex of the constraint cone. Define a cone boundary. The intersection of multiple cone boundaries forms a cone space, resulting in a recovery capacity constraint cone with the vertex of the constraint cone as the cone apex.

[0010] Based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, the system is layered and divided into multiple sequentially executed recovery stages, resulting in a phased recovery sequence, including: Calculate the Euclidean distance between each optimal projection point and the vertex of the recovery capacity constraint cone. Combined with the geometric topology of the recovery capacity constraint cone, calculate the curvature distribution of the constraint cone section where each optimal projection point is located. Then, divide the recovery capacity constraint cone radially into multiple nested curvature isosurface regions to obtain multiple recovery levels. For the curvature isosurface regions of each recovery level, a tangential recovery manifold is constructed. Based on the optimal projection point, the gradient direction of each distributed energy source on the tangential recovery manifold and the coupling strength with the radial contraction direction of the constraint cone are calculated. The distributed energy sources are sorted according to the coupling strength to determine the recovery execution order of each recovery level. The recovery levels are arranged in order from the outside to the inside according to the curvature isosurface, and the distributed energy sources arranged in the recovery execution order are connected in series to generate the phased recovery sequence. A recovery instruction corresponding to the phased recovery sequence is issued to each distributed energy source, and the recovery instruction carries a corresponding execution timing identifier.

[0011] A second aspect of the present invention provides a virtual power plant demand-side response and load management system, comprising: The data acquisition unit is used to acquire the demand-side response requests and historical response trajectories of each distributed energy source, wherein the demand-side response requests include the target adjustment power and the response duration. The load deviation calculation unit is used to perform time-frequency domain decomposition on the historical response trajectory, extract transient and steady-state components, and construct a dynamic response model that includes the current output power and the cumulative load deviation. Based on the dynamic response model, the target adjustment power is used as the control input, and the load deviation distribution at the end of the response duration is calculated by Kalman filtering. The decision optimization unit is used to construct a Markov decision process based on the load deviation distribution, define the combination of load deviation accumulation as a decision state set, combine it with the adjustable capacity constraint as a decision action set, calculate the reward function, and use the value iteration algorithm to solve the optimal power allocation strategy, generate a response execution plan and issue it. The load recovery analysis unit is used to collect the power change sequence of each distributed energy source after executing the response execution scheme, convert the power change sequence into a load recovery demand vector, construct a recovery capacity constraint cone with vertices based on the current available adjustment margin of the power grid, and solve for the optimal projection point of each load recovery demand vector in the constraint cone. The recovery scheduling unit is used to divide the recovery into layers based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, and to divide the recovery into multiple recovery stages that are executed sequentially, thus obtaining a phased recovery sequence.

[0012] A third aspect of the present invention, An electronic device is provided, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0013] Fourth aspect of the embodiments of the present invention, A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0014] The beneficial effects of this application are as follows: This invention acquires the demand-side response requests and historical response trajectories of each distributed energy source, performs time-frequency domain decomposition on the historical response trajectory to extract transient and steady-state components, and constructs a dynamic response model that includes the current output power and the cumulative load deviation. This model can accurately characterize the dynamic response characteristics of distributed energy sources, improve the accuracy of response behavior modeling, provide a reliable data foundation for subsequent power allocation decisions, and effectively enhance the accuracy and controllability of the virtual power plant's demand-side response.

[0015] This invention calculates load deviation distribution based on a dynamic response model using Kalman filtering, and constructs a Markov decision process to define the cumulative load deviation as a decision state set and the adjustable capacity constraint as a decision action set. It uses a value iteration algorithm to solve for the optimal power allocation strategy, achieving optimal decision-making under uncertain conditions. It can maximize the response effect while satisfying various constraints, significantly improving the load management efficiency and economy of virtual power plants.

[0016] This invention converts the power change sequence after the distributed energy source executes the response into a load restoration demand vector. Based on the current available adjustment margin of the power grid, it constructs a restoration capacity constraint cone with vertices to solve for the optimal projection point. Then, based on the distance between the optimal projection point and the vertex of the constraint cone, it divides the restoration into multiple restoration stages that are executed sequentially, realizing phased and orderly restoration. This avoids the impact on the power grid caused by simultaneous load restoration and ensures the safety and stability of power grid operation. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the virtual power plant demand-side response and load management method according to an embodiment of the present invention. Figure 2 A schematic diagram of the method for extracting transient and steady-state components and constructing a dynamic response model. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0020] Figure 1 This is a flowchart illustrating the virtual power plant demand-side response and load management method according to an embodiment of the present invention. Figure 1 As shown, the method includes: Obtain the demand-side response requests and historical response trajectories of each distributed energy source, wherein the demand-side response requests include the target adjustment power and the response duration; The historical response trajectory is decomposed in the time-frequency domain to extract transient and steady-state components, and a dynamic response model containing the current output power and the cumulative load deviation is constructed. Based on the dynamic response model, the target adjustment power is used as the control input, and the load deviation distribution at the end of the response duration is calculated by Kalman filtering. A Markov decision process is constructed based on the load deviation distribution. The cumulative load deviation is defined as the decision state set, and the adjustable capacity constraint is defined as the decision action set. The reward function is calculated, and the optimal power allocation strategy is solved using a value iteration algorithm. A response execution plan is generated and issued. Collect the power change sequence of each distributed energy source after executing the response execution scheme, convert the power change sequence into a load recovery demand vector, construct a recovery capacity constraint cone with vertices based on the current available adjustment margin of the power grid, and solve for the optimal projection point of each load recovery demand vector in the constraint cone; The system is divided into layers based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, and then into multiple recovery stages executed sequentially to obtain a phased recovery sequence.

[0021] Figure 2 This is a flowchart illustrating a method for extracting transient and steady-state components and constructing a dynamic response model. In one optional implementation, the historical response trajectory is decomposed in the time-frequency domain to extract transient and steady-state components, and a dynamic response model incorporating the current output power and the cumulative load deviation is constructed, including: Perform wavelet transform decomposition on the historical response trajectory to obtain a multi-scale time-frequency coefficient matrix, calculate the energy distribution of the coefficients at each scale, determine the transient interval, reconstruct the time-frequency coefficients by inverse wavelet transform to obtain the transient component, and reconstruct the time-frequency coefficients outside the transient interval by inverse wavelet transform to obtain the steady-state component. Calculate the transient response time constant of the transient component and establish the transient evolution equation of the current output power; calculate the power deviation sequence of the steady-state component and perform time integration to obtain the evolution data of the load deviation accumulation, fit the accumulation rate function and calculate the steady-state accumulation coefficient, and establish the steady-state evolution equation of the load deviation accumulation. The current output power and the cumulative load deviation are constructed into a state vector. The transient response time constant and the steady-state cumulative coefficient are embedded into the state transition matrix. The transient evolution equation and the steady-state evolution equation are discretized and combined into a dynamic response model in state space form.

[0022] When performing time-frequency domain decomposition of historical response trajectories, it is necessary to obtain power output trajectory data within a certain time window in the past. Taking an industrial load control scenario as an example, power output data recorded once per second over the past 24 hours is collected to form a data set containing 8.64 × 10⁻⁶ ohms. 4 The time series of sampling points records the complete process of the load gradually recovering from the initial state after experiencing a disturbance. The data range fluctuates from an initial power value of 150kW to a peak of 180kW and then gradually falls back to a stable value of 155kW.

[0023] A wavelet transform decomposition operation was performed on the historical response trajectory. Wavelet basis functions with good time-frequency localization characteristics were selected for multi-level decomposition, decomposing the original signal to level 8, with each level corresponding to a different frequency scale. During the decomposition process, each scale level generates a set of time-frequency coefficients, which reflect the energy distribution characteristics of the signal within a specific frequency range and time location. The high-frequency coefficients of levels 1 to 3 mainly capture the rapidly changing components in the signal, the mid-frequency coefficients of levels 4 to 6 reflect the moderately rapid power adjustment process, and the low-frequency coefficients of levels 7 and 8, as well as the final approximation coefficients, contain the slow changing trend of the signal.

[0024] When calculating the energy distribution of time-frequency coefficients at each scale, the sum of the squares of each element in the coefficient sequence of each layer is taken to obtain the total energy value of that scale layer. For a layer containing 8.64 × 10⁻⁶, the total energy value is calculated as follows: 4 The signal from the original sampling points, after the first layer decomposition, yields 4.32 × 10⁻⁶. 4 The energy of the detailed coefficients is calculated to be 2.85 × 10⁻⁶. 3 Unit; the second layer produces 2.16 × 10 4 The coefficient has an energy of 1.92 × 10⁻⁶. 3 Unit; Layer 3 produces 1.08 × 10 4 The coefficient has an energy of 1.45 × 10⁻⁶. 3 Unit, continuing to break down, the energy values ​​of layers 4 to 8 are 9.80 × 10⁸ respectively. 2 Unit, 6.20 × 10 2 Unit, 3.80×10 2 Unit, 2.10 × 10 2 Unit and 1.50×10 2 Units. By analyzing the energy distribution pattern with scale, it was found that the high-frequency coefficients in the first three layers account for approximately 42% of the total signal energy. These high-frequency components correspond to the rapid response stage after being disturbed.

[0025] When determining the transient interval, the amplitude changes of high-frequency coefficients are analyzed along the time axis. The absolute values ​​of the amplitudes of coefficients from layer 1 to layer 3 at corresponding time points are superimposed to form a comprehensive high-frequency activity intensity curve. This curve shows significant high-amplitude fluctuations in the first 600 seconds after the disturbance, with an amplitude range between 0.8 and 2.3. After 600 seconds, it rapidly decays to below 0.3 and remains stable. A threshold of 1.5 times the average amplitude, i.e., 0.45, is set, and the time period exceeding this threshold continuously is defined as the transient interval. After judgment, the period from the start of the disturbance to 720 seconds is defined as the transient response phase, during which the system exhibits obvious dynamic adjustment characteristics.

[0026] The time-frequency coefficients corresponding to the transient interval are extracted and reconstructed. Specifically, all detailed coefficients from layers 1 to 3 within the time range of 0 to 720 seconds are retained, while coefficients outside this time range are set to 0. Meanwhile, coefficients from layers 4 to 8 that contribute significantly to energy are retained, and coefficients with a cumulative contribution rate of over 85% are selected based on energy percentage. Wavelet inverse transform reconstruction is performed on these selected coefficients with time windows, reconstructing the signal layer by layer upwards. The reconstructed transient component time series shows that the power rapidly rises to a peak of 180kW in the early stages of the disturbance, and then gradually decreases to around 160kW within 720 seconds in an oscillating decay manner. The oscillation period is approximately 80 seconds, and the amplitude decreases by about 30% each time.

[0027] The time-frequency coefficients outside the transient interval are processed to extract the steady-state components. All scale layer coefficients after 720s are retained, and the low-frequency approximation coefficients of the 7th and 8th layers are also retained for the entire time range. Wavelet inverse transform is performed on these coefficients to reconstruct the steady-state component time series. This series shows a slow monotonic change trend, gradually adjusting from 160kW at 720s to 155kW at the end of 24h. The whole process does not have obvious oscillation characteristics and shows a smooth and gradual convergence.

[0028] When calculating the transient response time constant of the transient component, an exponential decay fitting method is used. The transient component data is represented as an exponential function of the initial amplitude multiplied by the natural logarithm base, with the exponential term including the ratio of time to the time constant. The transient component curve is fitted using the least squares method, and the time constant parameter that minimizes the sum of squared fitting errors is obtained. In this case, the fitted transient response time constant is 128 s, representing the characteristic time required for the transient response to decay to approximately 36.8% of its initial value. This time constant reflects the response rate of the fast adjustment mechanism; a smaller value indicates a faster transient adjustment of the system.

[0029] When establishing the transient evolution equation for the current output power, the power change rate is described as the negative of the ratio of the current power deviation from the target value to the transient response time constant. This equation shows that the speed at which the power approaches the target value is proportional to the current deviation. The proportionality coefficient is determined by the reciprocal of the time constant. Taking a time constant of 128s as an example, the proportionality coefficient is 7.8×10^-3, which means that the power adjustment per second is about 0.78% of the current deviation. When the power deviates from the target value by 10kW, the adjustment speed per second is about 78W.

[0030] When calculating the power deviation sequence of the steady-state component, the final stable value of 155kW is subtracted from the value of each sampling point in the steady-state component to obtain the power deviation at each time point. The deviation is 5kW at 720s, decreases to 3.2kW at 3600s, and reaches 1.08×10 4 The deviation at time s is 1.8 kW, down to 8.64 × 10 4 The deviation at time s is close to 0. A time integration operation is performed on this deviation sequence, using the trapezoidal integration method to accumulate the product of the deviation and the time interval in each time period. Integration begins at 720s, and the accumulated amount reaches 1.26 × 10⁻⁶ by 3600s. 4 kW·s, up to 1.08×10 4 The cumulative amount of s reached 6.84 × 10 4 kW·s, up to 8.64×10 4 The cumulative amount of s reached 2.28 × 10 5 kW·s, this cumulative amount reflects the overall degree of long-term deviation from the target state.

[0031] When fitting the cumulative rate function, the growth pattern of the cumulative load deviation over time was analyzed. The cumulative rate was obtained by differentiating the cumulative amount with respect to time, revealing a gradually decreasing trend. A linear fitting method was used, expressing the cumulative rate as the initial rate minus the product of the attenuation coefficient and time. Data fitting yielded an initial cumulative rate of 4.8 kW and an attenuation coefficient of 5.6 × 10⁻⁵ kW / s. The steady-state cumulative coefficient, defined as the ratio of the cumulative rate to the power deviation, was calculated to be 0.96, representing the cumulative increase per second per unit power deviation.

[0032] When establishing the steady-state evolution equation for the cumulative load deviation, the rate of change of the cumulative load is described as the product of the current power deviation and the steady-state cumulative coefficient. This equation shows that the growth rate of the cumulative load is proportional to the current power deviation of the system. When the power deviation is 3kW, the cumulative load increases by 2.88kW·s per second; when the deviation drops to 1kW, the cumulative load increases by only 0.96kW·s per second.

[0033] When constructing the state vector, the current output power and the cumulative load deviation are used as two state variables to form a two-dimensional vector. At a certain moment, the first component of the state vector is the power value of 162kW, and the second component is the cumulative deviation of 4.5×10. 4 kW·s. The state transition matrix describes the evolution of the state vector per unit time. The matrix embeds the reciprocal of the transient response time constant and the steady-state cumulative coefficient. The element in the first row and first column of the matrix is ​​the negative value of the reciprocal of the transient time constant, i.e., -7.8 × 10⁻⁶. -3 The element in the first column of the second row is the steady-state cumulative coefficient of 0.96, and the remaining elements are determined based on physical properties.

[0034] When discretizing the transient evolution equation and the steady-state evolution equation, the forward difference method is used to approximate the continuous time derivative as the ratio of the difference between state quantities at adjacent time steps to the time step size. If the discrete time step size is set to 1 second, then the power value at the next time step is equal to the current power value minus the product of the time step size, the reciprocal of the transient time constant, and the power deviation. The cumulative amount at the next time step is equal to the current cumulative amount plus the product of the time step size, the steady-state cumulative coefficient, and the power deviation.

[0035] When combined into a state-space dynamic response model, the discretized evolution equations are written in matrix-vector form. The state vector at the next time step is equal to the product of the state transition matrix and the current state vector, plus the input vector. This model fully describes the evolution trajectory of the state variables over time from any initial state under a given input. It can predict the output power and cumulative deviation at any future time, providing a dynamic prediction basis for load control strategy formulation. The model parameters are extracted from historical data and have clear physical meaning and numerical interpretability.

[0036] In one optional implementation, based on the dynamic response model, the target regulating power is used as the control input, and the load deviation distribution at the end of the response duration is calculated using Kalman filtering, including: The target regulation power is discretized according to the response duration to generate a control input sequence, which is substituted into the dynamic response model as the control vector. A process noise covariance matrix with a block diagonal structure is constructed, and the time update step of Kalman filtering is performed to calculate the prior state estimate and prior error covariance. Collect real-time power measurements of each distributed energy source, perform measurement update steps and calculate Kalman gain, correct the prior state estimate, and obtain the posterior state estimate and posterior error covariance. Based on the posterior error covariance, calculate the covariance submatrix corresponding to the cumulative load deviation, perform eigenvalue decomposition on the covariance submatrix, identify the eigenvector corresponding to the direction of the cumulative load deviation, construct a multidimensional confidence ellipsoid along the direction of the eigenvector, use its projection boundary on the dimension of the cumulative load deviation as the upper and lower bounds of the distribution, and combine it with the mean of the cumulative load deviation in the posterior state estimation to construct the load deviation distribution.

[0037] In this specific embodiment, the target regulating power input is used as the control variable and processed by the Kalman filter framework. Assuming the received target regulating power is 800kW and the response duration is 600s, the continuous control command needs to be converted into a discrete-time sequence. A time step of 10s is selected as the discretization interval, dividing the 600s duration into 60 time steps. Within each time step, the control input remains constant at 800kW, forming a control input sequence of length 60. This sequence is directly substituted into the established dynamic response model as a control vector, with each element of the control vector corresponding to the power regulation command for one time step.

[0038] The dynamic response model adopts a state-space representation. The state vector contains variables such as the response power, response rate, and cumulative load deviation of each distributed energy source. The dimension of the state vector depends on the number of distributed energy sources involved in regulation. Assuming there are 15 distributed energy devices, the state of each device is described by two variables: response power and response rate. In addition, the global cumulative load deviation is added, and the total dimension of the state vector is 31. The control vector is associated with the state vector through the state transition matrix. The state transition matrix describes how the control input drives the dynamic evolution of the power response of each device.

[0039] A process noise covariance matrix is ​​constructed to characterize the uncertainties in the dynamic evolution of the system. This matrix adopts a block-diagonal structure to separately describe the noise characteristics of different types of state variables. For the response power variables of 15 distributed energy devices, their process noise covariance is set as a diagonal block, and the power noise standard deviation of each device is set to 2% of the actual power. Specifically, if the current response power of a device is 50kW, the corresponding process noise variance is 1kW. 2 The process noise covariance of the response rate variable is also a diagonal block, and the noise standard deviation is set to 0.5 kW / s. The cumulative load deviation is used as the integral, and its process noise variance is set to 0.01 (kW·h). 2 The noise covariance matrix of the entire process is in a block diagonal form, with no cross-correlation between the blocks, and the matrix size is 31×31.

[0040] The time update step of the Kalman filter calculates the prior state estimate, starting from the initial time step, using the known posterior state estimate as the starting point. It is assumed that in the initial state estimate of the first time step, the response power of each device is 0 kW, the response rate is 0 kW / s, and the cumulative load deviation is 0 kW·h. The state transition matrix is ​​multiplied by the current posterior state estimate, and then added to the product of the control input matrix and the control vector to obtain the prior state estimate for the next time step. The elements of the state transition matrix are determined based on the dynamic response characteristics of the devices. For example, if the response time constant of a certain energy storage device is 5 s, the corresponding state transition coefficient is 0.8187. Matrix operations are used to obtain the prior state estimate at the end of the first time step, where the response power of each device reaches approximately 15% of the target regulation power.

[0041] The prior error covariance is calculated by propagating the current posterior error covariance through the state transition matrix, and then superimposing it with the process noise covariance matrix. Assuming the initial posterior error covariance matrix is ​​0.1 times the size of the identity matrix, indicating relatively low initial uncertainty, after the transformation of the state transition matrix and the injection of process noise, the elements of the prior error covariance matrix at the first time step increase, with the diagonal elements growing from approximately 0.1 to 0.15, reflecting the increased uncertainty during the prediction process. This matrix retains its symmetric positive definite property and has a size of 31×31.

[0042] A measurement vector is constructed by collecting real-time power measurements from each distributed energy source. A power measurement device is configured to sample the output power of 15 devices in real time at each time step. Assuming that at the end of the first time step, the measurement device has collected 15 power values ​​for each device, including 12kW, 8kW, and 15kW, etc., for a total of 15 values. The measurement vector is 15-dimensional and only includes directly observable response power variables, excluding non-measurable states such as response rate and cumulative load deviation. A 15×31 measurement matrix is ​​used to extract the observable portion from the complete state vector. The first 15 columns of this matrix correspond to the positions of response power variables as 1, and the remaining elements are 0.

[0043] The measurement noise covariance matrix describes the uncertainty of the measurement process and uses a diagonal structure. The standard deviation of the measurement noise for each power measurement device is set to 0.5 kW, corresponding to a measurement noise variance of 0.25 kW. 2 The measurement noise covariance matrix has a size of 15×15, with diagonal elements all having a value of 0.25 and off-diagonal elements having a value of 0, indicating that the noise of each measurement channel is independent of each other.

[0044] The Kalman gain matrix is ​​used to fuse prediction and measurement information. It is calculated using the prior error covariance, the measurement matrix, and the measurement noise covariance. Specifically, the measurement matrix is ​​multiplied by the prior error covariance matrix to obtain an intermediate matrix, which is then multiplied by the transpose of the measurement matrix and added to the measurement noise covariance matrix, forming a denominator matrix of size 15×15. This denominator matrix is ​​then inverted and multiplied by the product of the prior error covariance matrix and the transpose of the measurement matrix to obtain the Kalman gain matrix of size 31×15. The Kalman gain matrix typically ranges from 0 to 1. For example, a value of 0.6 indicates that the update of that state variable depends more on the measured value than the predicted value.

[0045] The posterior state estimate is obtained by correcting the prior state estimate. The difference between the measured vector and the predicted measured value obtained by applying the measured matrix to the prior state estimate is called the measurement residual. Assuming the actual measured power of a certain device is 12kW, the power predicted by the prior state estimate is 10kW, and the measurement residual is 2kW. Multiplying the Kalman gain matrix by the measurement residual vector yields the state correction, which is then added to the prior state estimate to obtain the posterior state estimate. In the corrected posterior state estimate, the response power of each device is closer to the actual measured value, and unmeasurable states such as response rate and cumulative load deviation are also indirectly corrected.

[0046] The posterior error covariance is obtained by adjusting the prior error covariance using the Kalman gain and the measurement matrix. The corrected matrix is ​​obtained by subtracting the product of the Kalman gain and the measurement matrix from the identity matrix. Multiplying this corrected matrix with the prior error covariance yields the posterior error covariance matrix. The diagonal elements of the corrected posterior error covariance are reduced, for example, from 0.15 to 0.08, indicating that the measurement update step reduces the uncertainty of the state estimation.

[0047] The time update and measurement update steps are executed recursively step by step until the 60th time step is completed, yielding the posterior state estimate and posterior error covariance at the end of the response duration, i.e., 600s. The covariance submatrix corresponding to the cumulative load deviation is extracted from the 31×31 posterior error covariance matrix. The cumulative load deviation is located at the 31st position in the state vector, and its corresponding covariance submatrix is ​​a scalar value, assumed to be 0.05 (kW·h). 2 .

[0048] Eigenvalue decomposition of the covariance submatrix identifies the main uncertainty direction. For the scalar case, the eigenvalue is the covariance value, 0.05 (kW·h). 2The eigenvectors are unit vectors. For multidimensional load deviation accumulation, the covariance submatrix is ​​a 2×2 or higher dimensional matrix. Eigenvalue decomposition yields multiple eigenvalues ​​and corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue indicates the main direction of change in the load deviation accumulation.

[0049] A multidimensional confidence ellipsoid is constructed to quantify the uncertainty range. A confidence level of 95% is selected, corresponding to a chi-square distribution quantile of 3.84. This quantile is multiplied by the square root of the eigenvalue to obtain the semi-axis length of the ellipsoid along the eigenvector direction, assuming a calculated value of 0.86 kW·h. The mean of the cumulative load deviation in the posterior state estimate is used as the center, assuming a mean of 2.5 kW·h. The confidence ellipsoid is constructed by extending the semi-axis length along the eigenvector direction.

[0050] The distribution boundary is obtained by projecting a multidimensional confidence ellipsoid onto the load deviation cumulative dimension. The projection operation extracts the coverage of the ellipsoid on the load deviation cumulative coordinate axis, resulting in the upper and lower bounds of the distribution. The lower bound is the mean minus the semi-axis length, calculated to be 1.64 kW·h. The upper bound is the mean plus the semi-axis length, calculated to be 3.36 kW·h. Combined with the mean of 2.5 kW·h, this fully describes the probability distribution characteristics of the load deviation cumulative at the end of the response duration. This distribution is presented in the form of a normal distribution, with a mean of 2.5 kW·h and a 95% confidence interval of 1.64 kW·h to 3.36 kW·h.

[0051] In one optional implementation, a Markov decision process is constructed based on the load deviation distribution. The cumulative load deviation is defined as a decision state set, and the adjustable capacity constraint is defined as a decision action set. A reward function is calculated, and a value iteration algorithm is used to solve for the optimal power allocation strategy. A response execution plan is then generated and issued, including: Based on the load deviation distribution, the expected value and variance of each distributed energy source are extracted, and a discretized state set of the cumulative load deviation is generated according to a preset step size. The Cartesian product of these states is then combined to construct the decision state set of the Markov decision process. Based on the adjustable capacity constraints of each distributed energy source, the target adjustable power is decomposed into power allocation increments, a feasible region that satisfies the total power balance and capacity constraints is constructed, and a set of decision actions is generated through grid sampling. Based on the superposition relationship between the cumulative load deviation and the power allocation increment under each decision state, the state transition probability matrix is ​​calculated using the variance of the deviation distribution, and a reward function is constructed that includes the square term of the adjustment power deviation, the weighted sum of squares of the load deviation, and the penalty term of the power allocation change rate. Based on the state transition probability and reward function, a value iteration algorithm is used to update the value function until convergence. According to the optimal value function, the power allocation increment scheme with the largest expected cumulative reward is selected and superimposed on the output power of each distributed energy source at the start of the response, which is then used as the response execution scheme and issued.

[0052] Before constructing the Markov decision process, key statistical characteristic parameters need to be extracted from the load deviation distribution. For three distributed energy devices in an industrial park, their load deviation data are recorded for consecutive time periods. The expected load deviation of the first device is -12kW with a variance of 25kW^2; the expected load deviation of the second device is +8kW with a variance of 16kW. 2 The expected value of the third device is -5kW, and the variance is 9kW. 2 The discretization step size is determined based on the variance. For the first device with a variance of 25, a step size of 5kW is set, generating discrete state points from -50kW to +50kW, for a total of 21 state values. The second device uses a step size of 4kW, generating discrete state points from -40kW to +40kW, for a total of 21 state values. The third device uses a step size of 3kW, generating discrete state points from -30kW to +30kW, for a total of 21 state values.

[0053] When performing the Cartesian product combination operation, the discrete state values ​​of the three devices are fully permuted and combined. In the specific operation, the -50kW state of the first device, the -40kW state of the second device, and the -30kW state of the third device form the first combined state. All combinations are traversed in turn to generate a decision state set containing 9261 combined states. Each combined state is represented by a triple. For example, the triple corresponding to the combined state number 5000 is (-10kW, +8kW, -6kW), which represents the current cumulative load deviation of the three devices.

[0054] When constructing the decision action set, it is necessary to clarify the adjustable capacity boundaries of each distributed energy source. The first device has an upward adjustable capacity of 30kW and a downward adjustable capacity of 20kW; the second device has an upward adjustable capacity of 25kW and a downward adjustable capacity of 15kW; and the third device has an upward adjustable capacity of 20kW and a downward adjustable capacity of 10kW. Assuming that the target power to be adjusted at the current response time is 50kW, this target power needs to be distributed across the three devices. The sampling step size for the power allocation increment is set to 5kW. The candidate value range for the power allocation increment is -20kW to +30kW for the first device, -15kW to +25kW for the second device, and -10kW to +20kW for the third device.

[0055] When establishing feasible region constraints, the sum of power allocation increments must equal the target regulating power of 50kW. A grid sampling method is used to traverse all candidate combinations. The combination of allocating 20kW to the first device, 20kW to the second device, and 10kW to the third device satisfies the total power balance constraint and is within the capacity constraints of each device; therefore, it is considered a feasible action. By exhaustively enumerating all combinations that satisfy the constraints, a decision action set containing 463 feasible power allocation schemes is generated. Each action is represented by a triplet for the power allocation increment; the triplet corresponding to action number 200 is (25kW, 15kW, 10kW).

[0056] When calculating the state transition probability matrix, the superposition effect of the cumulative load deviation and the power allocation increment is considered. In the current state, the cumulative load deviation of the first device is -10kW. After executing the action to allocate a 20kW power increment, the cumulative value becomes +10kW. Due to the presence of random disturbances, based on the variance of 25kW... 2 Calculate the transition probability distribution. Following a normal distribution centered at +10kW with a standard deviation of 5kW, the probability of transitioning to the +10kW state is 0.398, the probability of transitioning to the +15kW state is 0.242, the probability of transitioning to the +5kW state is 0.242, and the probability of transitioning to other adjacent states decreases sequentially. Calculate the individual transition probabilities for each of the three devices, and obtain the transition probability of the combined state by multiplying the probabilities.

[0057] The reward function is constructed from three core components: The first is the squared term for the power deviation, which measures the deviation between the actual adjustment effect and the target. When the sum of the power allocation increments of the three devices is 50kW but the target adjustment power is 52kW, the deviation is 2kW, and this term has a value of -4. The second is the weighted sum of squares term for load deviation, which penalizes states where the cumulative load deviation is too large. For example, the first device has a current cumulative load deviation of +15kW with a weighting coefficient of 0.6, contributing a penalty of -135; the second device has a cumulative load deviation of -8kW with a weighting coefficient of 0.8, contributing -51.2; and the third device has a cumulative load deviation of +6kW with a weighting coefficient of 1.0, contributing -36. The third is the penalty term for the power allocation change rate, which limits drastic changes in power allocation between adjacent time steps. For example, the first device's allocation increment was 15kW at the previous time step and is currently 25kW, with a change rate of 10kW and a penalty coefficient of 0.2, contributing a penalty of -20. The sum of these three components yields the immediate reward value for the current state action pair.

[0058] The value iteration algorithm approximates the optimal solution by repeatedly updating the value function. The value function for all states is initialized to 0, and a discount factor of 0.95 is set. In the first iteration, for all possible actions in state 5000, the weighted sum of the immediate reward and the value of the subsequent state is calculated. The immediate reward for action 200 in state 5000 is -200, and the expected value of the subsequent state is -150 (weighted by transition probabilities), discounted to -142.5. The expected cumulative reward for this action is -342.5. After traversing all possible actions, the action number with the largest expected cumulative reward and its corresponding value are selected to update the value function of state 5000. The algorithm is considered converged when the change in the value function of all states is less than 0.01 in two consecutive iterations.

[0059] Based on the converged optimal value function, the optimal action corresponding to each state is extracted. Assuming the current system is in state number 5000, the optimal action for this state is number 200, corresponding to power allocation increments of (25kW, 15kW, 10kW). The actual output powers of the three devices at the start of the response are obtained as 80kW, 60kW, and 40kW respectively. The power allocation increments are then added to the actual output powers, resulting in a response execution power of 105kW for the first device, 75kW for the second device, and 50kW for the third device. A response execution plan data packet containing the target power values ​​and execution timestamps of the three devices is generated and sent to the local controllers of each distributed energy source via the communication interface to execute power adjustment commands.

[0060] In one optional implementation, the power change sequence is converted into a load recovery demand vector, and a recovery capacity constraint cone with vertices is constructed based on the current available adjustment margin of the power grid. The optimal projection point of each load recovery demand vector within the constraint cone is then solved, including: The peak deviation of each distributed energy source is calculated based on the power change sequence, and the average recovery power demand is calculated in combination with the preset recovery time to generate a load recovery demand vector. The electrical distance between nodes is calculated based on the voltage measurement values ​​of the power grid bus nodes. Multiple voltage control zones are divided, and the reactive and active reserve capacity of each zone is obtained. The reactive reserve capacity is converted into the maximum load recovery power to form a zone reserve capacity vector. The remaining transmission margin is calculated based on the current transmission power and thermal stability limit of the tie line. A cross-zone power transmission constraint matrix is ​​established, and the negative values ​​of the zone reserve capacity vectors are used as vertices and the row vectors of the transmission constraint matrix are used as the cone surface normal vectors to generate the recovery capacity constraint cone. Calculate the line vector connecting each load recovery demand vector to the vertex of the recovery capacity constraint cone. Determine whether the load recovery demand vector is located inside the constraint cone by judging the sign of the inner product of the connecting line vector and the cone surface normal vector. For a load recovery demand vector located outside the constraint cone, an optimization problem is constructed with the square of the distance from the projection point to the original vector as the objective function. The optimal projection point is calculated by solving the Lagrange dual problem.

[0061] In the process of converting the power change sequence into a load restoration demand vector, real-time power data of each distributed energy access point is continuously monitored. For a certain photovoltaic inverter, the output power is 150kW during normal operation. When a voltage drop fault occurs, the output power drops to 20kW. The fluctuation of this power change sequence on the time axis is recorded. When calculating the peak deviation, the minimum value of 20kW in the power change sequence is extracted and the difference is calculated with the normal operating power of 150kW, resulting in a peak deviation of 130kW. According to the preset recovery time parameter, for example, setting the recovery time to 10s, the peak deviation of 130kW is divided by the recovery time of 10s to calculate the average restoration power demand of 13kW / s. For the 8 distributed energy access points in the distribution network, the average restoration power demand of each point is calculated as follows: 13kW / s, 9kW / s, 15kW / s, 11kW / s, 7kW / s, 12kW / s, 10kW / s, and 8kW / s. These 8 values ​​are arranged in order of access point number to form an 8-dimensional load restoration demand vector.

[0062] Voltage measurements at power grid bus nodes are collected in real time by voltage sensors deployed at each node. For example, in a distribution network with 20 bus nodes, node 1 has a voltage amplitude of 10.2 kV and a phase angle of 0°, while node 2 has a voltage amplitude of 10.1 kV and a phase angle of -2°. When calculating the electrical distance between nodes, the voltage phasors of the two nodes are used for vector difference calculation to obtain the voltage difference amplitude and phase angle. The voltage difference amplitude between node 1 and node 2 is 0.15 kV, and the phase angle difference is 2°. The voltage difference amplitude and phase angle difference are weighted and combined with weighting coefficients of 0.7 and 0.3 respectively, resulting in an electrical distance of 0.165 pu between node 1 and node 2. The electrical distances between all node pairs are calculated sequentially, constructing a 20×20 electrical distance matrix. Clustering is performed based on this matrix, setting an electrical distance threshold of 0.3 pu. Nodes with electrical distances less than this threshold are grouped into the same voltage control zone. This partitioning method divides the 20 nodes into 3 voltage control zones: Zone 1 contains nodes 1 to 7, Zone 2 contains nodes 8 to 14, and Zone 3 contains nodes 15 to 20.

[0063] When obtaining the reactive and active power reserve capacity of each zone, the rated capacity and current operating power of all grid-connected generating equipment in Zone 1 are queried. Zone 1 contains two synchronous condensers, each with a rated capacity of 5 Mvar, and current operating capacities of 3 Mvar and 2.5 Mvar respectively, resulting in a reactive power reserve capacity of 4.5 Mvar. This zone also contains one gas turbine with a rated active power of 10 MW and a current output power of 7 MW, resulting in an active power reserve capacity of 3 MW. When converting the reactive power reserve capacity to the maximum load recovery power, based on the typical load power factor of 0.9 for this zone, the ratio of reactive power to active power is calculated to be 0.484. Dividing the reactive power reserve capacity of 4.5 Mvar by this ratio yields a convertible active power of 9.3 MW. The converted active power of 9.3MW is superimposed with the original active reserve capacity of 3MW to obtain the total reserve capacity of 12.3MW for the first zone. The total reserve capacities of the second and third zones are calculated in the same way as 8.7MW and 10.5MW respectively. The three values ​​are arranged according to the zone number to form a 3D zone reserve capacity vector.

[0064] The calculation of the remaining transmission margin of tie lines involves the transmission status of lines connecting different voltage control zones. The current transmission power of the tie line between Zone 1 and Zone 2 is 4.2MW, and the thermal stability limit of this tie line is 8MW. Subtracting the current transmission power from the thermal stability limit yields a remaining transmission margin of 3.8MW. Similarly, the remaining transmission margin of the tie line between Zone 2 and Zone 3 is calculated to be 2.9MW, and the remaining transmission margin of the tie line between Zone 1 and Zone 3 is 4.1MW. When establishing the cross-zone power transmission constraint matrix, the number of rows in the matrix equals the number of tie lines (3), and the number of columns equals the number of zones (3). The first row of the matrix corresponds to the tie line between Zone 1 and Zone 2, with values ​​of 1, -1, and 0, indicating power flow from Zone 1 to Zone 2. The second row of the matrix corresponds to the connection line between the second and third partitions, with values ​​of 0, 1, and -1. The third row of the matrix corresponds to the connection line between the first and third partitions, with values ​​of 1, 0, and -1. By taking the negative values ​​of each component of the partition spare capacity vector, we obtain -12.3MW, -8.7MW, and -10.5MW. This point is used as the vertex of the recovery capacity constraint cone. Each row vector of the transmission constraint matrix is ​​used as the normal vector of the cone surface, defining the boundary surface of the constraint cone.

[0065] When determining whether the load restoration demand vector lies inside the constraint cone, a vector operation method is used. Given a load restoration demand vector with three components -9kW, -6kW, and -8kW, the differences between these components and the constraint cone vertex coordinates -12.3MW, -8.7MW, and -10.5MW are calculated, resulting in a connecting vector of 3.3MW, 2.7MW, and 2.5MW. Multiplying these connecting vectors by the corresponding components of the first cone surface normal vector (1, -1, 0) and summing the results, we get 3.3 - 2.7 = 0.6MW. 2 The inner product result is positive, indicating that the connecting vector and the normal vector point to the same half of space. Calculating the inner products of the connecting vector with the other two normal vectors yields -0.2MW^2 and 0.8MW respectively. 2 Because of the existence of a negative inner product, the load recovery demand vector is determined to be outside the constraint cone, and projection calculation is required.

[0066] When constructing the optimization problem, the three coordinate components of the projection point are set as variables to be solved, denoted as the first component, the second component, and the third component, respectively. The objective function is the squared distance from the projection point to the original load recovery demand vector, which is represented by the sum of the squared differences between each component of the projection point and the corresponding component of the original vector. The constraints include that the projection point must be located inside the constraint cone, and that the inner product of the line vector connecting the projection point to the cone vertex and the normal vectors of all cone surfaces must be non-negative. Lagrange multipliers are introduced to relax the constraints, and the number of multipliers is equal to the number of constraints, which is 3. By adding the constraint condition and the product of the Lagrange multipliers to the objective function, a Lagrange function is constructed. The partial derivatives of the Lagrange function with respect to the three components of the projection point and the three Lagrange multipliers are calculated respectively. Setting the partial derivatives to 0, a system of six equations is obtained. Solving this system of equations yields the three components of the optimal projection point as -10.8MW, -7.9MW, and -9.6MW. It is verified that the inner product of the line vector connecting the projection point to the vertex of the cone and all normal vectors is non-negative, satisfying the internal conditions of the constrained cone, and thus confirming it as the optimal projection point.

[0067] In one optional implementation, the remaining transmission margin is calculated based on the current transmission power and thermal stability limit of the tie line, and a cross-zone power transmission constraint matrix is ​​established. A recovery capacity constraint cone is generated using the negative values ​​of the zone reserve capacity vectors as vertices and the row vectors of the transmission constraint matrix as the cone surface normal vectors. This includes: Calculate the remaining transmission margin of each tie line based on the real-time transmission power and corresponding thermal stability limit of the tie lines between each voltage control zone. Based on the remaining transmission margin and the topological connection relationship of each voltage control zone, a cross-zone power transmission constraint matrix is ​​constructed to characterize the directional constraints of power flow between zones. Each row of the cross-zone power transmission constraint matrix corresponds to the power transmission boundary condition of a tie line, and the absolute value of the matrix element represents the power transmission capacity coefficient of the corresponding zone on the tie line. The negative values ​​of each component of the partition spare capacity vector are used to obtain the constraint cone vertex, which represents the recovery capability benchmark of each partition when there is no cross-partition power support. Extract each row vector of the cross-regional power transmission constraint matrix as a cone normal vector. Each cone normal vector is perpendicular to the hyperplane passing through the vertex of the constraint cone. Define a cone boundary. The intersection of multiple cone boundaries forms a cone space, resulting in a recovery capacity constraint cone with the vertex of the constraint cone as the cone apex.

[0068] In this specific embodiment, real-time operating data of the power system is collected, and the transmission power of the tie lines between each voltage control zone is monitored. In actual operation, it is assumed that a power system contains four voltage control zones, denoted as Zone A, Zone B, Zone C, and Zone D. Zone A and Zone B are connected by tie line L1, Zone A and Zone C are connected by tie line L2, Zone B and Zone D are connected by tie line L3, and Zone C and Zone D are connected by tie line L4. Power measurement devices deployed on each tie line collect the transmission power data of each tie line in real time. At a certain moment, the real-time transmission power of tie line L1 is 350MW, and its thermal stability limit is 500MW; the real-time transmission power of tie line L2 is 280MW, and its thermal stability limit is 450MW; the real-time transmission power of tie line L3 is 180MW, and its thermal stability limit is 300MW; and the real-time transmission power of tie line L4 is 220MW, and its thermal stability limit is 400MW.

[0069] For each tie line, its remaining transmission margin is calculated. This margin value is equal to the thermal stability limit value minus the current real-time transmission power value. The remaining transmission margin of tie line L1 is calculated as 500MW-350MW=150MW; the remaining transmission margin of tie line L2 is 450MW-280MW=170MW; the remaining transmission margin of tie line L3 is 300MW-180MW=120MW; and the remaining transmission margin of tie line L4 is 400MW-220MW=180MW. These remaining transmission margin values ​​represent the additional power capacity that each tie line can transmit without exceeding the thermal stability limit.

[0070] Based on the topological connections and remaining transmission margins between partitions, a cross-partition power transmission constraint matrix is ​​constructed. The number of rows in this matrix equals the number of tie lines, and the number of columns equals the number of partitions. Each row of the matrix corresponds to the power transmission boundary conditions for one tie line. For tie line L1 connecting partitions A and B, in the first row of the matrix, the elements in the corresponding column for partition A are filled with +1, the elements in the corresponding column for partition B are filled with -1, and the elements in the corresponding columns for the other partitions are filled with 0. This row represents the power transmission constraint from partition A to partition B. For tie line L2 connecting partitions A and C, in the second row of the matrix, the elements in the corresponding column for partition A are filled with +1, the elements in the corresponding column for partition C are filled with -1, and the elements in the corresponding columns for partitions B and D are filled with 0. For tie line L3 connecting partitions B and D, in the third row of the matrix, the elements in the corresponding column for partition B are filled with +1, the elements in the corresponding column for partition D are filled with -1, and the elements in the corresponding columns for partitions A and C are filled with 0. For the case where tie line L4 connects partitions C and D, in the fourth row of the matrix, the elements in the corresponding column for partition C are filled with +1, the elements in the corresponding column for partition D are filled with -1, and the elements in the corresponding columns for partitions A and B are filled with 0. The absolute values ​​of the matrix elements reflect the power transmission capacity coefficient of each partition on the corresponding tie line, and the positive and negative signs reflect the directionality of power flow.

[0071] Reserve capacity data for each zone is obtained through the power system energy management system. Assuming zone A has a reserve capacity of 200MW, zone B 150MW, zone C 180MW, and zone D 120MW, the zone reserve capacity vector represents the generation resources available for load restoration in each zone under the current state. Taking the negative value of each component of this vector yields the vertex of the constraint cone. The coordinate component corresponding to zone A is -200MW, zone B is -150MW, zone C is -180MW, and zone D is -120MW. This constraint cone vertex defines a specific location in 4D space, representing the baseline point of each zone's recovery capacity without cross-zone power support.

[0072] The row vectors of the cross-zone power transfer constraint matrix are extracted as the normal vectors of the cone surfaces. The first row vector of the matrix contains a vector with a component of +1 corresponding to zone A, a component of -1 corresponding to zone B, and a component of 0 corresponding to zones C and D. This vector serves as the normal vector of the first cone surface. The second row vector contains the same vector, which serves as the normal vector of the second cone surface. The third row vector contains the same vector, which serves as the normal vector of the third cone surface. The fourth row vector contains the same vector, which serves as the normal vector of the fourth cone surface.

[0073] Each cone normal vector defines a hyperplane boundary passing through the constrained cone vertex. This hyperplane is perpendicular to the corresponding normal vector, forming an interface in the cone space. The boundary conditions of the hyperplane corresponding to the connection line between partition A and partition B are determined by the fact that the difference between the recovery capacity of partition A and partition B does not exceed the remaining transmission margin of 150MW for connection line L1; the boundary conditions of the hyperplane corresponding to the connection line between partition A and partition C are determined by the fact that the difference between the recovery capacity of partition A and partition C does not exceed the remaining transmission margin of 170MW for connection line L2; the boundary conditions of the hyperplane corresponding to the connection line between partition B and partition D are determined by the fact that the difference between the recovery capacity of partition B and partition D does not exceed the remaining transmission margin of 120MW for connection line L3; and the boundary conditions of the hyperplane corresponding to the connection line between partition C and partition D are determined by the fact that the difference between the recovery capacity of partition C and partition D does not exceed the remaining transmission margin of 180MW for connection line L4.

[0074] Multiple hyperplane boundaries intersect in 4-dimensional space, and their intersection forms a cone-shaped spatial structure with the apex of the constraint cone as the cone's top. This cone-shaped space contains all combinations of partitioned recovery capacities that satisfy the tie-line transmission constraints. Any point inside the cone represents a feasible partitioned power support scheme that satisfies the transmission capacity limits of each tie-line without violating the thermal stability limit constraints.

[0075] The boundary of the cone is jointly determined by the remaining transmission margin of each tie line, and the shape and size of the cone reflect the cross-regional power transmission capacity of the power system under the current operating conditions. By searching for the optimal solution within this recovery capacity constraint cone, the load recovery amount of each region can be determined, thereby achieving coordinated recovery control of the power system.

[0076] In one optional implementation, the system is layered based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, and divided into multiple sequentially executed recovery stages to obtain a phased recovery sequence, including: Calculate the Euclidean distance between each optimal projection point and the vertex of the recovery capacity constraint cone. Combined with the geometric topology of the recovery capacity constraint cone, calculate the curvature distribution of the constraint cone section where each optimal projection point is located. Then, divide the recovery capacity constraint cone radially into multiple nested curvature isosurface regions to obtain multiple recovery levels. For the curvature isosurface regions of each recovery level, a tangential recovery manifold is constructed. Based on the optimal projection point, the gradient direction of each distributed energy source on the tangential recovery manifold and the coupling strength with the radial contraction direction of the constraint cone are calculated. The distributed energy sources are sorted according to the coupling strength to determine the recovery execution order of each recovery level. The recovery levels are arranged in order from the outside to the inside according to the curvature isosurface, and the distributed energy sources arranged in the recovery execution order are connected in series to generate the phased recovery sequence. A recovery instruction corresponding to the phased recovery sequence is issued to each distributed energy source, and the recovery instruction carries a corresponding execution timing identifier.

[0077] After obtaining the optimal projection points of each distributed energy source, a hierarchical recovery sequence generation process based on geometric topology is executed. For each optimal projection point, its 3D coordinate position in the recovery capacity constraint cone space is calculated. The vertex of the constraint cone represents the complete failure state of the system, and its coordinates are set as the origin. Taking a distribution network containing 20 distributed photovoltaic power generation units as an example, the optimal projection point coordinates of the first photovoltaic unit are obtained as follows: active power dimension value 86kW, reactive power dimension value 42kvar, and voltage dimension value 10.2kV. The spatial straight-line distance from this point to the origin is calculated. Using the distance calculation method between two points in 3D space, the coordinate differences in the three dimensions are squared and summed. Then, the square root of the sum is calculated, and the Euclidean distance from the optimal projection point to the vertex is found to be 96.8 units. The same calculation is performed on the optimal projection points of all 20 distributed photovoltaic units, obtaining distribution data with distance values ​​ranging from 75.3 to 158.6 units.

[0078] Based on the geometric and topological characteristics of the recovery capacity constraint cone, a cross-sectional analysis system is established inside the constraint cone. The constraint cone exhibits a conical structure expanding outward from the apex, and multiple cross-sections perpendicular to the cone axis are selected for curvature analysis. For the cross-section containing the optimal projection point, the local curvature of the boundary curve of this cross-section is calculated, and the curvature value is quantified by measuring the radius of the arc formed by three adjacent sampling points on the cross-section boundary. At the first cross-section location, a sampling point is taken at every 5° angle on the boundary curve, resulting in a total of 72 sampling points, and the curvature value at each sampling point is calculated. Due to the non-uniformity of the constraint cone, the curvature at different locations exhibits a differential distribution, and the recorded curvature values ​​range from 0.018 to 0.053 / unit length. Cross-sectional regions with similar curvature values ​​are clustered, and a curvature difference threshold of 0.005 / unit length is set, dividing the entire constraint cone space into 8 curvature isosurface regions. These regions exhibit a nested shell structure, with the recovery progress decreasing from the outer layer to the inner layer.

[0079] Eight curvature isosurface regions are defined as eight restoration levels. Each level contains several optimal projection points. The first restoration level is located in the outermost region and contains the five optimal projection points farthest from the vertex, corresponding to five distributed photovoltaic units, which have the lowest restoration difficulty. The second restoration level shrinks inward and contains four optimal projection points, and so on until the eighth restoration level is located in the innermost layer and contains only one optimal projection point. For the first restoration level, a tangential restoration manifold is constructed on the corresponding curvature isosurface. This manifold represents the restoration path space along the tangent direction of the isosurface while keeping the distance from the vertex essentially unchanged.

[0080] The local geometric characteristics of each optimal projection point on the tangential restoring manifold in the first restoration level are calculated. Taking photovoltaic unit number 03 in this level as an example, its optimal projection point coordinates are 152kW active power, 68kvar reactive power, and 10.8kV voltage. A local coordinate system of the tangential restoring manifold is constructed at this point. By calculating the normal vector of the manifold surface at this point, the tangential plane perpendicular to the radial direction is determined. Within the tangential plane, the steepest descent direction when the photovoltaic unit moves from the current projection point to a better operating state is analyzed. This direction points to the position with a larger capacity margin on the boundary of the constraint cone. By calculating the partial derivatives of the three power and voltage dimensions at this point, the gradient vector is obtained. This vector has a component of 8.3 in the active power dimension, a component of 3.6 in the reactive power dimension, and a component of 1.2 in the voltage dimension.

[0081] The spatial relationship between the gradient direction and the radial contraction direction of the constraint cone is analyzed. The radial contraction direction of the constraint cone points to the vertex, representing the decreasing trend of system capacity. The cosine of the angle between the gradient vector and the radial vector is calculated to quantify the coupling strength. The closer the two vector directions are, the higher the coupling strength, indicating that the recovery process of this distributed energy source is more sensitive to the overall system capacity contraction. The gradient vector of photovoltaic unit 03 has an angle of 38° with the radial vector, and the cosine of the angle is 0.788, indicating a strong coupling relationship. The coupling strength was calculated for all five photovoltaic units in the first recovery level, and the values ​​obtained were 0.823, 0.756, 0.788, 0.692, and 0.715, respectively. The units with higher coupling strength were sorted from highest to lowest, and the units with higher coupling strength were given priority for recovery, because the recovery operations of these units can more effectively promote the overall system to a safe state. The recovery execution order of the first recovery level was determined to be units 01, 03, 02, 05, and 04.

[0082] The same tangential manifold construction, gradient calculation, and coupling strength analysis process was performed on the remaining seven recovery levels. The second recovery level contains four photovoltaic units, with coupling strength calculations of 0.801, 0.769, 0.744, and 0.728, and the execution order was determined to be numbers 06, 08, 07, and 09. The third recovery level contains three photovoltaic units, with coupling strengths of 0.812, 0.776, and 0.751, and the execution order was numbers 11, 10, and 12. After completing the internal sorting of all eight recovery levels, these levels were arranged in spatial order from the outside to the inside according to the curvature isosurface, forming a complete phased recovery sequence. This sequence starts from the first recovery level and executes the second to eighth recovery levels sequentially, with distributed energy units connected in series within each level according to the determined execution order.

[0083] The complete phased recovery sequence is as follows: Phase 1 recovers 5 units numbered 01, 03, 02, 05, and 04; Phase 2 recovers 4 units numbered 06, 08, 07, and 09; Phase 3 recovers 3 units numbered 11, 10, and 12, and so on until Phase 8, which recovers unit number 20. A recovery instruction data package is generated for each distributed energy unit. The instruction content includes target operating parameters and execution timing identifiers. The recovery instruction for photovoltaic unit number 01 sets the target active power to 152kW, the target reactive power to 68kvar, and the target voltage to 10.8kV. Its execution timing identifier is Phase 1, position 1, indicating that it will be executed first in Phase 1. The execution timing identifier for photovoltaic unit number 03 is Phase 1, position 2, indicating that it will be executed immediately after unit number 01 is completed.

[0084] The distribution network communication network sends recovery commands to the local controllers of each distributed energy unit. The local controllers initiate recovery operations at the specified stage and sequence position according to the received execution timing identifier, and adjust the inverter output parameters to make the actual operating state close to the target parameter value, thereby realizing the orderly and phased recovery process of the distribution network.

[0085] This invention provides a virtual power plant demand-side response and load management system, comprising: The data acquisition unit is used to acquire the demand-side response requests and historical response trajectories of each distributed energy source, wherein the demand-side response requests include the target adjustment power and the response duration. The load deviation calculation unit is used to perform time-frequency domain decomposition on the historical response trajectory, extract transient and steady-state components, and construct a dynamic response model that includes the current output power and the cumulative load deviation. Based on the dynamic response model, the target adjustment power is used as the control input, and the load deviation distribution at the end of the response duration is calculated by Kalman filtering. The decision optimization unit is used to construct a Markov decision process based on the load deviation distribution, define the combination of load deviation accumulation as a decision state set, combine it with the adjustable capacity constraint as a decision action set, calculate the reward function, and use the value iteration algorithm to solve the optimal power allocation strategy, generate a response execution plan and issue it. The load recovery analysis unit is used to collect the power change sequence of each distributed energy source after executing the response execution scheme, convert the power change sequence into a load recovery demand vector, construct a recovery capacity constraint cone with vertices based on the current available adjustment margin of the power grid, and solve for the optimal projection point of each load recovery demand vector in the constraint cone. The recovery scheduling unit is used to divide the recovery into layers based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, and to divide the recovery into multiple recovery stages that are executed sequentially, thus obtaining a phased recovery sequence.

[0086] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0087] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0088] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0089] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A virtual power plant demand-side response and load management method, characterized in that, include: Obtain the demand-side response requests and historical response trajectories of each distributed energy source, wherein the demand-side response requests include the target adjustment power and the response duration; The historical response trajectory is decomposed in the time-frequency domain to extract transient and steady-state components, and a dynamic response model containing the current output power and the cumulative load deviation is constructed. Based on the dynamic response model, the target adjustment power is used as the control input, and the load deviation distribution at the end of the response duration is calculated by Kalman filtering. A Markov decision process is constructed based on the load deviation distribution. The cumulative load deviation is defined as the decision state set, and the adjustable capacity constraint is defined as the decision action set. The reward function is calculated, and the optimal power allocation strategy is solved using a value iteration algorithm. A response execution plan is generated and issued. Collect the power change sequence of each distributed energy source after executing the response execution scheme, convert the power change sequence into a load recovery demand vector, construct a recovery capacity constraint cone with vertices based on the current available adjustment margin of the power grid, and solve for the optimal projection point of each load recovery demand vector in the constraint cone; The system is divided into layers based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, and then into multiple recovery stages executed sequentially to obtain a phased recovery sequence.

2. The method according to claim 1, characterized in that, The historical response trajectory is decomposed in the time-frequency domain to extract transient and steady-state components, and a dynamic response model including the current output power and the cumulative load deviation is constructed, including: Perform wavelet transform decomposition on the historical response trajectory to obtain a multi-scale time-frequency coefficient matrix, calculate the energy distribution of the coefficients at each scale, determine the transient interval, reconstruct the time-frequency coefficients by inverse wavelet transform to obtain the transient component, and reconstruct the time-frequency coefficients outside the transient interval by inverse wavelet transform to obtain the steady-state component. Calculate the transient response time constant of the transient component and establish the transient evolution equation of the current output power; calculate the power deviation sequence of the steady-state component and perform time integration to obtain the evolution data of the load deviation accumulation, fit the accumulation rate function and calculate the steady-state accumulation coefficient, and establish the steady-state evolution equation of the load deviation accumulation. The current output power and the cumulative load deviation are constructed into a state vector. The transient response time constant and the steady-state cumulative coefficient are embedded into the state transition matrix. The transient evolution equation and the steady-state evolution equation are discretized and combined into a dynamic response model in state space form.

3. The method according to claim 1, characterized in that, Based on the dynamic response model, using the target regulating power as the control input, the load deviation distribution at the end of the response duration is calculated using Kalman filtering, including: The target regulation power is discretized according to the response duration to generate a control input sequence, which is substituted into the dynamic response model as the control vector. A process noise covariance matrix with a block diagonal structure is constructed, and the time update step of Kalman filtering is performed to calculate the prior state estimate and prior error covariance. Collect real-time power measurements of each distributed energy source, perform measurement update steps and calculate Kalman gain, correct the prior state estimate, and obtain the posterior state estimate and posterior error covariance. Based on the posterior error covariance, calculate the covariance submatrix corresponding to the cumulative load deviation, perform eigenvalue decomposition on the covariance submatrix, identify the eigenvector corresponding to the direction of the cumulative load deviation, construct a multidimensional confidence ellipsoid along the direction of the eigenvector, use its projection boundary on the dimension of the cumulative load deviation as the upper and lower bounds of the distribution, and combine it with the mean of the cumulative load deviation in the posterior state estimation to construct the load deviation distribution.

4. The method according to claim 1, characterized in that, A Markov decision process is constructed based on the load deviation distribution. The cumulative load deviation is defined as the decision state set, and the adjustable capacity constraint is defined as the decision action set. The reward function is calculated, and a value iteration algorithm is used to solve for the optimal power allocation strategy. A response execution plan is generated and issued, including: Based on the load deviation distribution, the expected value and variance of each distributed energy source are extracted, and a discretized state set of the cumulative load deviation is generated according to a preset step size. The state set of the Markov decision process is constructed by combining the Cartesian product of these states. Based on the adjustable capacity constraints of each distributed energy source, the target adjustable power is decomposed into power allocation increments, a feasible region that satisfies the total power balance and capacity constraints is constructed, and a set of decision actions is generated through grid sampling. Based on the superposition relationship between the cumulative load deviation and the power allocation increment under each decision state, the state transition probability matrix is ​​calculated using the variance of the deviation distribution, and a reward function is constructed that includes the square term of the adjustment power deviation, the weighted sum of squares of the load deviation, and the penalty term of the power allocation change rate. Based on the state transition probability and reward function, a value iteration algorithm is used to update the value function until convergence. According to the optimal value function, the power allocation increment scheme with the largest expected cumulative reward is selected and superimposed on the output power of each distributed energy source at the start of the response, which is then used as the response execution scheme and issued.

5. The method according to claim 1, characterized in that, The power change sequence is converted into a load recovery demand vector. A recovery capacity constraint cone with vertices is constructed based on the current available adjustment margin of the power grid. The optimal projection point of each load recovery demand vector within the constraint cone is then solved, including: The peak deviation of each distributed energy source is calculated based on the power change sequence, and the average recovery power demand is calculated in combination with the preset recovery time to generate a load recovery demand vector. The electrical distance between nodes is calculated based on the voltage measurement values ​​of the power grid bus nodes. Multiple voltage control zones are divided, and the reactive and active reserve capacity of each zone is obtained. The reactive reserve capacity is converted into the maximum load recovery power to form a zone reserve capacity vector. Calculate the remaining transmission margin based on the current transmission power and thermal stability limit of the tie line, establish a cross-zone power transmission constraint matrix, and generate a recovery capacity constraint cone with the negative values ​​of the zone reserve capacity vector as vertices and the row vectors of the transmission constraint matrix as the cone surface normal vectors. Calculate the line vector connecting each load recovery demand vector to the vertex of the recovery capacity constraint cone. Determine whether the load recovery demand vector is located inside the constraint cone by judging the sign of the inner product of the connecting line vector and the cone surface normal vector. For a load recovery demand vector located outside the constraint cone, an optimization problem is constructed with the square of the distance from the projection point to the original vector as the objective function. The optimal projection point is calculated by solving the Lagrange dual problem.

6. The method according to claim 5, characterized in that, Based on the current transmission power and thermal stability limit of the tie line, the remaining transmission margin is calculated, and a cross-regional power transmission constraint matrix is ​​established. Using the negative values ​​of the regional reserve capacity vectors as vertices and the row vectors of the transmission constraint matrix as the cone normal vectors, a recovery capacity constraint cone is generated, including: Calculate the remaining transmission margin of each tie line based on the real-time transmission power and corresponding thermal stability limit of the tie lines between each voltage control zone. Based on the remaining transmission margin and the topological connection relationship of each voltage control zone, a cross-zone power transmission constraint matrix is ​​constructed to characterize the directional constraints of power flow between zones. Each row of the cross-zone power transmission constraint matrix corresponds to the power transmission boundary condition of a tie line, and the absolute value of the matrix element represents the power transmission capacity coefficient of the corresponding zone on the tie line. The negative values ​​of each component of the partition spare capacity vector are used to obtain the constraint cone vertex, which represents the recovery capability benchmark of each partition when there is no cross-partition power support. Extract each row vector of the cross-regional power transmission constraint matrix as a cone normal vector. Each cone normal vector is perpendicular to the hyperplane passing through the vertex of the constraint cone. Define a cone boundary. The intersection of multiple cone boundaries forms a cone space, resulting in a recovery capacity constraint cone with the vertex of the constraint cone as the cone apex.

7. The method according to claim 1, characterized in that, Based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, the system is layered and divided into multiple sequentially executed recovery stages, resulting in a phased recovery sequence, including: Calculate the Euclidean distance between each optimal projection point and the vertex of the recovery capacity constraint cone. Combined with the geometric topology of the recovery capacity constraint cone, calculate the curvature distribution of the constraint cone section where each optimal projection point is located. Then, divide the recovery capacity constraint cone radially into multiple nested curvature isosurface regions to obtain multiple recovery levels. For the curvature isosurface regions of each recovery level, a tangential recovery manifold is constructed. Based on the optimal projection point, the gradient direction of each distributed energy source on the tangential recovery manifold and the coupling strength with the radial contraction direction of the constraint cone are calculated. The distributed energy sources are sorted according to the coupling strength to determine the recovery execution order of each recovery level. The recovery levels are arranged in order from the outside to the inside according to the curvature isosurface, and the distributed energy sources arranged in the recovery execution order are connected in series to generate the phased recovery sequence. A recovery instruction corresponding to the phased recovery sequence is issued to each distributed energy source, wherein each recovery instruction carries a corresponding execution timing identifier.

8. A virtual power plant demand-side response and load management system, used to implement the method as described in any one of claims 1-7, characterized in that, include: The data acquisition unit is used to acquire the demand-side response requests and historical response trajectories of each distributed energy source, wherein the demand-side response requests include the target adjustment power and the response duration. The load deviation calculation unit is used to perform time-frequency domain decomposition on the historical response trajectory, extract transient and steady-state components, and construct a dynamic response model that includes the current output power and the cumulative load deviation. Based on the dynamic response model, the target adjustment power is used as the control input, and the load deviation distribution at the end of the response duration is calculated by Kalman filtering. The decision optimization unit is used to construct a Markov decision process based on the load deviation distribution, define the combination of load deviation accumulation as a decision state set, combine it with the adjustable capacity constraint as a decision action set, calculate the reward function, and use the value iteration algorithm to solve the optimal power allocation strategy, generate a response execution plan and issue it. The load recovery analysis unit is used to collect the power change sequence of each distributed energy source after executing the response execution scheme, convert the power change sequence into a load recovery demand vector, construct a recovery capacity constraint cone with vertices based on the current available adjustment margin of the power grid, and solve for the optimal projection point of each load recovery demand vector in the constraint cone. The recovery scheduling unit is used to divide the recovery into layers based on the distance between the optimal projection point and the vertex of the recovery capacity constraint cone, and to divide the recovery into multiple recovery stages that are executed sequentially, thus obtaining a phased recovery sequence.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images