A Day-Ahead Scheduling of Large-Scale Thermostatically Controlled Loads Model Considering Second-Order Equivalent Thermal Parameters Model

Among various kinds of flexible loads resources in power system, the thermostatically controlled loads (TCLs) are considered as the excellent candidates for the demand response (DR) programs due to the heat/cool storage characteristics. Effective operation and coordination of TCLs in DR programs requires an accurate model to capture the aggregate power. One appealing approach is the equivalent energy storage (EES) method, which describes the change of aggregate power offered by the TCLs as the dynamic of an energy storage model with limits on the output power and the energy capacity. In this paper, an EES model with the second-order equivalent thermal parameters (second-order ETP) model is established for interpretation of TCLs dynamic characteristics. Considering the internal and external heat exchange of the TCLs, the relationships between the heat exchange power and energy storage of the TCLs EES model are redefined. In order to exert the actual potential of TCLs in day-ahead scheduling, the time-varying charging-discharging power and energy storage constraints of the TCLs EES model are developed. By this way, more feasible scheduling results can be obtained. The performance of the proposed method is verified by the testing results.


A. SETS AND INDICES
Maximum energy storage of aggregated TCLs P e,agg,max , P e,agg,min Upper and lower limit of power of aggregated TCLs a j , b j , c j Fuel cost coefficients of j th thermal power unit P Gmax,j , P Gmin,j Upper and lower limit of power generation of j th thermal power unit c SU,j , c SD,j Start-up and shut-down price of j th thermal power unit SU j , SD j Upper and lower limit of the ramp rate of j th thermal power unit c WC Price of wind power curtailment c DRup , c DRdown Price of DR with increased /decreased power consumption PL max,l Maximum power flow of the l th branch QG l,j , QW l,w , QD l,d Element of branches-generators incidence matrix, the brancheswind farm incidence matrix and the branches-load incidence matrix Outdoor temperature of i th TCL at time t T a,i (t) Indoor air temperature of i th TCL at time t T m,i (t) Indoor mass temperature of i th TCL at time t Q i (t) Cooling power of i th TCL at time t P ex_ao,i (t) Heat exchange power between indoor air and outdoor of i th TCL at time t P ex_am,i (t) Heat exchange power between indoor air and mass of i th TCL at time t E a,i (t) Energy storage of air of i th TCL at time t E m,i (t) Energy storage of indoor mass of i th TCL at time t E a,max,i Maximum energy storage of i th TCL P c,agg (t) Charging and discharging power of aggregated TCLs at time t P e,agg (t) Electric power of aggregated TCLs at time t P a,agg (t) Average power of aggregated TCLs at time t P c_a,agg (t) Charging and discharging power of indoor air at time t P c_m,agg (t) Charging and discharging power of indoor mass at time t P ex_ao,agg (t) Heat exchange power between indoor air and outdoor of aggregated TCLs at time t P ex_am,agg (t) Heat exchange power between indoor air and mass of aggregated TCLs at time t E a,agg (t), E m,agg (t) Energy storage of indoor air and mass of aggregated TCLs at time t E a,agg (t), Energy storage variation of E m,agg (t) indoor air and mass at time t E a,ave (t), E m,ave (t) Average energy storage of indoor air and mass of aggregated TCLs at time t T o,ave (t) Average outdoor temperature of aggregated TCLs at time t T max,ave (t) Average maximum indoor temperature of aggregated TCLs at time t F Total cost of power generation P G,j (t) Power output of j th thermal unit at time t C G,j (t) Fuel cost of j th thermal power unit at time t u G,j (t) Operating state of j th thermal power unit at time t x SU,j (t), x SD,j (t) Binary variable of start-up and shut-down of j th thermal power unit at time t C SU,j (t), C SD,j (t) Start-up and shut-down costs of j th thermal power unit at time t P W,w (t) Wind power output of w th wind power farm at time t P WC,w (t) Wind With the gradual large-scale construction of new energy power systems, many works have been reported on the energy systems scheduling in which different energy sources and storage options complemented by responsive electrical/thermal loads are taken into account [1]- [7], such as the optimal scheduling of the gas-electric integrated energy system (IES), the regional IES with electric vehicles (EVs) and air conditioning loads, the multicarrier energy supplies with a biogas-solar-wind hybrid renewable system, the IES with wind power and multi-type energy storage, and so on. It shall be noted that the anti-peak-regulation characteristic of new energy power is the main issue of these generation sources despite their various advantages. And energy storage systems and demand response (DR) programs can be regarded as tools to overcome the interruptible operation and volatility of new energy power, and provide energy owners flexibility to follow the day-ahead scheduling. The trend of further development and utilization of the demand side resources is getting increasing attention in power systems. The electric loads are encouraged to actively participate in the DR programs to achieve specific objectives such as cost reduction (CR), peak load shaving (PLS), and congestion management (CM). Many studies have been reported on the residential and commercial loads participating in the DR programs [8]. Among various kinds of flexible loads resources, the thermostatically controlled loads (TCLs) are considered as the excellent candidates for the DR programs due to the heat/cool storage characteristics [9], [10].
Due to the complex thermodynamic characteristic of TCLs, the potential of TCLs in DR programs is difficult to be estimated accurately, and the large number of aggregated TCLs also lead to obstacles in the power system operation. To make the DR programs with TCLs meet the requirements of day-head scheduling of the power systems, a variety of scheduling strategies have been proposed in Ref. [11]-Ref. [17]. The Markov transition matrix (MTM) is adopted to model the large-scale TCLs to accomplish PLS in Ref. [11]. Reference [12] studies the coordinated scheduling strategy for TCLs and generation units. The strategy in Ref. [13] could work well in meeting the increased or decreased scheduling requirements by narrowing down or shifting the comfort band of all TCLs, then this strategy may underestimate the potential of TCLs and decrease the advantage of TCLs participating in the scheduling. In Ref. [14], a decomposed scheduling framework for large-scale TCLs is proposed, which groups the TCLs based on the similarity of the TCL model parameters. Reference [15] establishes an aggregated TCLs model considering the heterogeneity. Reference [16] focuses on the potential of building-level TCLs, and develops the day-ahead scheduling strategy of a building aggregator. Besides, Reference [17] develops the human thermal comfort estimation module to satisfy the comfort of the customers in the scheduling.
Recently, the heat/cool storage characteristic of TCLs is received increasing attention. The storage of heat/cool can be equivalent to the energy storage, and the heating/cooling power can be equivalent to charging and discharging power. Further, the aggregated TCLs can consequently be viewed as an equivalent energy storage (EES) model [18]. Several EES models for TCLs are proposed in Ref. [19]-Ref. [23]. For getting the compatibility of the air conditionings (ACs) with the scheduling models, Ref. [19] establishes an equivalent thermal battery model to facilitate its participation in scheduling. Focusing on the heterogeneity parameters of TCLs and no short-cycling requirement, Ref. [20] proposes an improved EES model. The power of the EES models in Ref. [21] and Ref. [22] are calculated by the average power, resulting in fixed upper and lower limits of the power, and fail to reflect the time-varying characteristic of the power. Then, the time-varying heat exchange power is adopted to calculate the power in Ref. [23].
Notwithstanding the contributions of the above works, there exist limitations which need to be improved.
• The TCLs model that the EES constraints are not taken into account may be prone to overload of power [11]- [17].
• The average power EES model fails to reflect the time-varying characteristic of the power [21], 22].
• The deviation of first-order ETP model is hard to be eliminated, which may reduce the accuracy of scheduling with large-scale TCLs [23]. To fill these gaps, this paper introduces the second-order ETP model into the TCLs EES modeling. The heat exchange power is decoupled into the internal and external parts. The relationships between two parts of the heat exchange power and energy storage are developed, and the constraints of charging and discharging power, energy storage are redefined, so that the proposed EES model can make TCLs more accurately scheduled in the power system day-ahead scheduling, and beneficial to operation of the power system.
The rest of this paper is structured as follows: The day-ahead scheduling model considering DR with TCLs is formulated in Section II. The individual TCL model considering the second-order ETP model is described in Section III. Section IV introduces the proposed EES model considering the second-order ETP model. The case studies are developed in Section V. Finally, the conclusions are summarized in Section VI.

II. SCHEDULING MODEL CONSIDERING DR WITH TCLS
More flexibility of scheduling model is required for the power systems operation considering the increasingly penetration of renewable energy (e.g., wind power). However, thermal power units cannot always work at preferred operating points and follow loads in real-time, due to their operating economics and thermal efficiency issues. As the alternative resource for flexibility, the TCLs are considered for their distinguished capability to improve system operating efficiency and enhance system economic benefit. To obtain the optimal economy, the minimum cost of power generation is set as the goal to establish scheduling model [24].

A. OBJECTIVE FUNCTION CONSIDERING TCLs
The minimum cost of power generation is set as the objective function to establish the day-ahead scheduling model. Based on the data of wind power and load forecasting, the day-ahead scheduling model considering large amounts of TCLs can be established.
where F is the total cost of power generation, including operation cost, wind curtailment cost and compensation cost of TCLs. C G.j (t) is the fuel cost of thermal power units and P G.j (t) is the output power of j th thermal power unit at time t. a j , b j , c j are fuel cost coefficients of units. C SU,j (t) and C SD,j (t) are the costs of start-up and shut-down. x SU,j (t) and x SD,j (t) represent the station of on or off of j th thermal power unit at time t. c SU,j and c SD,j are the prices of start-up and shut-down. C WC,w (t) is the cost of wind power curtailment and P WC,w (t) is the curtailment of w th wind power farm at time t. c WC is the price of wind power curtailment. C DR,d (t) is the compensation cost of TCLs at time t. P DRup,d (t) and P DRdown,d (t) are the DR with increased power consumption and decreased power consumption of the d th load bus at time t, respectively.
where P G.j (t) is the output power of j th thermal power unit at time t. P W,w (t) and P WC,w (t) are the forecasted wind power and wind power curtailment of w th wind power farm at time t. L d (t) is the load of d th bus at time t. P DR,d (t) is the scheduled power of TCLs of d th bus at time t. If P DR,d (t) is positive, it means the load is shifted from time t to other time; if it is negative, it means the load is shifted from other time to time t.
2) START-UP AND SHUT-DOWN VARIABLE CONSTRAINT OF THERMAL POWER UNITS u G,j (t) = 1, if j th thermal power unitis on 0, if j th thermal power unit is off (11) where u G,j (t) is the state of j th thermal power unit at time t.

3) ON-OFF CONSTRAINT OF THERMAL POWER UNITS
where P Gmax,j and P Gmin,j are the upper and lower limits of the output power of j th thermal power unit, respectively.

5) RAMP RATE CONSTRAINT OF THERMAL POWER UNITS
where SU j and −SD j are the upper and lower limits of the ramp rate of j th thermal power unit, respectively.
where PL l (t) is the power flow of the l th branch at time t. QG l,j , QW l,w , QD l,d are the element of branches-generators incidence matrix, the element of branches-wind farm incidence matrix, and the element of branches-load incidence matrix, respectively. L d (t) is the power of the d th load bus at time t. PL max,l is the maximum value of the l th branch power flow. Based on the formulas above, the conventional day-ahead scheduling model considering DR with large amounts of TCLs is established. In order to make the scheduling more accurate, the TCLs EES model based on the second-order ETP models is introduced into the day-ahead scheduling model in this paper.

III. MODELING OF INDIVIDUAL TCL BASED ON SECOND-ORDER ETP MODEL
The individual TCL model is the basis for developing an aggregated TCLs model. Most of previous approaches adopt first-order differential equations for individual TCL ETP model shown in Fig. 1. In the first-order ETP model, the indoor mass temperature dynamics are ignored, which increases the uncertainty and deviation of modeling. For an aggregated TCLs model, large differences may be caused by the small deviation of individual TCL after multiple iterations. In practical applications, building materials and furnishing have a large heat capacity, and both the indoor air and mass temperature dynamics need to be considered at the same time. In this paper, the second-order ETP model shown in Fig. 2 is adopted to describe the thermal dynamics of each individual TCL [25].
According to Ref. [25] and Ref. [26], the second-order ETP model of individual TCL can be expressed as following: where T a (t) and T m (t) are the indoor air temperature and indoor mass temperature of TCL at time t, respectively. T o (t) is the outdoor temperature at time t. Subjected to the temperature set-point T set changing from 24.5 • C to 27.5 • C, the evolution of the T a (t) and T m (t) of a TCL is shown in Fig. 3. In Fig. 3, a much longer time is taken for T a (t) to increase from 26.5 • C to 27.5 • C during the time period [t 1 t 2 ] than during the period [t 3 t 4 ]. The increase of the T a (t) is slowed down significantly, due to the lower T m (t) during [t 1 t 2 ]. Neglecting T m (t), the time derivative of T a (t) would only depend on itself, and then the different duty cycles could not be observed. However, the entire trajectory of T a (t) shown in Fig. 3 cannot be accurately described by any first-order ETP model. This indicates that the critical importance of T m (t), especially for accurately capturing the transient response [26].
The second-order ETP model can be transformed into the state space form as follows: whereẋ(t) is the state vector for the TCL at time t, and the matrices A and B can be easily derived from equation (18) and equation (19). s(t) is the switching state of each individual TCL, which is typically regulated by a simple hysteretic controller based on a temperature dead-band, and s(t) can be defined as following: The TCL turns on when the indoor air temperature T a (t) reaches the upper boundary T set + δ/2, and the TCL turns off at the temperature T set − δ/2. The δ is the dead-band size.
There are heating and cooling operation modes for TCLs. However, only the cooling mode is considered in this paper. The relationship between the cooling power Q(t) and the electric power P(t) of a TCL can be formulated as: where η is the energy conversion efficiency of the TCL. Based on the equations (18)∼(28), the second-order ETP model of individual TCL is established.

IV. MODELING OF EQUIVALENT ENERGY STORAGE BASED ON AGGREGATED TCLS
According to the heat/cool storage characteristics of TCLs, the heating/cooling power can be equivalent to charging and discharging power of aggregated TCLs, and the storage of heat/cool can be equivalent to the energy storage. In this VOLUME 8, 2020 section, the energy storage models of aggregated TCLs are described by heterogeneous second-order ETP models.

A. HETEROGENEITY AND UNCERTAINTY OF TCL PARAMETERS
In order to characterize the aggregate dynamics of large amounts of TCLs, the easiest method is to model each individual TCL and then the aggregated charging and discharging power can be superimposed by individual TCL.
Due to the different of building structures, device types, and outdoor temperature, the parameters and the power consumption of a TCL population are heterogeneous [15]. Thus, the parameters of TCLs based on the second-order ETP models can be defined as and Q i (t) have time-varying characteristic. In this paper, the heterogeneity and uncertainty of the TCLs' parameters are assumed to be characterized by the normal distribution N (µ, χ 2 ), in which µ is mean and χ is standard deviation [27].

B. EQUIVALENT ENERGY STORAGE MODELS OF HETEROGENEOUS TCLs
In the day-ahead scheduling model, the charging and discharging power P c,agg (t) of TCLs EES model approximates to the scheduled power of TCLs P DR (t) [21]- [23], and the relationship between P DR (t) and P c,agg (t) can be expressed as following: If P c,agg (t) is positive, the EES model is charging, and DR is load-increasing; if P c,agg (t) is negative, the EES model is discharging, and DR is load-decreasing.
In Ref. [21], Ref. [22], P c,agg (t) is calculated directly by the average power P a,agg (t). P c,agg (t) = P e,agg (t) − P a,agg (t) (30) where the power of aggregated TCLs P e,agg (t) is calculated by the average cooling power Q ave (t).
The steady-state operation of TCLs is periodically on /off switching. The duty cycle of TCLs can be expressed by the on-state duration τ on and the off-state duration τ off . Thus, the average power P a,agg (t) is calculated by the duty cycle of steady-state operation as following: Obviously, according to equation (34), P a,agg (t) is a fixed value, and the time-varying power is ignored. Considering the accuracy of the model, Ref. [23] indicates that the indoor air temperature T a (t) changes with time, and the average power changes with the T a (t). The average power P a,agg (t) is replaced by the heat exchange power P ex,agg (t) to calculate the charging and discharging power P c,agg (t). P c,agg (t) = P e,agg (t) − P ex,agg (t) However, the heat exchange power P ex,agg (t) in equation (35) is only defined as a function of T a (t) in Ref. [23].
Combined with the second-order ETP model of individual TCL, indoor mass also has the heat exchange power. Thus, in order to make the scheduling more accurate, the indoor mass heat exchange power of heterogeneous aggregated TCLs is considered in this paper. The heat exchange power of the aggregated TCLs EES model P ex,agg (t) is decoupled into the heat exchange power between indoor air and outdoor P ex_ao,agg (t) and the heat exchange power between indoor air and mass P ex_am,agg (t). P ex_ao,agg (t) is the heat exchange power between the EES model and the outside. P ex_am,agg (t) is the heat exchange power inside the EES model. The relationship of the heat exchange power between indoor air, mass and outdoor is shown in Fig. 4. Then, the charging and discharging power P c,agg (t) can be expressed as the follows: • The charging and discharging power of the aggregated TCLs EES model is P c,agg (t): P c,agg (t) = P e,agg (t) − P ex_ao,agg (t) (36) • The charging and discharging power of indoor air is P c_a,agg (t): P c_a,agg (t) = P e,agg (t) − P ex_ao,agg (t) − P ex_am,agg (t) • The charging and discharging power of indoor mass is P c_m,agg (t): Further, from the micro level, based on the second-order ETP model, the individual EES model can be established.
where P ex_ao,i (t) is the heat exchange power from indoor air to outdoor of individual TCL. P ex_am,i (t) is the heat exchange power from indoor air to mass of individual TCL.
In the cooling mode, the energy at T max,i (t) is taken as the minimum of energy storage. For individual TCL, the energy storage of indoor air E a,i (t) and the energy storage of indoor mass E m,i (t) can be expressed as follows: By eliminating the variables T a,i (t) and T m,i (t) according to equations (43)∼(44), the relationships between P ex_ao,i (t), P ex_am,i (t) and E a,i (t), E m,i (t) are obtained as follows: Based on the heat exchange power of individual TCL, the heat exchange power P ex_ao,agg (t) and P ex_am,agg (t) of TCLs can be transformed into a equivalent form as follows: • Heat exchange power between indoor air and outdoor P ex_ao,agg (t).
• Heat exchange power between indoor air and mass P ex_am,agg (t).
where E a,ave (t), E m,ave (t), C a,ave , R a,ave , C m,ave , R m,ave , T o,ave (t), T max,ave (t) and η ave are the average values, respectively. E a,agg (t) is the energy storage of indoor air of aggregated TCLs at time t. E m,agg (t) is the energy storage of indoor mass of aggregated TCLs at time t.
TCLs have the characteristics of thermal energy storage, and the heat/cool variation of large amounts of TCLs during t can be expressed through the energy storage variation of air E a,agg (t) and the energy storage variation of indoor mass E m,agg (t), respectively. In the cooling mode, E a,agg (t) and E m,agg (t) can be expressed as follows: = P c_m,agg (t)dt ≈ P ex_am,agg (t) · t (50) According to equations (47)∼(48) and equations (49)∼(50), the recursion formulas between E a,agg (t+1) and E a,agg (t), E m,agg (t+1) and E m,agg (t) can be obtained.

C. CONSTRAINTS OF EQUIVALENT ENERGY STORAGE MODEL IN DAY-AHEAD SCHEDULING
In the day-ahead scheduling model with the TCLs EES model, the constraints include not only the conventional constraints of the scheduling model (7)∼(17), but also the relevant constraints of the TCLs EES model.

1) CHARGING AND DISCHARGING POWER CONSTRAINT OF TCLs
Each TCL may be in different states, on or off, so P e,agg (t) of TCLs has the upper and lower limits. P e,agg,min ≤ P e,agg (t) ≤ P e,agg,max where P e,agg,max is the upper limit of the total electric power and it can be obtained by setting all the TCLs on. Similarly, the lower limit P e,agg,min is obtained by setting all the TCLs off and it is usually zero. Further, combined with equation (37), the charging and discharging power constraint of TCLs participating in the scheduling can be expressed as following: P c_a,agg (t) ≥ −P ex_ao,agg (t)−P ex_am,agg (t) P c_a,agg (t) ≤ P e,agg,max −P ex_ao,agg (t)−P ex_am,agg (t) (54)

2) ENERGY STORAGE CONSTRAINT OF TCLs
In the cooling mode, the energy at T max,i (t) is taken as the minimum of energy storage. According to on equation (43) and equation (44), the maximum energy storage of individual TCL is the maximum energy storage of indoor air E a,max,i , and can be expressed as following: Then the maximum energy storage of large amounts of TCLs E a,agg,max can be obtained: According to equation (47) and equation (48), the heat exchange power of TCLs P ex,ao,agg (t) and P ex,am,agg (t) have the relationships with E a,agg (t) and E m,agg (t), respectively. Thus, the equality constraint between them can be obtained as equation (58) and equation (59).

4) TIME-COUPLING CONSTRAINT OF E a,agg (t ) AND E m,agg (t ).
According to equations (51)∼(52), the energy storage variation of indoor air E a,agg (t) and indoor mass E m,agg (t) have time-coupling characteristics. Then, combined with equations (47)∼(48), the time-coupling constraint of E a,agg (t) and E m,agg (t) can be obtained as equation (60) and equation (61).
The constraints used in the optimization problem are as follows: • Conventional constraints of the scheduling model, including the equations (7)∼(17).
• Constraints of the TCLs energy storage model, including the equations (53)∼(61). Along with the power system scheduling model introduced in Section II, the day-ahead scheduling model of large amounts of heterogeneous TCLs EES based on second-order ETP model is established.

V. CASE STUDIES
In the following cases, the proposed modeling approach of aggregated TCLs EES is validated against simulations on the 6-bus system and the 118-bus system [28]. In the two system, 50 thousand and 1 million TCLs working in cooling mode are considered, respectively. The forecasted load and wind power of summer are shown in Fig. 5.
The nominal values of model parameters referred from [23], [26], and [29] are listed in Table 1. Considering the heterogeneous parameter-distribution of massive TCLs, the parameters are set to random normal distribution N (µ, χ 2 ) with 0.1 relative standard deviation (RSD). The initial temperatures and operating states of TCLs is assumed to be stable in scheduling cases. The nonlinear problems of day-ahead scheduling are transformed into linear problems by quadratic  programming, and Matlab software with Matpower, Cplex, Yalmip toolboxes are used for solving the optimization problem.
To verify the effectiveness of the proposed modeling approach, different scheduling models are compared. The performance of them are separated into three parts for analysis, as follows: A. Performance of day-ahead scheduling. B. Response to load following signal. C. Verification on the 118-bus system.

A. PERFORMANCE OF DAY-AHEAD SCHEDULING
For comparison, four scheduling models are developed as follows: • Model 1 (Traditional Model): Only the charging and discharging power constraint of TCLs is considered and it is calculated by equation (30), in which the average power P a,agg (t) is calculated by equation (34).
• Model 3 (EESM_ETP_1 st ): The EES model of TCLs scheduling is established by the first-order ETP models [23]. In this model, the charging and discharging power P c,agg (t) is calculated through the heat exchange power P ex,agg (t), and P ex,agg (t) is only considers the heat exchange between indoor air and outdoor.
• Model 4 (EESM_ETP_2 nd ): The EES model of TCLs scheduling is established by the second-order ETP models. In this model, the charging and discharging power P c,agg (t) is calculated through two parts of heat exchange power P ex,ao,agg (t) (between indoor air and outdoor) and P ex,am,agg (t) (between indoor air and mass), according to equations (36)∼(38) and equations (47)∼(48). The EES of TCLs involved the two parts are E a,agg (t) and E m,agg (t), which follow the recursion formulas (51)∼(52). And the charging and discharging power constraint is redefined by equation (54). Besides, the ideal maximum energy storage of TCLs E a,agg,max is 120MWh calculated by equation (56). However, actually the indoor air temperature is hard to reach to the maximum or the minimum, so when considering the upper or lower limit of the energy storage, the energy storage constraint margin (M AI ) of ideal energy storage limits is set aside. In the simulations of this paper, M AI = 5%.
Then, based on the forecasted load and wind power, the performance of each scheduling model (Model 1∼4) in day-ahead scheduling of 6-bus system is tested. The scheduling results of each model are shown in Fig. 6∼Fig. 9.
From the performance of each model, the following observations can be obtained: 1) If the EES constraints are ignored, TCLs may not be sufficient to meet scheduling requirements.
From Fig. 6 (c), it can be seen that, although the charging and discharging power P c,agg (t) is constrained within the limits, the energy storage value E agg (t) exceeds the actual limits of energy storage, and the TCLs may not be sufficient to meet scheduling requirements.
2) The P c,agg (t) calculated by P a,agg (t) cannot reflect its time-varying characteristic.
From Fig. 7 (c), Fig. 8 (c) and Fig. 9 (c), it can be seen that, the Model 2, Model 3 and Model 4 all consider energy storage constraints in the day-ahead scheduling, and the energy storage value E agg (t) are all constrained within the limits. However, the upper and lower limits of E agg (t) in Model 2 are fixed values and do not reflect the time-varying characteristics of P c,agg (t).
3) The scheduled power of EESM_ETP_1 st and EESM_ETP_2 nd is more accurate.
From Fig. 6 (d), Fig. 7 (d), Fig. 8 (d) and Fig. 9 (d), it can be seen that, the scheduled power P DR,scheduling (t) of Model 1 and Model 2 obviously overestimate the response capacity of TCLs, and cannot reflect the charging and discharging characteristics of TCLs EES. From the power and distribution of DR, the scheduled power of Model 3 and Model 4 is more accurate.
4) Charging and discharging power constraint of EESM_ETP_2 nd is tighter and more effective.
From Fig. 8 and Fig. 9, it can be seen that, although the values of P c,agg (t) in Model 3 and Model 4 are almost the same (they all related to P ex,agg (t) between indoor air and outdoor), the limits of P c,agg (t) in Model 4 is closer to the maximum and minimum values of the P c,agg (t). And the time-varying characteristic of the limits of P c,agg (t) in Model 4 is more obvious than Model 3, indicating that P c,agg (t) is affected by the heat exchange between indoor air and mass.

B. RESPONSE TO LOAD FOLLOWING SIGNALS
In order to verify effectiveness of the proposed scheduling model (Model 4 EESM_ETP_2 nd ), load tracking control with PI controller is adopted, so that the actual aggregated power of TCLs is as close as possible to the scheduled power. Corresponding to the four scheduling models above, four control models are established for analysis and comparison. All the four control models consider the TCLs with second-order ETP models, and each individual TCL is modeled according to equations (20)∼ (26).
In addition, for comparison and analysis, one dynamic evaluation index of TCLs energy storage is introduced. The State-of-Charge (SOC) is an important evaluation index for energy storage [30]. The SOC of aggregated TCLs EES model at time t is defined as following: where E agg,control (t) is the energy storage of aggregated TCLs obtained by the tracking control.  Then, based on the load tracking control with PI controller, the response of load following signals of each control model are compared. The related results are shown in Fig. 10∼Fig. 13.
From the performance of each load tracking control model, the following observations can be obtained: 1) If the EES constraints are ignored, aggregated TCLs may be far from meeting the scheduling results. Compared with Fig. 11∼Fig. 13, the load tracking performance in Fig. 10 is the poorest. From Fig. 10, during 5h∼13h and 20h∼24h, the scheduling results of Model 1 (the EES is ignored) are hardly tracked by the load tracking control with aggregated TCLs, and the SOC quickly approaches the limits. The poor performance of load tracking control illustrates that most TCLs have lost the ability to participate in DR, and the actual potential of TCLs cannot meet the scheduling results.
2) The time-varying characteristic of P c,agg (t) is important for estimating the potential of TCLs.
From Fig. 11, during 7h∼13h, the load tracking control results obviously deviates from the scheduling results of Model 2 (the time-varying characteristic of P c,agg (t) is ignored), and the SOC is close to 100%. The P c,agg (t) is calculated with a fixed value P a,agg (t) in Model 2, and the potential of TCLs in the scheduling may be overestimated or underestimated in some periods.
3) The scheduling results of EESM_ETP_2 nd better match the actual capacity of aggregated TCLs.
Compared with Fig. 10∼Fig. 12, the load tracking performance in Fig. 13 is better. Although the performance load tracking of the Model 3 (both the EES constraints and the time-varying heat exchange power are taken into account) is a little better than Model 2, it still fails to fully tracked. The scheduling results of Model 4 can be successfully tracked in all periods, and the SOC is always controlled within the limits. It shows that the scheduling results based on the proposed model have higher accuracy.
Then, in order to clearly show the accuracy of the load tracking control results, the integrated square error (ISE) is adopted, and as shown in equation (63).

ISE =
T 0 P DR,control (t) − P DR,scheduling (t) 2 dt (63) where P DR,control (t) is the power obtained by the tracking control. P DR,scheduling (t) is the scheduled power of DR by day-ahead scheduling.   The results are shown in Table 2, from which we can see that ISE of Model 4 (EESM_ETP_2 nd ) is the smallest and far smaller than the other three models. This indicates that it is important to introduce the proposed EES model into the scheduling model with TCLs, which can improve the accuracy of the model.
Further, in order to compare the influence of the energy storage constraint margin (M AI ) in the scheduling model on  the load tracking control performance, ISE is also adopted to show the differences, as shown in Fig. 14.
From Fig. 14, when the M AI is smaller than 5%, ISE tends to be larger because the scheduling results are not tracked well. Besides, another reason is that the indoor air temperature cannot reach the maximum or the minimum and so that the energy storage cannot reach 100% or 0% too. When the M AI is higher than 5%, ISE tends to be smaller. However, at this case the potential of TCLs is not fully utilized because the energy storage is too constrained. Although the load tracking control with higher M AI performs well, the participation of TCLs in the scheduling will decrease a lot, resulting in decreasing the advantage of TCLs participating in the scheduling.

C. VERIFICATION ON THE 118-BUS SYSTEM
To further verify the feasibility of the proposed model, Model 4 (EESM_ETP_2 nd ) is scheduled on the 118-bus system. One million TCLs under the cooling mode are distributed proportionally to each load bus. The scheduling results and the load tracking control results are shown in Fig. 15 and Fig. 16, respectively. From Fig. 15 and Fig. 16, the aggregated TCLs established by the proposed model (Model 4 EESM_ETP_2 nd ) are scheduled successfully on the 118-bus system, and the accuracy of the scheduling results are verified through the load tracking control.
Moreover, the computational cost (processing time and memory requirement) of the methods are shown in Table 3. It can be seen that the computational cost of the proposed method is a little higher than other methods, due to the additional order of the ETP model. The processing time and the memory requirement is acceptable and competitive to other methods. Considering practical applications, the proposed method has low computational cost.

VI. CONCLUSION
In this paper, the second-order ETP model is taken into account in equivalent energy storage (EES) model and day-ahead scheduling model with TCLs, in order to make the scheduling more accurate and feasible. The main contributions of this paper can be summarized as follows: 1) An EES model of TCLs considering the second-order ETP model is established. Compared with Ref. [23], the proposed model can more accurately describe the actual dynamic of the TCLs time-varying heat exchange power.
2) The second-order ETP model of TCLs is taken into account in the day-ahead scheduling model, based on which the relationships between two part of the heat exchange power (internal and external of the TCLs EES model) and energy storage are developed, and the constraints of charging and discharging power, energy storage are redefined, so that more accurate and feasible scheduling results can be obtained. Through load tracking control, the scheduling results of the proposed model (EESM_ETP_2 nd ) can be accurately tracked. And the ISE of the proposed model is 10.25 (MW) 2 · h, which is far smaller than others. It shows that the scheduling results have higher accuracy and better feasibility when the second-order ETP model is considered.
Future work will focus on the accurate cost modeling of aggregated TCLs EES model with the second-order ETP model, and introduce it to the real-time retail electricity markets, in order to further taping the potential of TCLs in real-time conditions.