# Aggregation of building predictive energy flexibility in smart microgrid

Zhigang Zhang, Hervé Guéguen

► To cite this version:

Zhigang Zhang, Hervé Guéguen. Aggregation of building predictive energy flexibility in smart microgrid. International Journal of Electrical Power & Energy Systems, 2024, 160, pp.110073. 10.1016/j.ijepes.2024.110073. hal-04612757

HAL Id: hal-04612757

https://hal.science/hal-04612757v1

Submitted on 27 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Aggregation of building predictive energy flexibility in smart microgrid

Zhigang Zhang*, Hervé Guéguen

IETR, CentraleSupélec, 5 Avenue de la Boulée, Cesson-Sévigné, 35510, Ile-et-Vilaine, France

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Flexibility

Aggregation

Approximation

Building

Demand response

## ABSTRACT

Numerous demand-side resources possess significant potential to offer substantial flexible reserves and ancillary services, thereby enhancing the economic and technical performance of power systems within the demand response domain. A novel approach to aggregating the energy flexibility of buildings based on the virtual battery model is proposed in the smart microgrid, employing predictive control techniques while ensuring consumer privacy. The individual approximation process is implemented to extract the flexible capacity of local resources while protecting information privacy. Subsequently, the aggregator collects the flexibility information to develop an optimization strategy for maximum flexibility aggregation at the upper level, enabling the broader exploitation of global flexibility. The feasibility of power disaggregation is guaranteed according to the aggregation procedure. Furthermore, the comprehensive framework accommodates temporal preferences for flexibility requests from the aggregator. Numerical simulations are conducted to perform a comparative analysis with existing state-of-the-art methods. The results demonstrate that the proposed approach exhibits improved adaptability in characterizing the aggregate flexibility from buildings, effectively managing the heterogeneity of load characteristics and the variability in the collected flexibility information. Microgrids can realize cost savings and achieve dynamic flexibility utilization. This methodology can be extended and applied to diverse types of flexible resources.

## 1. Introduction

The increasing integration of renewable energy sources (RESs) and the significant rise in energy consumption present challenges for energy systems management, such as maintaining the power balance between supply and demand, ensuring effective operation, increasing energy use during peak periods, and managing reliability issues [1,2]. Compared to conventional fast-ramping generation units that would require significant investments in the power system, the energy flexibility can provide a promising low-cost alternative in the energy transition, which is predominantly realized through the integration of decentralized resources into demand response (DR) programs [3]. The rapid development of smart grid technologies has precipitated increased demand-side flexibility to provide operating reserves, which confers numerous advantages over supply-side flexibility, especially in improving the economic efficiency and reliability of the grid. Demand-side flexibility is predominantly achieved by integrating decentralized resources into demand response programs. The distribution system operator (DSO) often assumes the role on the grid side of communicating with DR service providers. In the smart grid, DR represents a favorable facet that enables consumers to modify the energy demand according to financial factors intelligently [4]. DR facilitates the active demand-side participation in grid operations, making better use of power network

capacities and resulting in improved power system resilience, economic efficiency as well as environmental effects [5]. DR measures are typically realized through the coordination of distributed energy resources (DERs), such as thermostatically controlled loads (TCLs), electric vehicles (EVs), and energy storage systems (ESSs). Among the various energy consumers, buildings constitute one of the largest sectors in terms of energy consumption. According to [6], buildings account for 40% of global energy consumption and 15% of total emissions in end-use sectors. Therefore, by altering the energy usage while preserving thermal comforts and convenience, buildings possess great potential to make substantial contributions towards providing energy flexibility for power grids.

In the academic literature, the definitions of demand-side energy flexibility primarily focus on three criteria: resource type, duration, and activation payments [7]. We quantify building flexibility as the technical ability of buildings to modulate their baseline consumption over a specific time frame [8]. The baseline level refers to the normal energy utilization of consumers without DR events, serving as a benchmark for assessing the building's capacity to provide flexibility [9]. Typically, individual flexibility capacity is characterized by the set of all admissible power profiles, which specify maximum upward and

<footnote>* Corresponding author.
E-mail addresses: zhigang.zhang@centralesupelec.fr (Z. Zhang), hervé.gueguen@centralesupelec.fr (H. Guéguen).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Schematic illustration of building thermal flexibility.</figcaption></figure>

downward power deviations from the baseline, as well as maximum energy (thermal) deviation.

Fig. 1 illustrates how a building's flexibility within its capacity can vary over different periods. The blue lines depict the baseline scenario's energy utilization and temperature development. Under the condition that the temperature remains within a comfortable range for occupants, as shown by the black lines in Fig. 1(a), the building can adjust its power demand freely. In Fig. 1(b), the orange line exemplifies a power increase relative to the baseline. The corresponding temperature profile, described by the orange curve in Fig. 1(a), naturally records higher values than the normal temperatures in the given period. We define such pattern of power augmentation as upward flexibility. Furthermore, constrained by the temperature bounds, the potential for offering flexibility is inherently limited and fluctuates over time. In the case depicted by the orange curves, the building exploits much of its flexibility capacity in the beginning stage, resulting in the temperature rapidly approaching its upper threshold. In contrast, the green curves illustrate an alternate scenario where the building gradually enhances its upward flexibility towards the end of the period, strategically utilizing upward capacity more extensively later on. Similarly, the power decrease illustrated by the red curves in Fig. 1 means the downward flexibility development.

The efficient DR management of potential distributed DERs in the utility grid represents a complex issue. The local DERs are often distinguished by their extensive heterogeneity, decentralized geographical distribution, and significant stochastic variability [10]. Meanwhile, in practical applications, incorporating distributed flexibility information directly into grid operations proves impracticable due to the elevated computation burden, complicated regulation frameworks, and the unpredictable errors in forecasts attributed to individual consumers [11]. To enhance the reliability of power system operations, it is imperative to develop effective collaborative coordination and control mechanisms for DERs including buildings [12]. Particularly, the employment of appropriate aggregation methodologies has come into view as a viable and promising strategy to address these challenges, with smart microgrids (MGs) recognized as vital assets in these implementations [13]. A microgrid is a small-scale power system, generally composed of distributed generation sources, energy storage, energy conversion equipment, and a set of loads [14]. To encourage flexible energy contribution, an incentive-based DR strategy is advocated, wherein local consumers are financially rewarded with payments for successfully adjusting their energy usage. In terms of energy aggregation, the microgrid operator (MGO) coordinates local-level flexible

loads and acts as the upper-level aggregator, serving as an intermediary between the grid and resources. Furthermore, the aggregation approach should characterize the feasible region for time-varying power, thus facilitating the appropriate integration and representation of DERs within the network.

In microgrid management applications, we prioritize online or real-time exploitation of building flexibility and dynamic discrete control strategies through decentralized approaches. Model predictive control (MPC) is recognized as a highly effective solution for this purpose. MPC is a receding-horizon scheme that repeatedly predicts future control signals, solves a constrained optimization problem using the updated current state and determines the optimal control sequence across a finite moving horizon [15]. MPC could well facilitate the quantification and exploitation of demand-side flexibility in microgrids, irrespective of the buildings' typology and classification (e.g., residential, educational, commercial, institutional) [16,17]. Within the predictive framework, the MGO serves as a microgrid energy management system (MEMS) that gathers data, including forecasts of solar generation, load consumption patterns, and updated information from the grid side.

To enhance clarity in describing the microgrid's aggregation framework, Fig. 2 depicts the generalized control framework into three stages: individual flexibility approximations, total flexibility aggregation, and flexibility disaggregation. Initially, the DSO presents the grid pricing scheme for the forthcoming predictive horizon to the microgrid. Concurrently, the MGO sets internal flexibility prices for all local consumers. Subsequently, local buildings establish the baseline and forecast their actual predictive flexibility, which may exhibit diverse characteristics and formulations. Therefore, an approximation process is essential to standardize the individual capacity information into a uniform format for control purposes. The orange diamond block in the diagram denotes this approximation measure of reformulating individual flexibility as stage 1. Depending on specific conditions and considerations, such as communication efficiency and privacy concerns, the approximation stage may be conducted by the aggregator or by the local flexible resources. Advanced techniques in information gathering, control, and communication enable the aggregator to consolidate these flexibility capacities into a simple and compact aggregation proposal during stage 2. This aggregated flexibility model is designed to optimize total flexibility potential and could be integrated into the grid's operations, simplifying system functionality by reducing the number of active nodes and centralizing communications with the aggregators [18]. This model is then applied to regulate the global flexibility demand. The process of determining this demand, represented by the green diamond block, varies based on the scenarios and may be executed by either the microgrid aggregator or the grid operator. Finally, flexibility disaggregation is conducted in stage 3 to decide the individual power decisions. To prevent scenarios where local buildings cannot implement the power assignments issued by the aggregator in the DR programs, we hold the principle that microgrid flexible energy delivery is outlined within contractual agreements among participants, commonly known as "contracted flexibility". The feasibility of energy disaggregation must be thoroughly considered throughout the aggregation process to ensure that the power decisions align with the commitments of all involved agents [19].

According to Fig. 2, two primary challenges arise in the aggregation problems. The first pertains to precise estimations to better capture the actual flexibility capacity of each DER, termed *individual flexibility approximation*; the second problem is related to exploring the total available flexibility from all DERs, referred to as *total flexibility aggregation*. The flexibility of various DERs can be defined and expressed in multiple forms, and it is crucial that the flexibility set information in the system adheres to a unanimously agreed format. Without sacrificing generality and considering regulation power as the input and the state of DER as the output, convex polytopes are widely employed to depict flexibility information precisely. The preferred representation of a convex polytope is expressed as the intersection of a finite number of closed
<figure><img src="image_5.png" /><figcaption>Fig. 2. Generalized predictive control framework of building flexibility aggregation in microgrid. Each building assesses the actual flexibility capacity. Next, Stage 1 processes the heterogeneous local flexibility information into a standardized format through individual approximations, implemented at the local or aggregator level. Subsequently, the microgrid calculates the aggregate flexibility model in Stage 2, which is used to ascertain the global flexibility demand. This demand may be determined by either the microgrid aggregator or the grid operator. In Stage 3, the aggregator dispatches power schedules to individual agents.</figcaption></figure>

half-spaces (referred to as the H-representation of polytope) [20]. The exact aggregate set is equivalent to the Minkowski sum of all the individual proposals, which is computationally challenging and intractable in general. Consequently, existing research predominantly focuses on utilizing simplified approximation models to estimate the Minkowski sum. The outer approximation methods tend to overestimate the total flexibility, potentially leading to impractical scenarios for final power decision-making. Therefore, despite flexibility sacrifices, we use conservative inner estimation methods for individual approximations and aggregation to ensure that all power decisions are feasible.

Numerous endeavors have been made to develop the aggregation model. The dominant geometry-based methods seek to evaluate the aggregate flexibility set by formulating the largest convex inner approximation of local resources, thereby simplifying the actual set's dimensionality. A distinct advantage is that the aggregate model is readily derived once the individual approximations are completed. Various fundamental geometries, often referred to as base sets or prototype sets, have been employed. Inner approximations employing a hyperbox (axis-aligned box) are noted for their conservatism but are highly efficient for computing the Minkowski sum through the straightforward addition of constraint limits [21]. Ellipsoidal approximations are another technique for individual flexibility capacity estimation [22,23]. Homothetic transformations, which involve the translation and dilation of a predefined base set, are widely employed due to their rapid computation of Minkowski sum [13,24]. Zonotope, a specific class of homothets (centrally symmetric polytopes), can also serve as the base set in researches such as [25,26]. Concurrently, many previous studies have leveraged the concept of the virtual battery (VB) [27,28]. This succinct model, defined by a few parameters, could effectively describe the characteristics of actual resource information.

Another significant category of aggregating flexibility involves non-deterministic approaches. These methods employ stochastic or random algorithms such as the Monte-Carlo sampling technique for the approximations [29]. Heuristic methods provide a practicable alternative for generating sampling points, but the effectiveness depends significantly on the chosen sampling size and technique [30]. The data-driven approaches typically necessitate substantial amounts of measurement data, which may be costly or unavailable [31,32]. Additionally, we could also identify methods designed for specific requirements and presented in distinctive forms, such as aggregation through the unions

of hyperboxes [33], and the inward-to-outward progressive expansion algorithm [34].

In addition to improving the characterization of flexibility capacity, our work identifies several research gaps associated with previous discussions on constructing aggregated energy flexibility. First, many aggregate methods are restricted in the application scope. For example, many flexibility aggregation methods are highly specialized for frequently studied loads, such as TCLs [27] and EVs [35,36]. Typically, TCLs are modeled with fixed constraint values, and the polytope matrices for EV-type load descriptions are identical, significantly simplifying the aggregation complexity. In contrast, buildings in the microgrid may exhibit vastly different properties, with power and energy constraints that are both unpredictable and time-variant. Accordingly, some of the latest methods to enhance aggregation quality may not be applied to scenarios involving building flexibility, such as illustrated in [37,38].

Second, it is generally presumed that aggregators have complete knowledge of local flexibility information without considering the protection of consumer privacy. DERs transmit their actual flexibility directly to the aggregator under this full-information assumption. It should be emphasized that the aggregator executes all individual flexibility approximations and the subsequent aggregation. However, building owners could perceive their actual capacities as confidential and be reluctant to provide them in DR programs. Suppose we continue this approach while incorporating privacy considerations, it becomes essential for each building to perform a preliminary local-level approximation to propose the flexibility set before it is further processed at the aggregator level. The apparent drawback following this plan is that two successive approximations may lead to a substantial degree of conservatism in the assessed flexibility potential. Therefore, incorporating measures to protect local privacy poses a challenge to the system.

Third, we also take into account the varying requirements for flexibility capacity related to time for microgrid management. Most methods assume that the significance of flexible contributions remains constant across the horizon, and a common objective in the optimizations is to regulate the polytope volume through indirect indices, such as scaling factor [24] or matrix determinant [37]. However, in practice, upon receiving real-time market power price schemes, the aggregator may prioritize certain periods over others, seeking enhanced flexibility utilization during some intervals to reduce operational costs. The temporal preferences introduce additional complexities in evaluating aggregation quality, as the volume metric no longer suits this
context. To the best of our knowledge, research that integrates temporal importance factors into aggregation methodologies is currently limited. An alternative criterion is required to accommodate these distinctions.

Motivated by the aforementioned challenges, this article proposes a hierarchical aggregation framework to exploit the predictive energy flexibility from heterogeneous buildings within microgrid operations. The main contributions are outlined as follows:

1. Maximum inner approximation optimization is formulated to characterize the individual building flexibility. The approximation treats resources uniformly as virtual batteries, effectively representing the geometric capacity of actual DER flexibility and reducing the problem complexity. The methodology captures the features of load flexibility and offers advantageous solutions in terms of the size of feasible regions. This flexibility approximation process is compatible with diverse DERs and applicable to those described by high-order dynamics.

2. The aggregation method is formulated to enhance the exploitation of total flexibility potential with lower conservatism. This process fully considers the heterogeneity of local resources. The feasibility of power disaggregation is consistently ensured.

3. The research addresses the privacy protection requirements of local consumers and incorporates the temporal preferences in flexibility capacity, which have received limited attention in the existing studies.

The rest of this article is organized as follows. Section 2 constructs the individual approximation process for building flexibility. Section 3 presents the aggregation optimization process of total flexibility capacity. Simulation experiments are performed to demonstrate the effectiveness of the proposed method in Section 4. Section 5 highlights the main conclusions and discusses future work.

Throughout the remainder of the paper, H-polytope is used for flexibility description, which is $P = \{u | Hu \le d\}$, where H and d are the matrix and vector of appropriate dimensions, respectively. Given our focus on convex shapes, the term 'convex' is assumed and not explicitly stated. Besides, we consider time as discretized and introduce the following notations: Instant k represents the current moment at which all predictive controllers are implemented; $T_k$ is the sampling time; $\mathcal{N} := \{0, 1, ..., N-1\}$ depicts the control horizon starting at time k, while N will be adequately selected based on the contractual specifications. Moreover, vectors or variable profiles are typically denoted in bold.

## 2. Individual flexibility approximation

This section aims to design the individual approximation process to exploit the building's energy flexibility capacity, while considering the local privacy protection and incorporating the temporal demands from the upper-level aggregator.

### 2.1. Representation of actual building flexibility

#### 2.1.1. Model setting and baseline determination

Accurate mathematical modeling is essential to ensure satisfactory prediction capabilities, especially in capturing the coupling characteristics of building power and heat. Various approaches are available for analyzing the thermal performance of buildings, such as steady-state models and dynamic models [39]. For this study, we have selected the commonly utilized second-order thermal network described in [40], which effectively characterizes the thermal capacity and resistance of buildings and facilitates the investigation of heat performance over long-time horizons.

We denote the thermal state variable of building as $x_{th}(i) = [\theta_{th}(i) \quad \theta_{w}(i)]^T \in \mathbb{R}^2$, including the internal temperature $\theta_{th}$ and the wall surface temperature $\theta_{w}$. The disturbances signal $\theta_{d}$ accounts for the ambient temperature. The heating power delivered to the

<figure><img src="image_6.png" /><figcaption>Fig. 3. Building nominal and enlarged thermal bounds.</figcaption></figure>

building is represented by $u_{th}$. The building's thermal dynamics are then discretized using the following equations:

$$x_{th}(i+1) = A_{th}x_{th}(i) + B_{th}^{u}u_{th}(i) + B_{th}^{d}\theta_{d}(i) \quad (1a)$$

$$y_{th}(i) = C_{th}x_{th}(i) \quad (1b)$$

where $A_{th}$, $B_{th}^{u}$, $B_{th}^{d}$ are constant matrices of appropriate dimensions, and $y_{th}$ equals the building's inner temperature.

The baseline estimation is a critical component in incentive-based DR programs to ensure fair and reasonable settlements for energy management [8]. Each building predicts its baseline consumption, $u_{th}^{bs}(k) = \{u_{th}^{bs}(i)\}_{i \in \mathcal{N}}$, prior to quantifying its subsequent flexibility capacity. In our approach, each building seeks to minimize energy cost and maintain thermal comfort, considering the nominal power price and potential discomfort penalties. The detailed methodology is developed in our previous work [41].

#### 2.1.2. Polytopic representation of actual building flexibility

The building could modulate its thermal behaviors, and the power variation could offer upward flexibility $u_{th}^{uw}$ or downward flexibility $u_{th}^{dw}$ to the MGO. The building's actual thermal power consumption $u_{th}$ is:

$$u_{th}(i) = \begin{cases} u_{th}^{bs}(i) + u_{th}^{uw}(i), & u_{th}^{f}(i) = u_{th}^{uw}(i) \ge 0 \\ u_{th}^{bs}(i) + u_{th}^{dw}(i), & u_{th}^{f}(i) = u_{th}^{dw}(i) \le 0 \end{cases} \quad (2)$$

where $u_{th}^{f}$ is used to describe the flexibility for both scenarios.

Because of the incentive rewards, it is beneficial for the building to enhance its flexibility capacity by enlarging the range of thermal comfort envelope, indicated by Eq. (3) and depicted in Fig. 3:

$$\underline{\theta}_{th}(i) \le \underline{\theta}_{th}^{n}(i) \le \theta_{th}(i) \le \bar{\theta}_{th}^{n}(i) \le \bar{\theta}_{th}(i) \quad (3)$$

where $\underline{\theta}_{th}^{n}$ and $\bar{\theta}_{th}^{n}$ represent the nominal thermal bounds; $\underline{\theta}_{th}$ and $\bar{\theta}_{th}$ are used as the expanded thermal bounds.

A representative case study for determining the optimal thermal envelope is illustrated in [42], in which an optimization problem is formulated considering the economic returns and consumer discomfort penalties. For the sake of simplicity, we utilize the maximally enlarged limits as the thermal bounds to quantify the building's flexibility. Consequently, any power profile that complies with the new thermal constraints can be considered as a valid flexibility proposal.

Within the MPC framework over the horizon $\mathcal{N}$, the actual thermal flexibility capacity, denoted as $P_{th}^{rl}$, is defined as the set of all admissible power profiles that comply with the following conditions: $\forall i \in \mathcal{N}$,

$$P_{th}^{rl} = \left\{ u_{th}^{f}(k) \begin{array}{l} x_{th}(i+1) = A_{th}x_{th}(i) + B_{th}^{u}u_{th}^{bs}(i) + B_{th}^{d}u_{th}^{dw}(i) \\ y_{th}(i) = C_{th}x_{th}(i) \\ \underline{u}_{th}(i) \le u_{th}^{bs}(i) + u_{th}^{f}(i) \le \bar{u}_{th}(i) \\ \underline{\theta}_{th}(i+1) \le \theta_{th}(i+1) \le \bar{\theta}_{th}(i+1) \end{array} \right\} \quad (4)$$
where $\mathbf{u}_{th}^f$ is the thermal flexibility vector.

For the following sections, we denote the thermal temperature profile, as well as the vectors representing its lower and upper bounds, as $\theta_{th}(k+1) = \{\theta_{th}(i+1)\}_{i \in \mathcal{N}}$, $\theta_{th}(k+1) = \{\theta_{th}(i+1)\}_{i \in \mathcal{N}}$, and $\bar{\theta}_{th}(k+1) = \{\bar{\theta}_{th}(i+1)\}_{i \in \mathcal{N}}$, respectively. We also represent the disturbance trajectory as $\theta_d(k) = \{\theta_d(i)\}_{i \in \mathcal{N}}$. The input power and the respective bound profiles are expressed in a similar manner as $\mathbf{u}_{th}(k) = \{\mathbf{u}_{th}(i)\}_{i \in \mathcal{N}}$, $\mathbf{u}_{th}(k) = \{\mathbf{u}_{th}(i)\}_{i \in \mathcal{N}}$ and $\bar{\mathbf{u}}_{th}(k) = \{\bar{\mathbf{u}}_{th}(i)\}_{i \in \mathcal{N}}$. By rewriting Eq. (1), the building thermal vector could be expressed as:

$$ \theta_{th}(k+1) = \bar{\mathbf{C}}_{th}^{x} x_{th}(k) + \bar{\mathbf{C}}_{th}^{u} \mathbf{u}_{th}^{bs}(k) + \bar{\mathbf{C}}_{th}^{u} \mathbf{u}_{th}^{f}(k) + \bar{\mathbf{C}}_{th}^{d} \theta_{d}(k) \quad (5) $$

By means of integrating $\theta_{th}(k+1)$ into Eq. (4), the H-polytopic form of $P_{th}^{rl}$ is shown as follows:

$$ P_{th}^{rl} = \left\{ \mathbf{u}_{th}^{f} \middle| H_{th}^{rl} \mathbf{u}_{th}^{f}(k) \leq d_{th}^{rl} \right\} \quad (6a) $$

$$ H_{th}^{rl} = [-I_N, I_N, -\bar{\mathbf{C}}_{th}^{u}, \bar{\mathbf{C}}_{th}^{u}]^T \quad (6b) $$

$$ d_{th}^{rl} = [-\bar{\mathbf{e}}_{th}^{rl,u}, \bar{\mathbf{e}}_{th}^{rl,u}, -\bar{\mathbf{e}}_{th}^{rl,x}, \bar{\mathbf{e}}_{th}^{rl,x}]^T \quad (6c) $$

where the values of the right-hand vector $d_{th}^{rl}$ are:

$$ \bar{\mathbf{e}}_{th}^{rl,u} = \mathbf{u}_{th}(k) - \mathbf{u}_{th}^{bs}(k) \quad (7a) $$

$$ \bar{\mathbf{e}}_{th}^{rl,u} = \bar{\mathbf{u}}_{th}(k) - \mathbf{u}_{th}^{bs}(k) \quad (7b) $$

$$ \bar{\mathbf{e}}_{th}^{rl,x} = \bar{\theta}_{th}(k+1) - \bar{\mathbf{C}}_{th}^{x} x_{th}(k) - \bar{\mathbf{C}}_{th}^{u} \mathbf{u}_{th}^{bs}(k) - \bar{\mathbf{C}}_{th}^{d} \theta_{d}(k) \quad (7c) $$

$$ \bar{\mathbf{e}}_{th}^{rl,x} = \bar{\theta}_{th}(k+1) - \bar{\mathbf{C}}_{th}^{x} x_{th}(k) - \bar{\mathbf{C}}_{th}^{u} \mathbf{u}_{th}^{bs}(k) - \bar{\mathbf{C}}_{th}^{d} \theta_{d}(k) \quad (7d) $$

Ultimately, the building's actual thermal flexibility capacity $P_{th}^{rl}$ is represented as an H-polytope, which is a full-dimensional convex polytope as long as not empty [43].

Apart from thermal energy, other types of resources, such as batteries, can also integrate and provide flexibility in a microgrid. Additionally, we can formulate the actual flexibility capacity of these resources into a convex polytope, following a similar procedure as described above. The process is not detailed for simplification.

## 2.2. Individual flexibility approximation optimization

The individual flexibility approximation process strives to preserve the building's flexibility capacity while ensuring privacy protection.

### 2.2.1. Virtual battery model

We use the concept of the virtual battery as the basis to implement individual flexibility approximations. This method is simple and defines the state of charge (SOC) as the energy state, while simulating the load dynamics using a set of standard parameters. One notable advantage of virtual battery is its adaptability to the flexibility characteristics of many types of DERs, making it an effective tool for promoting efficient energy utilization in microgrids. Furthermore, this method defines the set of feasible power regions with limits on state and power to meet real-time requirements [44]. For the predictive control period $\mathcal{N}$ beginning at time $k$, the set of feasible power signals $\mathbf{u}(k) = \{\mathbf{u}(i)\}_{i \in \mathcal{N}}$ is described by:

$$ P = \left\{ \mathbf{u}(k) \mid \begin{array}{l} x(i+1) = ax(i) + u(i), \forall i \in \mathcal{N} \\ \underline{u}(i) \le u(i) \le \bar{u}(i), \forall i \in \mathcal{N} \\ \underline{x}(i+1) \le x(i+1) \le \bar{x}(i+1), \forall i \in \mathcal{N} \end{array} \right\} \quad (8) $$

where $\mathbf{u}$ is the charge or discharge power from the battery; $x$ is the state to represent the battery energy; $a$ is the dissipation rate characterizing the dynamic property. Each VB model is fully specified by these parameters: $\Gamma_{VB} := (a, x(k), \underline{u}(k) = \{\underline{u}(i)\}_{i \in \mathcal{N}}, \bar{u}(k) = \{\bar{u}(i)\}_{i \in \mathcal{N}}, \underline{x}(k+1) = \{\underline{x}(i+1)\}_{i \in \mathcal{N}}, \bar{x}(k+1) = \{\bar{x}(i+1)\}_{i \in \mathcal{N}})$, in which $x(k)$ is the initial state; From Eq. (8), we denote the VB state profile $x(k)$ as

$$ x(k+1) = \bar{A}x(k) + \bar{B}\mathbf{u}(k) \quad (9) $$

where

$$ \bar{A} = [a, a^2, \dots, a^N]^T \in \mathbb{R}^N \quad (10a) $$

$$ \bar{B} = \begin{bmatrix} 1 \\ a & 1 \\ a^2 & a & 1 \\ \vdots & \vdots & \vdots \\ a^{N-1} & a^{N-2} & \dots & a & 1 \end{bmatrix} \in \mathbb{R}^{N \times N} \quad (10b) $$

We describe Eq. (8) into the H-polytopic form:

$$ P = \{ \mathbf{u}(k) \mid H\mathbf{u}(k) \le d \} \quad (11a) $$

$$ H = [-I_N, I_N, -\bar{B}, \bar{B}]^T \quad (11b) $$

$$ d = [-\underline{e}^u, \underline{e}^u, -\underline{e}^x, \underline{e}^x]^T \quad (11c) $$

where the right vector $d$ is defined by:

$$ \underline{e}^u = \underline{u}(k) \quad (12a) $$

$$ \bar{\underline{e}}^u = \bar{\underline{u}}(k) \quad (12b) $$

$$ \underline{e}^x = \underline{x}(k+1) - \bar{A}x(k) \quad (12c) $$

$$ \bar{\underline{e}}^x = \bar{\underline{x}}(k+1) - \bar{A}\underline{x}(k) \quad (12d) $$

The feasible region of the VB model is determined by the values of $a$ and $d$. Furthermore, we could assume throughout this paper that $d$ consists only of nonnegative values and consistently contains the origin point, which allows for zero power consumption for any flexible source in real scenarios.

### 2.2.2. Computation of individual flexibility approximation model

We investigate an individual flexibility approximation technique using the VB model, whose solution is denoted as $P$. The mathematical expression of our objectives involves maximizing the feasible region of the approximation model and ensuring that it is contained within the building's actual flexibility capacity $P^{rl}$.

The fundamental basis of our problem lies in establishing an optimization objective to characterize the size of the polytope region. The geometric volume serves as a direct and intuitive metric for quantifying the size of a convex polytope, as employed in many previous studies [22,24]. Unfortunately, the volume determination is computationally challenging of complexity #NP-hard [45]. Current methodologies offer numerical estimations, such as through sampling methods [46], yet these require substantial computational resources. As far as we know, a general and explicit algorithm to efficiently compute the polytope volume is unavailable.

Furthermore, as previously discussed, MGO may exhibit distinct expectations for flexibility utilization based on factors like anticipated future prices or other considerations. Certain time intervals might be prioritized to have more considerable flexibility potential. The time-dependent requirements can be inferred from the microgrid price plan, with the condition that buildings contribute a higher level of flexibility capacity during instants characterized by greater incentive rewards. To quantitatively accommodate this demand, each local building defines its time-weight coefficients: $\omega^{bd} = \{\omega_i^{bd}\}_{i \in \mathcal{N}}$. The values of $\omega^{bd}$ are determined based on the principle that moments with higher time weights should correspond to increased levels of flexibility incentives. Given these complexities, integrating volume to evaluate flexibility potential in our optimization endeavors is impractical, primarily due to the absence of a suitable estimation formula and its inability to capture temporal sensitivities.

For the aim of polytope region expansion, we introduce another effective measure offset difference, illustrated in Fig. 4 [26,47]. From Eq. (11), we know that $P$ is composed of $4N$ halfspaces, denoted by: $h_i, \mathbf{u}(k) \le d_i, i = 1, 2, \dots, 4N$. It is straightforward to find that increasing the size of $P$ in a specific direction $h_i$ generally results in a more extensive projection $\rho_i$ along this direction. We denote the
<figure><img src="image_7.png" /><figcaption>Fig. 4. Schematic illustration of offset difference. The outer polytope $P^{rl}$ represents the given flexibility capacity. $P$ is the inner VB approximation model. The offset difference along direction $h_i$ is defined as $|\rho_i^{rl} - \rho_i|$. The facets of the polytopes are degenerated into one-dimensional edges for the purpose of visualization.</figcaption></figure>

projection distance from the origin point of outer polytope $P^{rl}$ over the same direction as $\rho_i^{rl}$. Apparently, the greater the extension of inner $P$ along $h_i$, the smaller the difference between these two projections $|\rho_i^{rl} - \rho_i|$ becomes, which is the corresponding offset difference along the direction $h_i$. We describe the projection distances of the inner polytope $P$ along its halfspace normal vectors as $\rho = \{\rho_i\}_{i=1}^{4N}$. Likewise, the vector $\rho_i^{rl} = \{\rho_i^{rl}\}_{i=1}^{4N}$ symbolizes the projection distances of the outer polytope $P^{rl}$ along these directions.

Furthermore, $\rho$ is determined by the right-hand vector $d$ value from Eq. (11c). Accordingly, while meeting the inner containment relation, we strive to maximize the range of each $d_i$ to identify the largest approximation model. This can be accomplished by minimizing the values of the offset difference vector:

$$ \delta^{od} = \{\delta_i^{od}\}_{i=1}^{4N}, \delta_i^{od} = |\rho_i^{rl} - \rho_i| \quad (13) $$

To simplify the projection computation, we normalize the halfspace direction vectors $h_i$ as $\hat{h}_i$. Subsequently, the equivalent expression of VB model Eq. (11) is derived as Eq. (14), and it is straightforward to deduce that $\rho = \hat{d}$.

$$ P = \{u(k) | \hat{H}u(k) \le \hat{d}\} \quad (14a) $$

$$ \hat{H} = \{\hat{h}_i\}_{i=1}^{4N}, \hat{h}_i = \frac{h_i}{\|h_i\|} \quad (14b) $$

$$ \hat{d} = \{\hat{d}_i\}_{i=1}^{4N}, \hat{d}_i = \frac{d_i}{\|d_i\|} \quad (14c) $$

By incorporating time importance, the offset difference vector is assigned weights to provide a more comprehensive characterization of the flexibility capacity. Suppose MGO assigns greater significance to the ith instant, and it is expected that buildings will enhance the flexibility capacity of this moment by increasing the $m_i$th value of $\hat{d}$, where $m_i = i, N + i, 2N + i$, and $3N + i$. To facilitate expression simplicity, we define $\hat{\omega}^{bd} = \{\hat{\omega}_i^{bd}\}_{i=1}^{4N} = [\omega^{bd}, \omega_i^{bd}, \omega_i^{bd}, \omega_i^{bd}]^T$ to match the matrix size. Subsequently, the weighted offset difference vector is characterized as follows:

$$ \delta^{od,t} = \delta^{od} \odot \hat{\omega}^{bd} = \{\delta_i^{od,t}\}_{i=1}^{4N}, \delta_i^{od,t} = \hat{\omega}_i^{bd} \delta_i^{od} \quad (15) $$

where $\odot$ means element-wise product.

We now conclude that the feasible region of $P$ is fully expanded if the minimization of $\delta^{od,t}$ is realized. Notice that values of $\rho^{rl}$ could be solved through a series of linear programming. Different norms could be used to quantify the vector $\delta^{od,t}$ as the total weighted offset value, expressed as: $\Delta^{od,t} = \|\delta^{od,t}\|_{\ell}$, where $\ell \in \{1, 2, \infty\}$.

<figure><img src="image_8.png" /><figcaption>Fig. 5. Example of the two-dimensional results for building flexibility approximation using different dissipation values $a_1$ and $a_2$ in the VB model.</figcaption></figure>

The crucial requirement for the individual flexibility approximation process is to ensure $P \subseteq P^{rl}$. We use Farka's lemma to achieve the polytope containment condition, referring to [48]. Under the assumption of the predefined dissipation rate, we seek to obtain the largest VB-based polytope inscribed within the actual building flexibility capacity. By minimizing the weighted offset difference vector, the optimization of individual flexibility approximation is formulated as Eq. (16):

$$ \begin{array}{ll} \min_{\hat{d}, G} & \Delta^{od,t} = \|\delta^{od,t}\|_{\ell} \\ \text{s.t.} & G(i, j) \ge 0 \\ & G\hat{H} = H^{rl}, G\hat{d} \le d^{rl} \\ & 0 \le \hat{d} \le \rho^{rl} \\ & \text{Eq. (6), Eq. (7), Eq. (10)–Eq. (15)} \end{array} \quad (16) $$

The result also depends on the selection of the dissipate rate $a$. Fig. 5 presents a two-dimensional example to show the solutions using different values of $a$. In the VB model, power constraints are characterized as an axis-aligned box. The sloping hyperplanes of the VB state constraints truncate the box in the diagonal direction and generate new facets. It is evident that the quality of the estimation solution improves when the facets of the inner polytope and outer polytope demonstrate a higher level of parallelism. By experimenting with different values of the dissipate rate, the building could explore more feasible solutions and ultimately achieve the optimal approximation model by selecting the scenario with the minimal weighted offset difference.

Moreover, the approximation procedure outlined in this subsection exhibits good adaptability and can be extended to map the flexibility of other types of DERs independent of the resource model order. The mathematical expressions for these reformulations are straightforward and not detailed in this work.

### 3. Microgrid flexibility aggregation

This section is dedicated to developing the process of total flexibility aggregation at the microgrid level, which serves to optimize decision-making and facilitate communication with the grid. The aggregation process is designed to dynamically and efficiently manage flexible energy in real-time, capturing critical information such as load dynamics, constraints, and measurements of available flexibility while maintaining computational tractability. However, the heterogeneity in resource characteristics can lead to substantial flexibility loss in the aggregate solution [49]. Many studies utilize averaging techniques as a compromise, but these can lead to overly conservative outcomes. As the diversity in the shapes of the loads' flexibility polytopes increases,
the size of the feasible region for aggregate flexibility generally diminishes [50]. Consequently, it is crucial that the aggregation model is carefully proposed to take advantage of the heterogeneous features and conservatively exploit the total flexibility capacity.

We continue to employ the virtual battery as the foundation. Suppose the dissipate rate is $a_A$, then the aggregation model is represented as:

$$P_A = \{u_A(k) | H_A u_A(k) \le d_A\} \quad (17a)$$

$$H_A = [-I_N, I_N, -\tilde{B}_A, \tilde{B}_A]^T \quad (17b)$$

$$d_A = [-e_A^u, \tilde{e}_A^u, -e_A^x, \tilde{e}_A^x]^T \quad (17c)$$

where $\tilde{B}_A$ is composed of the power of $a_A$ and possesses the same form as $\tilde{B}$ in Eq. (10b). In the computation process of $d_A$, the characteristic parameter $a_A$ is assumed to be defined in advance.

### 3.1. Polytope containment for Minkowski sum of individual flexibility

The exact capacity of a collection of flexible resources is their Minkowski sum $P_{MS}$:

$$P_{MS} = \bigoplus_j P_j = \left\{ u_{MS}(k) \middle| u_{MS}(k) = \sum_j u_j(k), u_j(k) \in P_j \right\} \quad (18)$$

where $\oplus$ means the Minkowski sum computation. Specifically, $P_{MS}$ is a convex polytope if all the sets $P_j$ are defined by convex polytopes [51]. To ensure the disaggregation feasibility of any profile within the aggregated flexibility polytope $P_A$, it is essential that it remains a subset of $P_{MS}$.

We observe that the expression of Eq. (18) is highly abstract. The Minkowski sum of H-form polytopes is a difficult task with #NP-hard complexity and is seldom implemented in practical applications [52]. As far as we know, there is no efficient and widely applicable method available to provide a general and tractable solution to this problem [53]. The inherent variability of the parameters from the local flexibility information adds difficulty to its computation. Therefore, utilizing Farka's lemma to ensure inner aggregation is not feasible due to the unavailability of an explicit expression for $P_{MS}$.

Nevertheless, we shall consider the Minkowski sum from another perspective. The polytope $P_{sc}^C$ is formed to contain all the individual flexibility information:

$$P_{sc}^C = \{u_{sc}^C(k) | H_{sc}^C u_{sc}^C(k) \le d_{sc}^C\} \quad (19)$$

where $H_{sc}^C$ and $d_{sc}^C$ store the collected approximation models:

$$u_{sc}^C(k) = [u_1(k), u_2(k), \dots, u_{N_s}(k)]^T \in \mathbb{R}^{(N \times N_s)} \quad (20a)$$

$$H_{sc}^C = \operatorname{diag}(H_1, H_2, \dots, H_{N_s}) \in \mathbb{R}^{(4N \times N_s) \times (N \times N_s)} \quad (20b)$$

$$d_{sc}^C = [d_1, d_2, \dots, d_{N_s}]^T \in \mathbb{R}^{(4N \times N_s)} \quad (20c)$$

Then the Minkowski sum $P_{MS}$ of resources' flexibility is transformed into the projection of the high-dimensional polytope $P_{sc}^C$ [54,55], achieved by the following affine transformation:

$$P_{MS} = T_{sc}^{pj} P_{sc}^{C} \quad (21)$$

where the projection matrix is

$$T_{sc}^{pj} = [I_N \quad I_N \quad \dots \quad I_N] \in \mathbb{R}^{N \times (N \times N_s)} \quad (22)$$

in which $I_N \in \mathbb{R}^N$ represents the identity matrix.

The AH-polytope containment theorem is introduced as an advanced extension of Farka's lemma, serving for more general polytope containment problems [56]:

**Theorem 1 (AH-polytope Containment).** Given two polytopes:

$$\begin{align*} P_i &= \{u_i \in \mathbb{R}^{N_i} \mid H_i u_i \le d_i, H_i \in \mathbb{R}^{m_i \times N_i}\} \\ P_e &= \{u_e \in \mathbb{R}^{N_e} \mid H_e u_e \le d_e, H_e \in \mathbb{R}^{m_e \times N_e}\} \end{align*}$$

and two affine transformations:

$$\begin{align*} P_i^{AH} &= p_i + T_i P_i \\ P_e^{AH} &= p_e + T_e P_e \end{align*}$$

$$\text{where } p_i \in \mathbb{R}^{N_i^{AH}}, p_e \in \mathbb{R}^{N_e^{AH}}, T_i \in \mathbb{R}^{N_i^{AH} \times N_i}, T_e \in \mathbb{R}^{N_e^{AH} \times N_e}.$$

The inner containment relationship $P_i^{AH} \subseteq P_e^{AH}$ is guaranteed if the following conditions are valid:

$$\exists \Phi \in \mathbb{R}^{N_e \times N_i}, \beta \in \mathbb{R}^{N_e}, G \in \mathbb{R}^{m_e \times m_i}, \text{ so that:}$$

$$\begin{align} G(i, j) &\ge 0 \tag{23} \\ T_i &= T_e \Phi, p_e - p_i = T_e \beta \nonumber \\ GH_i &= H_e \Phi, G d_i \le d_e + H_e \beta \nonumber \end{align}$$

Based on Theorem 1, we have the following conditions to guarantee the aggregation conservatism $P_A \subseteq P_{MS}$:

$$\begin{align} \exists \Phi \in \mathbb{R}^{(N \times N_s) \times N}, \beta \in \mathbb{R}^{(N \times N_s)}, \nonumber \\ \exists G \in \mathbb{R}^{(4N \times N_s) \times 4N}, G(i, j) \ge 0, \nonumber \end{align} \quad (24) \\ T_{sc}^{pj} \Phi = I_N, T_{sc}^{pj} \beta = 0, \\ GH_A = H_{sc}^C \Phi, G d_A \le d_{sc}^C + H_{sc}^C \beta \quad (24)$$

### 3.2. Microgrid flexibility aggregation optimization

During the aggregation phase, the optimization objective of minimizing the offset difference between the projections of the aggregation polytope and the overall flexibility capacity remains the preferred metric. Besides, it is vital to incorporate the potential demand of aggregator temporal preference $\omega^A = \{\omega_i^A\}_{i \in N}$ into the problem formulation. The values of $\omega^A$ should be strategically designed following the same principles for incentivizing buildings, whereby time instants with higher incentive values should be assigned larger time weights.

We normalize the row vectors of matrix $H_A$ into $\hat{H}_A = \{\hat{h}_{A,i}\}_{i=1}^{4N}$ for the sake of projection computation. Subsequently, the aggregation model Eq. (17) is reformulated as $P_A = \{u_A(k) | \hat{H}_A u_A(k) \le \hat{d}_A\}$, where $\hat{d}_A = \{\hat{d}_{A,i}\}_{i=1}^{4N}$, $\hat{d}_{A,i} = \frac{d_{A,i}}{\|\hat{h}_{A,i}\|}$. We use $\rho_{MS} = \{\rho_{MS,j}\}_{j=1}^{4N}$ and $\rho_A = \{\rho_{A,i}\}_{i=1}^{4N}$ standing for the projection distances of $P_{MS}$ and $P_A$ over all the $\hat{h}_{A,i}$. The vector $\rho_{MS}$ is deemed as the sum of individual projection vectors $\rho_j$:

$$\rho_{MS} = \sum_j \rho_j, j \in S \quad (25)$$

where $\rho_j$ is solved through linear programmings according to each polytope $P_j$. The corresponding offset difference vector is thereafter defined as $\delta_A^{od} = \{\delta_{A,i}^{od}\}_{i=1}^{4N}$, $\delta_{A,i}^{od} = |\rho_{MS,j} - \hat{d}_{A,i}|$.

To ensure compatibility with the optimization variable size, the parameter $\hat{\omega}^A = \{\omega_i^A\}_{i=1}^{4N}$ is defined, which holds a similar structure to $\hat{\omega}^{hd}$ described in Section 2. The weighted offset difference vector for the aggregation process is:

$$\delta_A^{od,j} = \delta_A^{od} \odot \hat{\omega}^A = \{\delta_{A,i}^{od,j}\}_{i=1}^{4N}, \delta_{A,i}^{od,j} = \hat{\omega}_j^A \delta_{A,i}^{od} \quad (26)$$

The optimization for microgrid total flexibility aggregation is formulated as follows:

$$\begin{align} \min_{\hat{d}_A, \Phi, \beta, G} \quad & \Delta_A^{od,j} = \|\delta_A^{od,j}\|_\ell \nonumber \\ \text{s.t.} \quad & 0 \le \hat{d}_A \le \rho_{MS} \tag{27} \\ & \text{Eq. (19)–Eq. (22), Eq. (24)–Eq. (26)} \nonumber \end{align}$$

The predefined dissipate rate $a_A$ of the aggregation model is essential in determining the aggregate quality. Furthermore, certain tools are available to estimate the volume for limited-size and low-dimensional structures. Considering the typical operational horizons in
microgrid operation, we propose designing $a_A$ by the accumulation of the weighted dissipation rates from local flexibility information. The weights are assigned based on the volume of their respective polytope, thus enabling a compromise in negotiating the local contributions regarding flexibility capacity. Accordingly, the dissipation rate for the aggregate model is given by Eq. (28), where Vol($P_j$) denotes the polytope volume estimation for the proposed flexibility capacity of the building indexed by $j$.

$$a_A = \sum_j \frac{\text{Vol}(P_j)a_j}{\sum_j \text{Vol}(P_j)}, \forall j \in S \qquad (28)$$

### 3.3. Limitations

The inclusion of a bilinear term in the optimization Eq. (27) introduces complexities to ensure the convergence of the solution. Additionally, the size of the variable $G$ in optimization Eq. (24) increases with the population number $N_s$ and the predictive horizon $N$. The problem can potentially exceed the maximum capacity that the solvers can reliably handle. To mitigate this issue, we have restricted the length of the control horizon to prevent overloading. We have also implemented a multi-start optimization process in the computation, which significantly increases the likelihood of obtaining successful solutions.

### 4. Case studies

In this section, we substantiate the effectiveness of the proposed approach. Our objective is to better preserve the energy flexibility at both the local level of individual flexibility approximations and the upper level of total flexibility aggregation, while ensuring consumer privacy. We also examine the effects of incorporating varying preferences when the aggregator adjusts its temporal flexible regulation capacity to bolster global performance. All computational experiments were conducted in MATLAB 2022a. The modeling languages YALMIP [57] and CVX [58] were employed, and optimizations were resolved using either the linear programming solver GUROBI [59] or the nonlinear programming solver SNOPT [60]. The volume of the metric polytope serves as a criterion for the quality evaluation and comparison of polytopes. The MPT3 toolbox [61] and VolEsti [45] were utilized for estimating the polytope volume. For all norms involved in the optimizations, $\ell = 1$ was selected.

#### 4.1. Analysis of individual flexibility approximation results

We first devised a series of scenarios and compared the solutions for individual flexibility approximations derived through our approach with those achieved by geometry-based algorithms utilizing hyperbox, ellipsoid, and homothet models. The scalability of this analysis is confined to a single-time horizon. The hyperbox-based approximation was obtained using the optimization formulation described in [21]. The ellipsoid-based approximation was tackled by addressing the problem of finding a maximum inscribed ellipsoid (MVE) in [23]. With respect to the homothetic algorithm in [24], the virtual battery model is predetermined as the base set in the form of a unit parallelotope, considering the aggregator's limited access to actual resource flexibility information. Our approach adopts the same dissipation rate value $a = 0.8$ as used in the homothet battery. We assume a uniform temporal preference as other methodologies do not account for this aspect.

For the convenience of visualization, Fig. 6 illustrates the individual approximated results of the given flexibility polytope over the horizon $N = 2$. The experimental parameters of the outer polytope are detailed in Table 1, representing the actual building flexibility information during this period. It should be noted that instances are represented by $i = 0, 1$ and $2$. The numbers in square brackets indicate the parameter values at the corresponding time. For example, $\theta_{ih}(1, 2)$ denotes the minimum building temperature bound at moment $i = 1$ and $2$. As can

<figure><img src="image_9.png" /><figcaption>Fig. 6. Individual approximation results for horizon $N = 2$. The dark gray region represents the given polytope depicting the actual flexibility capacity. The approximated polytopes are plotted in blue.</figcaption></figure>

<table><caption>Table 1
Information of the building's actual flexibility.</caption><thead><tr><td>Parameter</td><td>Value (s)</td><td>Unit</td><td>Explanation</td></tr></thead><tbody><tr><td>θ<sub>ih</sub>(0)</td><td>24</td><td>°C</td><td>Initial temperature</td></tr><tr><td>u<sub>ih</sub><sup>max</sup>(0, 1)</td><td>[100; 100]</td><td>kW</td><td>Baseline power</td></tr><tr><td>u<sub>ih</sub><sup>min</sup>(0, 1)</td><td>[25; 0]</td><td>kW</td><td>Minimum power</td></tr><tr><td>u<sub>ih</sub><sup>max</sup>(0, 1)</td><td>[150; 175]</td><td>kW</td><td>Maximum power</td></tr><tr><td>θ<sub>ih</sub>(1, 2)</td><td>[23.5; 23]</td><td>°C</td><td>Minimum temperature</td></tr><tr><td>θ<sub>ih</sub><sup>max</sup>(1, 2)</td><td>[24.75; 24.5]</td><td>°C</td><td>Maximum temperature</td></tr></tbody></table>

<table><caption>Table 2
Comparison of individual flexibility approximation performance for $N = 2$.</caption><thead><tr><td>Approach</td><td>Volume ratio</td><td>Offset difference</td></tr></thead><tbody><tr><td>Hyperbox</td><td>0.46</td><td>NA</td></tr><tr><td>Ellipsoid</td><td>0.86</td><td>NA</td></tr><tr><td>Homothet</td><td>0.59</td><td>113</td></tr><tr><td>Our approach</td><td>0.96</td><td>16</td></tr></tbody></table>

be seen, Fig. 6 reveals a notable advantage of the proposed approach in terms of the covered area, indicating a greater potential for adaptability and enhanced estimation quality in the individual approximation process.

More specifically, Table 2 presents the geometric properties depicted in Fig. 6. The item volume ratio is calculated as the ratio of the volume of the computed polytope to that of the outer polytope. This metric directly quantifies the extent to which the approximation result accurately represents the size and magnitude of the given flexibility capacity. We also evaluate the offset difference between the homothetic scenario and our solution, as both approaches utilize the VB model and share an equal number of hyperplanes. The higher volume ratio and lower offset difference value suggest that our approach achieves greater accuracy and better utilization of flexibility, because the proposed estimation process thoroughly considers the characteristics of the given information. In this instance, the ellipsoid approach yields a larger set among the existing algorithms. However, we cannot conclude that it is definitely the superior geometric method under all operating conditions. It should be noted that the homothetic approximation is heavily dependent on the initial base set, which can lead to distinctive outcomes. Lastly, the hyperbox approximation approach offers the least flexibility capacity as it cannot model the coupling of power and state variables.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Comparison of the correlation between the normalized volume ratio and the horizon by the approximation approaches.</figcaption></figure>

To further investigate the scalability of the proposed approach, we conducted experiments across different time horizons (equal to polytope dimension $N$). High-dimensional regions are not visually perceptible. It is important to notice that the volume ratio, as shown in Table 2, decreases rapidly with increasing $N$. Consequently, to provide a fairer illustration, we introduce the evaluation criterion of normalized volume ratio to reflect the conservatism level of the solutions [47]. Eq. (29) demonstrates its formulation.

$$r_{Vol} = \sqrt[\frac{N}{2}]{\frac{Vol(P)}{Vol(P^{rl})}} \quad (29)$$

Fig. 7 depicts the correlation between the normalized volume ratio and the horizon. It is evident that the available flexibility from all approaches decreases as the horizon increases, which infers a common trend of increased loss in long-horizon applications. The proposed approach consistently yields solutions with larger volumes, effectively exploring the potential of flexibility with minimal conservatism. We could also observe the favorable performance of the ellipsoid algorithm in estimating individual polytopes. Furthermore, the invariably inferior performance of the hyperbox method regarding volume suggests that this method may lead to a substantial flexibility loss. Although the normalized volume ratio and horizon exhibit a trend similar to a negative linear correlation in this scenario, it is crucial to note that such a characteristic is not universally applicable. In fact, the performance of the normalized volume ratio by the proposed approach is heavily influenced by the properties of the outer polytope and VB parameters. Therefore, the relationship with the horizon can vary considerably under different conditions.

Meanwhile, we use the running time of the simulation as the metric to evaluate the computational costs of these approaches. Fig. 8 shows the correlation between the running time and the control horizon. As $N$ increases, all approaches, except the hyperbox algorithm, require more time to generate solutions. Notably, approximations based on the hyperbox and homothet demonstrate superior performance in this aspect, with both approaches taking less than approximately 0.2 s for all tested horizons. The ellipsoid approximation requires a longer time, which is a reasonable outcome considering the algorithm complexity in MVE problem. Our approach consumes the most time, attributable to the elevated complexity of the optimizations. This disadvantage should be considered when applying the proposed approach to large-scale applications. In applications, although a longer horizon can often provide greater flexibility and convenience for the system's energy management, it is necessary to choose an appropriate horizon length because of the corresponding increasing operation burden.

In conclusion, selecting the suitable approximation approach requires a comprehensive trade-off between complexity considerations and the requirements for flexibility capacity. The analyses demonstrate that the proposed approach excels in preserving the flexibility potential,

<figure><img src="image_11.png" /><figcaption>Fig. 8. Comparison of the correlation between the running time and the horizon by the approximation approaches.</figcaption></figure>

<table><caption>Table 3
Price scheme of grid power.</caption><thead><tr><th>Time (h)</th><th>0</th><th>5</th><th>13</th><th>13.5</th><th>14</th><th>16</th><th>20</th><th>20.5</th></tr></thead><tbody><tr><td>P<sup>g,in</sup> (€/kWh)</td><td>0.14</td><td>0.16</td><td>0.18</td><td>0.18</td><td>0.18</td><td>0.16</td><td>0.4</td><td>0.6</td></tr><tr><td>P<sup>g,out</sup> (€/kWh)</td><td>0.2</td><td>0.2</td><td>0.4</td><td>0.6</td><td>0.8</td><td>0.2</td><td>0.2</td><td>0.2</td></tr></tbody></table>

albeit at the expense of larger computation burden. We also deem the homothetic approach as a satisfactory geometric approach, as it produces relatively good results while consuming little time. Despite the potential substantial flexibility loss with the hyperbox approach, it remains widely used in many practical programs due to its simplicity, efficiency, and compatibility with various DERs [21].

## 4.2. Analysis of flexibility aggregation results in microgrid operation

Next, we demonstrate the effectiveness of our approach to explore the flexibility potential at the aggregator level. Numerical simulations of the microgrid based on predictive control are conducted within a 24-hour daytime scenario. The sampling time is fifteen minutes. The experimental microgrid configuration comprises six heterogeneous buildings and a solar energy generator, in which four buildings are capable of proposing thermal flexibility, and the battery of one building is activated. The battery has a capacity of 100 kW/100 kWh, with an efficiency of 0.8 for charging and 0.83 for discharging. The solar power generation is sourced from a summer day of the typical meteorological year weather file in Rennes, whose power curve is depicted in Fig. 9. The predictive horizon for price announcements and building baseline prediction is two hours, while the horizon for flexibility aggregation is defined as one hour.

The microgrid aims to achieve cost-effective operations and minimize economic expenses using the economic model predictive control (EMPC) strategy. Table 3 details the grid price scheme throughout the day, consisting of two distinct sets of prices. Specifically, $P^{g,in}$ represents the cost of power purchased from the grid, and $P^{g,out}$ indicates the expenditure associated with feeding energy to the grid. The microgrid internal price plan encompasses the profiles for the nominal price $P^n$ and upward/downward flexibility incentives $P^{uw}/P^{dw}$. These prices are set at fixed amounts of $P^n = 0.2€/kWh$, $P^{uw} = 0.02€/kWh$, $P^{dw} = 0.22€/kWh$ to simplify the analysis. The microgrid's involvement in the power market is not applicable in this context, and the injection of excess power into the grid results in extra costs rather than revenue for the microgrid's operation.

Based on the insights from the previous analysis, we consider that the homothet-based method exhibits the most favorable overall performance among the geometry-based approaches. This observation may also serve as an explanation for its extensive utilization in numerous academic research. Consequently, we conduct a comparative analysis
<figure><img src="image_12.png" /><figcaption>Fig. 9. Solar power generation profile.</figcaption></figure>

in which the microgrid implements each aggregation process by means of the homothetic approach and our proposed approach. Specifically, both are formulated on the foundation of the VB model. Given the requirement to protect local consumer privacy, the aggregator lacks access to the actual resource flexibility information. We assume the homothetic base set is predefined as a unit parallelotope, as the typically used averaging techniques to determine this set are not feasible. The influence of the temporal factor is omitted for a fair comparison.

To directly present the flexibility capacity in the aggregation models, Fig. 10 shows the volume of the aggregation polytopes obtained in each control horizon. We present the per-unit values calculated from the maximum volume. The proposed approach exhibits superior performance due to its comparatively larger feasible region of the aggregation model on a global scale. Furthermore, the extent of microgrid flexibility utilization is shown in Fig. 11.

During midnight and early morning hours, when the solar generator is inactive and the grid power costs are comparatively low, MGO refrains from activating flexibility requirements during this period. This correlates with the significant volume of aggregation models of both methods depicted in Fig. 10. As we can see, as the solar generator begins to ramp up power generation significantly after 8:00, the system transitions from a power deficit to a surplus. Therefore, MGO must consume the surplus energy for balancing purposes. According to Table 3, the upward flexibility incentive remains lower than the cost of power injection into the grid, prompting MGO to seek upward flexibility capacity from local resources rather than giving power into the grid. Fig. 10 highlights the advantage of our methodology in terms of the augmented volume of the aggregation model, which offers a greater possibility for microgrid regulation. The broader flexibility utilization by our approach in Fig. 11 also corroborates this point. The extensive application of upward flexibility lasts throughout most of the daytime hours until approximately 17:00, in line with the major period of solar generation. Likewise, the nighttime period after 20:00 corresponds to the peak price hours for purchasing grid electricity supply. Consequently, downward flexibility is utilized until midnight to mitigate the energy demand charge from the grid. Our approach continues to offer larger aggregation models and enables more downward flexibility.

From another perspective, Fig. 12 shows the power exchange between the microgrid and the grid. We observe the enhanced power exchange reduction by our approach. Overall, the comparisons regarding the region size of aggregation models and microgrid flexibility requests provide compelling evidence for the performance improvement of the proposed approach in exploring a greater capacity for microgrid energy management. Moreover, Fig. 13 displays the temperature profile of a sample building using our method, where the regulation by the aggregator strictly fulfills the thermal comfort requirements. This example validates the feasibility of disaggregation

<figure><img src="image_13.png" /><figcaption>Fig. 10. Comparison of the aggregation polytope volume at each control horizon.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Comparison of microgrid flexibility utilization. Positive/negative values indicate upward/downward flexibility.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Comparison of the microgrid power exchange with the grid. Positive/negative values represent power purchase from/injection to the grid.</figcaption></figure>

decisions. It is noteworthy that some consumers, such as the sample building, raise their minimum temperature threshold during night hours, leading to higher energy demand and reduced capacity for downward flexibility. This helps explain why buildings do not contribute as much to the overall downward flexibility as the upward scenario.

Table 4 summarizes the main expenses and revenues of microgrid operation through both aggregation approaches. Key expenses include flexibility incentive payments and power exchange costs with the grid, which are shown as positive values. Conversely, buildings are required to cover the costs associated with their consumption, which contribute
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Temperature profile of a sample building. The thermal limits are represented by the black dashed lines.</figcaption></figure>

<table><caption>**Table 4**
Comparison of microgrid economic performance (€) in detailed items.</caption><thead><tr><td>Item</td><td>Homothetic approach</td><td>Our approach</td></tr></thead><tbody><tr><td>Building energy consumption</td><td>-1367.6</td><td>-1402.7</td></tr><tr><td>Flexibility incentive</td><td>203.8</td><td>419.5</td></tr><tr><td>Grid power</td><td>2875.3</td><td>2520.1</td></tr><tr><td>Sum</td><td>1711.5</td><td>1536.9</td></tr></tbody></table>

to the microgrid's revenue and are indicated by negative values. Our approach offers incentive payments that are more than double those the homothetic process provides. Simultaneously, the reduced power with the grid implies that the microgrid has successfully contracted reliance on the external supply. Global economic savings are achieved through the enhanced utilization of cheaper flexibility from the increased capacity in the aggregation process to balance the microgrid power. The illustrative example demonstrates a 10.2% decrease in total costs, indicating effective resource flexibility exploitation. Various factors, such as the microgrid configuration and building characteristics, may influence the extent of the improvement. Nevertheless, it can be inferred that the proposed approach holds the potential to enhance the economic efficiency of the microgrid. This could be further augmented if the microgrid actively engages in the power market and provides promoted DR services.

### 4.3. Analysis of temporal preference for flexibility

The discussions above assume uniform importance across all time instances to facilitate fair comparison with other aggregation methods. As explained in the introduction, the microgrid operator might adjust the flexibility requirement according to time under specific scenarios. We provide illustrative examples to showcase the impact of temporal preferences on the aggregation solutions. We take a single predictive period $N$ of length 5, in which time instants are indexed from 1 to 5 by the notation $i$. Meanwhile, a building is chosen as an exemplar to analyze the influence of time weights on the local flexibility capacity.

Plotting the geometry in high-dimensional space is not feasible. As a result, we will present the right-vector value $d$ of Eq. (12) to further elucidate the details of polytopic solutions based on the VB model. We establish two constraint ranges: $\mathcal{E}^u = [\bar{e}^u, \bar{e}^u]$ shows the power limits of battery; the range $\mathcal{E}^x = [\bar{e}^x, \bar{e}^x]$ indicates the battery energy state capacity. The boundaries of $\mathcal{E}^u$ and $\mathcal{E}^x$ will be presented by time index to illustrate more directly the geometric properties of the battery polytope at each spatial dimension. Besides, the measurement unit for the boundary values of both ranges is standardized as one to facilitate the analysis of mathematical interpretation.

We conduct the first examination during the midday period with abundant solar generation. Hence, MGO primarily requires upward

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Comparison of the temporal preference influence on constraint ranges of the flexibility approximation model from a building over a daytime horizon.</figcaption></figure>

flexibility potential to consume the surplus power. The following two conditions are considered:

* Case 1: No temporal preference with even weights.

* Case 2: Incremental time importance with weights series $\omega^A = [0.05, 0.05, 0.1, 0.3, 0.5]$.

Case 2 incorporates the varying importance of time, in which more flexibility is required in the later moments of the given period. Accordingly, the incentives within this interval would be designed based on $\omega^A$ and subsequently communicated to each flexible resource.

From the example building, Fig. 14 displays the constraint ranges $\mathcal{E}^u$ and $\mathcal{E}^x$ of the flexibility approximation model of the building. The results of Case 1 with even weights and Case 2 with varying weights are plotted in blue and red, respectively. Each bar in Fig. 14(a) represents the interval $[\bar{e}^u(i), \bar{e}^u(i)]$, and the bars in Fig. 14(b) mean the interval $[\bar{e}^x(i), \bar{e}^x(i)]$. In Case 1, when temporal preferences are not considered, there is a smooth variation of $\mathcal{E}^u$ and $\mathcal{E}^x$. The evolution of these limits profiles suggests the corresponding model's relatively stable flexibility capacity at different moments. The implementation of different time weights in Case 2 reveals a significant decrease in $\mathcal{E}^x$ during the initial instants. In contrast, $\mathcal{E}^x$ exhibits larger distances at instants 4 and 5, which increases the possibility for flexibility regulation capacity of the aggregator at these instants. In both cases, the slight variation of $\mathcal{E}^u$ can be attributed to the dominant influence of energy state constraints in shaping the model polytope. Based on Fig. 14, we can infer that this building has effectively adapted its flexibility contribution by implementing distinct weights in the optimization to match the temporal demand from the aggregator.

The constraint ranges $\mathcal{E}^u$ and $\mathcal{E}^x$ of the aggregate models obtained in both cases are displayed in Fig. 15. On the one hand, the constraint ranges of even weights present larger values in the first three time slots, suggesting that the aggregation battery in Case 1 preserves more flexibility capacity during these moments. On the other hand, we observe that $\mathcal{E}^u$ and $\mathcal{E}^x$ by changing the weights adhere to the 'save-first-release-later' pattern. This pattern implies an increased flexibility capacity of the aggregation battery in Case 2 during the last two time slots. This result is consistent with our expectation to enhance the flexibility potential for periods characterized by elevated importance.
<figure><img src="image_18.png" /><figcaption>Fig. 15. Comparison of the temporal preference influence on constraint ranges of the flexibility aggregation result over a daytime horizon.</figcaption></figure>

We could also perceive that the introduction of varying weights in the aggregation model could shift the distribution of flexibility capacity.

Similarly, we present another experiment conducted at night, where the microgrid requires downward flexibility to reduce energy purchase expenses from the grid. We compare two cases: Case 3 considers uniform time weights, and Case 4 assumes varying time weights whose values remain unchanged as in Case 2. The constraint ranges of the individual flexibility of a building and the aggregation model are depicted in Fig. 16 and Fig. 17, respectively. It can be inferred that the optimizations implementing varying time weights also demonstrate increased flexibility capacity during moments of higher temporal importance values, and this finding aligns with the previous daytime analysis.

To further elucidate the potential variation in flexibility utilization based on the above aggregation models, we envisage the benefits coefficients to be $\Psi = \{\psi_i\}_{i=1,2,6,10}$, where $\psi_i$ denotes the comprehensive advantage earned by the aggregator for each unit of flexibility employed during the $i$th time interval. The economic effect is no longer dominant as in the previous microgrid operation. The benefit profile $\Psi$ is in accordance with the temporal preference $\omega^A$ since higher importance instants are assigned augmented benefits. It is crucial to emphasize that the values of $\Psi$ are determined through comprehensive considerations of multiple factors to enhance microgrid performance. The total benefit of aggregator flexibility utilization within the given horizon is defined as

$$J_{mg}^{be} = \sum_{i} \psi_{i} |u_{A}(i)|, i \in \mathcal{N} \qquad (30)$$

Taking Eq. (30) as the objection function into optimization with the daytime aggregation model presented in Fig. 15, the optimal profiles of microgrid upward flexibility utilization are shown in Fig. 18. We observe no flexibility in both profiles during the first three moments due to the low associated benefits. Moreover, when employing the aggregation model from Case 2, the amount of flexibility utilization surpasses that observed in Case 1. The total benefit incorporating the time weights amounts to approximately 10.1k, nearly double the value of 5.6k obtained without considering this factor. Similarly, utilizing the same benefit function, the optimal profiles of microgrid downward

<figure><img src="image_19.png" /><figcaption>Fig. 16. Comparison of the temporal preference influence on constraint ranges of the flexibility approximation model from a building over a night horizon.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. Comparison of temporal preference influence on flexibility aggregation model over a night horizon.</figcaption></figure>

flexibility utilization based on the aggregation models from Fig. 17 are displayed in Fig. 19. In Case 4, the total benefit is 3.98k, marking a 19.5% improvement over Case 3. This observation further proves the effectiveness of our approach in prioritizing model capacity according to temporal factors.

In conclusion, employing varying levels of temporal importance alters the geometric configuration of both the individual flexibility approximation model and the total flexibility aggregation model. Notably, the VB polytope exhibits an increased size along the dimensions correlated with higher weight values. This analysis illustrates the capability
<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Comparison of optimal microgrid flexibility utilization based on aggregation models considering temporal preferences over a daytime horizon.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Comparison of optimal microgrid flexibility utilization based on aggregation models considering temporal preferences over a night horizon.</figcaption></figure>

of the proposed approach to adjust the aggregate flexibility capacity in response to time-related requests. In addition, the time weights are exclusively used for flexibility estimation, and they do not affect the actual energy flexibility characterization.

## 5. Conclusion

The employment of aggregation is vital for the power system to manage flexible assets effectively. In this study, we have presented an approach for aggregating the predictive energy flexibility of buildings through optimizations. The virtual battery is employed as the reference model to describe the flexibility information. The individual approximation effectively preserves the flexibility potential of available resources and addresses local privacy concerns. The aggregation method unlocks the total flexibility capacity by thoroughly considering the heterogeneity of individual resources. Feasible dispatch decisions are ensured while respecting consumer energy constraints. Numerical experiments and comparative analyses with existing algorithms demonstrate that the proposed approach is superior in exploiting both individual and aggregate flexibility. During a day-long operation test, the microgrid achieves an overall cost reduction of around 10%, while the extent of improvement may vary depending on specific system configurations. Moreover, the study confirms that the aggregator's temporal preferences influence the flexibility utilization pattern. The implementation of this proposed approach in the microgrid has the potential to yield substantial improvements in economic and other associated benefits. This framework can accommodate a broad spectrum of flexible resources characterized by high-order models.

Future research endeavors will explore methods to improve the computational efficiency of optimizations, such as the clustering algorithm of load management. Moreover, given the increased computational burden observed with larger horizons, devising a trade-off strategy to determine the appropriate problem size could be meaningful. Additionally, we recognize the value of studying the synergy effects of aggregation.

## CRediT authorship contribution statement

**Zhigang Zhang:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Hervé Guéguen:** Writing – review & editing, Supervision, Resources, Methodology, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgment

This work was supported by the China Scholarship Council (CSC) [grant No. 202008610214].

## References

[1] Saha S, Saleem M, Roy T. Impact of high penetration of renewable energy sources on grid frequency behaviour. Int J Electr Power Energy Syst 2023;145:108701.

[2] Shuai J, Zhao Y, Wang Y, Cheng J. Renewable energy product competitiveness: Evidence from the United States, China and India. Energy 2022;249:123614.

[3] Dong J, Xue G, Dong M, Xu X. Energy-saving power generation dispatching in China: Regulations, pilot projects and policy recommendations—A review. Renew Sustain Energy Rev 2015;43:1285–300.

[4] El Geneidy R, Howard B. Contracted energy flexibility characteristics of communities: Analysis of a control strategy for demand response. Appl Energy 2020;263:114600.

[5] Sarker E, Halder P, Seyedmahmoudian M, Jamei E, Horan B, Mekhilef S, et al. Progress on the demand side management in smart grid and optimization approaches. Int J Energy Res 2021;45(1):36–64.

[6] IEA. World energy outlook 2022. France: IEA Paris; 2022.

[7] Degefa MZ, Sperstad IB, Søle H. Comprehensive classifications and characterizations of power system flexibility resources. Electr Power Syst Res 2021;194:107022.

[8] Lu X, Li K, Xu H, Wang F, Zhou Z, Zhang Y. Fundamentals and business model for resource aggregator of demand response in electricity markets. Energy 2020;204:117885.

[9] Wang F, Li K, Liu C, Mi Z, Shafie-Khah M, Catalão JP. Synchronous pattern matching principle-based residential demand response baseline estimation: Mechanism analysis and approach description. IEEE Trans Smart Grid 2018;9(6):6972–85.

[10] Yang S, Lao K-W, Chen Y, Hui H. Resilient distributed control against false data injection attacks for demand response. IEEE Trans Power Syst 2023.

[11] Gutiérrez-Alcaraz G, Tovar-Hernández J, Lu C-N. Effects of demand response programs on distribution system operation. Int J Electr Power Energy Syst 2016;74:230–7.

[12] Kundu S, Kalsi K, Backhaus S. Approximating flexibility in distributed energy resources: A geometric approach. In: 2018 power systems computation conference. IEEE; 2018, p. 1–7.

[13] Córdova S, Lorca Á, Olivares DE, et al. Aggregate modeling of thermostatically controlled loads for microgrid energy management systems. IEEE Trans Smart Grid 2023.

[14] Muhtadi A, Pandit D, Nguyen N, Mitra J. Distributed energy resources based microgrid: Review of architecture, control, and reliability. IEEE Trans Ind Appl 2021;57(3):2223–35.

[15] Mattingley J, Wang Y, Boyd S. Code generation for receding horizon control. In: 2010 IEEE international symposium on computer-aided control system design. IEEE; 2010, p. 985–92.
[16] Tian G, Sun QZ, Wang W. Real-time flexibility quantification of a building HVAC system for peak demand reduction. IEEE Trans Power Syst 2021;37(5):3862-74.

[17] Amato V, Knudsen M, Petersen S. Dual-zone economic model predictive control of residential space heating for demand response using a single heat meter. Energy Build 2023;281:112759.

[18] Sajjad IA, Chicco G, Napoli R. Definitions of demand flexibility for aggregate residential loads. IEEE Trans Smart Grid 2016;7(6):2633-43.

[19] Qdr Q. Benefits of demand response in electricity markets and recommendations for achieving them. 2006, US Dept. Energy, Washington, DC, USA, Tech. Rep 2006.

[20] Ziegler GM. Lectures on polytopes, vol. 152, Springer Science & Business Media; 2012.

[21] Chen X, Dall'Anese E, Zhao C, Li N. Aggregate power flexibility in unbalanced distribution systems. IEEE Trans Smart Grid 2019;11(1):258-69.

[22] Polymeneas E, Meliopoulos S. Aggregate modeling of distribution systems for multi-period OPF. In: 2016 power systems computation conference. IEEE; 2016, p. 1-8.

[23] Cui B, Zamzam A, Bernstein A. Network-cognizant time-coupled aggregate flexibility of distribution systems under uncertainties. In: 2021 American control conference. IEEE; 2021, p. 4178-83.

[24] Zhao L, Zhang W, Hao H, Kalsi K. A geometric approach to aggregate flexibility modeling of thermostatically controlled loads. IEEE Trans Power Syst 2017;32(6):4721-31.

[25] Müller FL, Sundström O, Szabó J, Lygeros J. Aggregation of energetic flexibility using zonotopes. In: 2015 54th IEEE conference on decision and control. IEEE; 2015, p. 6564-9.

[26] Müller FL, Szabó J, Sundström O, Lygeros J. Aggregation and disaggregation of energetic flexibility from distributed energy resources. IEEE Trans Smart Grid 2017;10(2):1205-14.

[27] Hao H, Sanandaji BM, Poolla K, Vincent TL. Aggregate flexibility of thermostatically controlled loads. IEEE Trans Power Syst 2014;30(1):189-98.

[28] Hughes JT, Domínguez-García AD, Poolla K. Identification of virtual battery models for flexible loads. IEEE Trans Power Syst 2016;31(6):4660-9.

[29] Gonzalez DM, Hachenberger J, Hinker J, Rewald F, Häger U, Rehtanz C, et al. Determination of the time-dependent flexibility of active distribution networks to control their TSO-dso interconnection power flow. In: 2018 power systems computation conference. IEEE; 2018, p. 1-8.

[30] Bhattarai BP, Myers KS, Bak-Jensen B, de Cerio Mendaza ID, Turk RJ, Gentle JP. Optimum aggregation of geographically distributed flexible resources in strategic smart-grid/microgrid locations. Int J Electr Power Energy Syst 2017;92:193-201.

[31] Wang K, Yin R, Yao L, Yao J, Yong T, Deforest N. A two-layer framework for quantifying demand response flexibility at bulk supply points. IEEE Trans Smart Grid 2016;9(4):3616-27.

[32] Meng W, Song D, Huang L, Chen X, Yang J, Dong M, et al. Distributed energy management of electric vehicle charging stations based on hierarchical pricing mechanism and aggregate feasible regions. Energy 2024;291:130332.

[33] Pan Z, Guo Q, Sun H. Feasible region method based integrated heat and electricity dispatch considering building thermal inertia. Appl Energy 2017;192:395-407.

[34] Liu X, Li Y, Wang L, Tang J, Qiu H, Berizzi A, et al. Dynamic aggregation strategy for a virtual power plant to improve flexible regulation ability. Energy 2024;131261.

[35] Wang M, Mu Y, Shi Q, Jia H, Li F. Electric vehicle aggregator modeling and control for frequency regulation considering progressive state recovery. IEEE Trans Smart Grid 2020;11(5):4176-89.

[36] Zhang M, Xu Y, Shi X, Guo Q. A fast polytope-based approach for aggregating large-scale electric vehicles in the joint market under uncertainty. IEEE Trans Smart Grid 2023.

[37] Taha FA, Vincent T, Bitar E. An efficient method for quantifying the aggregate flexibility of plug-in electric vehicle populations. 2022, arXiv preprint arXiv: 2207.07067.

[38] Wen Y, Hu Z, You S, Duan X. Aggregate feasible region of DERs: Exact formulation and approximate models. IEEE Trans Smart Grid 2022;13(6):4405-23.

[39] Cui B, Fan C, Munk J, Mao N, Xiao F, Dong J, Kuruganti T. A hybrid building thermal modeling approach for predicting temperatures in typical, detached, two-story houses. Appl Energy 2019;236:101-16.

[40] Boodi A, Beddiar K, Amirat Y, Benbouzid M. Building thermal-network models: a comparative analysis, recommendations, and perspectives. Energies 2022;15(4):1328.

[41] Zhang Z, Gueguen H. Bilevel optimization based on building dynamic flexibility capacity in microgrid. IFAC-Papersonline 2022;55(9):274-9.

[42] Chen Y, Guéguen H, Wang X. Managing flexibility of power consumption of smart buildings on microgrid. IFAC-PapersOnLine 2019;52(4):383-8.

[43] Henk M, Richter-Gebert J, Ziegler GM. Basic properties of convex polytopes. In: Handbook of discrete and computational geometry. Chapman and Hall/CRC; 2017, p. 383-413.

[44] Wang J, Huang S, Wu D, Lu N. Operating a commercial building HVAC load as a virtual battery through airflow control. IEEE Trans Sustain Energy 2020;12(1):158-68.

[45] Emiris ΙΖ, Fisikopoulos V. Practical polytope volume approximation. ACM Trans Math Softw 2018;44(4):1-21.

[46] Emiris ΙΖ, Fisikopoulos V. Efficient random-walk methods for approximating polytope volume. In: Proceedings of the thirtieth annual symposium on computational geometry. 2014, p. 318-27.

[47] Jin B, Maasoumy M, Nuzzo P, Sangiovanni-Vincentelli A. Online computation of polytopic flexibility models for demand shifting applications. In: 2017 13th IEEE conference on automation science and engineering. IEEE; 2017, p. 900-5.

[48] Eaves BC, Freund RM. Optimal scaling of balls and polyhedra. Math Program 1982;23(1):138-47.

[49] Huang X, Dong P, He X, Liu M, Deng W. Research on aggregation flexibility method of AC/DC distribution network considering flexibility balance and its application. IEEE Syst J 2023.

[50] Nazir MS, Hiskens IA, Bernstein A, Dall'Anese E. Inner approximation of minkowski sums: A union-based approach and applications to aggregated energy resources. In: 2018 IEEE conference on decision and control. IEEE; 2018, p. 5708-15.

[51] Tiwary HR. On the hardness of computing intersection, union and minkowski sum of polytopes. Discrete Comput Geom 2008;40:469-79.

[52] Song Y, Xia M, Chen Q. The robust synchronization control scheme for flexible resources considering the stochastic and delay response process. Appl Energy 2023;343:121224.

[53] Xu W-R, Bebiano N, Chen G-L. A reduction algorithm for reconstructing periodic Jacobi matrices in Minkowski spaces. Appl Math Comput 2022;419:126853.

[54] Yi Z, Xu Y, Gu W, Yang L, Sun H. Aggregate operation model for numerous small-capacity distributed energy resources considering uncertainty. IEEE Trans Smart Grid 2021;12(5):4208-24.

[55] Wen Y, Hu Z, Liu L. Aggregate temporally coupled power flexibility of DERs considering distribution system security constraints. IEEE Trans Power Syst 2022.

[56] Sadraddini S, Tedrake R. Linear encodings for polytope containment problems. In: 2019 IEEE 58th conference on decision and control. IEEE; 2019, p. 4367-72.

[57] Löfberg J. YALMIP : A toolbox for modeling and optimization in MATLAB. In: In proceedings of the CACSD conference. 2004.

[58] Grant M, Boyd S. CVX: Matlab software for disciplined convex programming, version 2.1. 2014, http://cvxr.com/cvx.

[59] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual. 2023, URL https://www.gurobi.com.

[60] Gill PE, Murray W, Saunders MA. SNOPT: An SQP algorithm for large-scale constrained optimization. SIAM Rev 2005;47(1):99-131.

[61] Herceg M, Kvasnica M, Jones C, Morari M. Multi-Parametric Toolbox 3.0. In: Proc. of the European control conference. 2013, p. 502-10, http://control.ee.ethz.ch/~mpt.