Aggregate Control Strategy for Thermostatically Controlled Loads with Demand Response

The improvement of intelligent appliances provides the basis for the demand response (DR) of residential loads. Thermostatically controlled loads (TCLs) are one of the most important DR resources and are characterized by a large load and a high degree of control. Due to its distribution characteristic, the aggregation of TCLs and their control are key issues in implementing the load control for the DR. In this study, we focus on air conditioning loads as an example of TCLs and propose a simple and transferable aggregate model by establishing a virtual house model, which accurately captures the aggregate flexibility. The deviation of the aggregate model is analyzed for the model evaluation. An air conditioning DR control scheme is proposed based on the aggregate model; it has the advantage of simple implementation and convenient control for the individual units. Simulations are performed in Gridlab-D to evaluate the accuracy and effectiveness of the proposed model and control method.


Introduction
Power consumers participate in the control and operation of power systems through the demand response (DR), which results in considerable benefits, such as improving the operation efficiency of the power system, reducing the cost of electricity, and even increasing the social benefits [1,2]. With the rapid development of advanced metering infrastructures and bidirectional communication, residential consumers can easily take part in the DR [3]. Thermostatically controlled loads (TCLs) such as air conditioning loads are one of the most important DR resources. The DR of TCLs can be applied in frequency regulations, capacity reserves, ancillary services, etc. Current research on TCLs and its DR has mainly focused on aggregate models and DR schemes.
The purpose of aggregate models of TCLs is to describe the probability density evolution of a population of TCLs. Many mathematical models have been used for this purpose such as regression models [4], physics-based models [5], stochastic Fokker-Planck diffusion models [6], partial differential equation (PDE) models [7,8], and Markov chain models [9,10]. These models are used to calculate the relationship between the control variables of TCLs and the change in the load level.
An additional research focus is the achievement of the DR. Lu et al. [11] presented a queueing model to analyze the DR based on price aggregate TCLs. Callaway [12] regulated aggregated homogeneous TCLs by changing the setpoint temperature. Bashash and Fathy [7] achieved the DR for homogeneous TCLs by setting changing rate of the setpoint temperature. Using on/off control signals,

Modeling of Air Conditioning Load
In this section, the equivalent thermal parameters (ETP) model, which defines the temperature dynamics of the air conditioner is presented. To develop an aggregate model for air conditioners, we propose a virtual house model to approximate the relationship between power consumption and average indoor temperature. Based on this model, a DR scheme is put forward. It is worth mentioning that the virtual house model is only used for an approximate analysis and the deviation is presented in Section 5. A. Physical model of air conditioning loads The components of the air conditioning system include a compressor, condenser, fan, and other parts. The compressor cools the liquid by compressing the refrigerant (e.g., Freon); an endothermic reaction occurs and the liquid changes its state at atmospheric pressure [17]. Similarly, an air conditioning system can produce heat [18]. The deadband is an important parameter to consider in a thermostat [19]. An air conditioning thermostat is used to control the indoor temperature and is used to adjust the operating time constant of the air conditioning equipment to ensure that the indoor temperature remains within a given range.
The heat flow between the inside of the house and the outside is a key factor affecting the air conditioning load. For a house, the following types of heat exchanges occur: (1) heat conduction through facades, roofs, and windows, (2) heat exchange due to indoor-outdoor airflow, (3) solar heat irradiating to the house directly (sun), and (4) heat generated by electrical equipment (internal gains) [20].
The thermal flow is shown in Figure 1. As shown in Figure 1, we can obtain the expressions of the air temperature node and the wall temperature node: Equation (1) [20] is solved for TM, differentiated with respect to time to provide dTM/dt, and both are substituted into Equation (2) to create a second-order linear differential equation in TA of the form: The temperature, thermal mass, and heat flow in physics are equivalent to the voltage, capacitor, and current flow in an electric circuit. Therefore, a circuit relationship equivalent to the thermal relationship can be developed, as shown in Figure 2. In practice, this circuit is always overdamped. In other words, the temperature and load exponentially decay and approach steady-state (not oscillatory) conditions. As shown in Figure 1, we can obtain the expressions of the air temperature node and the wall temperature node: Equation (1) [20] is solved for T M , differentiated with respect to time to provide dT M /dt, and both are substituted into Equation (2) to create a second-order linear differential equation in T A of the form: where we can obtain an expression for T A : where The temperature, thermal mass, and heat flow in physics are equivalent to the voltage, capacitor, and current flow in an electric circuit. Therefore, a circuit relationship equivalent to the thermal relationship can be developed, as shown in Figure 2. In practice, this circuit is always overdamped. In other words, the temperature and load exponentially decay and approach steady-state (not oscillatory) conditions.  [7] or other special cases [21]. As mentioned above, many studies on aggregate models of TCLs have been conducted. Different aggregation methods have been used in these aggregate models depending on whether the goal was to improve the accuracy or an easy implementation of the load control.
From the perspective of economy and feasibility, in aggregate models, the independent control of a temperature controlled appliance (TCA) should be avoided. The requirement of independent control of each TCA means that control devices and instant communication devices have to be installed in each TCA. In fact, aggregate control of TCAs is often implemented in residential houses, shopping malls, dormitories, and hotels. In these situations, there is a strong consistency between every TCL. The scheme achieves DR of the aggregated TCLs using simple control signals. This type of scheme is proposed and is described next.
The aggregate model for achieving DR of the aggregated TCLs using simple control signals is developed by creating a virtual house that includes all TCAs. The virtual house ignores the state parameters of the aggregated units and calculates the total loads over a period by using only the initial state quantities and characteristic parameter values of the aggregated units. In other words, we propose a virtual house, in which the heat exchange with the outdoors is equivalent to the sum of the heat exchange of the aggregated units.
The virtual house is created using the following steps. A virtual house consisting of two houses is used as an example.
1. Two houses have their own initial indoor temperatures and setpoints. 2. We assume that there is a pipe with a switch between the two houses. There is no heat exchange between the pipe and the outside. If the switch is on, the pipe creates infinite heat exchange between the two houses. It is worth mentioning that this pipe does not exist in reality and is only used to understand the analysis. 3. If the switch is on, the air temperatures inside the two houses are the same in the next unit of time. The two houses are equivalent to one virtual house and the area of the virtual house is the sum of the area of the two houses. The heat exchange coefficient with the outside world is equal to the sum of the heat exchange coefficients of the two houses. Although the coefficient of heat exchange is the same, the heat exchange between the virtual house and the outside is not strictly equal to the heat exchange between the two houses and the outside because of the indoor temperature difference. 4. In Gridlab-D, this progress means establishing a house whose area is equal to the sum of the area of two small houses and modifying the corresponding heat exchange parameters to make them equivalent, as shown in Figure 3. In the circuit equivalent model, the resistance and capacitance are parallel, as shown in Figure 3. . The virtual house model In TCLs, aggregate load models are used to study the effects of DR control [7] or other special cases [21]. As mentioned above, many studies on aggregate models of TCLs have been conducted. Different aggregation methods have been used in these aggregate models depending on whether the goal was to improve the accuracy or an easy implementation of the load control.
From the perspective of economy and feasibility, in aggregate models, the independent control of a temperature controlled appliance (TCA) should be avoided. The requirement of independent control of each TCA means that control devices and instant communication devices have to be installed in each TCA. In fact, aggregate control of TCAs is often implemented in residential houses, shopping malls, dormitories, and hotels. In these situations, there is a strong consistency between every TCL. The scheme achieves DR of the aggregated TCLs using simple control signals. This type of scheme is proposed and is described next.
The aggregate model for achieving DR of the aggregated TCLs using simple control signals is developed by creating a virtual house that includes all TCAs. The virtual house ignores the state parameters of the aggregated units and calculates the total loads over a period by using only the initial state quantities and characteristic parameter values of the aggregated units. In other words, we propose a virtual house, in which the heat exchange with the outdoors is equivalent to the sum of the heat exchange of the aggregated units.
The virtual house is created using the following steps. A virtual house consisting of two houses is used as an example.

1.
Two houses have their own initial indoor temperatures and setpoints.

2.
We assume that there is a pipe with a switch between the two houses. There is no heat exchange between the pipe and the outside. If the switch is on, the pipe creates infinite heat exchange between the two houses. It is worth mentioning that this pipe does not exist in reality and is only used to understand the analysis.

3.
If the switch is on, the air temperatures inside the two houses are the same in the next unit of time. The two houses are equivalent to one virtual house and the area of the virtual house is the sum of the area of the two houses. The heat exchange coefficient with the outside world is equal to the sum of the heat exchange coefficients of the two houses. Although the coefficient of heat exchange is the same, the heat exchange between the virtual house and the outside is not strictly equal to the heat exchange between the two houses and the outside because of the indoor temperature difference.

4.
In Gridlab-D, this progress means establishing a house whose area is equal to the sum of the area of two small houses and modifying the corresponding heat exchange parameters to make them equivalent, as shown in Figure 3.
In the circuit equivalent model, the resistance and capacitance are parallel, as shown in Figure 3. The virtual house model is expected to reflect the relationship between the average indoor temperature of the aggregate and the total load, which means that the heat exchange between the virtual house and the outside equals the sum of the heat exchange between each house and the outside. This needs to be verified by mathematical derivation.

Aggregate Model validity analysis
Two scenarios are created and analyzed to improve the accuracy of the proposed aggregate (virtual house) model and the effectiveness of the aggregated load control.
Scenario one: Homogeneous model The homogeneous model consists of houses with the same areas, structures, and building materials. Assuming that there are house A1 and house A2, according to the ETP model, the following equation is obtained.
In the homogeneous model, a1 = a2 = … = an, the same applies to b, c, and d. Equation (5) can be simplified as.
Equation (6) has the form of a second-order linear differential equation for TA (Equation (3)). Therefore, the virtual house is equivalent to all houses. Similarly, a virtual house can be equivalent to several houses with the same structure but different area.
In this case, the initial temperature of the equivalent house can be calculated using the law of energy conservation. At the initial time step, the thermal mass of the virtual house should be equal to the sum of the thermal masses of all houses, as shown in Equation (7): where c is the specific heat capacity of air, ρ is the air density, h is the height of the house, and n is the number of floors of the house. Assuming the height of each house is the same and there is only one floor, the initial temperature of the equivalent house can be calculated as follows. The virtual house model is expected to reflect the relationship between the average indoor temperature of the aggregate and the total load, which means that the heat exchange between the virtual house and the outside equals the sum of the heat exchange between each house and the outside. This needs to be verified by mathematical derivation.

Aggregate Model Validity Analysis
Two scenarios are created and analyzed to improve the accuracy of the proposed aggregate (virtual house) model and the effectiveness of the aggregated load control.
Scenario one: Homogeneous model The homogeneous model consists of houses with the same areas, structures, and building materials. Assuming that there are house A1 and house A2, according to the ETP model, the following equation is obtained.
In the homogeneous model, a 1 = a 2 = . . . = a n , the same applies to b, c, and d. Equation (5) can be simplified as.
Equation (6) has the form of a second-order linear differential equation for T A (Equation (3)). Therefore, the virtual house is equivalent to all houses. Similarly, a virtual house can be equivalent to several houses with the same structure but different area.
In this case, the initial temperature of the equivalent house can be calculated using the law of energy conservation. At the initial time step, the thermal mass of the virtual house should be equal to the sum of the thermal masses of all houses, as shown in Equation (7): Energies 2019, 12, 683 where c is the specific heat capacity of air, is the air density, h is the height of the house, and n is the number of floors of the house. Assuming the height of each house is the same and there is only one floor, the initial temperature of the equivalent house can be calculated as follows.
For the homogeneous model, the heat exchange of the virtual house is equivalent to the heat exchange between the small houses and the outside. In this case, the average temperature calculated by the virtual house model is equal to the actual value and no error is introduced in the load control.
Scenario two: Heterogeneous model When the areas and structures of the two rooms are different, for the heterogeneous model, Equation (5) cannot be expressed as the form of a second-order linear differential equation for T A (Equation (3)). This means that there will be an error in the proposed aggregate model. To verify the applicability of the aggregate model, in this case, an error analysis is required and the effects of the error on our control strategy have to be considered.
To analyze the error of the heterogeneous loads model, it is necessary to introduce the model of the houses in Gridlab-D. In Gridlab-D, the specific relationship between the building structure, building materials, and area are shown in the following equation. Table 1 shows that the structure and the heat transfer coefficient of the house are related to the area. Therefore, the heterogeneous model can be established by selecting different floor areas of the houses in Gridlab-D. In this study, an aggregate of 2000 houses is used and different areas are evaluated to observe the error of the proposed aggregate model for the heterogeneous conditions. Table 1. Relationship between the house structure and area Table 2 lists the areas of the different houses. Rows one to five represent the areas of the houses that meet the uniform distributions. Rows six to eight represent the areas of the houses that meet the Gaussian distribution. Rows nine to twelve represent the areas of the houses that meet a different distribution range. The average area of the houses in all cases is 2000 square feet. In these cases, the indoor average temperature error caused by the heterogeneity of the house is shown in Figure 4. It is evident that there is no direct relationship between the degree of the heterogeneity of the houses and the size of the error. The difference in the average temperature never exceeds 0.08 degrees Fahrenheit in these cases.  Table 2 lists the areas of the different houses. Rows one to five represent the areas of the houses that meet the uniform distributions. Rows six to eight represent the areas of the houses that meet the Gaussian distribution. Rows nine to twelve represent the areas of the houses that meet a different distribution range. The average area of the houses in all cases is 2000 square feet. In these cases, the indoor average temperature error caused by the heterogeneity of the house is shown in Figure 4. It is evident that there is no direct relationship between the degree of the heterogeneity of the houses and the size of the error. The difference in the average temperature never exceeds 0.08 degrees Fahrenheit in these cases.  In fact, many studies [10,22] have shown that homogeneous load populations often result in strong oscillations, whereas well-diversified loads result in natural damping and a more stable aggregated response.  In fact, many studies [10,22] have shown that homogeneous load populations often result in strong oscillations, whereas well-diversified loads result in natural damping and a more stable aggregated response. Figure 5 shows a comparison of the aggregated responses of 2000 homogeneous and heterogeneous TCLs for a heat setpoint change from 23.3 • C to 24.4 • C. The power curves of 2000 houses with an average area of 185.8 m 2 are shown. U100 and U200 represent the uniform distribution of the house area with a variance of (2 × 100) 2 /12 and (2 × 200) 2 /12 respectively. For the heterogeneous model, the suppression of load oscillations is strong and the higher the degree of dispersion of the aggregated units, the stronger the ability to suppress the oscillations is.  Table 2 lists the areas of the different houses. Rows one to five represent the areas of the houses that meet the uniform distributions. Rows six to eight represent the areas of the houses that meet the Gaussian distribution. Rows nine to twelve represent the areas of the houses that meet a different distribution range. The average area of the houses in all cases is 2000 square feet. In these cases, the indoor average temperature error caused by the heterogeneity of the house is shown in Figure 4. It is evident that there is no direct relationship between the degree of the heterogeneity of the houses and the size of the error. The difference in the average temperature never exceeds 0.08 degrees Fahrenheit in these cases.  In fact, many studies [10,22] have shown that homogeneous load populations often result in strong oscillations, whereas well-diversified loads result in natural damping and a more stable aggregated response. Figure 5    This temperature error will eventually affect the accuracy of the load control under DR and even lead to load bound and oscillations. In the following section, the effect of the temperature error on the load control will be analyzed. For an air conditioning system, the indoor air temperature can be determined using Equation (3):

Symbol Calculation
where Figure 6 shows the process of raising the setpoint temperature from T 1 to T 2 at time t 3 for a single air conditioning system. It is observed that the air conditioner is in the rated power (in the running state) during the period of temperature increase and it is in a suspended state during the period of temperature decrease, of which the power is zero. Therefore, for a single air conditioning unit, the probability p when the air conditioning unit is in the rated power running state can be expressed as the period of the temperature increase/the period of one cycle. It can be calculated using the following equation: Energies 2019, 12 FOR PEER REVIEW 8 This temperature error will eventually affect the accuracy of the load control under DR and even lead to load bound and oscillations. In the following section, the effect of the temperature error on the load control will be analyzed. For an air conditioning system, the indoor air temperature can be determined using Equation (3):  Figure 6 shows the process of raising the setpoint temperature from T1 to T2 at time t3 for a single air conditioning system. It is observed that the air conditioner is in the rated power (in the running state) during the period of temperature increase and it is in a suspended state during the period of temperature decrease, of which the power is zero. Therefore, for a single air conditioning unit, the probability p when the air conditioning unit is in the rated power running state can be expressed as the period of the temperature increase/the period of one cycle. It can be calculated using the following equation: Figure 6. Indoor air temperature for a single air conditioning system. As mentioned above, the operating state of the air conditioning system is divided into the running period and the suspending period. In each period, the indoor air temperature can be obtained by Equation (9). The time required for the room temperature to reach a certain temperature value can be obtained by Equation (10).
Therefore, t2, t3, and t4 can be solved by Equation (10). When the aggregated unit is in the steady state, the value NA of the number of air conditioning units at the rated power running condition is At any time t, the probability that the air conditioning unit i is in the rated power running state is pi(t) and the probability that the air conditioning unit i is in the suspended state is 1-pi(t). Therefore, the power of the aggregate is expected to satisfy Equation (11). As mentioned above, the operating state of the air conditioning system is divided into the running period and the suspending period. In each period, the indoor air temperature can be obtained by Equation (9). The time required for the room temperature to reach a certain temperature value can be obtained by Equation (10).
Therefore, t 2 , t 3 , and t 4 can be solved by Equation (10). When the aggregated unit is in the steady state, the value N A of the number of air conditioning units at the rated power running condition is At any time t, the probability that the air conditioning unit i is in the rated power running state is p i (t) and the probability that the air conditioning unit i is in the suspended state is 1-pi(t). Therefore, the power of the aggregate is expected to satisfy Equation (11).
where p i (t) = If the aggregated unit exits the load reduction stage in the next unit of time, the setpoint temperature will be restored. In the heating mode, this can be divided into the following two situations: 1.
If ∆T ≥ 2 × deadband, then all aggregated units will operate at the rated power in the next unit of time, which means P = 1.

2.
If ∆T ≤ 2 × deadband and if the air conditioner is in the running state, the state will not change; if it is in the suspended state, the next state will change according to the indoor temperature.
The specific determination method is as follows. If T A > T set + ∆T + deadband, the pause operation state will occur in the next unit of time. If T A > T set + ∆T + deadband, the running state will occur in the next unit of time. In these two cases, the probability that the air conditioner is in the running state after increasing the setpoint temperature can be calculated using Equation (12). where In the case of 2000 air conditioning units, the expected instantaneous power of the aggregated unit when the ∆T setpoint temperature is increased is defined in Equation (13).
The relationship between the instantaneous power and the change in the setpoint temperature is shown in Figure 7.
If the aggregated unit exits the load reduction stage in the next unit of time, the setpoint temperature will be restored. In the heating mode, this can be divided into the following two situations: 1. If ∆T ≥ 2 × deadband, then all aggregated units will operate at the rated power in the next unit of time, which means P = 1. 2. If ∆T ≤ 2 × deadband and if the air conditioner is in the running state, the state will not change; if it is in the suspended state, the next state will change according to the indoor temperature. The specific determination method is as follows. If > + ∆T + deadband, the pause operation state will occur in the next unit of time. If ≤ + ∆T + deadband, the running state will occur in the next unit of time. In these two cases, the probability that the air conditioner is in the running state after increasing the setpoint temperature can be calculated using Equation (12).
In the case of 2000 air conditioning units, the expected instantaneous power of the aggregated unit when the ∆T setpoint temperature is increased is defined in Equation (13).
The relationship between the instantaneous power and the change in the setpoint temperature is shown in Figure 7. It is observed that in the heating mode, an increase in the setpoint temperature of the aggregate will cause an increase in the load. The range of values of the load limit of the virtual house model can be obtained using the relationship shown in Figure 7. Figure 8 shows that the load rebound of the aggregated unit contains the detailed models of 2,000 air conditioning systems created in Gridlab-D; an increase in the temperature settings at a given moment can be determined. The rebound peak load shown in Figure 8 is in agreement with the relationship shown in Figure 7. However, the rebound peak load does not occur at the moment of changing the setpoint temperature; there is a short period of increase (shown as the red dotted line in Figure 8  It is observed that in the heating mode, an increase in the setpoint temperature of the aggregate will cause an increase in the load. The range of values of the load limit of the virtual house model can be obtained using the relationship shown in Figure 7. Figure 8 shows that the load rebound of the aggregated unit contains the detailed models of 2000 air conditioning systems created in Gridlab-D; an increase in the temperature settings at a given moment can be determined. The rebound peak load shown in Figure 8 is in agreement with the relationship shown in Figure 7. However, the rebound peak load does not occur at the moment of changing the setpoint temperature; there is a short period of increase (shown as the red dotted line in Figure 8) because of the cycle time of the air conditioning compressor. This occurs because the compressor requires some time to return to the running state from the pause state. This time is called the cycle time.
Energies 2019, 12 FOR PEER REVIEW 10 conditioning compressor. This occurs because the compressor requires some time to return to the running state from the pause state. This time is called the cycle time. In this case, the initial temperature of the equivalent house can also be calculated using the law of energy conservation.

Load Control Scheme
The proposed load control scheme is divided into two stages, i.e., the load reduction stage and temperature recovery stage. The purpose of the load reduction stage is to achieve load control and the purpose of the temperature recovery stage is to restore the room temperature while suppressing the load rebound and oscillation.
In the load reduction stage, depending on the requirements of the load reduction, the average indoor temperature at the end of load reduction is first calculated by the aggregate model. Then the rate of decrease in the average temperature can be calculated. The rate of decrease will be broadcasted to all TCAs. Finally, according to the minimum control time (△t), the set temperature difference (△T) before and after can be calculated and implemented. Figure 9 shows the DR process of an air conditioning system. The green background represents the load reduction state. As shown in the figure, as the set temperature decreases, the time required for the temperature to increase in each heating cycle of the air conditioner decreases. In addition, it takes longer for the indoor temperature to decrease; therefore, the air conditioning power decreases at the load reduction stage. The load can be controlled by setting the value △T based on the aggregate model.  In this case, the initial temperature of the equivalent house can also be calculated using the law of energy conservation.

Load Control Scheme
The proposed load control scheme is divided into two stages, i.e., the load reduction stage and temperature recovery stage. The purpose of the load reduction stage is to achieve load control and the purpose of the temperature recovery stage is to restore the room temperature while suppressing the load rebound and oscillation.
In the load reduction stage, depending on the requirements of the load reduction, the average indoor temperature at the end of load reduction is first calculated by the aggregate model. Then the rate of decrease in the average temperature can be calculated. The rate of decrease will be broadcasted to all TCAs. Finally, according to the minimum control time ( t), the set temperature difference ( T) before and after can be calculated and implemented. Figure 9 shows the DR process of an air conditioning system. The green background represents the load reduction state. As shown in the figure, as the set temperature decreases, the time required for the temperature to increase in each heating cycle of the air conditioner decreases. In addition, it takes longer for the indoor temperature to decrease; therefore, the air conditioning power decreases at the load reduction stage. The load can be controlled by setting the value T based on the aggregate model.
Energies 2019, 12 FOR PEER REVIEW 10 conditioning compressor. This occurs because the compressor requires some time to return to the running state from the pause state. This time is called the cycle time. In this case, the initial temperature of the equivalent house can also be calculated using the law of energy conservation.

Load Control Scheme
The proposed load control scheme is divided into two stages, i.e., the load reduction stage and temperature recovery stage. The purpose of the load reduction stage is to achieve load control and the purpose of the temperature recovery stage is to restore the room temperature while suppressing the load rebound and oscillation.
In the load reduction stage, depending on the requirements of the load reduction, the average indoor temperature at the end of load reduction is first calculated by the aggregate model. Then the rate of decrease in the average temperature can be calculated. The rate of decrease will be broadcasted to all TCAs. Finally, according to the minimum control time (△t), the set temperature difference (△T) before and after can be calculated and implemented. Figure 9 shows the DR process of an air conditioning system. The green background represents the load reduction state. As shown in the figure, as the set temperature decreases, the time required for the temperature to increase in each heating cycle of the air conditioner decreases. In addition, it takes longer for the indoor temperature to decrease; therefore, the air conditioning power decreases at the load reduction stage. The load can be controlled by setting the value △T based on the aggregate model.  In the temperature recovery stage, the total load required to restore all the air conditioning systems to their original room temperature can be calculated by the aggregate model. Then, according to the rebound load limit value, the minimum time required (Tmin) to exit the DR can be obtained. Subsequently, the set temperature of the air conditioner can be calculated by Equation (14). ∆T is heterogeneous because the minimum control time ( t) is different in each air conditioner. As shown in Figure 9, the red background represents the temperature recovery state.
The control diagram is shown in Figure 10. We monitor the indoor temperature and total load and when the values exceed the alarm values, we repeat the process. The proposed DR scheme achieves load control using infrequent communications and few control signals, whereas the queuing method requires frequent communication with and control of each air conditioner; the higher the requirement for load control accuracy, the stricter the demand for communication and control is.
Energies 2019, 12 FOR PEER REVIEW 11 In the temperature recovery stage, the total load required to restore all the air conditioning systems to their original room temperature can be calculated by the aggregate model. Then, according to the rebound load limit value, the minimum time required (Tmin) to exit the DR can be obtained. Subsequently, the set temperature of the air conditioner can be calculated by Equation (14).
T ′ Δ is heterogeneous because the minimum control time (△t) is different in each air conditioner. As shown in Figure 9, the red background represents the temperature recovery state.
The control diagram is shown in Figure 10. We monitor the indoor temperature and total load and when the values exceed the alarm values, we repeat the process. The proposed DR scheme achieves load control using infrequent communications and few control signals, whereas the queuing method requires frequent communication with and control of each air conditioner; the higher the requirement for load control accuracy, the stricter the demand for communication and control is.

Results and Discussion
To verify the applicability of the proposed aggregate model and the effectiveness of the DR scheme, models of 2000 air conditioning units are created in Gridlab-D. These air conditioning units are in the heating mode when the outdoor temperature is set to 10 °C.
Three scenarios are simulated and discussed, namely the homogeneous model and the heterogeneous model. In the third model, each air conditioning unit has a different setpoint

Results and Discussion
To verify the applicability of the proposed aggregate model and the effectiveness of the DR scheme, models of 2000 air conditioning units are created in Gridlab-D. These air conditioning units are in the heating mode when the outdoor temperature is set to 10 • C.
Three scenarios are simulated and discussed, namely the homogeneous model and the heterogeneous model. In the third model, each air conditioning unit has a different setpoint temperature based on the heterogeneous model. The detailed models are created in Gridlab-D and the virtual house is created by the proposed aggregate method. The aggregated system receives a DR command that reduces the load by 50 percent from 18:30 to 19:00. The time period for exiting the DR is 19:00 to 20:00. The parameters of the virtual house are solved by the proposed aggregate model.
The specific parameters of the three scenarios are shown in Table 3. The first scenario represents the homogeneous model, in which the aggregated area of the houses is 185.8 m 2 . The second scenario represents the heterogeneous model with the same setpoint temperature; the aggregated area of the houses meets a uniform distribution and is 167-204 m 2 . The third scenario represents the heterogeneous model with the different setpoint temperatures; the aggregated area of the houses also meets the uniform distribution and is 167-204 m 2 . The initial setpoint temperature of the aggregated houses meets the uniform distribution at 23.3-24.4 • C. The average temperature is calculated for the virtual house model and broadcast to each air conditioning unit at a time interval of 5 min. The load-average temperature curve of the aggregated units under a DR scheme is shown in Figure 11; (a) represents the first scenario, (b) represents the second scenario, (c) represents the third scenario, and (d) represents the queuing method. Figure 12 shows the indoor temperature curves of 200 units randomly selected from the aggregated unit.  The specific parameters of the three scenarios are shown in Table 3. The first scenario represents the homogeneous model, in which the aggregated area of the houses is 185.8 m 2 . The second scenario represents the heterogeneous model with the same setpoint temperature; the aggregated area of the houses meets a uniform distribution and is 167-204 m 2 . The third scenario represents the heterogeneous model with the different setpoint temperatures; the aggregated area of the houses also meets the uniform distribution and is 167-204 m 2 . The initial setpoint temperature of the aggregated houses meets the uniform distribution at 23.3-24.4 °C. The average temperature is calculated for the virtual house model and broadcast to each air conditioning unit at a time interval of 5 min. The load-average temperature curve of the aggregated units under a DR scheme is shown in Figure 11; (a) represents the first scenario, (b) represents the second scenario, (c) represents the third scenario, and (d) represents the queuing method. Figure 12 shows the indoor temperature curves of 200 units randomly selected from the aggregated unit.
It can be seen in Figure 11a that the total load quickly decreased to 50% and the load was maintained for half an hour; the average temperature of the aggregate unit fell by more than 0.5 degree centigrade. Compared with Figure 11a, a load rebound is observed at 20:00 in Figure 11b because the proposed virtual house model has an error for heterogeneous houses. As can be seen from Figure 11c, different setpoint temperatures do not affect the load control but suppress the load rebound. In Figure 11d, the air conditioning units are divided into 10 groups for queuing and the performance of the load control accuracy or load oscillation suppression is lower than for the proposed DR scheme.  Figure 12 shows the load curve of the aggregated under the DR schemes. It is evident that the proposed DR scheme exhibits a good load control for the heterogeneous model and the homogeneous model. For a large number of heterogeneous models, the proposed aggregate model is applicable and achieves accurate load control and suppression of load rebound.
As shown in Figure 13, the proposed DR scheme not only reflects the equality of consumers by It can be seen in Figure 11a that the total load quickly decreased to 50% and the load was maintained for half an hour; the average temperature of the aggregate unit fell by more than 0.5 degree centigrade. Compared with Figure 11a, a load rebound is observed at 20:00 in Figure 11b because the proposed virtual house model has an error for heterogeneous houses. As can be seen from Figure 11c, different setpoint temperatures do not affect the load control but suppress the load rebound. In Figure 11d, the air conditioning units are divided into 10 groups for queuing and the performance of the load control accuracy or load oscillation suppression is lower than for the proposed DR scheme. Figure 12 shows the load curve of the aggregated under the DR schemes. It is evident that the proposed DR scheme exhibits a good load control for the heterogeneous model and the homogeneous model. For a large number of heterogeneous models, the proposed aggregate model is applicable and achieves accurate load control and suppression of load rebound.
As shown in Figure 13, the proposed DR scheme not only reflects the equality of consumers by being able to participate in the DR (unlike in Figure 13d, each consumer in Figure 13a, Figure 13b, and Figure 13c participates in the load control) but also retains the different temperature settings of the consumers (in Figure 13c, consumers with different temperature settings retain their temperature setting preferences).  Figure 12 shows the load curve of the aggregated under the DR schemes. It is evident that the proposed DR scheme exhibits a good load control for the heterogeneous model and the homogeneous model. For a large number of heterogeneous models, the proposed aggregate model is applicable and achieves accurate load control and suppression of load rebound.
As shown in Figure 13, the proposed DR scheme not only reflects the equality of consumers by being able to participate in the DR (unlike in Figure 13d, each consumer in Figure 13a, Figure 13b, and Figure 13c participates in the load control) but also retains the different temperature settings of the consumers (in Figure 13c, consumers with different temperature settings retain their temperature setting preferences).

Conclusions
In this paper, an aggregate model for air conditioning loads and their DR control strategy has been proposed. To validate the aggregate model, a detailed model of 2000 air conditioners was established in Gridlab-D to conduct comparisons. The simulation results showed that the