Optimization Strategy for Demand Response Participation of Heterogeneous Flexible Load Aggregators Based on Inner Approximation Method

Abstract

Flexible load resources on the demand side are characterized by diverse types, small capacity, and wide distribution. Since these load resources cannot be individually invoked at the system level, it is becoming increasingly important to consider how to establish heterogeneous flexible load aggregates and utilize their adjustable characteristics. In response to the above problems, this paper establishes mechanistic models for three specific heterogeneous load resources, uses the basic homothetic polytope to approximate the original feasible region of loads, and aggregates the feasible region of heterogeneous load resources using the Minkowski sum, while ensuring the aggregation speed and accuracy. Finally, for scenarios of system load shedding and high wind and solar penetration, a demand response strategy on the load side is introduced, and the role of this method in system power supply guarantee and new energy consumption is illustrated through a numerical example.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Accompanied by the continuous increase in the proportion of new energy generation, traditional generation-side resources have difficulty in dealing with the problem that the output of wind and solar power fluctuates strongly due to environmental factors [1, 2], posing challenges to the security and stability of the power system. The electricity consumption of flexible demand-side resources such as air conditioners (ACs), electric vehicles (EVs), and data centers (DCs) continues to rise. For instance, in some regions of China, ACs already account for over 50% of total electricity demand during summer. According to the Energy and AI Report released by the International Energy Agency (IEA), in 2024, DCs consumed approximately 1.5% of global electricity demand. By 2030, their electricity consumption is projected to more than double, reaching around 945 TWh [3]. The IEA estimates that by 2030, DCs will contribute to about one-tenth of the global growth in electricity demand, though their share remains lower than that of industrial motors, ACs, or EVs. The flexibility potential of these demand-side resources is determined by their inherent mechanisms, reflected in dimensions such as adjustable capacity ranges [4, 5], response characteristics [6], and reliability boundaries [7]. Notably, flexible loads like EVs and DCs exhibit bidirectional regulation capabilities. They can increase consumption during off-peak hours to absorb renewable energy and reduce consumption during peak hours to alleviate grid stress, thereby enhancing the power system’s flexibility. Price or incentive signals can guide flexible load resources to change their power consumption behaviors through demand response, thereby compensating for the system’s inability to effectively balance power and energy due to the stochastic fluctuations of new energy sources [8, 9]. This has become a new challenge for the power system. Therefore, how to accurately characterize the adjustable capabilities of heterogeneous load resources [10] and how to effectively guide flexible load resources to participate in demand response [11, 12] have become focal issues in academic research.

Existing studies often focus on a single flexible load, combining their flexibility in electricity consumption for applications such as virtual power plant (VPP) economic dispatch, new energy consumption, and demand response markets. The adjustable power range of aggregated cluster loads is typically obtained by summing the power ranges of similar loads within a region using the boundary Minkowski sum. However, this method effectively enlarges the feasible region. For a more precise description of the feasible region of cluster loads, most existing methods adopt a bottom-up approach. Specifically, individual load mechanistic models are first established, and the feasible region of flexible loads is represented using a vertex expression [13] or a polytope [14]; the Minkowski sum [15] is then performed in Euclidean space. However, as the time scale and the number of aggregated loads increase, the summation process faces nondeterministic polynomial (NP-hard) problems [16], significantly increasing the computational complexity. To address these challenges, researchers have proposed approximate methods for the feasible region, sacrificing part of the feasible region to ensure efficient computation of the aggregated model. For instance, Reference [17] introduces the volume box approximation method, which uses square boxes in Euclidean space to approximate the original feasible region of loads. This method is simple in principle and efficient in computation but suffers from significant flexibility loss. Reference [18] employs an inner approximation Minkowski sum aggregation method based on affine transformations of a convex polytope to quantify the adjustable capabilities of diverse flexible loads. However, this method has a complex solution model and high computational complexity. Reference [19] adopts the Zonotope method, which uses symmetric polyhedrons in space to approximate the feasible region. However, as the time scale increases, the selection of limited directional vectors restricts its extension to only a limited number of directions, resulting in a loss of load flexibility. Reference [20] adopts the basic homothetic polytope method, which like affine transformations of convex polytope, uses basic homothetic polytope to approximate the original feasible region. Compared to fully symmetric polyhedrons, it has better adaptability, theoretically higher approximation accuracy, a simpler solution model, and faster computation. Therefore, this paper selects the basic homothetic polytope to approximate the load feasible region. Reference [21] takes hydrogen hybrid EVs as the research object and studies the low-carbon environmental-economic dispatch issue of their participation in demand response. Reference [22] models the inverter air conditioner (IAC) in VPP and proposes a strategy to coordinate a vast number of IACs with heterogeneous parameters and states to participate in demand response, providing regulation services for VPP. The above references only consider the homogeneous flexible loads participating in demand response. However, VPP contains a large number of heterogeneous flexible loads with huge adjustable capabilities. Aggregation of a single type cannot tap the greater regulation potential on the demand side. Therefore, comprehensively considering multiple heterogeneous loads is of great significance for dispatching operation control in the power system.

Reference [23] models the flexible loads in the industrial park as transferable loads and proposes a coordinated optimization scheduling strategy involving power generation, loads, and energy storage resources to achieve the economy and low-carbon characteristics of the industrial park. Reference [24] takes into account three types of loads in the park, namely, shiftable loads, transferable loads, and curtailable loads, establishes flexible load models, and reduces the cost during the operation of the power system by flexibly scheduling these loads. The abovementioned research failed to take into account the inherent regulatory characteristics of the load and predetermined the adjustable capacity range of all three load categories (shiftable, transferable, and curtailable), indicating that there is a lot of room for improvement in the method. However, they only classify the loads into three categories (shiftable, transferable, and curtailable) and lack credibility by predefining the adjustable capacity ranges of the three types of loads. Reference [25] considers ACs and EVs as aggregates and proposes an effective demand response strategy for residential areas according to time-of-use electricity prices to reduce the electricity burden of households and greenhouse gas emissions into the atmosphere. Reference [26] explores the adjustable capabilities of ACs and EVs in microgrids and reduces the operation costs of microgrids by dispatching the adjustable capabilities of flexible loads through distribution network operators. Most of the existing researches focus on using EVs and ACs as flexible load, lacking consideration of other loads such as DCs. Meanwhile, when modeling, the coupling of the time series power and energy constraints of loads is not specifically considered, which expands the original feasible region of loads and indirectly affects the effect and accuracy of loads aggregation.

To address the above issues, this paper selects EVs, ACs, and DCs as the research objects of heterogeneous flexible loads. This study develops a universal mechanistic model for flexible loads by holistically integrating the operational principles of the demand-side resources and their inherent power–energy coupling characteristics across temporal consumption profiles. The model leverages half-plane expressions to mathematically delineate the power adjustable domains of individual flexible loads, employs the basic homothetic polytope approximation to derive computationally tractable solutions for these domains, and utilizes Minkowski sum operations to aggregate heterogeneous load clusters into a unified feasible region model—thereby systematically quantifying the flexibility capabilities of load aggregators. In parallel, a hierarchical demand response strategy is formulated to incentivize participation of heterogeneous load clusters in demand response, while rigorously accounting for their operational constraints and energy requirements. Finally, through case analysis, it verifies that the proposed heterogeneous load aggregate model and demand response strategy enhances the economic efficiency and cleanliness in the VPP.

2. Characterization of Adjustable Characteristics of Heterogeneous Load Resources

2.1. EV Model

In recent years, EVs have gradually entered households, serving as an important flexibility resource. Characterized by wide distribution and small individual capacity, EVs are modeled in this paper using the concept of virtual batteries [27]. The operational characteristics of EVs, including charging and discharging states, are constructed as follows:$$\begin{array}{}\text{(1)}& P_{\max }^{\mathrm{d}}\le P_t^{\text{EV}}\le P_{\max }^{\mathrm{c}},\text{(2)}& E_{\min }^{\text{EV}}\le E_t^{\text{EV}}\le E_{\max }^{\text{EV}},\text{(3)}& E_{t+1}^{\text{EV}}=1-\nu E_t^{\text{EV}}+{\eta }^{\text{EV}}{P}_t^{\text{EV}}\mathrm{\Delta }t,\end{array}$$where $P_t^{\text{EV}}$ represents the charging and discharging power at time t; $P_{\max }^{\mathrm{d}}$ denotes the maximum discharging power of EVs; $P_{\max }^{\mathrm{c}}$ denotes the maximum charging power of EVs; $E_t^{\text{EV}}$ represents the energy state of EVs at time t; $E_{\min }^{\text{EV}}$ and $E_{\max }^{\text{EV}}$ represent the boundary conditions of EV’s energy state; $\nu$ represents the charging and discharging loss factor; and ${\eta }^{\text{EV}}$ represent the charging and discharging efficiency.

2.2. AC Model

According to the relevant report of the National Bureau of Statistics, Chinese AC electricity consumption ranks first in the world, making it an important type of adjustable load. This paper utilizes the thermodynamic model of AC to depict its operating mechanism [28]:$$\begin{array}{}\text{(4)}& \frac{\mathrm{d}{T}_t^{\text{in}}}{\mathrm{d}t}=\frac{T_t^{\text{out}}-T_t^{\text{in}}}{R_{\text{in}}{C}_{\text{in}}}-\frac{D_t^{\text{AC}}}{C_{\text{in}}},\text{(5)}& D_t^{\text{AC}}={\eta }^{\text{AC}}{P}_t^{\text{AC}},\end{array}$$where $T_t^{\text{in}}$ and $T_t^{\text{out}}$ represent the indoor temperature and outdoor temperature at time t, respectively; $C_{\text{in}}$ represents the air thermal capacity; ${\eta }^{\text{AC}}$ represents the working efficiency of the AC; $R_{\text{in}}$ represents the thermal resistance of the room; $P_t^{\text{AC}}$ represents the AC’s power at time t; and $D_t^{\text{AC}}$ represents the cooling capacity of the AC.

Solving the first-order linear nonhomogeneous differential equation of equation (4) yields the expression between the indoor environmental temperature and the AC load at adjacent time steps, as shown in equation (6). Meanwhile, the constraints on indoor temperature and AC power are provided:$$\begin{array}{}\text{(6)}& T_{t+1}^{\text{in}}=T_t^{\text{in}}{e}^{-\mathrm{\Delta }t/R_{\text{in}}{C}_{\text{in}}}+{\eta }^{\text{AC}}{R}_{\text{in}}{P}_t^{\text{AC}}+T_t^{\text{out}}1-e^{-\mathrm{\Delta }t/R_{\text{in}}{C}_{\text{in}}},\text{(7)}& T_{\min }^{\text{in}}\le T_t^{\text{in}}\le T_{\max }^{\text{in}},\text{(8)}& 0\le P_t^{\text{AC}}\le P_{\max }^{\text{AC}},\end{array}$$where $T_{\min }^{\text{in}}$ and $T_{\max }^{\text{in}}$ represent the upper and lower boundaries of indoor temperature that the human body can tolerate; $P_{\max }^{\text{AC}}$ represents the maximum power of AC; and $\mathrm{\Delta }t$ represents the time step.

2.3. DC Model

New-generation DCs [29] have a wide range of application scenarios in society, such as real-time navigation, large AI models, data storage, and downloading. According to relevant data, the annual power consumption of DCs accounts for 2% of China’s total electricity consumption. With the construction of a smart nation, the demand for computing power is increasing. Relevant institutions predict that by 2030, the annual electricity consumption of national intelligent computing centers alone will account for 5%–10% of the total electricity consumption of society.

The majority of DC energy consumption comes from IT equipment. IT equipment processes load data dispatched by front-end servers, which is typically divided into two categories: delay-tolerant loads and latency-sensitive loads [30]. Delay-tolerant loads refer to tasks that can be processed by IT equipment within a certain time frame, such as cloud storage downloads and data training, which do not require real-time processing. Latency-sensitive loads consist of real-time tasks such as web browsing and real-time navigation, which have high latency sensitivity. Existing research demonstrates that the energy consumption of DCs is closely correlated with parameters such as workload volume, data scale, and processing duration [31]. As workload increases, server power consumption rises accordingly, accompanied by greater heat dissipation, which in turn drives up cooling system power usage, ultimately leading to a significant increase in the DCs’ total power consumption. To simplify the modeling process, this study employs the power usage effectiveness metric [32] to quantify the relationship between server power consumption and the data DCs’ total electricity usage. Furthermore, based on the linear energy consumption modeling approach, we establish an association model between workload processing volume and server power consumption, providing a more intuitive representation of how DCs’ electricity usage dynamically varies with operational workload:$$\begin{array}{}\text{(9)}& L_{t+1}^{\text{DC}}=L_t^{\text{DC}}+m_t^{\mathrm{\alpha }}+m_t^{\mathrm{\beta }}\mathrm{\Delta }t,\text{(10)}& P_t^{\text{IT}}=\frac{N\mathrm{·}{P}^{\text{idl}}+m_t^{\mathrm{\alpha }}+m_t^{\mathrm{\beta }}{P}^{\text{peak}}-P^{\text{idl}}}{{\lambda }^{\text{DC}}},\text{(11)}& P_t^{\text{DC}}=1+{\eta }^{\text{DC}}{P}_t^{\text{IT}},\text{(12)}& P_t^{\text{DC},\min }\le P_t^{\text{DC}}\le P_t^{\text{DC},\max },\text{(13)}& L_t^{\text{DC},\min }\le L_t^{\text{DC}}\le L_t^{\text{DC},\max },\end{array}$$where $m_t^{\mathrm{\alpha }}$ and $m_t^{\mathrm{\beta }}$ represent the number of delay-tolerant loads and latency-sensitive loads allocated to the server at time t, respectively; $L_t^{\text{DC}}$ represents the total number of loads allocated to the server at time t; $P^{\text{idl}}$ and $P^{\text{peak}}$ represent the energy consumption of the server in idle and running states, respectively; $N$ represents the number of servers owned by the DCs; ${\lambda }^{\text{DC}}$ represents load processing efficiency, i.e., the amount of load each server can process per unit time; $P_t^{\text{IT}}$ and $P_t^{\text{DC}}$ represent the power consumption of IT equipment and DCs, respectively; ${\eta }^{\text{DC}}$ represents the energy efficiency index of the DCs, which is the ratio of the electricity consumed by the DCs to that consumed by IT equipment; $P_t^{\text{DC},\min }$ and $P_t^{\text{DC},\max }$ represent the boundaries for power consumption of the DCs at time t, respectively; and $L_t^{\text{DC},\min }$ and $L_t^{\text{DC},\max }$ represent the boundaries of the number of loads that the DCs can handle at time t.

3. Aggregation Method for Feasible Region of Multielement Heterogeneous Loads

3.1. Approach for Building the Feasible Region of Load Resources

The existing approaches for depicting the feasible region of a single load consist of vertex expression and half-plane expression, as detailed below:$$\begin{array}{}\text{(14)}& \mathbf{P}={\displaystyle \sum }_{i=1}^I{\alpha }^i{\mathit{\upsilon }}^i{\mathit{\upsilon }}^i\in {\mathfrak{R}}^n,{\alpha }^i\mathfrak{\in }\mathfrak{R},{\alpha }^i\ge 0,{\displaystyle \sum }_i^I{\alpha }^i=1,\text{(15)}& \mathbf{P}={\mathit{p}}_i\in {\mathfrak{R}}^{T\times 1}{\mathit{A}}_i{\mathit{p}}_i\le {\mathit{\beta }}_i,\end{array}$$where $\mathbf{P}$ represents the adjustable feasible region of load resources; ${\mathit{\upsilon }}^i$ represents the $i-\text{th}$ vertex; ${\alpha }^i$ represents the nonnegative coefficients of the convex hull; ${\mathit{p}}_i$ represents the power vector of the $k-\text{th}$ load; and ${\mathit{A}}_i$ and ${\mathit{\beta }}_i$ represent the coefficient matrices of the $k-\text{th}$ load.

As the time scale increases, when using the vertex expression method, the computational difficulty increases exponentially, and the solution accuracy also decreases with the increase of numerical calculation errors. Compared with the vertex expression, the plane expression has good scalability. When the time scale increases, only the constraints at the corresponding moments need to be added to the original expression. Therefore, this paper chooses the half-plane expression to represent the feasible region of loads.

This paper first unifies the mechanism models of various load resources given in Section 1 into a general mechanistic model with the same structure as equation (16). Then, the half-plane expression is used to represent the feasible region of the power of various loads. When constructing a specific general mechanism model for flexible loads, for DCs loads, since the delay-tolerant loads can be processed with a delay within a certain time period, the maximum amount of loads to be processed at each moment is related to the original load curve, and the load processing amount between adjacent moments directly affects the power boundary of the DCs. Considering this point, unlike other papers that give fixed power and “energy” boundaries, considering the coupling of power and energy constraints in the time sequence of loads, time-varying boundary conditions are given, which improves the accuracy of representing the feasible region of the original loads.$$\begin{array}{}\text{(16)}& \begin{array}{c}{P}_t^{\min }\le P_t\le P_t^{\max },\\ E_t^{\min }\le E_t\le E_t^{\max },\\ E_{t+1}=\alpha E_t+\eta P_t\mathrm{\Delta }t+{\beta }_t,\end{array}\end{array}$$where $P_t$ represents the power consumption of the flexible load at time t; $P_t^{\min }$ and $P_t^{\max }$ represent the lower and upper bounds of the flexible load’s power consumption at time t, respectively; $E_t$ represents the energy state of the flexible load at time t; $E_t^{\min }$ and $E_t^{\max }$ represent the minimum and maximum energy limits at time t; and $\alpha$, $\eta$, ${\beta }_t$ are dynamic characteristic parameters governing the load’s temporal behaviors.

3.2. Load Feasible Region Aggregation Method Based on Basic Homothetic Polytope Method

The Minkowski sum of a polyhedron formed by multiple point sets is defined as in equation (17). In simple terms, the Minkowski sum process involves sequentially arranging and correspondingly adding the point sets of each set. The resulting new vertex set represents the power feasible region of the aggregated load resources.$$\begin{array}{}\text{(17)}& \mathbf{S}n+\cdots +m=c_i+\cdots +c_j\forall i\in 1,n,\dots \forall j\in 1,m,\end{array}$$where $c_i$ and $c_j$ represent point sets and $n$ and $m$ represent the number of points in each point set, respectively.

In the context of massive load scenarios, load resources exhibit diverse operating characteristics and significant parameter differences. If the polyhedrons representing individual load resources are aggregated using the vertex method for Minkowski sum, the computational space complexity becomes high, the solution time is long, and the problem becomes an NP-hard problem. Therefore, this paper adopts the idea of transforming the original power feasible region of load resources using a set of representation forms with the same structure but different parameters to achieve an inner approximation of the feasible region for massive heterogeneous load resources, thereby realizing the aggregation function.

By sequentially performing benchmark feasible region determination, scaling and displacement transformations of the benchmark feasible region for inner approximation of the original feasible region, and Minkowski sum of the approximated feasible region, the aggregated power feasible region of heterogeneous load resources is obtained.

This paper adopts the averaging method to set the benchmark feasible region, as follows:$$\begin{array}{c}\text{(18)}& {\mathit{a}}_0=\frac{{\sum }_{i=1}^n{\mathit{a}}_i}{n},{\mathit{b}}_0=\frac{{\sum }_{i=1}^n{\mathit{b}}_i}{n},\end{array}$$where ${\mathit{a}}_0$ and ${\mathit{b}}_0$ represent the coefficient matrices of the selected benchmark feasible region and $n$ represents the number of loads to be aggregated.

The expression for the power feasible region after scaling and displacement transformation of the benchmark feasible region is as follows:$$\begin{array}{}\text{(19)}& {\tilde{\mathbf{\Omega }}}_i={\mathit{p}}_i\in {\mathfrak{R}}^{T\times 1}{\mathit{a}}_0{\mathit{p}}_i-{\mathit{\psi }}_i\le {\phi }_i{\mathit{b}}_0,\end{array}$$where ${\tilde{\mathbf{\Omega }}}_i$ represents the $i-\text{th}$ load power feasible region after inner approximation and ${\mathit{\psi }}_{\mathit{i}}$ and ${\phi }_i$ represent the translation coefficient and scaling coefficient matrix of the $i-\text{th}$ benchmark feasible region, respectively.

The approximate solution process within the feasible region can be transformed into the problem of solving for the maximum translation coefficient and scaling coefficient within the original feasible region range. This paper adopts the method proposed in Reference [20]. The model solution is as follows:$$\begin{array}{}\text{(20)}& \underset{{\phi }_{i,}{\mathit{\psi }}_i}{\max }{\phi }_i,\text{(21)}& {\tilde{\mathbf{\Omega }}}_i={\phi }_i{\tilde{\mathbf{\Omega }}}_0+{\mathit{\psi }}_i\in {\mathbf{\Omega }}_i,\text{(22)}& {\phi }_i\ge 0.\end{array}$$

Transform the above solution model into an LP problem for solution by introducing auxiliary variables, ${\lambda }_i=1/{\phi }_i$, ${\mathit{U}}_i$, and ${\mathit{\mu }}_i=-{\lambda }_i{\mathit{\psi }}_i$:$$\begin{array}{}\text{(23)}& \underset{{\lambda }_i,{\mathit{\mu }}_i,{\mathit{U}}_i}{\min }{\lambda }_i,\text{(24)}& {\mathit{U}}_i{\mathit{a}}_0={\mathit{a}}_i,\text{(25)}& {\mathit{U}}_i,{\lambda }_i\ge 0,\text{(26)}& {\mathit{U}}_i{\mathit{b}}_0\le {\mathit{b}}_i{\lambda }_i+{\mathit{a}}_i{\mathit{\mu }}_i.\end{array}$$

Since the feasible region of the aggregated load resources after inner approximation processing has the same structure but different parameters, the aggregation of the feasible region of massive heterogeneous loads can be directly achieved by equation (27) through the Minkowski sum of the scaling and translation coefficient matrices of the loads to be aggregated. Thus, the overall feasible region of the aggregate can be expressed as equation (28):$$\begin{array}{}\text{(27)}& \begin{array}{c}{\phi }^{\text{AGG}}={\displaystyle \sum }_{i=1}^m{\phi }_i,\\ {\mathit{\psi }}^{\text{AGG}}={\displaystyle \sum }_{i=1}^m{\mathit{\psi }}_i,\end{array}\text{(28)}& {\tilde{\mathbf{\Omega }}}^{\text{AGG}}={\mathit{P}}^{\text{AGG}}\in {\mathfrak{R}}^{T\times 1}{\mathit{a}}_0{\mathit{P}}^{\text{AGG}}-{\mathit{\psi }}^{\text{AGG}}\le {\phi }^{\text{AGG}}{\mathit{b}}_0,\end{array}$$where ${\phi }^{\text{AGG}}$ and ${\mathit{\psi }}^{\text{AGG}}$, respectively, represent the scaling coefficient and translation coefficient of the aggregate; ${\tilde{\mathbf{\Omega }}}^{\text{AGG}}$ represents the feasible region of inner-approximated aggregate; and ${\mathit{P}}^{\text{AGG}}$ represents the power vector of the aggregate.

3.3. Multivariate Heterogeneous Load Aggregation Method

Although homogeneous load resources have similar operational modes, differences in their initial operating parameters lead to variations in the feasible region among different individuals. Using a single benchmark feasible region to approximate the original feasible regions of each load within a massive cluster of homogeneous loads will result in significant errors. Therefore, this paper adopts a progressive approach, namely, the concept of “clustering before aggregation,” to solve for the feasible region of homogeneous load resources. Load resource initial operating parameters are selected as load characteristics, and homogeneous load resources are grouped through the K-means algorithm. For each cluster of homogeneous loads with similar characteristics, the original feasible region approximation algorithm and Minkowski summation method proposed in Section 2.2 are applied for aggregation, resulting in the feasible region of each cluster aggregate as shown in equation (28). The aforementioned method is then used to aggregate the feasible regions of each cluster aggregate into the feasible region of the overall homogeneous load aggregate resource. This method can largely avoid the issue where selecting a single benchmark feasible region cannot effectively represent the feasible regions of a vast number of regional load resources due to differences in load resource power and “capacity.”

To facilitate readers’ understanding, Algorithm 1 presents the pseudocode of the heterogeneous load aggregation method. Additionally, Algorithm 2 provides the pseudocode for the heterogeneous load clustering method.

Algorithm 1: Heterogeneous load aggregate algorithm (pseudocode).

  Input: $\mathit{X}={\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_i,\dots ,{\mathit{x}}_m$-operational parameter of loads; M-number of loads; K-the number of clusters for loads by K-means algorithm; $N^i$-number of loads in cluster i.

  Output: ${\tilde{\mathbf{\Omega }}}^{\text{AGG}}$-aggregate feasible region.

  for m ⟵ 1 to M do

   Construct the feasible region based on mechanistic constraints with the same structure as equation (16);

   Obtain the coefficient matrices of the feasible region expression $\mathit{A}={\mathit{a}}_1,{\mathit{a}}_2,\dots ,{\mathit{a}}_i,\dots ,{\mathit{a}}_m$ and $\mathit{B}={\mathit{b}}_1,{\mathit{b}}_2,\dots ,{\mathit{b}}_i,\dots ,{\mathit{b}}_m$;

  end for

  Assign the coefficient matrices of different loads into distinct clusters generated by the K-means algorithm.

  for k ⟵ 1 to K do

   Calculate the coefficient matrices of the selected benchmark feasible region as equation (18).

   for n ⟵ 1 to $N^i$do

    Obtain the values of the optimization variables ${\lambda }_n$, ${\mu }_n$, and ${\mathit{U}}_n$ by solving the optimization problem with equation (23) as the objective function, subject to constraints equations (24)–(26).

    Obtain the translation coefficient and scaling coefficient matrix ${\phi }_n$ ⟵ $1/{\lambda }_n$, ${\mathit{\psi }}_n$ ⟵ $-{\mathit{\mu }}_n/{\lambda }_n$;

   end for

   Calculate the translation coefficient and scaling coefficient of difference aggregates as equation (27);

   Obtain the coefficient matrices that represent the feasible region of difference aggregates based on equation (28);

   Obtain the values of the optimization variables ${\lambda }_k$, ${\mu }_k$, and ${\mathit{U}}_k$ by solving the optimization problem with equation (23) as the objective function, subject to constraints equations (24)–(26).

   Obtain the translation coefficient and scaling coefficient matrix ${\phi }_k$ ⟵ $1/{\lambda }_k$, ${\mathit{\psi }}_k$ ⟵ $-{\mathit{\mu }}_k/{\lambda }_k$;

  end for

  Calculate the translation coefficient and scaling coefficient of the aggregate as equation (27);

  Derive the feasible region of the inner-approximated aggregate formulation ${\tilde{\mathbf{\Omega }}}^{\text{AGG}}$ by equation (28);

  Return${\tilde{\mathbf{\Omega }}}^{\text{AGG}}$

Algorithm 2: K-means algorithm (pseudocode).

  Input: $\mathit{X}={\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_i,\dots ,{\mathit{x}}_m$-operational parameter of loads; M-number of loads; K-the number of clusters for loads; $T$-max iterations; $\epsilon$-early stopping criterion; $D,$-distance metric.

  Output: $\mathit{C}$-cluster assignments of loads;

   $\mathit{C}⟵\phi$, $t⟵0$, $T⟵100$, $K⟵0$, $\epsilon ⟵{10}^{-5}$

  Determine the optimal number of clusters K for the K-means algorithm using the elbow method on dataset.

  While$2\le N\le \sqrt{M}$do

    Randomly select n samples from the dataset $\mathit{X}$ as the initial mean vectors ${\mathit{\mu }}_1,{\mathit{\mu }}_2,\dots {\mathit{\mu }}_n$.

  end while

  While t < T do

    for i ⟵ 1 to M do

     for j ⟵ 1 to K do

      Calculate the distance between the ith data point and the jth cluster center $d_{ij}$

     If $d_{ij}=\mathrm{Min}D{\mathit{x}}_i,{\mathit{\mu }}_i$ then

      ${\mathit{x}}_i\le {\mathit{C}}_j$

     end for

    end for

  for k ⟵ 1 to K do

   Calculate the jth cluster center $C_j$

   Calculate the new mean vectors: ${\mathit{\mu }}_i^{\prime }=1/{\mathit{C}}_i{\sum }_{\mathit{x}\in {\mathit{C}}_i}\mathit{x}$

   if ${\mathit{\mu }}_i^{\prime }-{\mathit{u}}_i\le \epsilon$ then

    ${\mathit{\mu }}_i^{\prime }={\mathit{u}}_i$

   end if

   $t+=1$

   end for

  end while

  Return$\mathit{C}$

4. Day-Ahead Dispatching Model Accounting the Contribution of Multivariate Heterogeneous Loads in Demand Response

4.1. Day-Ahead Dispatching Model

In order to minimize the total operation cost of VPP while meeting the power and energy balance, the dispatching model optimizes the output of various generating units and explores the optimal economic dispatch plan. The objective function includes the operating cost of thermal generation, penalty for curtailment of new energy, and the payment for demand response, as shown in equation (29):$$\begin{array}{}\text{(29)}& \min {\displaystyle \sum }_t^T{\displaystyle \sum }_{g\in {\mathrm{\Omega }}_{\mathrm{G}}}{F}_g{P}_{g,t}+S{U}_{g,t}+S{D}_{g,t}+{\displaystyle \sum }_{e\in {\mathrm{\Omega }}_{\mathrm{E}}}{\sigma }_e{P}_{e,t}+{\displaystyle \sum }_{r\in {\mathrm{\Omega }}_{\mathrm{R}}}{\varsigma }_r{P}_{r,t}^{\text{curt}}+{\displaystyle \sum }_{l\in {\mathrm{\Omega }}_{\mathrm{L}}}{\tau }_l{P}_{l,t}^{\text{curt}}+{\displaystyle \sum }_{d\in {\mathrm{\Omega }}_{\mathrm{D}}}{F}_d{P}_{d,t}\mathrm{\Delta }t,\end{array}$$where $F_g{P}_{g,t}$ represents the cost of coal consumption for thermal power units; $\sigma$, $\varsigma$, and $\tau$, respectively, represent the ESS cost coefficient, curtailment penalty coefficient for new energy, and curtailment penalty coefficient for loads; $F_d{P}_{d,t}$ represents the demand response cost; $S{U}_{u,t}$ and $S{D}_{u,t}$, respectively, represent the start-up cost and shutdown cost of thermal power units; ${\mathrm{\Omega }}_{\mathrm{G}}$ represents set of thermal power units; ${\mathrm{\Omega }}_{\mathrm{R}}$ represents set of new energy generating units; ${\mathrm{\Omega }}_{\mathrm{L}}$ represents set of load nodes; ${\mathrm{\Omega }}_{\mathrm{D}}$ represents set of various types of flexible loads; ${\mathrm{\Omega }}_{\mathrm{E}}$ represents set of energy storage units; $P_{r,t}^{\text{curt}}$ represents the curtailed power of the $r-th$ new energy generating unit; $P_{l,t}^{\text{curt}}$ represents the load curtailment of the $l-th$ node; and $P_{d,t}$ represents the power of the $d-th$ flexible load.

The constraint conditions of the day-ahead scheduling model include the operation constraints of thermal power generating units, the operation constraints of WT and PV generation, the operation constraints of ESS, and the power balance constraints, etc.

The traditional day-ahead scheduling model, when taking into account the participation of load-side resources in demand response, often simply models the load-side resources as shiftable loads, transferable loads, and curtailable loads to simplify computational complexity, and the solution speed of the scheduling model can be improved. However, this mode does not fully consider the adjustment characteristics of various flexible loads, which affects the accuracy and reliability of the dispatch results. This paper models the load-side resources as the multitype heterogeneous load aggregation model described in Section 3, taking full account of the adjustment characteristics of various flexible loads. The flexible load feasible region calculated by equation (28) is used to determine the adjustable boundary of the load aggregator. Meanwhile, incentive-based demand response policies are introduced to enhance the economy and flexibility of the VPP system while reducing carbon emissions.

4.2. Load-Side Regulation Strategy

This paper dynamically adjusts the demand response incentive electricity price based on the regulation amount of flexible loads and the actual load regulation amount to guide loads to participate in demand response and increase the participation degree of load resources. The specific strategy function is as follows:$$\begin{array}{}\text{(30)}& {\displaystyle \sum }_t{P}_{d,t}^0={\displaystyle \sum }_t{P}_{d,t},d\in {\mathrm{\Omega }}_{\mathrm{D}},\text{(31)}& \mathrm{\Delta }{P}_{d,t}=P_{d,t}-P_{d,t}^0,d\in {\mathrm{\Omega }}_{\mathrm{D}},\text{(32)}& P_{d,t}^{\min }\le P_{d,t}\le P_{d,t}^{\max },d\in {\mathrm{\Omega }}_{\mathrm{D}},\end{array}$$where $P_{d,t}^0$ represents the original load curve of the $d-\text{th}$ type of load; $\mathrm{\Delta }{P}_{d,t}$ represents the $d-\text{th}$ type of load regulation amount at time t; $P_{d,t}$ represents the load curve of the $d-\text{th}$ type of load after demand response; and $P_{d,t}^{\min }$ and $P_{d,t}^{\max }$ represent the boundary of the flexible load feasible region obtained through equation (28).

This paper adopts a user incentive electricity price formulation strategy of “benchmark value + floating up and down”:$$\begin{array}{}\text{(33)}& {\lambda }_{d,t}={\lambda }^0+k\mathrm{·}\mathrm{\Delta }{P}_{d,t},d\in {\mathrm{\Omega }}_{\mathrm{D}},\text{(34)}& {\lambda }^{\text{DR},\min }\le {\lambda }_{d,t}\le {\lambda }^{\text{DR},\max },d\in {\mathrm{\Omega }}_{\mathrm{D}},\end{array}$$where ${\lambda }_{d,t}$ represents the incentive price of the $d-\text{th}$ type of load at time t; ${\lambda }^0$ represents the benchmark value of the incentive electricity price; and $k$ represents the price coefficient.

The optimization objectives of the day-ahead dispatch model are modified in the way described below:$$\begin{array}{}\text{(35)}& \min {\displaystyle \sum }_t^T{\displaystyle \sum }_{g\in {\mathrm{\Omega }}_{\mathrm{G}}}{F}_g{P}_{g,t}+S{U}_{g,t}+S{D}_{g,t}+{\displaystyle \sum }_{e\in {\mathrm{\Omega }}_{\mathrm{E}}}{\sigma }_e{P}_{e,t}+{\displaystyle \sum }_{r\in {\mathrm{\Omega }}_{\mathrm{R}}}{\varsigma }_r{P}_{r,t}^{\text{curt}}+{\displaystyle \sum }_{l\in {\mathrm{\Omega }}_{\mathrm{L}}}{\tau }_l{P}_{l,t}^{\text{curt}}+{\displaystyle \sum }_{d\in {\mathrm{\Omega }}_{\mathrm{D}}}{\lambda }_{d,t}\mathrm{\Delta }{P}_{d,t}\mathrm{\Delta }t.\end{array}$$

4.3. Algorithm Flowchart

The day-ahead scheduling process considering the participation of diverse heterogeneous loads in demand response is shown in Figure 1. The process consists of three parts. The first part is the stage of characterizing the feasible region of similar loads and cluster processing. The feasible regions of similar loads, which are widely distributed and have large operational differences, are transformed into half-plane expression through equations (18) and (19).

[figure(s) omitted; refer to PDF]

The second part is the cluster aggregation processing stage of similar loads. Firstly, the similar loads in the region are grouped through the K-means algorithm. The approximate feasible region of the loads in each cluster after K-means processing is calculated through equations (23)–(28). Then, the aggregates of the load clusters are further aggregated to obtain the feasible region of the aggregated resources of massive similar loads.

The third part is the day-ahead scheduling stage. The VPP side collects the feasible region of the aggregated resources of flexible loads in the region, combines the next-day load forecast curve, aims to minimize its own cost as the optimization objective, formulates generation and incentive strategies, and stimulates the electricity flexibility of multivariate heterogeneous loads, enabling not only the transmission of electrical energy but also the conveyance of incentive prices and information about the adjustable capacity of loads between the VPP side and the load side.

5. Case Study

5.1. Parameter Settings of Equipment in the Case

This paper conducts a case analysis using the improved PCUC-6 bus system [33]. The PCUC-6 bus system utilizes the actual renewable energy output scenarios and load data of a certain location in China. In order to adapt to the actual installed capacity of renewable energy and load conditions, the thermal power units plugged in Buses 1, 2, and 6 have been improved. A 200-MW PV generation and a 250-MW WT generation are installed at Bus 1, and energy storage is only plugged in Bus 1. The load buses include 1000 EVs, 1000 ACs, and 10 small-scale DCs. The original load curves of EVs and ACs are referenced from dataset [34], and the original load demand of DCs is simulated based on reference [35]. The detailed data are provided in dataset [36], and it provides key data such as the day-ahead prediction curves in two scenarios of system supply guarantee and high generation of renewable energy, the arrival of DC loads, equipment parameter settings such as flexible loads, thermal power units, energy storage units, the original load curves, and the topology structure of the PCUC-6 six-node system. Among them, the data load in DC are divided into two categories: delay-tolerant and delay-sensitive: Normal distribution sampling was adopted for the delay-tolerant load with time series migration capability; for delay-sensitive loads, a day is divided into four periods according to the human daily routine, and Monte Carlo random sampling is conducted respectively.

5.2. Analysis on the Aggregation Effect of Heterogeneous Loads

This paper sets the volume box approximation method and the Zonotope approximation method as comparative examples to confirm the advantage of the adopted method with respect to approximation accuracy. Since the power feasible region has the characteristics of high dimensionality and being challenging to observe within the time dimension, only the approximation effects of the power feasible region in adjacent time periods are presented.

As can be seen from Figure 2, the horizontal and vertical coordinates, respectively, represent the load power consumption in adjacent time periods. According to the flexible load mechanism model, the power outputs between adjacent time periods are mutually coupled. The light blue area represents the original feasible region of the corresponding flexible load, the region enclosed by the red solid line denotes the approximate feasible region obtained using our proposed method, and the areas bounded by black dashed lines and blue solid lines show the load feasible regions approximated using the volume box approximation method and the Zonotope approximation method, respectively.

[figure(s) omitted; refer to PDF]

The volume box approximation method assumes that the approximate feasible region is an inscribed cube. This method does not take into account the actual shape of the feasible region and cannot represent the energy constraint of the load as shown in equation (16), resulting in a significant loss of the feasible region. However, due to the original power feasible region not being symmetric, compared with the basic homothetic polytope method adopted in this paper, there is still a relatively large loss of feasible region.

The area ratios of three types of loads under three inner approximation methods were calculated to compare their approximation performance. As demonstrated in Table 1, the proposed method in this paper achieves significantly better approximation results compared to both the volume box approximation method and the Zonotope approximation method.

Table 1

Comparative evaluation of flexible load adjustable domains under different strategies.

Strategy	EV	DC	AC
Basic homothetic polytope	0.9447	0.8698	0.9363
Volume box approximation	0.6412	0.8319	0.8867
Zonotope	0.8224	0.5840	0.7254

This paper considers the aggregation of massive heterogeneous loads in the VPP, and it is necessary to verify the accuracy of the aggregation of feasible region of various loads with indicator.

This paper develops an accuracy evaluation framework that systematically compares inner-approximated feasible regions with their original counterparts by computing directional width ratios along multiple sampled unit normal vectors spanning the feasible regions, where each ratio quantifies the approximation fidelity in a specific direction and the collective set of ratios serves as a comprehensive metric for assessing overall approximation quality. The calculation methods are as follows (formulas (36)–(38) referenced from [20]):$$\begin{array}{c}\text{(36)}& {\mathrm{\Lambda }}_{\text{load}}=\frac{1}{S}{\displaystyle \sum }_{s=1}^S\frac{d_s^{\text{apr}}}{d_s^{\mathrm{p}}},\text{(37)}& d_s^{\mathrm{p}}=\frac{\underset{\mathit{p}}{\max }{\mathit{\omega }}_s\mathit{p}-\epsilon -\underset{\mathit{p}}{\min }{\mathit{\omega }}_s\mathit{p}-\epsilon }{{{\mathit{\omega }}_s}_2},\mathrm{s}.\mathrm{t}.\mathit{a}\mathit{p}\le \mathit{b}.\\ \text{(38)}& d_s^{\text{apr},\text{zonotope}}=2{\mathit{\omega }}_s\mathit{G}{\mathit{\beta }}_{\max },\text{(39)}& d_s^{\text{apr},\text{box}}=\frac{{\mathit{\omega }}_s\mathit{u}-\mathit{l}}{{{\mathit{\omega }}_s}_2},\end{array}$$where ${\mathrm{\Lambda }}_{\text{load}}$ denotes the accuracy metric; $S$ represents the number of unit normal vectors sampled; ${\omega }_s$ represents the normal vector; $d_s^{\text{apr}}$ is the width of the approximated adjustable domain along the sth normal vector direction; $d_s^{\mathrm{p}}$ is the width of the original adjustable domain along the same direction; $d_s^{\text{apr},\text{zonotope}}$ and $d_s^{\text{apr},\text{box}}$ denote the adjustable domain widths derived from the Zonotope method and box approximation method, respectively; $\mathit{p}$ is the high-dimensional power vector; $\epsilon$ signifies an unbounded value; $\mathit{G}$, ${\mathit{\beta }}_{\max }$ represent the generator matrix and the scaling length along the ith generator direction; and $\mathit{u}$ and $\mathit{l}$ are the lower and upper bounds of the box domain.

This study compares the approximation accuracy metrics of flexible loads across three temporal dimensions under multiple spatial scales, as detailed in Table 2. The results demonstrate that our proposed method maintains significant precision advantages over both the volume box approximation and Zonotope approximation methods at all time scales. As temporal dimensionality increases, approximation complexity grows accordingly. While conventional methods exhibit rapid accuracy degradation, our approach maintains high fidelity even in 24-h high-dimensional spaces, particularly for EV and DC load feasible regions. Unlike symmetrical polytope-based approximations, our method utilizes asymmetric polyhedral approximations to better conform to the original nonsymmetric feasible regions in high-dimensional spaces. This geometric adaptability enables tighter approximation boundaries and superior precision.

Table 2

Comparison of load accuracy metrics under different time horizons.

Strategy	Load type	Average accuracy metrics
		6 h	12 h	24 h
Basic homothetic polytope	AC	0.9677	0.9595	0.8853
	EV	0.9620	0.9493	0.9580
	DC	0.9335	0.9483	0.9571
Volume box approximation	AC	0.8855	0.8361	0.8223
	EV	0.4545	0.4212	0.3857
	DC	0.5565	0.5218	0.4964
Zonotope	AC	0.7855	0.7117	0.6381
	EV	0.8455	0.7552	0.6615
	DC	0.5488	0.4819	0.4053

This paper presents a computational time comparison for approximating the feasible regions of 100 EVs. As observed in Figure 3, the volume box approximation method achieves the shortest computation time when approximating high-dimensional temporal feasible regions, followed by the Zonotope approximation method, while our method requires relatively longer processing time. The reason for this is that among the three types of approximation methods, the volume box approximation method has a low modeling complexity. However, the model of the method used in this paper is complex, but the overall time is within an acceptable range. Meanwhile, considering the approximation accuracy, the method used in this paper has advantages.

[figure(s) omitted; refer to PDF]

5.3. Optimized Scheduling Considering the Participation of Heterogeneous Load Aggregates in Demand Response

This paper constructs two cases. In Case 1, the influence of the participation of heterogeneous load aggregates in demand response on the operation of the VPP is examined when the VPP is in a load-shedding situation; in Case 2, the impact of the participation of heterogeneous load aggregates in demand response on the economy and cleanliness of the VPP under the situation of high generation of wind and solar resources is explored. Against the background of the two cases, two scheduling optimization strategies are designed. Strategy 1 involves not taking into account plenty of flexible loads in the system scheduling optimization; Strategy 2 is that plenty of flexible loads take part in the demand response regulation.

5.3.1. Large-Scale Load Shedding in the VPP

It can be observed from Figure A1 that the scenario set in this paper experiences load-shedding phenomena from 17:00 to 23:00. In Strategy 1, the energy storage devices are charged intensively during the high-generation period of new energy from 11:00 to 13:00, as shown in Figure 4. The load-shedding cost set in this paper is 10,000 yuan/MW. The cost caused by load-shedding in the system is huge as shown in Table 3. Therefore, the energy storage devices perform discharging operations from 17:00 to 19:00, alleviating the occurrence of load-shedding phenomena during this period. However, load shedding still occurs from 19:00 to 23:00.

[figure(s) omitted; refer to PDF]

Table 3

Comparison of optimization results under load shedding.

Scenario	Generation cost (CHY)	Start-up costs (CHY)	Charge–discharge cost (CHY)	Load curtailment cost (CHY)	Penalties for wind and solar curtailment (CHY)	New energy consumption rate (%)	Demand response cost (CHY)	Total cost (CHY)
No adjustment	61,512.805	23,300	12,236	1,820,275.44	0	99.999	0	1,894,024.245
Participate in regulation	65,210.417	16,300	14,031.2137	50,783.323	0	99.999	17,646.40	149,940.140

Compared with Strategy 1, the energy storage devices are charged intensively from 11:00 to 14:00, and at the same time, the charging amount is increased. Meanwhile, new energy units are in high demand. Due to the wind and solar curtailment penalties set in the case study, flexible loads increase electricity consumption during this period to enhance the utilization rate of new energy. Since the regulation ability of a large number of heterogeneous loads during the period from 17:00 to 18:00 is sufficient to cope with the load-shedding problem caused by excessive system load demand and insufficient supply on the source side, the load aggregates reduce the electricity consumption during this period and delay the energy storage discharging period to 19:00–23:00, avoiding the occurrence of load-shedding phenomena.

Due to Strategy 1, a demand response strategy is adopted on the load side in the system. According to Figure 5, the system dispatches flexible load resources to reduce the load when the system is under great pressure to supply guarantee and increase the load at other times simultaneously. Flexible loads can not only meet their own electricity demands but also reduce the risk of load shedding in the power system, ultimately achieving source–load interaction. Although the system dispatching process bears the demand response cost, it avoids high load-shedding costs and simultaneously enhances the stability of system operation.

[figure(s) omitted; refer to PDF]

5.3.2. High-Generation Situation of PV and WT

It can be observed from Figure A2 that in the scenario set in this paper, the phenomena of wind curtailment and photovoltaic curtailment occur from 01:00 to 16:00. In Strategy 1, the wind curtailment volume is 521.73 MW and the photovoltaic curtailment volume is 449.713 MW. Compared with Strategy 1, in Strategy 2, the wind curtailment volume is 576.251 MW and the photovoltaic curtailment volume is 321.617 MW, reducing the total curtailment volume by 73.585 MW. In Strategy 2, due to the increase in electricity consumption by flexible loads from 00:00 to 16:00, the curtailment volumes of wind and photovoltaic in the system are reduced, as shown in Figure 6. Based on the case study, it is found that in Strategy 2, the wind curtailment volume will be reduced from 00:00 to 08:00. From 09:00 to 14:00, even if the wind curtailment phenomenon still occurs, the reduction of the photovoltaic curtailment volume in the system will be prioritized. This phenomenon occurs because the PV generation is at a high level in this period. Theoretically, the photovoltaic curtailment volume is greater than the wind curtailment volume in the same period. Due to the limited ability of load aggregates to take part in demand response and the direct consumption of photovoltaic energy by energy storage devices, one type of power generation equipment is preferentially selected for consumption.

[figure(s) omitted; refer to PDF]

In the selected scenario, the original load increases from 17:00 to 23:00. At the same time, all types of loads need to meet the power balance constraint within the scheduling cycle. Therefore, Strategy 2 is compared with Strategy 1, and the loads within the aggregated load system reduce electricity consumption from 17:00 to 23:00 through demand response strategies, as shown in Figure 7.

[figure(s) omitted; refer to PDF]

During this period, the thermal power generation decreases. Through the experimental results, it is observed that one of the three thermal power generating units in Strategy 2 does not start within 24 h, reducing the generation cost. Finally, from Table 4, it can be observed that adopting the demand response strategy with the participation of load aggregates can cope with the phenomena of wind and photovoltaic curtailment, enhance the consumption of WT power and PV power, and improve the economic efficiency of system operation.

Table 4

Comparison of optimization results under the condition of high wind and solar incidence.

Scenario	Generation cost (CHY)	Start-up costs (CHY)	Charge–discharge cost (CHY)	Load curtailment cost (CHY)	New energy consumption rate (%)	Penalties for wind and solar curtailment (CHY)	Demand response cost (CHY)	Total cost (CHY)
No adjustment	21,382.901	23,300	12,236	0	76.748	485,726.335	0	542,645.236
Participate in regulation	17,981.864	16,300	12,236	0	78.983	448,933.999	5272.314	500,724.177

6. Conclusions

This paper proposes a multielement heterogeneous load aggregation model and its control method considering DCs, ACs, and EVs. Firstly, a general mechanism model of multielement flexible loads is established, considering the time series power and “energy” coupling characteristics of DCs and constructing the original feasible region of multielement loads. Secondly, using the theory of homogeneous polytopes, the original feasible region of a single load is approximated. The feasible region of the load aggregator is obtained through the method of Minkowski sum. Finally, the minimization of the VPP operation cost is set as the optimization objective, and an operation framework for multielement flexible loads participating in demand response is constructed and verified by case studies under two scenarios: system load shedding and peak generation of wind and solar power.

The current research has limitations. The adjustable capability of flexible loads has irrational characteristics. Therefore, user preferences should be analyzed in depth, a user utility function should be constructed, and the adjustable capability of loads should be more accurately characterized to provide a stronger reference for optimal scheduling.

Disclosure

There are no persons or third-party services involved in the research or manuscript preparation who are not listed as an author and have not been acknowledged.

Author Contributions

Chenzheng Xie participated in the development of the methodology, modeling, software, writing – original draft, and investigation. Rui Zhu participated in the development of software, formal analysis, and visualization. Xinge Xu participated in the data collection, data curation, and the result verification in this manuscript. Yingyun Sun participated in the writing – review and editing, and supervision of the manuscript.

Funding

This work was funded by the National Natural Science Foundation of China (No. U24B2077).

Acknowledgments

No AI software was used, either in whole or in part, during the preparation of this manuscript.

Appendix

Appendix A

[figure(s) omitted; refer to PDF]

References

[1] F. Fang, S. Yu, X. Xin, "Data-Driven-Based Stochastic Robust Optimization for a Virtual Power Plant With Multiple Uncertainties," IEEE Transactions on Power Systems, vol. 37 no. 1, pp. 456-466, DOI: 10.1109/TPWRS.2021.3091879, 2022.

[2] S. Yu, F. Fang, Y. Liu, J. Liu, "Uncertainties of Virtual Power Plant: Problems and Countermeasures," Applied Energy, vol. 239, pp. 454-470, DOI: 10.1016/j.apenergy.2019.01.224, 2019.

[3] IEA, Energy and AI, 2025. https://www.iea.org/reports/energy-and-ai

[4] N. Sadeghianpourhamami, N. Refa, M. Strobbe, C. Develder, "Quantitive Analysis of Electric Vehicle Flexibility: A data-driven Approach," International Journal of Electrical Power & Energy Systems, vol. 95, pp. 451-462, DOI: 10.1016/j.ijepes.2017.09.007, 2018.

[5] J. Cui, M. Gao, X. Zhou, "Demand Response Method Considering Multiple Types of Flexible Loads in Industrial Parks," Engineering Applications of Artificial Intelligence, vol. 122,DOI: 10.1016/j.engappai.2023.106060, 2023.

[6] T. Liu, P. Wang, Q. Peng, M. Zhang, T. Wang, J. Meng, "Operation-Area-Constrained Adaptive Primary Frequency Support Strategy for Electric Vehicle Clusters," Journal of Modern Power Systems and Clean Energy, vol. 11 no. 4, pp. 1982-1994, DOI: 10.35833/MPCE.2023.000233, 2023.

[7] Z. Yi, Y. Xu, W. Gu, W. Wu, "A Multi-Time-Scale Economic Scheduling Strategy for Virtual Power Plant Based on Deferrable Loads Aggregation and Disaggregation," IEEE Transactions on Sustainable Energy, vol. 11 no. 3, pp. 1332-1346, DOI: 10.1109/TSTE.2019.2924936, 2020.

[8] J. Tong, W. Liu, J. Mao, M. Ying, "Role and Development of Thermal Power Units in New Power Systems," IEEE Journal of Radio Frequency Identification, vol. 6, pp. 837-841, DOI: 10.1109/JRFID.2022.3205465, 2022.

[9] R. Zhang, J. Yu, "An Orderly Power Utilization Method for New Urban Power Grids Facing Severe Electricity Shortages," Journal of Modern Power Systems and Clean Energy, vol. 12 no. 6, pp. 1710-1723, DOI: 10.35833/MPCE.2023.000874, 2024.

[10] J. Tu, M. Zhou, H. Cui, F. Li, "An Equivalent Aggregated Model of Large-Scale Flexible Loads for Load Scheduling," IEEE Access, vol. 7, pp. 143431-143444, DOI: 10.1109/ACCESS.2019.2944233, 2019.

[11] J. Wu, W. Zhou, W. Zhong, J. Liu, "Multi-Energy Demand Response Management in Energy Internet: A Stackelberg Game Approach," Chinese Journal of Electronics, vol. 28 no. 3, pp. 640-644, DOI: 10.1049/cje.2019.03.011, 2019.

[12] S. Zheng, Y. Sun, B. Qi, B. Li, "Incentive-Based Integrated Demand Response Considering S and C Effect in Demand Side With Incomplete Information," IEEE Transactions on Smart Grid, vol. 13 no. 6, pp. 4465-4482, DOI: 10.1109/TSG.2022.3149959, 2022.

[13] J. K. Scott, D. M. Raimondo, G. R. Marseglia, R. D. Braatz, "Constrained Zonotopes: A New Tool for Set-Based Estimation and Fault Detection," Automatica, vol. 69, pp. 126-136, DOI: 10.1016/j.automatica.2016.02.036, 2016.

[14] M. Althoff, O. Stursberg, M. Buss, "Computing Reachable Sets of Hybrid Systems Using a Combination of Zonotopes and Polytopes," Nonlinear Analysis: Hybrid Systems, vol. 4 no. 2, pp. 233-249, DOI: 10.1016/j.nahs.2009.03.009, 2010.

[15] E. Öztürk, T. Faulwasser, K. Worthmann, M. PreißInger, K. Rheinberger, "Alleviating the Curse of Dimensionality in Minkowski Sum Approximations of Storage Flexibility," IEEE Transactions on Smart Grid, vol. 15 no. 6, pp. 5733-5743, DOI: 10.1109/TSG.2024.3420156, 2024.

[16] J.-R. Jiang, C.-W. Chu, "Classifying and Benchmarking Quantum Annealing Algorithms Based on Quadratic Unconstrained Binary Optimization for Solving NP-Hard Problems," IEEE Access, vol. 11, pp. 104165-104178, DOI: 10.1109/ACCESS.2023.3318206, 2023.

[17] F. L. Müller, O. Sundström, J. Szabó, "Aggregation of Energetic Flexibility Using Zonotopes," pp. 6564-6569, DOI: 10.1109/CDC.2015.7403253, .

[18] J. Wang, M. Liu, J. Chen, "Research on Multi-Resources Aggregation Methods of User-Side Flexibility Resources Based on Affine Transformations of Convex Polytope," Power System Technology, vol. 68,DOI: 10.13335/j.1000-3673.pst.2024.1595, 2024.

[19] F. L. Muller, J. Szabo, O. Sundstrom, J. Lygeros, "Aggregation and Disaggregation of Energetic Flexibility From Distributed Energy Resources," IEEE Transactions on Smart Grid, vol. 10 no. 2, pp. 1205-1214, DOI: 10.1109/TSG.2017.2761439, 2019.

[20] Z. Yi, Y. Xu, W. Gu, L. Yang, H. Sun, "Aggregate Operation Model for Numerous Small-Capacity Distributed Energy Resources Considering Uncertainty," IEEE Transactions on Smart Grid, vol. 12 no. 5, pp. 4208-4224, DOI: 10.1109/TSG.2021.3085885, 2021.

[21] B. Moarref, S. M. Abedi, S. M. Modaresi, "Optimal Scheduling of Power System Integrated With Hydrogen Vehicles and Flexible Resources: a Hybrid Uncertainty Management Method," International Journal of Hydrogen Energy, vol. 100, pp. 658-667, DOI: 10.1016/j.ijhydene.2024.12.042, 2025.

[22] L. Dong, Q. Wu, J. Hong, Z. Wang, S. Fan, G. He, "An Adaptive Decentralized Regulation Strategy for the Cluster With Massive Inverter Air Conditionings," Applied Energy, vol. 330,DOI: 10.1016/j.apenergy.2022.120304, 2023.

[23] Y. Zhang, H. Zhao, Z. Qi, B. Li, "A Two-Stage Low-Carbon Economic Coordinated Dispatching Model for Generation-Load-Storage Resources Considering Flexible Supply-Demand Balance," Applied Energy, vol. 373,DOI: 10.1016/j.apenergy.2024.123981, 2024.

[24] Z. Zhu, M. Wang, Z. Xing, Y. Liu, S. Chen, "Optimal Configuration of Power/Thermal Energy Storage for a Park-Integrated Energy System Considering Flexible Load," Energies, vol. 16 no. 18,DOI: 10.3390/en16186424, 2023.

[25] F. Pallonetto, S. Oxizidis, F. Milano, D. Finn, "The Effect of Time-of-Use Tariffs on the Demand Response Flexibility of an All-Electric Smart-Grid-Ready Dwelling," Energy and Buildings, vol. 128, pp. 56-67, DOI: 10.1016/j.enbuild.2016.06.041, 2016.

[26] M. Bagheritabar, S. M. Hakimi, G. Derakhshan, A. Rezaee Jordehi, "A Three-Stage Optimization Framework for Unlocking demand-side Flexibility in Highly Renewable Electricity Grids," Energy, vol. 320,DOI: 10.1016/j.energy.2025.135158, 2025.

[27] M. Ebrahimi, M. Shafie-khah, H. Laaksonen, "Uncertainty-Observed Virtual Battery Model for an Electric Vehicle Parking Lot Enabling Charger-Sharing Modelling," Journal of Energy Storage, vol. 89,DOI: 10.1016/j.est.2024.111578, 2024.

[28] J. T. Hughes, A. D. Domínguez-García, K. Poolla, "Identification of Virtual Battery Models for Flexible Loads," IEEE Transactions on Power Systems, vol. 31, pp. 4660-4669, DOI: 10.1109/TPWRS.2015.2505645, 2016.

[29] J. Wan, Y. Duan, X. Gui, C. Liu, L. Li, Z. Ma, "SafeCool: Safe and Energy-Efficient Cooling Management in Data Centers With Model-Based Reinforcement Learning," IEEE Transactions on Emerging Topics in Computational Intelligence, vol. 7 no. 6, pp. 1621-1635, DOI: 10.1109/TETCI.2023.3234545, 2023.

[30] Y. Cao, M. Cheng, S. Zhang, "Data-Driven Flexibility Assessment for Internet Data Center Towards Periodic Batch Workloads," Applied Energy, vol. 324,DOI: 10.1016/j.apenergy.2022.119665, 2022.

[31] F. Chen, J. Grundy, Y. Yang, J.-G. Schneider, Q. He, "Experimental Analysis of Task-based Energy Consumption in Cloud Computing Systems," Proceedings of the 4th ACM/SPEC International Conference on Performance Engineering (ICPE ‘13), pp. 295-306, DOI: 10.1145/2479871.2479911, .

[32] M. S. Hossain, S. Rahaman, A.-L. Kor, K. Andersson, C. Pattinson, "A Belief Rule Based Expert System for Datacenter PUE Prediction Under Uncertainty," IEEE Transactions on Sustainable Computing, vol. 2, pp. 140-153, DOI: 10.1109/TSUSC.2017.2697768, 2017.

[33] Illinois Institute of Technology, IIT. Index of Data Illinois Institute of Technology, 2016. https://motor.ece.iit.edu/data/

[34] K. Ganz, A. Guminski, M. Kolb, S. von Roon, "Normalized Industrial Electrical Load Profiles (Germany)," The FfE Open Data Portal, 2021. https://opendata.ffe.de/dataset/normalized-industrial-electrical-load-profiles-germany/

[35] Z. Ding, L. Xie, Y. Lu, P. Wang, S. Xia, "Emission-Aware Stochastic Resource Planning Scheme for Data Center Microgrid Considering Batch Workload Scheduling and Risk Management," IEEE Transactions on Industry Applications, vol. 54 no. 6, pp. 5599-5608, DOI: 10.1109/TIA.2018.2851516, 2018.

[36] C. Xie, "Optimization Scheduling Method for Demand Response Considering Heterogeneous Load Aggregators," OSCHINA, . https://gitee.com/xie-chenzheng/optimization-scheduling-method-for-demand-response-considering-heterogeneous-load-aggregators