Virtual power plant adjustable load aggregation and real-time response regulation method and system

Abstract

The application provides a virtual power plant adjustable load aggregation and real-time response regulation method and system, relates to the technical field of virtual power plants, and comprises the following steps: constructing and decomposing a load state space model, aggregating to form a group equivalent representation, calculating the deviation of the group equivalent representation from a dispatching instruction and reversely correcting a model coupling coefficient, and finally performing simulation prediction and generating a dispatching decision based on the updated equivalent representation. The method realizes accurate aggregation and dynamic regulation of distributed flexible loads, and improves the real-time performance and accuracy of the virtual power plant in responding to the dispatching instruction.

Description

Technical Field

[0001] This invention relates to virtual power plant technology, and more particularly to a method and system for adjustable load aggregation and real-time response control in virtual power plants. Background Technology

[0002] With the transformation of energy structure and the intelligent development of power systems, virtual power plants, as a technological form that aggregates distributed energy resources to participate in grid regulation, are at their core the efficient aggregation and precise regulation of massive, heterogeneous, and distributed flexible load devices. In existing technologies, the aggregation and response regulation of adjustable loads by virtual power plants typically rely on pre-set aggregation models and fixed regulation strategies.

[0003] Current conventional practices commonly employ static aggregation methods based on load characteristic classification or statistical regularities. These methods often analyze historical data to group load devices with similar physical characteristics or response patterns into one category, and use fixed equivalent parameters (such as equivalent capacity and equivalent response time) to characterize the entire load group. During real-time control, the virtual power plant's central controller, based on received grid dispatch instructions, directly issues unified control signals to the aggregated load group, or performs load switching control according to preset priority rules. This method focuses on describing the macroscopic characteristics of the load group and assumes that the aggregation model remains stable during control.

[0004] However, the aforementioned conventional approaches have significant drawbacks. Because the dynamic behavior of load equipment is influenced by multiple factors, including its internal physical constraints, external environment, and user habits, its response characteristics are not static. Static aggregation models struggle to accurately depict the dynamic evolution of load groups during the response process, especially when there are discrepancies between dispatch commands and the actual operating state of the loads. Model mismatch can lead to decreased control accuracy and even cause oscillations in the load response or inadequate command execution. Furthermore, existing methods typically treat the load group as a whole for "black box" control, lacking a fine-grained perception and dynamic correction mechanism for changes in the state of individual loads and their inter-load coupling relationships. This results in insufficient adaptability and robustness of the control strategy, making it difficult to cope with the complex and ever-changing real-time power grid control requirements. Summary of the Invention

[0005] The present invention provides a method and system for adjustable load aggregation and real-time response control of virtual power plants, which can solve the problems in the prior art.

[0006] A first aspect of the present invention provides a method for adjustable load aggregation and real-time response control of a virtual power plant, comprising:

[0007] Real-time operation data and historical response data of multiple distributed flexible load devices are obtained. Based on the physical characteristic constraints and energy conversion relationships of the load devices reflected in the real-time operation data and the historical response data, a state-space model describing the dynamic behavior of a single load is constructed. The state-space model is decomposed into a fast response sub-model and a slow drift sub-model through a multi-scale time decoupling mechanism.

[0008] Extract the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combine the topological connection relationship of each load device in the power distribution network. By weighted aggregation of the state variables, multiple individual loads are aggregated into a group equivalent representation.

[0009] The deviation between the output state of the group equivalent representation and the virtual power plant dispatch command is calculated. The deviation is then propagated in reverse to the load state space model parameters of each individual unit constituting the group equivalent representation. The coupling coefficient between the fast response sub-model and the slow drift sub-model is dynamically corrected based on the deviation.

[0010] The modified fast response sub-model and slow drift sub-model are re-aggregated into an updated group equivalent representation. Based on the updated group equivalent representation, the flexible load group is simulated and predicted. Based on the simulation prediction results, scheduling decision support information is generated and sent to the corresponding load equipment for execution.

[0011] Based on the constraints of physical characteristics and energy conversion of load equipment reflected in the real-time operating data and the historical response data, a state-space model describing the dynamic behavior of individual loads is constructed, including:

[0012] Extract the power change curve, energy storage status change curve, and control input sequence of the load equipment from the real-time operating data and the historical response data;

[0013] Based on the power peak value, upper and lower limits of energy capacity and slope of state change reflected in the power change curve and the energy storage state change curve, the physical characteristic constraints of the load equipment are identified, and the value boundaries and rate of change limits of the state variables are determined based on the physical characteristic constraints.

[0014] Based on the energy conversion relationship, a coupling equation is established between the state variables. The value boundary and the rate of change limit are used as constraints of the coupling equation, and combined with the coupling equation to form a system of differential equations.

[0015] The system of differential equations is discretized in time, transforming the continuous-time system of differential equations into a discrete-time state transition matrix form.

[0016] The control input sequence is used as the input to the state transition matrix, and the power change curve is used as the expected output of the state transition matrix. By minimizing the error between the predicted output of the state transition matrix and the power change curve, a state-space model is constructed.

[0017] Based on the energy conversion relationship, a coupling equation is established between the state variables. The boundary values ​​and the rate of change limits are used as constraints on the coupling equation. The coupling equation is combined with the following to form a system of differential equations:

[0018] Based on the energy flow path inside the load equipment described in the energy conversion relationship, the causal transmission chain between the energy storage state variable, the power input state variable and the power output state variable is identified, and based on the conversion delay characteristics of each state variable in the causal transmission chain, a differential coupling equation describing the time evolution of the state variable is established.

[0019] In the differential coupling equation, a nonlinear correction term characterizing energy conversion loss is introduced, the value boundary is expressed as the magnitude constraint function of the state variable, the rate of change constraint is expressed as the slope constraint function of the time derivative of the state variable, and the magnitude constraint function and the slope constraint function are embedded as penalty terms into the objective functional corresponding to the differential coupling equation.

[0020] By using the variational method to find the extremum conditions of the objective functional containing the penalty term, the objective functional is transformed into a system of differential equations in Euler-Lagrange form.

[0021] Extracting the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combining the topological connection relationship of each load device in the power distribution network, the state variables are weighted and aggregated to form a group equivalent representation, including:

[0022] Extract the fast state variables of each load device from the fast response sub-model, and extract the slow state variables of each load device from the slow drift sub-model;

[0023] Based on the topological connection relationship of each load device in the power distribution network, a topological adjacency matrix reflecting the electrical coupling strength between load devices is constructed, and the contribution weight of each load device to the group response characteristics is calculated based on the topological adjacency matrix.

[0024] The fast state variables are weighted and summed according to the contribution weight to obtain the fast equivalent state of the group, and the slow state variables are weighted and summed according to the contribution weight to obtain the slow equivalent state of the group.

[0025] The fast equivalent state of the group and the slow equivalent state of the group are combined into a group equivalent state vector, and a group equivalent representation is constructed based on the group equivalent state vector.

[0026] Calculating the deviation between the output state of the group equivalent representation and the virtual power plant dispatch command, and propagating the deviation backward to the individual load state space model parameters constituting the group equivalent representation includes:

[0027] Obtain the desired power trajectory and desired response time contained in the virtual power plant dispatch instruction, calculate the power deviation between the output state of the group equivalent representation and the desired power trajectory, and calculate the time delay deviation between the actual response time and the desired response time.

[0028] Construct an inverse mapping link from the output state of the group equivalent representation to the parameters of the state space model of each individual load. Based on the sensitivity coefficient of each individual load to the output state of the group in the inverse mapping link, decompose the power deviation and the time delay deviation according to the sensitivity coefficient to obtain the parameter correction amount allocated to the state space model of each individual load.

[0029] The parameter correction is applied to the system matrix and control matrix in the state-space model of the corresponding unit load, and the state-space model parameters are updated.

[0030] Based on the updated equivalent representation of the group, simulation prediction is performed on the flexible load group. Scheduling decision support information is generated based on the simulation prediction results and sent to the corresponding load equipment for execution, including:

[0031] Using the sequence of scheduling instructions for multiple scenarios within a future time window as input, and based on the updated group equivalent representation, we perform multi-trajectory parallel simulation and deduction to obtain the set of predicted response trajectories under different scenarios.

[0032] Extract the group response envelope surface corresponding to each scenario from the predicted response trajectory set, and determine the robust adjustable domain of the flexible load group within the future time window by solving the intersection region of the group response envelope surface;

[0033] A multi-objective optimization model is constructed within the robust adjustable domain, with the objectives of minimizing power grid dispatching costs and minimizing load equipment response deviations. The Pareto optimal dispatching strategy set is obtained by solving the multi-objective optimization model, and the current optimal strategy is selected from the Pareto optimal dispatching strategy set based on the real-time operating status of the power grid.

[0034] The current optimal strategy is decomposed into the target power trajectory and execution timing of each individual load, and combined with the physical constraint boundaries of each individual load, scheduling decision support information containing trajectory tracking targets and constraint tolerance ranges is generated and sent to the corresponding load devices for execution through the communication interface.

[0035] Extracting the population response envelope surface corresponding to each scenario from the predicted response trajectory set, and determining the robust adjustable domain of the flexible load population within the future time window by solving the intersection region of the population response envelope surface, includes:

[0036] Extract the extreme points of the group response power at each time under each scenario from the predicted response trajectory set, and connect the extreme points of the group response power in the power-time coordinate space to form the response trajectory boundary curves of each scenario.

[0037] The response trajectory boundary curves of each scenario are subjected to convex hull processing, and the power-time point set of the response trajectory boundary curve envelope is mapped into a convex polyhedron to obtain the group response envelope surface corresponding to each scenario.

[0038] The group response envelopes of multiple scenarios are geometrically intersected in the power-time coordinate space. The intersection interval of the group response envelopes in the power dimension is solved time by time to obtain the power range that all scenarios can reach at each time.

[0039] The power ranges at each moment are combined according to the time series to form a power reachable set covering future time windows, and the upper and lower bounds of the power reachable set are used as the boundaries of the robust adjustable domain.

[0040] A second aspect of the present invention provides a virtual power plant adjustable load aggregation and real-time response control system, comprising:

[0041] The data acquisition unit is used to acquire real-time operating data and historical response data of multiple distributed flexible load devices. Based on the physical characteristic constraints and energy conversion relationships of the load devices reflected in the real-time operating data and the historical response data, a state-space model describing the dynamic behavior of a single load is constructed. The state-space model is decomposed into a fast response sub-model and a slow drift sub-model through a multi-scale time decoupling mechanism.

[0042] The model building unit is used to extract the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combine the topological connection relationship of each load device in the power distribution network to aggregate multiple individual loads into a group equivalent representation by weighted aggregation of the state variables;

[0043] The state aggregation unit is used to calculate the deviation between the output state of the group equivalent representation and the virtual power plant dispatch command, and to propagate the deviation to the state space model parameters of each individual load constituting the group equivalent representation, and to dynamically correct the coupling coefficient between the fast response sub-model and the slow drift sub-model based on the deviation.

[0044] The deviation correction unit is used to re-aggregate the corrected fast response sub-model and the slow drift sub-model into an updated group equivalent representation, perform simulation prediction on the flexible load group based on the updated group equivalent representation, generate scheduling decision support information based on the simulation prediction results, and send it to the corresponding load equipment for execution.

[0045] A third aspect of the present invention provides an electronic device, comprising:

[0046] processor;

[0047] Memory used to store processor-executable instructions;

[0048] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0049] 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.

[0050] The beneficial effects of this application are as follows:

[0051] This method significantly improves the aggregation accuracy and real-time control capability of virtual power plants for distributed flexible loads. By constructing a state-space model based on physical characteristics and energy conversion relationships, and using a multi-scale time decoupling mechanism to decompose it into fast and slow sub-models, it can accurately characterize the dynamic behavior of load equipment at different time scales, effectively separating the load's fast response capability from the slow energy accumulation process, thus laying a model foundation for refined control.

[0052] Weighted aggregation of state variables based on load equipment topology connections results in a group equivalent representation that not only reflects the overall dynamics of the aggregate but also retains its correlation with the physical connection to the power grid. This gives the aggregation model clear physical meaning and network interpretability, facilitating interaction with the power grid dispatching system. By calculating the deviation between the group output and dispatching commands and propagating it back to the individual model parameters to dynamically correct the coupling coefficients between sub-models, adaptive online calibration of the aggregation model is achieved. This enables real-time tracking of load equipment characteristic changes and external environmental disturbances, ensuring the accuracy of the aggregation model under all operating conditions.

[0053] By using the updated aggregation model for simulation prediction and generating scheduling decisions, the formulation of control instructions becomes forward-looking and predictive, and the effectiveness of different scheduling strategies can be evaluated. This allows for the issuance of safer and more economical control instructions, directly improving the agility, reliability, and economy of virtual power plants in responding to grid demands. Attached Figure Description

[0054] Figure 1This is a flowchart illustrating the adjustable load aggregation and real-time response control method for virtual power plants according to an embodiment of the present invention.

[0055] Figure 2 This is a scheduling flowchart for multi-scenario group response capability prediction and evaluation in an embodiment of the present invention. Detailed Implementation

[0056] 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.

[0057] 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.

[0058] Figure 1 This is a flowchart illustrating the adjustable load aggregation and real-time response control method for virtual power plants according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0059] Real-time operation data and historical response data of multiple distributed flexible load devices are obtained. Based on the physical characteristic constraints and energy conversion relationships of the load devices reflected in the real-time operation data and the historical response data, a state-space model describing the dynamic behavior of a single load is constructed. The state-space model is decomposed into a fast response sub-model and a slow drift sub-model through a multi-scale time decoupling mechanism.

[0060] Extract the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combine the topological connection relationship of each load device in the power distribution network. By weighted aggregation of the state variables, multiple individual loads are aggregated into a group equivalent representation.

[0061] The deviation between the output state of the group equivalent representation and the virtual power plant dispatch command is calculated. The deviation is then propagated in reverse to the load state space model parameters of each individual unit constituting the group equivalent representation. The coupling coefficient between the fast response sub-model and the slow drift sub-model is dynamically corrected based on the deviation.

[0062] The modified fast response sub-model and slow drift sub-model are re-aggregated into an updated group equivalent representation. Based on the updated group equivalent representation, the flexible load group is simulated and predicted. Based on the simulation prediction results, scheduling decision support information is generated and sent to the corresponding load equipment for execution.

[0063] In one optional implementation, based on the physical characteristic constraints and energy conversion relationships of the load equipment reflected in the real-time operating data and the historical response data, a state-space model describing the dynamic behavior of a single load is constructed, including:

[0064] Extract the power change curve, energy storage status change curve, and control input sequence of the load equipment from the real-time operating data and the historical response data;

[0065] Based on the power peak value, upper and lower limits of energy capacity and slope of state change reflected in the power change curve and the energy storage state change curve, the physical characteristic constraints of the load equipment are identified, and the value boundaries and rate of change limits of the state variables are determined based on the physical characteristic constraints.

[0066] Based on the energy conversion relationship, a coupling equation is established between the state variables. The value boundary and the rate of change limit are used as constraints of the coupling equation, and combined with the coupling equation to form a system of differential equations.

[0067] The system of differential equations is discretized in time, transforming the continuous-time system of differential equations into a discrete-time state transition matrix form.

[0068] The control input sequence is used as the input to the state transition matrix, and the power change curve is used as the expected output of the state transition matrix. By minimizing the error between the predicted output of the state transition matrix and the power change curve, a state-space model is constructed.

[0069] In establishing a state-space model for distributed flexible load equipment, real-time operational data streams containing timestamps, power measurements, and temperature sensor readings are acquired from the data acquisition system. Simultaneously, response records for the past 30 days stored in the historical database are retrieved. Taking an energy storage air conditioner as an example, the power change curve records the power fluctuation characteristics of the compressor during start-up and shutdown. The curve shows that the power jumps from 0.2kW to 3.5kW at startup, maintains a steady-state power of 2.8kW, and decreases to 0.15kW in standby mode after shutdown. The energy storage state change curve reflects the evolution of indoor temperature. When the compressor is working, the temperature decreases at a rate of 0.3℃ per minute, and after shutdown, it recovers at a rate of 0.12℃ per minute due to environmental heat exchange. The control input sequence records the trigger time and duration of scheduling commands, including timestamps for start-up commands, power adjustment commands, and shutdown commands.

[0070] The peak parameters extracted from the power change curve indicate that the rated power of the device is 3.5kW. The temperature control range identified from 18℃ to 26℃ is determined from the energy storage state change curve. These values ​​constitute the core boundary of the physical characteristic constraints. The state change slope is obtained by numerical differentiation calculation of the curve. The maximum slope during the temperature drop phase is -0.35℃ / min, and the maximum slope during the temperature rise phase is +0.15℃ / min. These slopes limit the rate of change constraints of the corresponding state variables. The power state variable P(t) is set to a range of [0.15, 3.5]kW, and the temperature state variable T(t) is set to a range of [18, 26]℃. The rate of temperature change is constrained within the interval [-0.35, 0.15]℃ / min.

[0071] The energy conversion relationship is manifested as the conversion of electrical energy into heat energy, and a coupled equation is established based on the first law of thermodynamics. Changes in indoor temperature are influenced by three factors: compressor cooling power, ambient temperature, and heat conduction through the building envelope. This can be expressed as the derivative of temperature with respect to time equals the product of cooling power and energy efficiency ratio divided by heat capacity, minus the product of temperature difference and thermal conductivity divided by heat capacity. Changes in power state are driven by control commands. When a start-up command is received, the power rises exponentially to the rated value; when a shutdown command is triggered, the power decays exponentially to standby power consumption. The temperature boundary values ​​are transformed into inequality constraints. The rate of change limit is expressed as These constraints, together with the coupled equations, constitute a system of differential equations containing two differential equations and four inequality constraints.

[0072] The time discretization process employs the Euler forward difference method with a sampling period of 1 minute, approximating the temperature derivative as... ,in Minutes. The discretized equation is rewritten as follows: , where f represents the coupling function on the right and u_k is the control input. The discrete form of the power state is: The coefficients α and β are determined by fitting the start-stop characteristics of the power curve. The discrete equations described above are organized into a matrix form, and a state vector is defined. Input vector u_k, establish state transition relationship The output equation is as follows The output vector y_k corresponds to the measured power value.

[0073] During model parameter calibration, the control input sequence from historical response data is used as the input sample of u_k, and the corresponding power measurement value is used as the desired output. The predicted output is calculated using a state-space model. The error function is constructed as the sum of squares of the differences between the predicted and expected outputs at all times. The gradient descent method is used to iteratively adjust the elements of matrices A, B, and C until the error function converges to its minimum. After 150 iterations, the root mean square error of the model on the validation set is reduced to 0.08 kW, and the goodness of fit between the predicted output curve and the actual power curve reaches 0.96. The state-space model with completed parameter calibration accurately reflects the dynamic response characteristics of the air conditioning load when receiving different control commands, providing a precise mathematical representation of the individual load for subsequent multi-scale time decoupling and group aggregation.

[0074] In one optional implementation, a coupling equation is established between the state variables based on the energy conversion relationship, and the value boundaries and the rate of change limits are used as constraints on the coupling equation. The coupling equation is then combined with these constraints to form a system of differential equations, including:

[0075] Based on the energy flow path inside the load equipment described in the energy conversion relationship, the causal transmission chain between the energy storage state variable, the power input state variable and the power output state variable is identified, and based on the conversion delay characteristics of each state variable in the causal transmission chain, a differential coupling equation describing the time evolution of the state variable is established.

[0076] In the differential coupling equation, a nonlinear correction term characterizing energy conversion loss is introduced, the value boundary is expressed as the magnitude constraint function of the state variable, the rate of change constraint is expressed as the slope constraint function of the time derivative of the state variable, and the magnitude constraint function and the slope constraint function are embedded as penalty terms into the objective functional corresponding to the differential coupling equation.

[0077] By using the variational method to find the extremum conditions of the objective functional containing the penalty term, the objective functional is transformed into a system of differential equations in Euler-Lagrange form.

[0078] When establishing coupling equations between state variables based on energy conversion relationships, it is necessary to clarify the energy flow path within the load equipment. Taking an electric water heater as an example, electrical energy input is converted into heat energy by the heating element and stored in the water tank. The water tank temperature serves as the energy storage state variable, the electrical power input as the power input state variable, and the heat loss power as the power output state variable. In an air conditioning system, electrical energy drives the compressor to do work, and the refrigerant circulation realizes heat transfer. The indoor temperature serves as the energy storage state variable, the compressor power as the power input state variable, and the heat exchange power through the building envelope as the power output state variable. Identifying the causal transmission chain between these state variables requires tracing the physical process of energy flow. Changes in the power input state variable directly affect the accumulation rate of the energy storage state variable, and the current value of the energy storage state variable determines the release intensity of the power output state variable.

[0079] The transition time delay characteristics of each state variable in a causal transmission chain reflect the physical inertia of the energy conversion process. The water tank temperature of an electric water heater responds to the electrical power input with a heat capacity delay, typically ranging from 300 to 2000 seconds, depending on the ratio of tank volume to heating power. The indoor temperature of an air conditioner responds to the compressor power with a combined effect of room heat capacity and heat exchanger heat transfer delay, typically ranging from 5 to 30 minutes. Establishing differential coupled equations describing the time evolution of state variables requires expressing the law of conservation of energy in differential form. The derivative of the energy storage state variable with respect to time equals the power input state variable minus the power output state variable, divided by the equivalent heat capacity parameter. The unit of the equivalent heat capacity parameter is joules per unit of state variable. For temperature-type state variables, the equivalent heat capacity of an electric water heater is the specific heat capacity of water multiplied by the mass of water; the equivalent heat capacity of an air-conditioned room is the specific heat capacity of air multiplied by the mass of air plus the heat capacity of the building envelope. The power output state variable is usually expressed as the difference between the energy storage state variable and the environmental reference value multiplied by the heat transfer coefficient. The unit of the heat transfer coefficient is watts per unit state variable. The heat transfer coefficient of an electric water heater reflects the heat preservation performance of the water tank, with a typical value of 0.5 to 3 watts per degree Celsius. The heat transfer coefficient of an air-conditioned room reflects the overall thermal resistance of the building envelope, with a typical value of 50 to 300 watts per degree Celsius.

[0080] Introducing a nonlinear correction term characterizing energy conversion losses into the differential coupling equations is to accurately depict efficiency reduction in actual physical processes. The efficiency of the heating element in an electric water heater is not constant; it decreases when the power input is low or the water temperature is close to the boiling point. The nonlinear correction term is expressed as the power input state variable multiplied by an efficiency function dependent on the energy storage state variable. This efficiency function can be approximated using a piecewise linear approach or an exponential decay form. The cooling efficiency of an air conditioning system is significantly affected by the indoor-outdoor temperature difference. The nonlinear correction term is expressed as the power input state variable multiplied by an energy efficiency ratio function dependent on the energy storage state variable and the difference between the ambient temperature and the ambient temperature. The energy efficiency ratio function decreases as the temperature difference increases, with a typical decay slope of 0.05 to 0.15 per degree Celsius. The boundary conditions are expressed as amplitude constraint functions for the state variables. For temperature-type state variables, the water temperature of an electric water heater must meet the following conditions: a lower bound not lower than the ambient temperature plus 5 degrees Celsius and an upper bound not higher than the safety valve opening temperature minus 10 degrees Celsius. The room temperature of an air-conditioned room must meet the following conditions: a lower bound not lower than the set temperature minus 2 degrees Celsius and an upper bound not higher than the set temperature plus 2 degrees Celsius. The rate of change constraint is expressed as the slope constraint function of the derivative of the state variable with respect to time. The rate of change of water temperature in electric water heaters is limited by the ratio of heating power to heat capacity, with a typical upper limit of 0.02 to 0.1 degrees Celsius per second. The rate of change of room temperature in air-conditioned rooms is limited by the ratio of cooling power to room heat capacity, with a typical upper limit of 0.005 to 0.03 degrees Celsius per second.

[0081] When embedding amplitude and slope constraint functions as penalty terms into the objective functional corresponding to the differential coupling equation, a penalty function reflecting the degree of constraint violation needs to be constructed. The penalty term for amplitude constraints uses a weighted sum of the squares of the portions of the state variables exceeding their upper and lower bounds. The weighting coefficients reflect the rigidity of the constraint; larger weights (e.g., 1000 to 10000) are set for safety constraints, while smaller weights (e.g., 10 to 100) are set for comfort constraints. The penalty term for slope constraints uses a weighted sum of the squares of the portions of the time derivative of the state variables exceeding their upper and lower bounds. The weighting coefficients reflect the strength of the dynamic response constraint; larger weights (e.g., 500 to 5000) are set for equipment protection constraints, while smaller weights (e.g., 5 to 50) are set for user experience constraints. The objective functional is expressed as the integral of the time evolution trajectory of the state variables over a future time window. The integrand includes a quadratic term describing operating costs, a quadratic term describing tracking errors, and the aforementioned penalty term. The upper and lower bounds of the integral correspond to the start and end times of the time window, with typical window lengths ranging from 15 minutes to 4 hours.

[0082] When finding the extremum condition of the objective functional including the penalty term using the variational method, it is necessary to apply a small perturbation to the trajectory of the state variable and calculate the first-order variation of the objective functional. The necessary condition for the first-order variation to be zero derives the Euler-Lagrange equation. The left side of this equation contains the partial derivative of the integrand with respect to the state variable minus the partial derivative of the integrand with respect to the time of the state variable, and the right side is zero. The partial derivative of the integrand with respect to the state variable includes the linear coefficients of the running cost term and the tracking error term with respect to the state variable, as well as the nonlinear derivative of the penalty term with respect to the state variable. The nonlinear derivative is zero when the state variable is inside the constraint boundary; when it is near the constraint boundary, it is twice the difference between the state variable and the boundary value multiplied by the weighting coefficient. The partial derivative of the integrand with respect to the time derivative of the state variables includes the nonlinear derivative of the penalty term with respect to the time derivative. This derivative is 0 when the time derivative is inside the slope constraint, and is twice the difference between the time derivative and the constraint value multiplied by the weighting coefficient when it is near the slope constraint boundary. Differentiating this partial derivative with respect to time yields the term with respect to the second-order time derivative of the state variables. The differential coupling equations are combined with the Euler-Lagrange equations to form a system of differential equations. This system contains first-order and second-order time derivative equations for the state variables. The order of the system is equal to the number of state variables multiplied by 2, typically 2 to 6 for a single load device.

[0083] In one optional implementation, the state variables of each load device in the fast response sub-model and the slow drift sub-model are extracted. Combining the topological connections of each load device in the distribution network, the state variables are weighted and aggregated to form a group equivalent representation, including:

[0084] Extract the fast state variables of each load device from the fast response sub-model, and extract the slow state variables of each load device from the slow drift sub-model;

[0085] Based on the topological connection relationship of each load device in the power distribution network, a topological adjacency matrix reflecting the electrical coupling strength between load devices is constructed, and the contribution weight of each load device to the group response characteristics is calculated based on the topological adjacency matrix.

[0086] The fast state variables are weighted and summed according to the contribution weight to obtain the fast equivalent state of the group, and the slow state variables are weighted and summed according to the contribution weight to obtain the slow equivalent state of the group.

[0087] The fast equivalent state of the group and the slow equivalent state of the group are combined into a group equivalent state vector, and a group equivalent representation is constructed based on the group equivalent state vector.

[0088] After obtaining the state-space model of the individual load and completing multi-scale time decoupling, it is necessary to perform aggregate modeling of the distributed flexible load devices. For the fast response sub-model, fast state variables of each load device are extracted. These variables describe the power response characteristics of the load device within a time scale of seconds to minutes, such as the start-stop state of an air conditioner compressor and the rate of change of heating power of an electric water heater. For the slow drift sub-model, slow state variables reflecting the thermal capacity characteristics or energy storage state of the load device are extracted, such as the state of charge of an energy storage device and the temperature drift caused by the thermal inertia of the building. The change period of these variables is usually on the order of tens of minutes to hours.

[0089] After extracting the state variables, a topological adjacency matrix reflecting the electrical coupling relationships between load devices is constructed based on the actual topology of the distribution network. The elements of this matrix reflect the electrical distance or network impedance between two load devices. Load devices located in the same distribution area have stronger electrical coupling, and their corresponding matrix elements have larger values; load devices on different feeders have weaker coupling, and their corresponding element values ​​are close to zero. The contribution weight of each load device to the group response characteristics is calculated by performing eigenvalue decomposition or graph-based centrality analysis on the topological adjacency matrix. This weight comprehensively considers the installed capacity of the load device, its importance in the network location, and its potential for coordinated response with other devices.

[0090] After weighting, the fast state variables of each load device are summed according to their corresponding contribution weights. If there are N load devices in the system, the fast state variable of the i-th device is... If the corresponding weight is w_i, then the group's fast equivalent state is represented as the sum of the products of the fast state variables of each device and their weights. Similarly, the slow state variables of each load device are weighted and summed with equal weights to obtain the group's slow equivalent state. This weighted aggregation method ensures that the equivalent representation at the group level accurately reflects the actual contribution of each individual load.

[0091] The calculated fast and slow equivalent states of the group are concatenated column-wise to form a group equivalent state vector. This vector has a dimension much lower than the total dimension of all individual load state variables, achieving a dimensionality reduction mapping from a high-dimensional discrete state space to a low-dimensional group state space. Based on this group equivalent state vector, a group equivalent representation model capable of describing the dynamic response characteristics of the entire flexible load group is constructed. This model retains the dual-timescale characteristics of fast response and slow drift, while condensing the complex interactions of hundreds or even thousands of individual loads into a unified equivalent mathematical expression, providing a computational foundation for subsequent deviation correction and scheduling decisions.

[0092] In one optional implementation, calculating the deviation between the output state of the group equivalent representation and the virtual power plant dispatch command, and propagating the deviation backward to the state-space model parameters of each individual load constituting the group equivalent representation includes:

[0093] Obtain the desired power trajectory and desired response time contained in the virtual power plant dispatch instruction, calculate the power deviation between the output state of the group equivalent representation and the desired power trajectory, and calculate the time delay deviation between the actual response time and the desired response time.

[0094] Construct an inverse mapping link from the output state of the group equivalent representation to the parameters of the state space model of each individual load. Based on the sensitivity coefficient of each individual load to the output state of the group in the inverse mapping link, decompose the power deviation and the time delay deviation according to the sensitivity coefficient to obtain the parameter correction amount allocated to the state space model of each individual load.

[0095] The parameter correction is applied to the system matrix and control matrix in the state-space model of the corresponding unit load, and the state-space model parameters are updated.

[0096] The deviation between the output state of the equivalent representation of the load group and the virtual power plant dispatch instructions requires extracting key information from the virtual power plant dispatch instructions. The virtual power plant dispatch instructions are issued in the form of structured data messages, which include an instruction identifier field, a timestamp field, a desired power trajectory field, and a desired response time field. The desired power trajectory field is stored as a time-power value pair sequence, with a time resolution configurable from 5 to 60 seconds, and power values ​​in kilowatts with precision retained to two decimal places. The desired response time field is stored as the time interval from the moment the instruction is issued until the load group reaches the desired power, in seconds with precision retained to integers. Obtaining the desired power trajectory requires parsing the time-power value pair sequence and storing the power values ​​at each moment in the sequence into an array of length N. Typical values ​​for N are 12 to 720, corresponding to a trajectory window of 1 minute to 1 hour.

[0097] The output state of the group equivalent representation is calculated in real time by the group equivalent representation model. The output state includes two components: group aggregated power and group response delay. Group aggregated power represents the total active power of the flexible load group at the current moment, in kilowatts. Group response delay represents the time interval from receiving the dispatch command to the flexible load group reaching its steady-state power, in seconds. Calculating the power deviation between the output state of the group equivalent representation and the expected power trajectory requires aligning the time base. The time series of the group aggregated power and the time series of the expected power trajectory are aligned according to timestamps. After alignment, the difference between the group aggregated power and the expected power is calculated moment by moment, and the difference is stored in a power deviation array of length N. The positive or negative sign of the power deviation reflects whether the actual power is higher or lower than the expected power. A typical deviation range is -500 kW to +500 kW, corresponding to a group size of 50 to 200 individual loads. Calculating the delay deviation between the actual response time and the expected response time requires extracting the actual response time from the group response delay. Subtracting the expected response time from the actual response time yields the delay deviation, with a typical deviation range of -30 seconds to +60 seconds.

[0098] Constructing an inverse mapping link from the group's equivalent representation output state to the state-space model parameters of each individual load requires clarifying the dependency relationship between the group's output state and the individual model parameters. The group's aggregate power equals the sum of the output powers of each individual load. The output power of each individual load is calculated from the output equation of its state-space model, which includes the product of the output matrix and the state variables. The group's response delay is determined by the weighted average of the response delays of each individual load, with the weighting coefficient proportional to the rated power of each individual load. The state-space model parameters include a system matrix and a control matrix. The elements of the system matrix reflect the coupling relationships between state variables, while the elements of the control matrix reflect the strength of the control input's effect on the state variables. The inverse mapping link is expressed as the partial derivatives of the power deviation with respect to the elements of the system matrix and the control matrix, as well as the partial derivatives of the delay deviation with respect to the elements of the system matrix and the control matrix.

[0099] Based on the sensitivity coefficients of each individual load in the reverse mapping link to the group output state, the partial derivatives need to be calculated. The sensitivity coefficient reflects the magnitude of the change in the group output state caused by a small change in the individual load model parameters; a larger sensitivity coefficient indicates a more significant impact of that individual load on the group output state. The sensitivity coefficient is calculated using a numerical differentiation method. A perturbation of 0.01 times the absolute value of the element is applied to the system matrix element. The difference in the group output state before and after the perturbation is divided by the perturbation magnitude to obtain an approximate value of the sensitivity coefficient. For power deviation, the sensitivity coefficient of each individual load is obtained by chain differentiation of the output equation. For time delay deviation, the sensitivity coefficient of each individual load is equal to the ratio of the rated power of that individual load to the total rated power of the group multiplied by the partial derivative of the response time delay with respect to the system matrix element.

[0100] Decomposing power deviation and time delay deviation according to sensitivity coefficients requires solving a system of linear equations. The deviation value at each time step in the power deviation array is expressed as the sum of the products of the correction amount of each individual load parameter and the corresponding sensitivity coefficient. The time delay deviation is also expressed as the sum of the products of the correction amount of each individual load parameter and the corresponding sensitivity coefficient. The coefficient matrix of this system of linear equations consists of the sensitivity coefficients of each individual load at each time step. The number of rows in the matrix is ​​equal to the length N of the power deviation array plus 1, and the number of columns is equal to the number of individual loads M multiplied by the dimension D of the state-space model parameters. For a second-order state-space model, D is 6 to 8. The least squares method is used to solve the system of linear equations. The least squares solution is obtained by performing singular value decomposition on the sensitivity coefficient matrix. Singular value decomposition expresses the sensitivity coefficient matrix as the product of the left singular matrix, the singular value diagonal matrix, and the transpose of the right singular matrix. The parameter correction amount is equal to the right singular matrix multiplied by the inverse of the singular value diagonal matrix multiplied by the transpose of the left singular matrix multiplied by the deviation vector. When calculating the inverse of a singular value diagonal matrix, singular values ​​smaller than a threshold need to be truncated. The threshold is set to 0.001 to 0.01 times the maximum singular value.

[0101] After obtaining the parameter corrections allocated to the state-space models of each individual load, these corrections need to be applied to the system matrix and control matrix in the corresponding individual load state-space models. The parameter corrections are vectors of length D, with the first part corresponding to the corrections for elements in the system matrix and the second part corresponding to the corrections for elements in the control matrix. The system matrix elements are arranged in row-major order as a one-dimensional vector, and the control matrix elements are also arranged in row-major order as a one-dimensional vector. When applying the parameter corrections to the system matrix, the correction is read element-by-element and added to the corresponding element of the system matrix. The updated system matrix element equals the original element plus the correction. Similarly, when applying the parameter corrections to the control matrix, the correction is read element-by-element and added to the corresponding element of the control matrix. The updated control matrix element equals the original element plus the correction. After updating the state-space model parameters, a rationality check is required. The real parts of the eigenvalues ​​of the system matrix must be negative to ensure stability, and the absolute values ​​of the control matrix elements must be within the physically realizable range. If the updated parameters do not meet the rationality constraints, the parameter correction amount is scaled, with the scaling factor decreasing in increments of 0.1 from 1 until the updated parameters meet the rationality constraints.

[0102] In one optional implementation, the process of simulating and predicting the flexible load group based on the updated group equivalent representation, generating scheduling decision support information based on the simulation prediction results, and sending it to the corresponding load equipment for execution includes:

[0103] Using the sequence of scheduling instructions for multiple scenarios within a future time window as input, and based on the updated group equivalent representation, we perform multi-trajectory parallel simulation and deduction to obtain the set of predicted response trajectories under different scenarios.

[0104] Extract the group response envelope surface corresponding to each scenario from the predicted response trajectory set, and determine the robust adjustable domain of the flexible load group within the future time window by solving the intersection region of the group response envelope surface;

[0105] A multi-objective optimization model is constructed within the robust adjustable domain, with the objectives of minimizing power grid dispatching costs and minimizing load equipment response deviations. The Pareto optimal dispatching strategy set is obtained by solving the multi-objective optimization model, and the current optimal strategy is selected from the Pareto optimal dispatching strategy set based on the real-time operating status of the power grid.

[0106] The current optimal strategy is decomposed into the target power trajectory and execution timing of each individual load, and combined with the physical constraint boundaries of each individual load, scheduling decision support information containing trajectory tracking targets and constraint tolerance ranges is generated and sent to the corresponding load devices for execution through the communication interface.

[0107] Simulation prediction of flexible load groups based on the updated group equivalent representation requires a sequence of multi-scenario dispatch instructions within a future time window as input. The future time window length is set from 15 minutes to 4 hours, with a typical value of 1 hour. The multi-scenario dispatch instruction sequence contains 3 to 10 dispatch instructions, with a typical number of scenarios being 5. Each scenario reflects the uncertainty of grid dispatch demand. When performing multi-trajectory parallel simulation based on the updated group equivalent representation, the group equivalent representation model is called for each scenario. The state evolution trajectory is calculated by solving the group equivalent differential equation system, and the numerical sequence of the group aggregated power changing over time is output. The multi-trajectory parallel simulation adopts multi-threaded concurrent execution, with each scenario assigned to an independent thread, resulting in a set of predicted response trajectories under different scenarios.

[0108] Extracting the collective response envelope for each scenario from the predicted response trajectory set requires geometrization of the trajectories. Each scenario trajectory is represented as a curve in the power-time coordinate space, and the collective response envelope is defined as the smallest convex polyhedron containing this curve. Convex hull transformation employs an incremental construction method, checking point by point whether a new point lies outside the constructed convex hull. If so, it is added to the convex hull, and the boundary vertex set is updated; the typical number of vertices ranges from 10 to 100. Geometric intersection is used to solve for the intersection region of the collective response envelopes, traversing the boundary vertices of each envelope and determining whether each vertex lies simultaneously inside all envelopes. This determination is made using a system of linear inequalities; if all inequalities are satisfied after substituting the vertex coordinates into the system, the vertex lies inside. The typical number of boundary vertices in the intersection region ranges from 5 to 50, and the robust adjustable domain is expressed as a set of inequality constraints related to time and power.

[0109] Constructing a multi-objective optimization model within the robustly adjustable domain requires defining optimization variables and objective functions. The optimization variable is the power dispatch trajectory, expressed as a sequence of time-power numerical pairs. The objective function for minimizing grid dispatch cost is defined as the integral of the product of the power dispatch trajectory and the electricity price curve, while the objective function for minimizing load equipment response deviation is defined as the square integral of the difference between the power dispatch trajectory and the desired power trajectory. Constraints include inequality constraints within the robustly adjustable domain and power change rate limits, typically with an upper bound of 0.5 kW to 5 kW per second. A non-dominated sorting genetic algorithm is used to solve the multi-objective optimization model, with a population size of 50 to 200 and an iteration number of 50 to 500, resulting in a Pareto-optimal dispatch strategy set containing 20 to 100 strategies. The current optimal strategy is selected based on the real-time operating status of the grid. The selection criterion is defined as the weighted sum of the two objective functions, with the weight coefficients dynamically adjusted according to the grid load level. When the load level is higher than 90% of the rated load, the dispatch cost weight is 0.8 and the response deviation weight is 0.2.

[0110] The current optimal strategy is decomposed into a mapping relationship between the target power trajectory of each individual load and its execution time sequence. The group power scheduling trajectory is equal to the sum of the power scheduling trajectories of each individual load. The target power trajectory of each individual load is obtained by solving a linear allocation problem. The objective function of the linear allocation problem is the sum of the squares of the differences between the power trajectory of each individual load and its historical average power trajectory. The constraint condition is that the sum of the trajectories equals the group trajectory and satisfies the upper and lower bounds of the individual power. The linear allocation problem is solved using a quadratic programming solver with the interior-point method, and the iterative convergence accuracy is set to 0.0001 to 0.01. The execution time sequence of each individual load is defined as the time sequence of power value changes. Scheduling decision support information is generated by combining the physical constraint boundaries of each individual load. The physical constraint boundaries include the upper and lower bounds of power and the upper and lower bounds of the rate of change. The trajectory tracking target is the target power trajectory of each individual load, and the constraint tolerance range is the margin between the physical constraint boundary and the target power trajectory, with a typical margin range of 0.1 kW to 1 kW. The scheduling decision support information is encapsulated in the form of structured data messages. The messages include individual load identifiers, target power trajectories, execution timing, and constraint tolerance range fields, and are sent to the corresponding load devices for execution through the communication interface.

[0111] In one optional implementation, extracting the group response envelope surface corresponding to each scenario from the predicted response trajectory set, and determining the robust adjustable domain of the flexible load group within the future time window by solving the intersection region of the group response envelope surface includes:

[0112] Extract the extreme points of the group response power at each time under each scenario from the predicted response trajectory set, and connect the extreme points of the group response power in the power-time coordinate space to form the response trajectory boundary curves of each scenario.

[0113] The response trajectory boundary curves of each scenario are subjected to convex hull processing, and the power-time point set of the response trajectory boundary curve envelope is mapped into a convex polyhedron to obtain the group response envelope surface corresponding to each scenario.

[0114] The group response envelopes of multiple scenarios are geometrically intersected in the power-time coordinate space. The intersection interval of the group response envelopes in the power dimension is solved time by time to obtain the power range that all scenarios can reach at each time.

[0115] The power ranges at each moment are combined according to the time series to form a power reachable set covering future time windows, and the upper and lower bounds of the power reachable set are used as the boundaries of the robust adjustable domain.

[0116] like Figure 2 As shown, the method includes:

[0117] After obtaining the predicted response trajectory set for multiple scenarios, it is necessary to extract the actual control boundary of the flexible load group from the trajectory data. For the time-series response trajectory generated for each scenario, the total power value of the group corresponding to each sampling time is located in the power-time two-dimensional coordinate space. By traversing the power data points of all times under a certain scenario, the maximum and minimum sampling points of the trajectory in the power dimension of that scenario are identified. These extreme points are connected sequentially in time order to form the response trajectory boundary curve of that scenario. This boundary curve essentially outlines the outer contour of the load group's response capability under that scenario.

[0118] Because the start-up and shutdown switching of actual load equipment and power ramping are affected by physical constraints, the boundary curves formed by simply connecting discrete extreme points exhibit a non-convex structure, leading to numerical instability in subsequent intersection calculations. Therefore, a convex hull transformation algorithm is used to process the boundary curves of each scenario. In practice, all power-time sampling points within the boundary curve envelope are considered as a spatial point set, and a minimal convex polyhedron containing all point sets is constructed in two-dimensional space using the convex hull calculation algorithm. The projection of this convex polyhedron onto the power-time plane is the convexized group response envelope surface, whose vertices are composed of points from the original extreme points located at the convex hull boundary, while the internal points are covered by the envelope surface.

[0119] After convexifying the envelope surfaces of each scenario, geometric intersection operations are performed on the convex polyhedra corresponding to multiple scenarios. In practical implementation, the time axis is discretized into several time-slices. For each time-slice, a cross-section perpendicular to the time axis is established. The intersection of this cross-section with the envelope surfaces of each scenario yields several power intervals. The minimum value of the upper bound of the power intervals for all scenarios at that time-slice cross-section and the maximum value of the lower bound of the power intervals for all scenarios are extracted. The closed interval formed by these two values ​​is the power intersection range that all scenarios can reach at that time. This intersection range ensures that, under any considered uncertainty scenario, the load group can provide regulation capability within this interval.

[0120] By traversing all time slices within the future time window and calculating the power intersection interval for each time point, these discrete intervals are arranged and combined in chronological order to form a time-series-based power reachable set. Each element of this set corresponds to the upper and lower limits of the power adjustment range achievable by the flexible load group at a specific time. The time series formed by the upper power bounds of each time point in this set is extracted as the upper boundary curve of the robust adjustable domain, and the time series formed by the lower power bounds of each time point is extracted as the lower boundary curve. The area enclosed by the upper and lower boundary curves represents the range of load regulation capabilities that the virtual power plant can safely invoke when facing various uncertain scenarios. Any dispatch command within this range can achieve an effective response and execution from the load group under all considered scenarios. Through this robust adjustable domain, the dispatch decision module can formulate real-time control strategies while ensuring response reliability.

[0121] A second aspect of the present invention provides a virtual power plant adjustable load aggregation and real-time response control system, comprising:

[0122] The data acquisition unit is used to acquire real-time operating data and historical response data of multiple distributed flexible load devices. Based on the physical characteristic constraints and energy conversion relationships of the load devices reflected in the real-time operating data and the historical response data, a state-space model describing the dynamic behavior of a single load is constructed. The state-space model is decomposed into a fast response sub-model and a slow drift sub-model through a multi-scale time decoupling mechanism.

[0123] The model building unit is used to extract the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combine the topological connection relationship of each load device in the power distribution network to aggregate multiple individual loads into a group equivalent representation by weighted aggregation of the state variables;

[0124] The state aggregation unit is used to calculate the deviation between the output state of the group equivalent representation and the virtual power plant dispatch command, and to propagate the deviation to the state space model parameters of each individual load constituting the group equivalent representation, and to dynamically correct the coupling coefficient between the fast response sub-model and the slow drift sub-model based on the deviation.

[0125] The deviation correction unit is used to re-aggregate the corrected fast response sub-model and the slow drift sub-model into an updated group equivalent representation, perform simulation prediction on the flexible load group based on the updated group equivalent representation, generate scheduling decision support information based on the simulation prediction results, and send it to the corresponding load equipment for execution.

[0126] A third aspect of the present invention provides an electronic device, comprising:

[0127] processor;

[0128] Memory used to store processor-executable instructions;

[0129] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0130] 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.

[0131] 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.

[0132] 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 therein. Such 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 adjustable load aggregation and real-time response control method, characterized in that, include: Acquire real-time operation data and historical response data of multiple distributed flexible load devices. Based on the physical characteristic constraints and energy conversion relationships of the load devices reflected in the real-time operation data and the historical response data, construct a state-space model describing the dynamic behavior of a single load, including: extracting the power change curve, energy storage state change curve and control input sequence of the load devices from the real-time operation data and the historical response data. Based on the power peak value, upper and lower limits of energy capacity and slope of state change reflected in the power change curve and the energy storage state change curve, the physical characteristic constraints of the load equipment are identified, and the value boundaries and rate of change limits of the state variables are determined based on the physical characteristic constraints. Based on the energy conversion relationship, a coupling equation is established between the state variables. The value boundary and the rate of change limit are used as constraints of the coupling equation, and combined with the coupling equation to form a system of differential equations. The system of differential equations is discretized in time, transforming the continuous-time system of differential equations into a discrete-time state transition matrix form. The control input sequence is used as the input of the state transition matrix form, and the power change curve is used as the expected output of the state transition matrix form. By minimizing the error between the predicted output of the state transition matrix form and the power change curve, a state space model is constructed. The state space model is then decomposed into a fast response sub-model and a slow drift sub-model through a multi-scale time decoupling mechanism. Extract the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combine the topological connection relationship of each load device in the power distribution network. By weighted aggregation of the state variables, multiple individual loads are aggregated into a group equivalent representation. The deviation between the output state of the group equivalent representation and the virtual power plant dispatch command is calculated. The deviation is then propagated in reverse to the load state space model parameters of each individual unit constituting the group equivalent representation. The coupling coefficient between the fast response sub-model and the slow drift sub-model is dynamically corrected based on the deviation. The modified fast response sub-model and slow drift sub-model are re-aggregated into an updated group equivalent representation. Based on the updated group equivalent representation, the flexible load group is simulated and predicted. Based on the simulation prediction results, scheduling decision support information is generated and sent to the corresponding load equipment for execution.

2. The method according to claim 1, characterized in that, Based on the energy conversion relationship, a coupling equation is established between the state variables. The boundary values ​​and the rate of change limits are used as constraints on the coupling equation. The coupling equation is combined with the following to form a system of differential equations: Based on the energy flow path inside the load equipment described in the energy conversion relationship, the causal transmission chain between the energy storage state variable, the power input state variable and the power output state variable is identified, and based on the conversion delay characteristics of each state variable in the causal transmission chain, a differential coupling equation describing the time evolution of the state variable is established. In the differential coupling equation, a nonlinear correction term characterizing energy conversion loss is introduced, the value boundary is expressed as the magnitude constraint function of the state variable, the rate of change constraint is expressed as the slope constraint function of the time derivative of the state variable, and the magnitude constraint function and the slope constraint function are embedded as penalty terms into the objective functional corresponding to the differential coupling equation. By using the variational method to find the extremum conditions of the objective functional containing the penalty term, the objective functional is transformed into a system of differential equations in Euler-Lagrange form.

3. The method according to claim 1, characterized in that, Extracting the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combining the topological connection relationship of each load device in the power distribution network, the state variables are weighted and aggregated to form a group equivalent representation, including: Extract the fast state variables of each load device from the fast response sub-model, and extract the slow state variables of each load device from the slow drift sub-model; Based on the topological connection relationship of each load device in the power distribution network, a topological adjacency matrix reflecting the electrical coupling strength between load devices is constructed, and the contribution weight of each load device to the group response characteristics is calculated based on the topological adjacency matrix. The fast state variables are weighted and summed according to the contribution weight to obtain the fast equivalent state of the group, and the slow state variables are weighted and summed according to the contribution weight to obtain the slow equivalent state of the group. The fast equivalent state of the group and the slow equivalent state of the group are combined into a group equivalent state vector, and a group equivalent representation is constructed based on the group equivalent state vector.

4. The method according to claim 1, characterized in that, Calculating the deviation between the output state of the group equivalent representation and the virtual power plant dispatch command, and propagating the deviation backward to the individual load state space model parameters constituting the group equivalent representation includes: Obtain the desired power trajectory and desired response time contained in the virtual power plant dispatch instruction, calculate the power deviation between the output state of the group equivalent representation and the desired power trajectory, and calculate the time delay deviation between the actual response time and the desired response time. Construct an inverse mapping link from the output state of the group equivalent representation to the parameters of the state space model of each individual load. Based on the sensitivity coefficient of each individual load to the output state of the group in the inverse mapping link, decompose the power deviation and the time delay deviation according to the sensitivity coefficient to obtain the parameter correction amount allocated to the state space model of each individual load. The parameter correction is applied to the system matrix and control matrix in the state-space model of the corresponding unit load, and the state-space model parameters are updated.

5. The method according to claim 1, characterized in that, Based on the updated equivalent representation of the group, simulation prediction is performed on the flexible load group. Scheduling decision support information is generated based on the simulation prediction results and sent to the corresponding load equipment for execution, including: Using the sequence of scheduling instructions for multiple scenarios within a future time window as input, and based on the updated group equivalent representation, we perform multi-trajectory parallel simulation and deduction to obtain the set of predicted response trajectories under different scenarios. Extract the group response envelope surface corresponding to each scenario from the predicted response trajectory set, and determine the robust adjustable domain of the flexible load group within the future time window by solving the intersection region of the group response envelope surface; A multi-objective optimization model is constructed within the robust adjustable domain, with the objectives of minimizing power grid dispatching costs and minimizing load equipment response deviations. The Pareto optimal dispatching strategy set is obtained by solving the multi-objective optimization model, and the current optimal strategy is selected from the Pareto optimal dispatching strategy set based on the real-time operating status of the power grid. The current optimal strategy is decomposed into the target power trajectory and execution timing of each individual load, and combined with the physical constraint boundaries of each individual load, scheduling decision support information containing trajectory tracking targets and constraint tolerance ranges is generated and sent to the corresponding load devices for execution through the communication interface.

6. The method according to claim 5, characterized in that, Extracting the population response envelope surface corresponding to each scenario from the predicted response trajectory set, and determining the robust adjustable domain of the flexible load population within the future time window by solving the intersection region of the population response envelope surface, includes: Extract the extreme points of the group response power at each time under each scenario from the predicted response trajectory set, and connect the extreme points of the group response power in the power-time coordinate space to form the response trajectory boundary curves of each scenario. The response trajectory boundary curves of each scenario are subjected to convex hull processing, and the power-time point set of the response trajectory boundary curve envelope is mapped into a convex polyhedron to obtain the group response envelope surface corresponding to each scenario. The group response envelopes of multiple scenarios are geometrically intersected in the power-time coordinate space. The intersection interval of the group response envelopes in the power dimension is solved time by time to obtain the power range that all scenarios can reach at each time. The power ranges at each moment are combined according to the time series to form a power reachable set covering future time windows, and the upper and lower bounds of the power reachable set are used as the boundaries of the robust adjustable domain.

7. A virtual power plant adjustable load aggregation and real-time response control system, used to implement the method of any one of claims 1-6, characterized in that, include: The data acquisition unit is used to acquire real-time operating data and historical response data of multiple distributed flexible load devices. Based on the physical characteristic constraints and energy conversion relationships of the load devices reflected in the real-time operating data and the historical response data, a state-space model describing the dynamic behavior of a single load is constructed, including: extracting the power change curve, energy storage state change curve and control input sequence of the load devices from the real-time operating data and the historical response data. Based on the power peak value, upper and lower limits of energy capacity and slope of state change reflected in the power change curve and the energy storage state change curve, the physical characteristic constraints of the load equipment are identified, and the value boundaries and rate of change limits of the state variables are determined based on the physical characteristic constraints. Based on the energy conversion relationship, a coupling equation is established between the state variables. The value boundary and the rate of change limit are used as constraints of the coupling equation, and combined with the coupling equation to form a system of differential equations. The system of differential equations is discretized in time, transforming the continuous-time system of differential equations into a discrete-time state transition matrix form. The control input sequence is used as the input of the state transition matrix form, and the power change curve is used as the expected output of the state transition matrix form. By minimizing the error between the predicted output of the state transition matrix form and the power change curve, a state space model is constructed. The state space model is then decomposed into a fast response sub-model and a slow drift sub-model through a multi-scale time decoupling mechanism. The model building unit is used to extract the state variables of each load device in the fast response sub-model and the slow drift sub-model, and combine the topological connection relationship of each load device in the power distribution network to aggregate multiple individual loads into a group equivalent representation by weighted aggregation of the state variables; The state aggregation unit is used to calculate the deviation between the output state of the group equivalent representation and the virtual power plant dispatch command, and to propagate the deviation to the state space model parameters of each individual load constituting the group equivalent representation, and to dynamically correct the coupling coefficient between the fast response sub-model and the slow drift sub-model based on the deviation. The deviation correction unit is used to re-aggregate the corrected fast response sub-model and the slow drift sub-model into an updated group equivalent representation, perform simulation prediction on the flexible load group based on the updated group equivalent representation, generate scheduling decision support information based on the simulation prediction results, and send it to the corresponding load equipment for execution.

8. 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 6.

9. 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 6.

Images