Electrothermal coupling modeling of air-conditioning-building systems and aggregated direct load control dispatch tracking for office parks

1 Introduction

Against the backdrop of accelerating climate change and the global energy transition, China has committed to achieving peak carbon emissions before 2030 and carbon neutrality by 2060 (Hou et al., 2022). Buildings represent a large share of end-use energy consumption, with air-conditioning systems alone accounting for roughly 40%–60% of total building energy demand (Wang Q. et al., 2024). At the same time, the large-scale integration of renewable energy has increased peak-shaving and frequency regulation pressures on power systems, making it necessary to exploit flexible resources on the demand side (Huang et al., 2021). In this context, building air-conditioning loads—which offer considerable aggregate capacity, wide geographical distribution, and fast response capability—have become attractive candidates for grid-interactive applications (Wang H. et al., 2025).

Air-conditioning load modeling and optimal control have received sustained research attention. Modeling approaches fall broadly into three categories: physics-based mechanistic models, data-driven methods, and hybrid formulations (Yasaman et al., 2023). Physics-based models derive room thermal balance equations from the first law of thermodynamics, predicting indoor temperature by characterizing heat transfer through building envelopes (exterior walls, roofs, and windows), internal heat gains, and HVAC processes (Liu et al., 2016). Early work introduced simplified equivalent thermal parameter models that reduce complex building thermal dynamics to resistance-capacitance networks, achieving substantial model order reduction at acceptable accuracy (Goyal and Barooah, 2012). This framework was later extended to two-resistance, two-capacitance configurations capable of capturing dynamic thermal response (Ejenakevwe et al., 2025). While such models offer clear physical interpretability, parameter acquisition remains difficult and computational demands can be high. A further limitation is that these formulations focus primarily on thermodynamic processes, with relatively little attention paid to the electrical characteristics of air-conditioning equipment.

The spread of smart grid technologies and demand response programs has brought growing interest to the electrothermal coupling behavior of air-conditioning loads. Power-temperature coupling models for variable-frequency units have been developed to examine how compressor frequency modulation affects cooling performance (Guo et al., 2024). Dynamic COP models that account for wall temperature effects provide a basis for load aggregation (Anakottapary et al., 2026). From an aggregate control perspective, large populations of air-conditioning units have been treated as virtual energy storage systems, leading to state-queue-based aggregation models (Zhou et al., 2019). Markov chain descriptions of on/off state transitions across unit ensembles have supported probabilistic load forecasting (Austin et al., 2023). These contributions, though valuable, share several persistent limitations. Most existing models address isolated aspects of the problem rather than systematically tracing the full coupling chain from electrical power input through heating/cooling output to indoor temperature. Unit-level and aggregate-level model development also tends to proceed independently, with insufficient regard for aggregability, so accuracy degrades during the aggregation step. Model-driven optimal control methods suited to practical office park-scale implementation remain underdeveloped.

On the control side, hierarchical architectures have attracted attention for their scalability. Model predictive control-based methods for building HVAC systems have shown effectiveness in demand response settings (Wang X. et al., 2025). Two-tier optimization frameworks built around aggregator-user interactions achieve resource allocation through iterative negotiation (Zhu et al., 2026). Reinforcement learning applied to HVAC control has pointed to the potential of data-driven policy optimization (Liu et al., 2024). A shared weakness of these hierarchical methods is their reliance on simplified load representations that do not fully use the physical insight available from electrothermal coupling models, which constrains optimization performance under complex operating conditions.

The architectural scaling of HVAC loads for demand response has been extensively advanced through hierarchical coordination. For instance, Reference (Liu et al., 2026) developed a two-tier hierarchical optimal control framework that aggregates multi-zone building thermal dynamics via state-space order reduction to track grid dispatch targets. While (Wang T. H. et al., 2024) successfully mitigates computational complexity at the tracking layer, it relies on a linear thermal baseline that does not fully account for variable thermodynamic efficiencies (dynamic COP) or envelope thermal storage effects during rapid transient periods. The present work complements this lineage by embedding these physical realities directly into an aggregation-oriented design.

To precisely delineate the contributions of the present work relative to the closest prior studies, Table 1 provides a structured comparison across six dimensions: thermal model order, COP modeling, aggregability design, aggregation method, control architecture, and end-to-end chain closure from electrical input to temperature output.

TABLE 1. Comparison with representative prior works.

First, previous RC-based thermal models (Goyal and Barooah, 2012; Ejenakevwe et al., 2025) were developed for the simulation or calibration of individual buildings, and their design did not account for subsequent aggregation requirements. When such models are aggregated, their accuracy decreases. In this study, we treated aggregability as a primary design requirement from the outset. The 2R2C model structure was selected because its parameters can be linearly superimposed in weighted combinations and can be directly mapped to the eigenvectors used in the clustering stage. None of the compared references embodied this “aggregation-oriented design” philosophy, which is the key structural reason why the loss of aggregation accuracy in the case studies was limited to less than 8%. Second, References (Guo et al., 2024) and (Wang X. et al., 2025) discuss MPC control of HVAC loads, but employ a first-order thermal model that ignores the thermal storage of the building envelope, which is precisely the primary cause of the time lag between power regulation and temperature response. Although References (Goyal and Barooah, 2012) and (Ejenakevwe et al., 2025) capture the dynamic characteristics of the building envelope using a 2R2C model, they do not incorporate a dynamic COP—their heat input is treated as an exogenous variable rather than a function of operating conditions. References (Anakottapary et al., 2026)– (Austin et al., 2023) model COP or total load behavior but fail to construct a complete chain from electrical power input through dynamic heat transfer to indoor temperature evolution. This study systematically closes this chain within a single state-space model—electric power input → dynamic COP correction → heat output → 2R2C building envelope and aerodynamics → indoor temperature—and the model is applicable to parameter identification, aggregation, and real-time optimization. Third, although References (Zhou et al., 2019) and (Zhu et al., 2026) performed load aggregation, they employed state queues or rule-based methods and failed to utilize the identified physical model parameters. This study employs a two-stage clustering method—first performing coarse-grained attribute classification, followed by applying K-means clustering based on silhouette coefficient selection in the parameter space—to fully utilize the complete parameter identification results from the unit-level model, thereby forming thermodynamically homogeneous clustering units. Subsequently, capacity-weighted parameter equivalence ensures that buildings with larger installed capacities, which dominate the aggregated power response, are more accurately represented in the reduced-order office park model.

The post-2023 literature has further consolidated the integration of HVAC loads into virtual power plant (VPP) frameworks. Recent work has explored distributionally-robust optimization for hierarchical VPP-HVAC dispatch (Liu et al., 2024), incentive-elastic aggregation of air-conditioning loads (Yang et al., 2026), and cross-zone meta-modeling for single-building MPC at scale (Raman et al., 2021). These contributions, together with the broader cluster-level scheduling literature (Ge et al., 2024), converge on a shared recognition: the bottleneck is no longer whether HVAC loads can be aggregated, but how to do so with sufficient model fidelity and aggregability for the resulting dispatch decisions to be physically realizable. The present work positions itself within this trajectory by treating aggregability as a structural design criterion rather than as a post hoc consequence of model reduction.

In light of the research gaps identified above and summarized in Table 1, this paper proposes an integrated electro-thermal coupling modeling and hierarchical cooperative optimization framework for aggregated scheduling of building HVAC loads. The innovation of this study does not lie in any single component—RC thermodynamic modeling, COP correction functions, cluster analysis, and hierarchical control have all been discussed in previous literature—but rather in three specific comprehensive contributions: (i) a 2R2C electro-thermal coupling model designed for aggregation from the outset, which closes the complete power-temperature chain through a dynamic, empirically calibrated COP; (ii) A hierarchical clustering and capacity-weighted equivalence strategy that utilizes identified model parameters to form thermodynamically homogeneous aggregation units within an acceptable range of accuracy loss; and (iii) A two-tier hierarchical control architecture that propagates the physical fidelity of unit-level models upward to office park-level multi-objective optimization, rather than relying on simplified surrogate models at the higher level. A simulation case study of an office park comprising 15 buildings validates the effectiveness of this methodology. The present work is complementary to recent VPP-oriented contributions: where (Liu et al., 2023) addresses uncertainty through distributionally-robust optimization at the dispatch layer, this work absorbs uncertainty through the structural feedback path of the lower-level MPC, and supplies a physically-grounded aggregation primitive that such robust frameworks could in turn build upon.

2 Methods

2.1 Overall framework

A typical office park consists of multiple buildings, each served by independent or centralized HVAC systems. During heating or cooling seasons, these systems draw electrical energy to produce thermal output, which is delivered to indoor spaces through air distribution networks. Indoor thermal conditions are jointly governed by heat transfer through the building envelope, internal heat gains, and fresh air infiltration. This cascaded process—from electrical energy consumption through thermal generation to indoor temperature variation—forms the physical basis of air-conditioning-building electrothermal coupling.

From a grid dispatch perspective, the aggregate air-conditioning load of a office park can be treated as a flexible demand-side resource. The dispatch center issues power regulation commands to the office park according to system operating requirements, and the office park scheduling layer must then allocate power across individual buildings so that the aggregate office park load tracks the dispatch target while indoor temperatures remain within comfort bounds. Two competing objectives are involved: minimizing the deviation between actual office park power and the dispatch command, and keeping indoor temperatures within prescribed ranges. Building thermal inertia introduces a measurable lag between power adjustment and temperature response, a characteristic that provides temporal flexibility for scheduling but at the same time makes control strategy design more demanding.

To address these challenges, a three-tier progressive technical framework is constructed, consisting of a unit modeling layer, an aggregation extension layer, and an optimization control layer, as illustrated in Figure 1.

The unit modeling layer establishes electrothermal coupling energy efficiency models for buildings of different types, characterizing the dynamic relationship between electrical power input and indoor temperature output. These unit-level models serve as the basis on which subsequent aggregation is built. The aggregation extension layer then performs hierarchical clustering based on building characteristics, applying parameter equivalence to groups of similar buildings to produce representative aggregated models. This step reduces state variable dimensionality while keeping accuracy loss within acceptable bounds. The optimization control layer adopts a hierarchical architecture in which office park-level dispatch and building-level control operate at separate tiers. At the upper level, multi-objective optimization problems are solved using the aggregated models to determine power allocation across building categories. At the lower level, predictive control is implemented on the basis of unit electrothermal coupling models to execute specific regulation actions. Information exchange between the two levels forms a closed-loop feedback mechanism that links dispatch decisions to physical responses.

2.2 Unit-level electrothermal coupling energy efficiency model

The electrothermal coupling characteristics of air-conditioning-building systems manifest in the cascaded transfer relationships among electrical energy consumption, thermal output generation, and temperature variation. This section first establishes room thermal balance equations, then introduces the dynamic COP model for air-conditioning systems, and finally presents a unified expression for the cascaded transfer function.

2.2.1 Room thermal balance equation

For a representative office room, indoor air can be approximated as a lumped-parameter system whose thermal balance relationship is expressed by Equation 1:where C in denotes indoor air thermal capacitance (kJ/°C); T in t represents indoor temperature at time t (°C); Q ac t is the heat delivered to the room by the air-conditioning system (kW), with positive values indicating heating and negative values cooling; Q wall t represents heat transfer through the building envelope (kW); Q int t) denotes internal heat source dissipation (kW); and Q vent t accounts for heat exchange due to fresh air infiltration (kW).

C in d T in t d t = Q ac t + Q wall t + Q int t + Q vent t (1)

The heat transfer process through building envelopes employs an equivalent thermal resistance-capacitance network model. This representation treats the envelope as a single node with lumped thermal capacitance, exchanging heat with indoor and outdoor environments through internal and external thermal resistances. The thermal balance equation for the envelope node is given by Equation 2:where T out t denotes outdoor ambient temperature (°C); T wall t represents the equivalent envelope temperature (°C); R out and R in are the external and internal equivalent thermal resistances (°C/kW), respectively; and C wall is the equivalent envelope thermal capacitance (kJ/°C). The first term on the right-hand side of Equation 2 represents heat transfer from the outdoor environment to the envelope through the external thermal resistance, while the second term represents heat transfer from the envelope to the indoor space through the internal thermal resistance. The difference between these two terms governs the rate of envelope temperature change.

C wall d T wall t d t = T out t − T wall t R out − T wall t − T in t R in (2)

Based on the thermal balance relationship in Equation 2, the heat flux transmitted through the envelope to the indoor space can be expressed as Equation 3:

Q wall t = T wall t − T in t R in (3)

This equation describes the convective heat transfer process between the envelope interior surface and indoor air. When envelope temperature exceeds indoor temperature, heat flows from the structure toward the interior (Q wall > 0); the reverse applies otherwise.

Heat exchange induced by fresh air infiltration is expressed by Equation 4:where ρ denotes air density (kg/m³); c ρ is the specific heat capacity of air at constant pressure (kJ/(kg·°C)); and V vent represents the fresh air ventilation rate (m³/s).

Q vent t = ρ c ρ V vent T out t − T in t (4)

2.2.2 Dynamic coefficient of performance model

Air-conditioning systems convert electrical energy into heating or cooling capacity, and this conversion efficiency is expressed through the COP. Thermodynamically, such a system operates as a reverse Carnot cycle, whose theoretical COP is bounded by the temperatures of the heat source and heat sink (Feidt, 2020). For heating mode operation, the theoretical Carnot COP is given by Equation 5:where T in and T out are indoor and outdoor absolute temperatures (K). Practical air-conditioning systems, however, achieve actual dynamic COP values substantially below theoretical limits due to heat transfer temperature differences, compressor efficiency losses, and piping resistances. Accordingly, the actual COP is expressed as the product of rated COP and multiple correction factors, as shown in Equation 6:where COP nom represents the energy efficiency ratio under rated conditions; f T Δ T is the temperature difference correction function characterizing the influence of indoor-outdoor temperature differential Δ T on COP; and f LR λ LR is the load ratio correction function reflecting the effect of load ratio λ LR on COP. Each factor's influence on COP is assumed to be relatively independent and separable.

COP Carnot = T in T in − T out (5)

COP t = COP nom · f T Δ T · f LR λ LR (6)

The temperature difference correction function f T Δ T can be derived from Carnot efficiency considerations (Hu and Xiao, 2018). Taylor expansion of Equation 5 about the rated temperature difference Δ T nom, retaining first-order terms, yields Equation 7:

COP Carnot Δ T ≈ COP Carnot Δ T nom + ∂ CO P Carnot ∂ Δ T | Δ T nom Δ T − Δ T nom (7)

Differentiation of Equation 5 and substitution produces Equation 8:

∂ COP Carnot ∂ Δ T = − T in T in − T out 2 ≈ − COP Carnot 2 T in (8)

Substituting Equation 8 into Equation 7 and normalizing yields the linearized temperature difference correction function given by Equation 9:where k T is the temperature difference influence coefficient (1/°C), representing the relative COP variation per unit temperature difference change. In practical applications, k T can be determined through regression on equipment performance test data. The first-order Taylor expansion in Equation 9 is formally accurate only within an infinitesimally small neighborhood of point Δ T nom; as the operating point deviates from the rated conditions, its accuracy gradually decreases. However, the fitting of point k T, combined with the piecewise linearization process described in Supplementary Appendix A, mitigates linearization errors, ensuring that the linearization error in f Δ T remains within an acceptable range under the winter operating conditions of the case study.

f T Δ T = 1 − k T Δ T − Δ T nom (9)

The load ratio correction function f LR λ LR employs a quadratic polynomial fit as given by Equation 10:

f LR λ LR = a 0 + a 1 λ LR + a 2 λ LR 2 λ LR = P ac t / P rated (10)

Where a 0, a 1, and a 2 are fitting coefficients satisfying a 0 + a 1 + a 2 = 1; P ac t is the air-conditioning system input electrical power (kW); and P rated is the rated air-conditioning power.

Based on the dynamic COP model, the relationship between air-conditioning thermal output and input electrical power is expressed by Equation 11:

Q ac t = P ac t · COP t (11)

Substituting Equations 6, 9, 10 into Equation 11 yields Equation 12:which comprehensively describes the coupling relationships among electrical power, indoor-outdoor temperature differential, load ratio, and thermal output, establishing the foundation for subsequent state-space model development.

Q ac t = P ac t · COP nom · 1 − k T Δ T − Δ T nom · a 0 + a 1 λ LR + a 2 λ LR 2 (12)

2.2.3 Cascaded transfer function

Integrating the thermal balance equations with the dynamic COP model, this section derives a state-space model describing the complete coupling relationship from electrical power through heating capacity to indoor temperature.

First, substituting Equation 11 into the room thermal balance Equation 1, combined with Equations 3, 4, yields the dynamic equation for indoor temperature given by Equation 13.

C in d T in d t = P ac · COP t + T wall − T in R in + ρ c p V vent T out − T in + Q int (13)

Rearranging Equation 13 by moving state-variable-related terms to the right-hand side produces Equation 14.

d T in d t = − 1 C in 1 R in + ρ c p V vent T in + 1 C in R in T wall + COP t C in P ac + ρ c p V vent C in T out + 1 C in Q int (14)

The envelope thermal balance Equation 2 undergoes similar rearrangement to yield Equation 15.

d T wall d t = 1 C wall R in T in − 1 C wall 1 R in + 1 R out T wall + 1 C wall R out T out (15)

Defining the state vector as x t = T in t , T wall t T, input vector as u t = P ac t , T out t , Q int t T, and output y t = T in t, Equations 14, 15 can be expressed in standard state-space form as Equations 16, 17.

x ˙ t = A x t + B t u t (16)

y t = C x t (17)

Based on the coefficients in Equations 14, 15, the system matrix A is given by Equation 18:where a 11 represents the negative feedback effect of indoor temperature on its own variation, a 21 represents the positive influence of envelope temperature on indoor temperature, a_21 represents the positive influence of indoor temperature on envelope temperature, and a 22 characterizes the thermal dissipation properties of the envelope. The eigenvalues of matrix A determine the system time constants and are typically negative real numbers, ensuring asymptotic stability.

A = a 11 a 12 a 21 a 22 = − 1 C in 1 R in + ρ c p V vent 1 C in R in 1 C wall R in − 1 C wall 1 R in + 1 R out (18)

The input matrix B t is given by Equation 19:where b 11 t = COP t / C in represents the influence coefficient of electrical power on indoor temperature; since COP t depends on indoor-outdoor temperature differential and load ratio, this element exhibits time-varying characteristics. Term b 12 represents the influence of outdoor temperature on indoor conditions through fresh air infiltration; b 13 represents the direct influence of internal heat sources; and b 22 represents the influence of outdoor temperature on the envelope through its external surface.

B t = b 11 t b 12 b 13 b 21 b 22 b 23 = COP t C in ρ c p V vent C in 1 C in 0 1 C wall R out 0 (19)

The output matrix C is given by Equation 20:

C = 1 0 (20)

Owing to the time-varying nature of COP t, input matrix B t varies dynamically with operating conditions, rendering the system nonlinear. For subsequent optimization problem solutions, piecewise linearization treatment can be adopted, with specific procedures detailed in Supplementary Appendix A.

The model structure is illustrated in Figure 2, which shows the cascaded relationships among electrical input, thermodynamic processes, and temperature output. Electrical power is first converted to heating capacity through the dynamic COP, and this thermal output then acts on the room heat balance together with other heat sources to determine how indoor temperature evolves over time. The thermal storage effect of the building envelope introduces a temporal buffer into this process, with the equivalent envelope temperature serving as a state variable that carries the accumulated influence of prior thermal inputs.

2.2.4 Model parameter identification

This paper employs a data-driven parameter identification method that utilizes historical data recorded by building energy management systems for model calibration.

Required historical data includes: indoor temperature T in k sampled at 5–15 min intervals; outdoor temperature T out k; air-conditioning electrical power P ac k; air-conditioning operating status (on/off, setpoint temperature, operating mode, etc.); and internal load data (lighting power, equipment power, occupancy counts, etc.) for estimating internal heat sources Q int k. The data was collected from three representative buildings, each located in a different aggregation unit, over a continuous 60-day period during the winter heating season. Each building contributed approximately 17,280 records at 5-min intervals (60 days × 24 h × 12 readings per hour), resulting in a total of approximately 51,840 raw data records from the three buildings.

The continuous state-space model of Equations 16, 17 is discretized using the forward Euler method to obtain the discrete-time model given by Equations 21, 22: where the discretized matrices are A d = I + Δ t · A and B d = Δ t · B, with Δ t denoting the sampling period.

x k + 1 = A d x k + B d u k (21)

y k = C x k (22)

The parameter vector to be identified is defined by Equation 23.

θ = C in , C wall , R in , R out , COP nom , k T , a 0 , a 1 , a 2 T (23)

A least-squares criterion constructs the objective function as Equation 24:where N is the number of historical data samples; T in meas k is the measured indoor temperature at sampling instant k; and T in pred k , θ is the predicted value based on the model in Equations 21, 22. In this experiment, the 60-day dataset was divided chronologically: the first 45 days were used for parameter identification (training set), and the remaining 15 days were used to validate the model independently (test set).

min θ J θ = ∑ k = 1 N T in meas k − T in pred k , θ 2 (24)

Parameter identification must satisfy physical feasibility constraints given by Equations 25–29.

C in , ⁡ min ≤ C in ≤ C in , ⁡ max (25)

C wall , ⁡ min ≤ C wall ≤ C wall , ⁡ max (26)

R in , ⁡ min ≤ R in ≤ R in , ⁡ max , R out , ⁡ min ≤ R out ≤ R out , ⁡ max (27)

COP nom , ⁡ min ≤ COP nom ≤ COP nom , ⁡ max (28)

k T > 0 , a 0 + a 1 + a 2 = 1 (29)

Constraint boundaries can be estimated based on building type and equipment specifications. For instance, for typical office buildings: C in ∈ 50 , 200 kJ/°C (depending on room volume), R in ∈ 0.005 , 0.05 °C/kW, and COP nom ∈ 2.5 , 5.0.

To ensure model stability, eigenvalues of system matrix A must be negative real numbers, satisfying Equation 30.

a 11 + a 22 < 0 , a 11 · a 22 − a 12 · a 21 > 0 (30)

Considering parameter coupling and identifiability issues, a two-stage identification strategy is adopted, dividing the nine parameters into two groups for sequential identification.

Stage One: Thermal resistance-capacitance parameter identification. Temperature decay data from air-conditioning shutdown periods (P ac = 0) are selected for identification. Equation 21 then simplifies to Equation 31.

T in k + 1 T wall k + 1 = A d T in k T wall k + ρ c p V vent Δ t C in R in Δ t C wall R out T out k + Δ t C in 0 Q int k (31)

Parameters to be identified are θ 1 = C in , C wall , R in , R out T. Since T wall is typically not directly measurable, it can be treated as a state variable, employing extended Kalman filtering to jointly estimate states and parameters (Ziqing et al., 2022). The identification procedure comprises Steps 1.1 through 1.4:

Step 1.1: Filter periods when the air conditioner is off and the room is unoccupied to extract the temperature decay curve.

Step 1.2: Set initial parameter values θ 1 0 based on building drawings or empirical estimates.

Step 1.3: Solve Equation 32 using an optimization algorithm (such as a genetic algorithm):

θ 1 * = arg min θ 1 ∑ k ∈ K off T in meas k − T in pred k , θ 2 (32)

where K off denotes the set of sampling points during air conditioner shutdown periods.

Step 1.4: Verify the physical plausibility of the identification results τ = C in R eq by checking whether the thermal time constant lies within a reasonable range.

Stage Two: COP-related parameter identification. With thermal resistance-capacitance parameters fixed from Stage One, air-conditioning operating period data are utilized to identify COP model parameters θ 2 = COP nom , k T , a 0 , a 1 , a 2 T. The identification procedure comprises Steps 2.1 through 2.5:

Step 2.1: Identify periods of stable air conditioner operation characterized by minimal power fluctuations and gradual temperature changes.

Step 2.2: Based on Equations 11, 12, an implicit COP expression is constructed as Equation 33.

COP k = C in T in k + 1 − T in k / Δ t − Q wall k − Q vent k − Q int k P ac k (33)

Stage One parameters are used to calculate Q wall and Q vent, obtaining the COP observation sequence COP obs k for each sampling point.

Step 2.3: The corresponding temperature differences Δ T k and load ratios λ LR k are extracted to establish the regression model in Equation 34:

COP obs k = COP nom · 1 − k T Δ T k − Δ T nom · a 0 + a 1 λ LR k + a 2 λ LR 2 k + ϵ k (34)

where ϵ k denotes observation noise.

Step 2.4: Nonlinear least squares solves for parameters as in Equation 35, subject to constraint Equation 29.

θ 2 * = arg min θ 2 ∑ k ∈ K on COP obs k − COP pred k , θ 2 2 (35)

Step 2.5: If Stage Two results cause overall model prediction error to increase, iteration returns to Stage One for parameter fine-tuning until global error convergence.

2.3 Office park-level aggregated model construction

Buildings within a office park typically exhibit distributional patterns by type—for instance, by function (office buildings, conference centers, commercial facilities) or by envelope thermal characteristics (high-insulation versus standard buildings). This paper adopts a hierarchical clustering method implementing a two-stage strategy of coarse classification followed by fine clustering to achieve effective partitioning of building populations.

Stage One: Attribute-based coarse classification. Coarse classification partitions the M buildings within a office park into L basic categories based on fundamental physical attributes and functional characteristics. Classification criteria include: (1) Building functional type: office, conference center, restaurant, commercial facility, etc., as different functions exhibit markedly different usage patterns and load characteristics. (2) Structural form: high-rise/mid-rise/low-rise, frame/shear wall construction, etc., influencing envelope thermal performance. (3) Envelope grade: classified as high-insulation/standard-insulation/low-insulation based on exterior wall insulation thickness, window-to-wall ratio, glazing type, etc. (4) HVAC system type: variable refrigerant flow, air-cooled chiller, ground-source heat pump, etc., as different systems exhibit substantially different COP characteristics.

The attribute label vector for building i is defined by Equation 36:where l i func is the functional type label (e.g., 1-office, 2-conference, 3-dining), l i struct is the structural form label, l i envelope is the envelope grade label, and l i HVAC is the HVAC system type label.

l i = l i func , l i struct , l i envelope , l i HVAC T (36)

Coarse classification rules employ combinatorial matching: buildings with identical attribute labels are assigned to the same basic category. Let C l denote the building set of basic category l (l = 1 , 2 , ⋯ , L), as expressed by Equation 37:where l l is the standard attribute label for category l. The purpose of coarse classification is to ensure buildings within the same category possess comparable basic physical characteristics, establishing the foundation for subsequent parameter clustering.

C l = i : l i = l l , ∪ l = 1 L C l = 1 , 2 , ⋯ , M (37)

Stage Two: Parameter-based fine clustering. Within each basic category, fine clustering proceeds based on model parameter similarity. For buildings in basic category C l, the parameter feature vector is defined by Equation 38,where A i is building air-conditioned area (m²); V i is air-conditioned volume (m³); U i is the envelope overall heat transfer coefficient (W/(m²·°C)), calculable as U i = 1 / R in , i + R out , i; COP nom , i is rated air-conditioning COP; and τ i = C in , i · R eq , i is building thermal time constant (h), with R eq , i = R in , i · R out , i / R in , i + R out , i being the equivalent thermal resistance. These parameters collectively reflect building thermodynamic response characteristics.

f i = A i , V i , U i , COP nom , i , τ i T , i ∈ C l (38)

Given differing units and magnitudes among parameters, normalization processing is first performed according to Equation 39:where μ m and σ m are the mean and standard deviation of parameter m within category C l, and f ∼ i , m is the normalized feature value.

f ∼ i , m = f i , m − μ m σ m , m = 1 , 2 , ⋯ , 5 (39)

Within the normalized feature space, K-means clustering is applied to buildings in category C l (Ashani et al., 2022). If category C l is to be partitioned into K l sub-clusters, the optimization objective is given by Equation 40:where S l , j is the building set of sub-cluster j within basic category C l; f i is the normalized feature vector of building i; and μ l , j is the feature centroid vector of sub-cluster S l , j.

min S l , 1 , ⋯ , S l , K l ∑ j = 1 K l ∑ i ∈ S l , j f i − μ l , j 2 (40)

The number of sub-clusters K l is determined using the silhouette coefficient method. The silhouette coefficient for building i is defined by Equation 41,where a i is the average distance from building i to other buildings within the same cluster, and b i is the average distance from building i to the nearest neighboring cluster. The silhouette coefficient s i ∈ − 1 , 1, with larger values indicating better clustering quality. The K l value maximizing the average silhouette coefficient s ¯ = 1 C l ∑ i ∈ C l s i is selected.

s i = b i − a i max a i , b i (41)

Step 2.5 (Outlier rejection): After K-means clustering, the silhouette coefficient s i is recomputed for each building. Buildings with s i < 0 are flagged as outliers — they sit closer to a neighboring cluster than to their assigned one — and are either reassigned to the nearest cluster or, if persistent, isolated as singleton aggregation units. This step prevents data-quality anomalies (e.g., sensor drift, atypical retrofit) from contaminating cluster centroids and propagating into the capacity-weighted equivalent parameters.

Integration of two-stage clustering results yields a final partition of the office park into K = ∑ l = 1 L K l aggregation units upon completion of fine clustering across all basic categories. Each aggregation unit S j (j = 1 , 2 , ⋯ , K) contains buildings sharing both identical basic attribute labels and similar thermodynamic parameter characteristics. This hierarchical strategy holds certain advantages over directly clustering all buildings together. The coarse classification step ensures that buildings within the same aggregation unit are physically comparable, preventing buildings with substantially different functions or system types from being grouped together. The fine clustering step then draws on the precise parameter identification results to keep aggregation error under control. Because the two stages are executed independently, the granularity of classification can be adjusted to suit the specific conditions of a given office park.

For buildings within cluster S j, model parameters must be equivalenced to representative aggregated parameters. This paper adopts weighted averaging with weights corresponding to each building's air-conditioning installed capacity, as expressed by Equation 42,where θ ¯ j is the equivalent parameter of cluster j; P rated , i is the rated air-conditioning power of building i (kW); and θ i is the corresponding individual parameter. Using installed capacity as weights is justified because buildings with larger capacity exert more pronounced influence on aggregate office park power, warranting more adequate representation of their model characteristics in the aggregated model.

θ ¯ j = ∑ i ∈ S j P rated , i · θ i ∑ i ∈ S j P rated , i (42)

The aggregated state-space model is expressed by Equations 43, 44, where overlined variables denote post-aggregation equivalent quantities. The aggregated model state x ¯ j represents the average thermodynamic state of all buildings within the cluster, with the power component in input u ¯ j being the total air-conditioning power of that cluster.

x ¯ ˙ j t = A ¯ j x ¯ j t + B ¯ j u ¯ j t (43)

y ¯ j t = C ¯ j x ¯ j t (44)

2.4 Hierarchical control architecture and upper-level optimization model

Given the scale of a typical office park and the need for real-time computation, a two-tier hierarchical control architecture is adopted. The upper tier serves as the office park dispatch layer, where grid regulation commands are received, multi-objective optimization problems are solved using the aggregated models, and power allocation schemes are generated for each aggregation unit. The lower tier operates as the building execution layer, where unit-level electrothermal coupling models are used to implement predictive control, translating the allocated power into specific instructions for individual buildings and carrying out temperature regulation accordingly. Coordination between the two tiers is maintained through periodic information exchange: power commands and temperature bounds are passed downward from the upper tier, while actual power consumption and measured temperature states are fed back upward from the lower tier.

Decision variables for upper-level optimization are power sequences for each aggregation unit over the scheduling horizon. With a scheduling horizon of H time periods and K aggregation units, decision variable dimensionality equals K × H. Optimization objectives encompass both dispatch tracking and comfort assurance, converted to a single-objective function through weighted summation as Equation 45,where P j , h is allocated power for aggregation unit j in period h (kW); ω 1 and ω 2 are weight coefficients satisfying ω 1 + ω 2 = 1, adjusting the relative importance between objectives; J track is the dispatch tracking deviation term; and J comfort is the comfort deviation term.

min P j , h ⁡ J = ω 1 J track + ω 2 J comfort (45)

The dispatch tracking deviation term is defined as the weighted sum of squared deviations between aggregate office park power and dispatch commands, given by Equation 46:where P ref , h is the dispatch power command for period h (kW). This term drives optimization results toward minimizing deviation between aggregate office park power and dispatch targets.

J track = ∑ h = 1 H ∑ j = 1 K P j , h − P ref , h 2 (46)

The comfort deviation term employs a penalty formulation for temperature excursions beyond comfort bounds, given by Equation 47:where n j is the number of buildings in aggregation unit j; T ¯ j , h is the predicted average indoor temperature of unit j in period h; and g · is the temperature penalty function. The penalty function adopts a piecewise quadratic form as Equation 48,where T min and T max are the lower and upper comfort temperature limits (°C).

J comfort = ∑ h = 1 H ∑ j = 1 K n j · g T ¯ j , h (47)

g T = T − T max 2 , T > T max 0 , T min ≤ T ≤ T max T min − T 2 , T < T min (48)

Power bound constraints given by Equation 49:where P j , ⁡ min and P j , ⁡ max are power lower and upper limits for aggregation unit j (kW).

P j , ⁡ min ≤ P j , h ≤ P j , ⁡ max , ∀ j , h (49)

Power ramp rate constraints given by Equation 50:where Δ P j , ⁡ max is the maximum power rate of change for aggregation unit j (kW/period), employed to prevent drastic power fluctuations that could impact equipment and the grid.

P j , h − P j , h − 1 ≤ Δ P j , ⁡ max , ∀ j , h (50)

Temperature dynamic constraints given by Equation 51, representing the temperature dynamic relationship described by the aggregated model,where f · is the discretized form of Equations 43, 44. This constraint links optimization decisions with temperature evolution processes, ensuring physical feasibility of optimization results. The above optimization problem constitutes a quadratically-constrained quadratic program solvable using sequential quadratic programming or interior point methods (Aslam et al., 2024; Wang T. H. et al., 2024).

T ¯ j , h + 1 = f T ¯ j , h , P j , h , T out , h , ∀ j , h (51)

The weight coefficients ω 1 and ω 2 in the upper-level objective function (Equation 45) represent the aggregator's operational posture between grid compliance and contract comfort preservation. In practical deployment scenarios where a requested grid profile is physically infeasible due to the exhaustion of building envelope thermal storage, the quadratic nature of the penalty function (Equation 48) ensures that comfort preservation strictly dominates the objective. Rather than violating hard indoor boundaries, the framework execution layer automatically curtails its tracking participation, providing a physical boundary that can feed upward into market-level negotiation stages (such as OpenADR event acceptance). Furthermore, while a 120-period horizon is demonstrated, the framework is designed to run in an intra-day rolling horizon manner to easily absorb short-notice notification updates and weather forecast errors.

A note on disturbance handling: the lower-level MPC operates on measured indoor temperature with a 5-min rolling horizon. Any disturbance not represented in the upper-level aggregated model — including solar gain forecast errors, transient occupancy fluctuations, infiltration variability, or sensor noise — manifests as a deviation between predicted and measured T in at the next sampling instant and is corrected within the subsequent control cycle. This feedback structure is by design: the upper layer carries physical fidelity sufficient for grid-side dispatch, while the lower layer carries robustness against unmodeled local disturbances. The empirical robustness boundaries of this architecture are quantified in Section 3.2.4.

2.5 Lower-level predictive control

The core task of lower-level control is appropriately allocating received upper-level power commands to individual buildings and implementing specific regulation. For building i, given current state x i t 0 and the input sequence u i t t = t 0 t 0 + T p over the prediction horizon, predicted indoor temperature values can be obtained by solving the state-space model as Equations 52, 53, where T p is prediction horizon length; Δ t is prediction step size; and e A i Δ t is the matrix exponential of the state transition matrix. Practical computations typically employ Euler or Runge-Kutta methods for numerical solution.

x i t + Δ t = e A i Δ t x i t + ∫ 0 Δ t e A i τ B i u i t + Δ t − τ d τ (52)

T in , i pred t + Δ t = C i x i t + Δ t (53)

Based on temperature prediction results, a prediction-assisted rule-based control strategy is adopted to refine upper-level allocated aggregate power to individual buildings. First, regulation priority is computed for each building according to temperature deviation magnitude, given by Equation 54,where π i t is the regulation priority of building i at time t, with positive values indicating need for increased cooling/reduced heating and negative values indicating need for increased heating/reduced cooling. Power allocation to buildings within the aggregation unit then follows principles expressed by Equations 55, 56, where P i 0 t is the baseline power required for building i to maintain current temperature, and Δ P i t is the adjustment relative to baseline power.

π i t = T in , i pred t + T p − T max , T in , i pred t + T p > T max T min − T in , i pred t + T p , T in , i pred t + T p < T min 0 , (54)

P i t = P i 0 t + Δ P i t (55)

Δ P i t = π i t ∑ k ∈ S j π k t · P ¯ j t − ∑ k ∈ S j P k 0 t (56)

3 Results and discussion

3.1 Case overview

Simulation validation is conducted using a technology office park located in an eastern Chinese city. The office park has a total building area of approximately 120,000 square meters, comprising 15 office buildings, each equipped with variable refrigerant flow air-conditioning systems. The aggregate installed cooling and heating capacity of these systems is approximately 4200 kW. An energy management system is in place, collecting air-conditioning power consumption, indoor temperature, and other operational data at 5-min sampling intervals. The simulation was performed on MATLAB R2024b running on a personal computer equipped with an Intel Core i7-14700 KF processor (3.40 GHz, 20 cores), 64 GB of RAM, and a 1 TB NVMe SSD.

Based on architectural characteristics, the 15 buildings were divided into two basic categories during the coarse classification stage: the first category consists of 11 office buildings, and the second category consists of 4 conference centers and dining facilities. K-means-based fine-grained clustering was performed on the basis of this coarse classification. Setting the number of candidate sub-clusters as K l ∈ 2 , 3 , 4, the calculated average silhouette coefficient s ¯ is shown in Table 2. It can be seen that s ¯ reaches its maximum when K l = 2, indicating that dividing these 11 office buildings into two sub-clusters is most appropriate. The two sub-clusters correspond to high-rise buildings with large floor areas, low thermal resistance, and high rated COP (Unit A), and mid-rise buildings with moderate thermal and structural parameters (Unit B). The four buildings in the second category have extremely similar thermal parameters and are retained as a single aggregated unit (Unit C). Therefore, the final office park-level classification consists of three aggregation units (K l = 3). Table 3 presents the main parameters of each aggregated unit, while the simulation scenarios and optimized control parameters are shown in Supplementary Appendix B, Supplementary Appendix Table B1.

TABLE 2. Quantitative Selection of Cluster Counts (K) based on Average Silhouette Coefficients.

TABLE 3. Aggregated thermophysical and operational parameters for each aggregation unit (K = 3).

Outdoor temperature employs local meteorological data, with daytime temperatures varying between −2 °C and 6 °C. The comfort temperature range is set at 20–24 °C. Grid dispatch commands employ time-varying power curves simulating actual demand response scenarios, including three regulation periods: morning peak shaving (8:30–10:00), midday valley filling (12:00–13:30), and evening peak shaving (16:00–17:30).

Three comparative schemes are established to validate the proposed method: Scheme 1 (Proposed Method) adopts the proposed electrothermal coupling modeling and hierarchical coordinated optimization strategy, with upper-level multi-objective optimization based on aggregated models and lower-level predictive control based on unit models. Scheme 2 (Simplified Model Method) employs first-order simplified thermodynamic models (neglecting envelope thermal storage effects) for upper-level optimization with conventional PID temperature control at the lower level. Scheme 3 (No-Optimization Baseline Profiling Reference) forgoes office park-level optimization scheduling, with each building operating independently and air-conditioning systems executing only local temperature setpoint-based on/off control. This scheme provides a synthetic but realistic estimate of the park's natural, uncoordinated load profile, serving as the required baseline to gauge the absolute participation level and flexibility volume achieved by Scheme 1. To benchmark the proposed framework against state-of-the-art alternatives, two additional schemes are introduced. Scheme 4 implements a deep deterministic policy gradient controller with state space and action space trained on the same 60-day operational dataset under a reward function combining dispatch tracking and comfort penalty. Scheme 5 implements a simplified two-stage distributionally-robust optimization framework, using a Wasserstein ambiguity set of radius epsilon = 0.1 over 50 calibrated scenarios. Both schemes operate on the same aggregated model and dispatch commands as Scheme 1 to ensure controlled comparison.

To benchmark the proposed framework against state-of-the-art alternatives, two additional schemes are introduced. Scheme 4 implements a deep deterministic policy gradient controller with state space (T i n , j, T o u t, P r e f; P t o t a l) and action space (Δ P 1 , … , Δ P k), trained on the same 60-day operational dataset under a reward function combining dispatch tracking and comfort penalty. Scheme 5 implements a simplified two-stage distributionally-robust optimization framework following, using a Wasserstein ambiguity set of radius ε = 0.1 over 50 calibrated scenarios. Both schemes operate on the same aggregated model and dispatch commands as Scheme 1 to ensure controlled comparison.

Evaluation metrics include: (1) Dispatch tracking metric, employing normalized root mean square error (NRMSE) as Equation 57; (2) Comfort assurance metric (CAM), employing temperature violation duration ratio as Equation 58; and (3) Modeling accuracy metric (MAM), employing indoor temperature prediction root mean square error (RMSE) as Equation 59.

NRMSE track = 1 H ∑ h = 1 H P total , h − P ref , h 2 P ref , mean × 100 % (57)

VT comfort = ∑ i , h I T i , h ∉ T min , T max M × H × 100 % (58)

RMSE temp = 1 N ∑ k = 1 N T in pred k − T in meas k 2 (59)

To ensure absolute mathematical rigor and experimental reproducibility, the advanced baseline methods are configured as follows:

• Scheme 4 (Deep Deterministic Policy Gradient - DDPG): The state space S comprises four continuous variables: (T in avg, T out; P grid t). The action space A represents the continuous total power allocation vector for the aggregation units. The reward function R t is formulated as the negative sum of tracking deviations and comfort penalties, matching our objective function (Equation 45). Actor and Critic networks both utilize a 3-layer MLP structure (128-64–32 neurons) with ReLU activations. Hyperparameters include a learning rate α a = 10 − 4 for the Actor and α c = 10 − 3 for the Critic; a discount factor γ = 0.95; a replay buffer size of 10^–5; and a batch size of 64. The model was trained over 2000 episodes until moving-average reward convergence.
• Scheme 5 (Distributionally Robust Optimization - DRO): The operational and weather uncertainties are modeled using a data-driven Wasserstein ambiguity set centered around 50 calibrated historical scenarios of outdoor temperatures and internal loads. The Wasserstein radius is set precisely to ε = 0.1 to represent a robust 90% confidence envelope over empirical distributions. The mathematical structure is formulated as a two-stage min-max-min optimization problem. The inner maximization searches for the worst-case probability distribution within the Wasserstein ball, while the outer optimization minimizes the aggregate dispatch tracking error and comfort violation costs under this adversarial distribution. It was solved using CVX coupled with MOSEK via linear decision rules (LDR) to ensure tractability.

3.2 Results analysis

3.2.1 Modeling accuracy verification

Prediction accuracy of the proposed electrothermal coupling model is first verified. Three representative buildings, one from each aggregation unit, are selected; following parameter identification using historical operational data, model predictions are compared against measured data. Figure 3 displays indoor temperature prediction curves for a typical day.

As observable from Figure 3, model prediction curves for all three buildings exhibit good agreement with measured data, with predictions accurately capturing dynamic temperature variation trends. Specifically, the Unit A representative building achieves temperature prediction RMSE of 0.246 °C, Unit B achieves 0.096 °C, and Unit C achieves 0.0315 °C, satisfying engineering application accuracy requirements. Temporal distribution analysis of prediction errors reveals that errors concentrate primarily in two types of periods: first, during early-morning air-conditioning startup (8:00–9:00), when indoor temperature rises rapidly from overnight lows and the model exhibits slight deviation in characterizing this transient process; second, around midday (11:00–13:00), when frequent indoor occupant activity causes substantial internal heat source fluctuations that interfere with model predictions.

Figure 4 compares temperature predictions from unit model ensembles versus aggregated models under identical input conditions. Solid lines represent weighted averages of all unit model predictions within each aggregation unit, dashed lines represent corresponding aggregated model predictions, and shaded regions indicate the distribution range of unit model predictions. The aggregated model prediction curves closely track unit model weighted average curves, with deviations not exceeding 0.3 °C during most periods, validating the effectiveness of the weighted parameter equivalence method. Unit A exhibits the narrowest shaded region, indicating highest thermodynamic characteristic consistency among buildings in that cluster and best aggregation performance. Near 9:00 and 14:00, deviations between aggregated model and unit average slightly increase due to amplified differences in building response speeds during transient processes, though deviations remain within acceptable bounds.

For aggregated models, accuracy loss relative to unit models merits particular attention. Aggregation accuracy loss ε j is defined as indicated in Equation 60:

ε j = 1 S j ∑ i ∈ S j T in , i pred − T ¯ in , j pred 2 T in , i pred 2 (60)

In the formula, S j represents the number of buildings in cluster j; T in , i pred represents the temperature forecast for building i using the individual model; T ¯ in , j pred represents the temperature forecast using the aggregated model. If ε j exceeds the set threshold, the number of clusters must be increased or the clustering strategy adjusted with calculation results showing: Unit A aggregation accuracy loss of 6.2%, Unit B of 7.8%, and Unit C of 5.5%, with three-unit average loss of 6.5%. This implies that employing aggregated models in place of unit model ensembles for prediction increases temperature prediction error by less than 8%, within engineering-acceptable bounds.

3.2.2 Optimization control effectiveness verification

Figure 5 presents comparative aggregate office park power tracking results for the three schemes over the simulation day.

Differences in dispatch tracking performance across the three schemes are visible in Figure 5. The aggregate office park power curve under Scheme 1 follows the dispatch command curve most closely, maintaining acceptable tracking across all three regulation periods. Scheme 2 shows some tracking tendency, but response lag and overshoot are apparent, with deviations becoming larger at the boundaries between regulation periods. Under Scheme 3, the power curve is largely unresponsive to the dispatch commands, being governed instead by indoor temperature control requirements, and deviates from the dispatch target in an irregular manner throughout the simulation period.

Table 4 presents quantified performance comparisons across the three schemes. In terms of dispatch tracking, the proposed method achieves an NRMSE of 3.4%, which represents a 52.78% reduction relative to the simplified model method and a 86.18% reduction relative to the no-optimization baseline. This improvement can be attributed to the electrothermal coupling model's accurate characterization of the dynamic relationship between electrical power input and indoor temperature, which allows the upper-level optimizer to anticipate the consequences of control actions on future power and temperature states before committing to a decision. The average response delay metric is consistent with this interpretation: the proposed method achieves a delay of 3.2 min, compared with 8.7 min for the simplified model method. The simplified model, by neglecting envelope thermal storage effects, introduces systematic error into temperature change predictions, which in turn causes the timing of control actions to be imprecise.

TABLE 4. Performance metrics comparison of each solution.

On the comfort side, the proposed method yields a temperature violation duration ratio of 2.1%, 67% lower than that of the simplified model method. The no-optimization baseline produces a lower violation ratio than the simplified model method—a result that at first appears counterintuitive, but in fact points to an inherent weakness of the simplified model approach. When dispatch tracking is pursued with an inaccurate model, temperature constraint predictions become unreliable, and the resulting control actions can generate more comfort violations than would occur without any active dispatch coordination. The no-optimization baseline, by contrast, does not attempt dispatch tracking and instead allows each building to regulate temperature independently, which is sufficient to satisfy local comfort requirements. The gap between the simplified model method and the no-optimization baseline on this metric suggests that model accuracy is not merely a matter of prediction quality, but directly conditions whether comfort constraints can be respected under coordinated dispatch.

To verify the symmetric applicability of the proposed framework under cooling-mode operation, the simulation was repeated for a typical summer day (outdoor temperature 28–35 °C, comfort range 24–28 °C, dispatch commands inverted: morning valley filling and afternoon peak shaving). The aggregated model parameters were re-identified using a separate 60-day summer dataset from the same park. Under this scenario the proposed method achieved NRMSE = 3.7% and CAM = 2.4%, consistent with the winter case (3.4% and 2.1%, respectively). This supports the climate/season-agnostic nature of the framework: adaptation to a new operating regime is achieved through re-identification of the same model structure, not through structural modification of the methodology itself.

3.2.3 Synthetic extended park validation

To probe the framework's behavior beyond the 15-building case, a 45-building synthetic park is constructed via statistical replication of the calibrated parameter distribution (σ = 10% multiplicative perturbation per parameter, sampled with replacement from the empirical distribution of the 15 identified buildings). All five schemes are evaluated on the extended scenario. Table 4 reports aggregate results. The proposed method maintains NRMSE = 3.7% and CAM = 2.4% on the 45-building park (vs. 3.4%/2.1% on the 15-building case). Scheme 4 performance degrades more sharply with park size due to action space dimensionality (NRMSE = 5.2%), while Scheme 5 tracks closely (NRMSE = 4.1%). This indicates that the structural advantage of the proposed framework — physical aggregation rather than learned aggregation — is more pronounced at scale.

3.2.4 Solver performance and scalability

Table 5 summarizes the solver performance metrics for the two schemes under the baseline scenario.

TABLE 5. Solver performance metrics (K = 3; H = 120).

Compared to the simplified model approach, the proposed method takes slightly longer to solve, which is attributed to the additional temperature dynamics constraints imposed by the second-order aggregation model. This model introduces an envelope temperature as an additional state variable for each aggregation unit, increasing the number of constraints from 720 to 1,086. Nevertheless, for a full 10-h scheduling day, the total actual runtime is 46.8 s, which falls entirely within the 5-min scheduling cycle specified by the energy management system, confirming its real-time feasibility in the current scenario.

Table 6 reports the upper-level per-period solve time and total wall-clock time as a function of aggregation unit count K, with H = 120 fixed. As shown in Table 5, solve time grows approximately as O K 1.6 over the tested range K ∈ 3 , 5 , 10 , 20 , 30, consistent with the polynomial complexity of interior-point methods for QCQPs. At K = 10, the per-period solve time remains below 1.2 seconds—well within the 5-min dispatch cycle. Even at K = 30, the per-period solve time of 9.12 s remains comfortably below the 300-s real-time boundary, indicating that the proposed framework scales gracefully to office parks considerably larger than the present case study. For office parks requiring K > 30 aggregation units, however, decomposition-based approaches such as ADMM would be a natural extension to maintain real-time tractability with margin, and this is identified as a direction for future work.

TABLE 6. Exploratory statistics.

3.2.5 Sensitivity, error margins, and structural robustness analysis

To quantify how modeling error margins propagate into real-time dispatch decisions, two distinct uncertainty boundaries must be examined: static aggregation error margins and dynamic unmodeled disturbances.

As established in Section 3.2.1, the capacity-weighted parameter equivalence induces a bounded structural error margin, yielding an average aggregation accuracy loss of 6.5% (≤8% across all clusters). This modeling margin represents the physical trade-off required to collapse a high-dimensional multi-building state-space into a low-dimensional aggregated model suitable for real-time upper-level coordination.

Dynamic, unmodeled uncertainties—such as transient internal heat gains from occupancy spikes or solar irradiance deviations—are handled structurally through the lower-level MPC's 5-min rolling horizon loop. Any mismatch between the aggregated upper-level prediction and the local reality manifests as a localized temperature deviation at the next sampling step. The lower-level prioritized adjustment rule (Equations 54–56) then realigns individual building baseline power allocations to dampen these deviations within a single control cycle. The statistical resilience of this decoupled design is confirmed by the 45-building synthetic experiment. Even when a sigma = 10% multiplicative perturbation is forced onto the physical parameters of the extended building fleet, the optimization metrics remain remarkably stable (with NRMSE bounding at 3.7% and CAM at 2.4%). This demonstrates that the physical fidelity preserved during the aggregation stage, coupled with the fast feedback stabilization at the execution stage, gives the framework structural immunity against typical engineering error margins without requiring explicit robust optimization.

4 Conclusion

Addressing the refined modeling and large-scale scheduling requirements for office park air-conditioning loads participating in demand response, this paper has proposed an electrothermal coupling modeling and hierarchical coordinated optimization methodology oriented toward aggregate dispatch of air-conditioning-building systems. The principal contributions and findings are as follows:

• A unit-level electrothermal coupling model encompassing envelope thermal storage effects and dynamic air-conditioning energy efficiency characteristics has been established, achieving indoor temperature prediction RMSE controlled within 0.35 °C.
• A office park-level aggregated model construction method based on aggregability principles has been proposed. While ensuring aggregation accuracy loss below 8%, this method effectively reduces optimization problem scale, rendering real-time office park-level optimization scheduling feasible.
• A hierarchical coordinated optimization architecture integrating upper-level multi-objective optimization with lower-level predictive control has been designed. Simulation validation on an office park demonstrates that, compared to conventional simplified model methods, the proposed approach reduces dispatch tracking deviation by 42%, decreases comfort violation duration by 67%, and achieves more stable power fluctuation characteristics.

In most cases, office parks involve multiple building owners, tenants, and energy contracts. The next research objective is to determine how to achieve integrated scheduling for office parks in a multi-stakeholder environment.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg.2026.1844532/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Muhammad Waseem, Maynooth University, Ireland

Reviewed by: Michael Short, Teesside University, United Kingdom

Ahmed Salman, Northern Technical University, Iraq