Aggregating Large-Scale Generalized Energy Storages to Participate in Energy Market and Regulation Market

This paper proposes a concept of generalized energy storage (GES) to facilitate the integration of large-scale heterogeneous flexible resources with electric/thermal energy storage capacity to participate in multiple markets. First, a generalized state variable referred to as degree of satisfaction (DoS) is defined, and dynamic models with a unified form are derived for different types of GESs. Second, a real-time market-based coordination framework is proposed to facilitate control, and ensure user privacy and device security. Demand curves of different GESs are then developed based on DoS to express their demand urgencies as well as flexibilities. Furthermore, a low-dimensional aggregate dynamic model of a GES cluster is derived thanks to the DoS-equality control feature provided by the design of demand curve. At last, an optimization model for a large-scale GESs to participate in both the energy market and regulation market is established based on the aggregate model. Simulations results demonstrate that the optimization algorithm could effectively reduce the total cost of an aggregator. Additionally, the proposed coordination method has high tracking accuracy and could well satisfy users' diversified power demand.


Introduction
The increasing integration of distributed energy resources (DERs) poses threat on stability and reliability of power system operation. With the development of smart grid technologies, the control of flexible loads, e.g., electric energy storage (EES), thermostatically control load (TCL) and electric vehicle (EV), has become a promising research field to address the problem brought by DERs [1,2]. Due to the great variety, large scale, wide distribution and small individual capacity of flexible loads, load aggregators (LAs) are required to aggregate and dispatch such loads and provide flexible services to power grid. The coordination strategies have been widely studied. For example, a three-step EV charging algorithm is presented in [3] to minimize charging cost. In [4], a priority-stack-based method is put forward for fixed-frequency air-conditioner (FFA) fleet to track automatic generation control signal. Literature [5] establishes physical model of inverter air-conditioners (IVAs) and proposes a hierarchical control framework. A unified state model for EVs and TCLs is developed in [2]. In [1], a generalized battery model is presented to describe the operational characteristics of commercial/residential building HVAC and energy storage system, based on which optimal control methods are proposed to provide various services.
However, most of these methods have the following disadvantages. First, most approaches can only be applied to a certain type of load [3][4][5]. An LA should provide different interfaces to integrate different types of resources, which makes it difficult for the LA to fully utilize flexibility of various flexible loads and increase the control cost. Second, many literatures use centralized control method based on direct load control [1,2,4]. An LA should collect detailed parameters of all controlled loads and specify the response power of each load, which has heavy computational burden and arXiv:1901.03941v1 [cs.SY] 13 Jan 2019 communication traffic, and may lead to privacy issues and device security problems. Such methods may not be suitable for coordination of large-scale flexible loads. Some papers develop unified models to coordinate different types of resources [1,2], while they depend on centralized control as they mainly focus on individual behaviour rather than aggregate performance.
The objective of LAs is to provide flexible services to power grid. The method in [4] could provide reliable regulation services. Literature [2] demonstrates its effectiveness in power fluctuation smoothing. In [3,5], LAs participate in optimal dispatch and achieve benefits from energy markets. However, the above researches only consider a single market. It is pointed in [6] that since different services have different requirements, participating in multiple markets helps to better utilize the control flexibility of resources and obtain more benefits. To stack multiple services, energy storage systems are controlled in [7] to simultaneously participate in N-1 contingency requirement, voltage management and frequency regulation. An optimal control method is presented in [8] for energy storages to provide grid services including energy arbitrage, balancing service, capacity value and distribution system deferral, and outage mitigation. [9] and [10] introduce a Markov decision process model and a stochastic control method respectively to co-optimize battery storage for energy arbitrage and frequency regulation. However, these papers only consider energy storages. Although literature [1] considers coordination of various flexible loads to provide multiple grid services, it has to solve some centralized optimization problems, which may lead to a high computational cost.
To address the above problems, a unified modelling method and coordination strategy for generalized energy storages (GESs) is proposed. It has unified information interface and control method, and requires a relatively low communication and computation cost, which is helpful for an LA to conduct coordination control over large-scale GESs. The contributions of our work are threefold: 1. Dynamic models with a unified form are developed for heterogenous GESs based on a generalized state variable referred to as degree of satisfaction (DoS). 2. A real-time coordination framework based on the market equilibrium mechanism is presented to allocate aggregate power to individual GESs. General demand curves are constructed under the framework to achieve equal degree of satisfaction and meet users' diversified requirements. 3. A low-dimensional aggregate dynamic model for a GES cluster is derived and used in an optimization model for an LA to participate in both the energy market and regulation market.
The rest of this paper is organized as follows. Section 2 introduces the dynamic model of different GESs. In section 3, the market-based real-time coordination framework is proposed and demand curve construction methods for GESs are presented. Section 4 derives the aggregate dynamic model for a GES cluster. In section 5, the optimization problem that considers both energy and regulation markets is given. Section 6 shows the simulation results which demonstrate the effectiveness of our method. Finally, section 7 summarizes our contribution and future work.

Dynamic Models for GESs
This paper focuses on four typical types of GESs, i.e., EES, EV, inverter air-conditioner (IVA) and fixed-frequency air-conditioner (FFA) since they account for a large share of flexible resources in the demand side. These resources are called GES for that they all can store energy, i.e., electric energy or cold/thermal energy, thus their power consumption could be adjusted without affecting the user satisfaction. Meanwhile, they have similar dynamic characteristics, which makes it possible to establish a unified physical model for them. For example, there have been a few literatures that develop battery modelling method for the thermal energy storage loads. However, most researches only focus on FFAs [11][12][13]. A battery-type reduced-order model for building HVAC system is proposed in [14], but it only considers the dynamics of individual load. Literature [5] studies the battery modelling for individual IVA and aggregated IVA cluster, while the aggregating strategy would lead to the result that the IVA cluster's regulation ability could not be fully utilized. This paper will put forward a unified dynamic model for different types of GESs, and study both individual model (introduced in this section) and aggregate model (detailed in section.4).

Degree of Satisfaction, DoS
Before establishing the dynamic models, a dimensionless state variable referred to as degree of satisfaction (DoS) is defined for GESs, with the following purposes: 1. It could be used to measure the user satisfaction. The range of DoS is set to [-1,1], and the closer DoS is to 0, the higher user satisfaction is. 2. DoS could reflect a GES's state of energy: DoS equalling 0 indicates that the stored energy is at the expected level, while DoS close to ±1 means the stored energy is near the allowed range. 3. Since a GES can deviate from its ideal state (DoS=0) to provide services, the DoS can be used to quantify its current flexility, i.e., a DoS value close to 0 implies a high flexility reservation. 4. Finally, as DoS is a generalized index, it can be used to establish a unified dynamic model for various GESs.

Electric Energy Storage (EES)
Ignoring the charge/discharge efficiency, the dynamic model of an EES is given by where E EES i,k and P EES i,k denote the electric energy and power of EES i at time k, respectively (P EES i,k > 0 when EES is charging); ∆t is the control cycle.
Variable S is hereafter used to denote DoS. Definition of EES's DoS is given by Substituting Eq.(2) and SOC EES = E EES /C EES into Eq.(1), the following dynamic model can be derived: where C EES i denotes the nominal capacity of EES i.

Electric Vehicle (EV)
The physical model of an EV is where E EV i,k and P EV i,k denote the electric energy and power of EV i at time k, respectively; η EV i,k is the charge efficiency.
According to the current charger technology, we focus on the prevailing EV type, which operates at two discrete states: idle and charging with a fixed rate [15,16]. Let E tar,i denote the target energy of EV i at user-specified departure time t dep,i , and E in,i denote the energy at the time EV i is connected into power grid t in,i . The average power required to charge EV i to E tar,i can be calculated by If EV i charges at P req,i , the expected energy profile is An EV can provide its flexibility by deviating from E exp . Referring to the hysteretic model in [17], DoS of EV is defined as where C EV i denotes the nominal capcity of EV i; r i denotes the energy deadband, which limits the error between E tar,i and the actual energy at t dep,i within ±C EV i × r i %. Combined with Eq.(5)-(7), Eq.(4) can be transformed into

Inverter Air-conditioner (IVA)
Without loss of generality, cooling air-conditioners are studied in this paper. The thermal dynamic process is modelled by a first-order differential equation [4,5] aṡ where T a,i,k and T o,i,k denote the indoor air temperature and outdoor temperature at time k, respectively; Q IVA i,k is the heat rate of IVA i; the thermal parameter a i = 1/(R th,i C th,i ).
The analytical solution of Eq.(9) in recursive form is The electrical model of an IVA adopts the simplified linear model in [5], which is given by where P IVA i denotes the electric power of IVA i; f i is the operation frequency of the compressor; p 1 , p 2 , q 1 and q 2 are coefficients.
DoS of an IVA is defined as where T set,i denotes the setpoint; T dev,i denotes the allowed temperature deviation. Combined with Eq.(11) and Eq.(12), Eq.(10) can be transformed into where α i = e −a i ∆t , β i = q 1 /(a i C th,i p 1 ), γ i = (p 1 q 2 − p 2 q 1 )/(a i C th,i p 1 ).

Fixed-Frequency Air-conditioner (FFA)
The thermal model of an FFA can also be established by Eq.(9). Its electrical model is given by where COP i is the coefficient of performance. FFA is ON/OFF controlled load, so its electric power P FFA i equals to nominal power P FFA N,i when it is ON and P FFA i = 0 when it is OFF.
FFA's definition of DoS is identical to that of IVA, which is given by Eq. (12). Combined with Eq. (12) and Eq. (14), Eq.(10) can be transformed to where β i = R th,i × COP i .

Unified Dynamic Model
The dynamic models of the above GESs can now be represented in a unified form: where m i,1 , m i,2 , m i,3 are coefficients of the unified dynamic model. Mappings between DoS and the original variables of different GESs are summarized in Fig.1. This paper would study a real-time coordination method of GESs, which could make DoS of individual GESs (or average DoS of a cluster of GESs) approximately equal. Such control characteristic is referred to as DoS-equality control in this paper and will be detailed in the next section.

DoS-equality Control Based on Market Equilibrium Mechanism
Market equilibrium mechanism [18] is introduced in this paper to coordinate large-scale GESs. Main stages of the coordination method are: (1) Bidding: Each GES expresses its urgency and flexibility by constructing a demand curve. The demand curve is denoted by d i (λ) in this paper, which is a non-increasing function.
(2) Aggregating and clearing: LA collects demand curves from all GESs and forms the aggregate where N is the number of controlled GESs; Assume the aggregate target power is P tar , then LA can calculate the clearing price by λ * = D −1 (P tar ).
(3) Disaggregating: LA broadcasts λ * to all GESs. Each GES responds to λ * locally according to its demand curve. The response power of GES i can be obtained by P res,i = d i (λ * ).
Following the above steps, the aggregate target power can be allocated among the GESs, thus realizing accurate power tracking. It should be noted that: (1) The target power P tar depends on applications. In this paper, it will be determined by an optimization problem to be discussed in the next section; (2) The clearing price λ * is only a control signal rather than a price signal. It is dimensionless and its range is set to [-1,1]. Therefore, it is called "virtual price" in this paper.
The proposed market-based coordination method highlights the following advantages: 1. It improves the autonomy of the GES. Each GES can convert its private information, e.g., user preferences, current adjustable range and security constraints, into demand curve. Since demand curves of all GESs have a unified form, it can shield the differences among various GESs and effectively protect user privacy. Besides, the LA does not have permission to directly control GES, which improves device security. 2. It simplifies the control of the LA. An LA does not need to specify each GES's type and is able to coordinates various GESs via an identical signal, i.e., the virtual price signal λ * , which significantly reduces control complexity and requirement of communication bandwidth.
In addition, the proposed method can realize the DoS-equality control to obtain the following advantages: 1. GESs could have same degree of user satisfaction regardless of the resource type or capacity, which ensures control fairness. In addition, since DoS reflects a GES's state of energy, the DoS-equality control could avoid some GESs going beyond their adjustable range prematurely, thus better utilizing the regulation ability of a GES cluster. 2. The unique DoS of a GES cluster can be a state variable to derive an aggregate dynamic model, making it possible to treat the whole GES cluster as a virtual storage, which will be detailed in section 4.
According to the operation characteristic, the GES can be further classified into two types: GES operating at continuous power (CP-GES) and GES operating at discrete power with discrete states (DP-GES). A CP-GES, e.g., an EES or an IVA, is able to keep its DoS at the desired value by adjusting its operating power, while for a DP-GES, e.g., an EV, an FFA or an electric heater, its DoS generally fluctuates within the allowed range. Besides, the state switching frequency of a DP-GES should generally be controlled to prolong the device's lifetime.
The key of the proposed coordination method is the construction of demand curves for different GESs, which will be introduced in the following subsections.

Demand Curve
Let S i denote the current DoS of GES i. The following construction principle of demand curve is established in this paper: where P CONST,i denotes the electric power required to maintain the current S i over a control cycle. Eq. (17) is used to realize the DoS-equality control feature. For explanation, the clearing price λ * is assumed to be constant. When λ * < S i , the response power P res,i = d i (λ * ) is higher than P CONST,i , leading to the decrease of S i ; when λ * > S i , P res,i is lower than P CONST,i , leading to the increase of S i ; when λ * = S i , P res,i equals P CONST,i , keeping S i unchanged. Therefore, S i values of all GESs following the principle in Eq.(17) will approach the control signal λ * .
Based on Eq. (17), the demand curve of CP-GES is constructed as shown in Fig.2. Note again that the price herein is virtual and is limited between -1 and 1. The demand curve consists of 5 key points. The anchor point A(P CONST , S) could satisfy the condition in Eq. (17). Points B(P SAT,MI N , 1) and C(P SAT,MAX , −1) are used to keep the DoS within limits when responding to any clearing price λ * . Thus, P SAT,MI N and P SAT,MAX are minimum and maximum power that would not make DoS go beyond limit over a certain period of time t p (which is set 5min in this paper), without considering operation constraints. Points D and E which lies on line AB and line AC respectively are further introduced to guarantee the response power would not exceed operational constraints. Thus, P OPT,MI N and P OPT,MAX are minimum and maximum power that a GES could operate at in current control cycle. The calculation method of the above characteristic power will be detailed in the next section.

Characteristic Power
The EES and IVA are typical CP-GESs. They can both adopt the demand curve shown in Fig.2. An EES's characteristic power can be calculated by where P EES N denotes the nominal charging/discharging power; η EES char and η EES disc are charge/discharge efficiency, respectively; E EES min and E EES max are the minimum and maximum allowed energy, respectively; E EES t is the stored energy at time t. IVA's characteristic power can be calculated by where P IVA min and P IVA max denotes the minimum and maximum power; the function g IVA P (T tar , t p ) obtains the electric power that makes the indoor air temperature change from current value to the target temperature T tar over a period of time t p . The derivation of g IVA P (T tar , t p ) is detailed in Appendix A.

Demand Curve
A DP-GES typically has two states, i.e., ON and OFF states. The proposed demand curve for a DP-GES is illustrated in Fig.3(a), where P ON denotes the operating power when the GES is ON, and S is transformed from its DoS value by offsetting and then normalizing: 8 of 20 To explain the principle, the aggregate demand curve of a cluster of DP-GESs is shown in Fig.3(b), which is obtained by sorting DP-GESs in descending order of S . As can be seen, the proposed bidding strategy can divide DP-GESs into two groups according to their operation states, i.e., an ON group in the upper half-plane and an OFF group in the lower half-plane, and achieve the following purposes: (1) For DP-GESs in the same group, a DP-GES's S reflects its power consumption priority. The higher S is, the higher probability to maintain or switch to ON state is, and vice versa.
(2) A DP-GES in the ON group always has a higher S than that in the OFF group, which gives high priority for DP-GESs to maintain their current states, thus avoiding frequent switching. To explain the change rule of a single DP-GES's DoS, the clearing price λ * is assumed to be constant over a period of time. Note that the state change of a DP-GES may be triggered either by its S value, or by S to ensure comfort. The trajectory of S and S at two different λ * , i.e., λ * 1 > 0 and λ * 2 < 0, are illustrated in Fig.4(a) and Fig.4(b), and the following laws can be found: when λ * > 0, S ranges in [-1+2λ * ,1]; when λ * < 0, S ranges in [-1,1+2λ * ], which means that the DoS value fluctuates within a symmetrical range around λ * . Therefore, if DoS values of DP-GESs are assumed to be uniformly distributed in such range, the average DoS of the cluster equals λ * .

Characteristic Power
EV and FFA are typical DP-GESs. They can both adopt the demand curve in Fig.3(a). An EV's characteristic power is where P EV N denotes the EV's nominal charging power. An FFA's characteristic power is

Locked State
In certain control cycles, a GES may get into the locked state, which means it should maintain its current operating state and power. There are two typical situations: (1) To reduce mechanical wear and protect the device [19,20], DP-GESs, such as FFA and EV, should satisfy lockout time constraints before they switch state; (2) Restricted by the device capability, response cycle of some GESs may be longer than the real-time control cycle (which is 10s in this paper).
The lockout mechanism can be easily realized in this paper, as a GES can simply submit the following demand curve during the lockout time: which means the response power of GES i maintains its current operating power P i,t for any λ.

Aggregate Dynamic Model of a GES Cluster
One significant advantage of our method is that the aggregate dynamic model of a large-scale GES cluster can be easily derived thanks to the DoS-equality control feature, as well as the unified dynamic model for different GESs defined in Eq. (16).
For CP-GESs, add their dynamic models together: where Ω c represents the set of CP-GESs. Under the DoS-equality control, DoS of all CP-GESs becomes equal. Let S agg−c denote DoS of the CP-GES cluster, and P agg−c denote the aggregate power. The aggregate dynamic model of CP-GESs can then be derived as For DP-GESs, to facilitate analysis, we assume coefficients of their dynamic models to be equal first. Add their dynamic models together, and we obtain ∑ i∈Ω d where Ω d represents the set of DP-GESs. The instantaneous DoS value of each DP-GES is a random variable. Denote the average DoS of DP-GESs as S agg−d , and denote the aggregate power as P agg−d . Then Eq.(26) can be written as: where |Ω d | denotes the number of DP-GESs. Under the DoS-equality control, S agg−c of CP-GESs and S agg−d of DP-GESs tend to be equal (both equal λ * ), thus can both be denoted by S agg . Therefore, Eq.(25) and Eq.(27) can be added together: where N = |Ω c | + |Ω d | denotes the total number of GESs; P agg,k is the aggregate power. Eq.(28) can be represented in a compact form as P agg,k = M 1,k S agg,k+1 + M 2,k S agg,k + M 3,k .
The above derivation is based on the assumption that all DP-GESs have equal model coefficients. For those with different coefficients, we can first divide them into different groups according to their coefficients, then aggregate each group using Eq. (27), and finally form the aggregate dynamic model of all the GESs (including CP-GESs and DP-GESs) by Eq. (28).
The aggregate dynamic model of heterogeneous GESs proposed above has the advantage that the LA could obtain the aggregate model easily by adding model coefficients of all GESs with no need to identify GES's type, thus having a low computational cost. Furthermore, from the control aspect, the low-dimensional aggregate model greatly reduces the complexity of the optimization problem, which will be discussed in the next section.

Optimal Multi-Market Flexibility Allocation
An LA can aggregate flexibility of large-scale GESs to provide multiple services to power grid. For example, it can schedule an optimal power consumption profile according to the electricity price of energy market [1,5,21], as the profile P sch in Fig.5. LA can provide other ancillary services at the same time to obtain higher benefits [1,21], e.g., responding to the regulation signal P reg in Fig.5.
The regulation signal can be decomposed into the low frequency part denoted by regA, and the high frequency part denoted by regD [22]. Literature [23] analyses the regulation signal of a certain power grid, and finds that the high frequency part could account for up to 30%. In this paper, the LA responds to the regD signal considering the following facts: First, the regD signal has zero-mean over a period of time [22], which could significantly reduce requirements for the capacity of GESs; Second, the extra energy introduced by the regD signal is close to 0, thus having little impact on electricity bills; Third, if regulation payments are determined by the performance-based policy used in PJM regulation market [24], the LA could obtain high benefits.
Considering both the energy and regulation markets, the target power of an LA is given by where P sch is the hourly scheduled power, which determines the bill paid to the energy market; regD is the regulation signal normalized to [0,1]; C reg denotes the contracted regulation capacity. Let the optimization cycle be 1 hour, then the LA solves the following convex optimization problem in the nth cycle to allocate its flexibility to the two markets: (ω score µ cap,k C reg,k capacity payments +ω score µ mile,kωmile C reg,k mileage payments ) + f s (S agg,k ) s.t.P sch,k = M 1,k S agg,k+1 + M 2,k S agg,k + M 3,k , ∀k − 1 ≤ S agg,k ≤ 1, ∀k P sch,k + C reg,k ≤ P max agg,k , ∀k P sch,k − C reg,k ≥ P min agg,k , ∀k C reg,k ≥ 0, ∀k, where µ ele,k , µ cap,k , µ mile,k denote respectively the electricity price, regulation capacity price and regulation mileage price in the kth cycle;ω score denotes the statistical value of the regulation performance score defined by PJM [24], andω mile denotes the statistical value of the regulation mileage [24]; P max agg,k and P min agg,k are the maximum and minimum power of the GES cluster, which are calculated in every optimization cycle as P max agg = ∑ N i=1 P max,i and P min agg = ∑ N i=1 P min,i , where Note that EESs are able to discharge, thus the LA may sell electricity to power grid at this time. The purchase price and sale price are assumed to be equal in this paper. The first constraint in Eq.(31) is the aggregate dynamic model of the GES cluster, the second constraint ensures the user satisfaction, and the last three ones constrain the regulation capacity. Thanks to the established aggregate dynamic model, the scale of the GES cluster does not affect the computational complexity of the optimization problem.
To improve user satisfaction, Eq.(31) includes a penalty term f s (S agg,k ), which is defined as where ω scale is a proportionality coefficient, which is assigned 0.1 here; µ avg ele is the daily average electricity price.

Three-Layer Control Structure
To sum up, this paper develops four models, i.e., a unified dynamic model and a unified demand model for individual GES, as well as an aggregate dynamic model and an optimization model for an LA to approximate dynamics of a GES cluster and participate in multi-markets. These four models are organized in a a three-layer control structure, as illustrated in Fig.6. (1) Rolling Optimization Layer At the beginning of current optimization cycle n, the LA gathers model coefficients m i,1 ∼ m i, 3 from GESs and update the aggregate dynamic model according to Eq. (28). Meanwhile, current DoS values of all GESs are collected to update the initial recursive value of S agg,n bŷ The optimization problem in Eq.(31) can then be solved, and the optimal scheduled power sequence [P sch,n , ..., P sch,24 ] as well as the regulation capacity sequence [C reg,n , ..., C reg,24 ] are obtained. Implement only the first elements of these optimal sequences, i.e., P sch,n and C reg,n in the current optimization cycle, and repeat the above steps each hour. It's well known that this idea of rolling optimization comes from the model predictive control (MPC) [19], which is adopted here to update the aggregate dynamic model iteratively, and to consider constraints in future time slots explicitly.
(2) Real-time Coordination Layer In each control cycle (10s in this paper), the LA receives the regulation signal regD from the control center, and then calculates the real-time target power P tar according to Eq.(30). The virtual market is cleared according to section 3, and the clearing price λ * is broadcast to each GES.
(3) GES Autonomy Layer In each optimization cycle (1h), each GES updates its DoS, model coefficients in Eq. (16) and power constraints in Eq.(32), and then reports these information to the LA. In each control cycle (10s), each GES reports its flexibility and responds to the clearing price both through the demand curve.
It is worth mentioning that, since the LA interacts with different GESs through a unified set of information, i.e., DoS, model coefficients, power constraints, demand curve and virtual price, the method in this paper supports a flexible tree-like structure. For example, a local concentrator can be deployed in an community to pre-aggregate the information. Therefore, the method has high scalability and is suitable for wide-area coordination of large-scale GESs.

Simulation Settings
The simulation cases are based on a residential community system. The simulation lasts for 24h. Electricity price µ ele , regulation capacity price µ cap and mileage price µ mile adopt the data in literature [5] and [25], and their profile are illustrated in Fig.7. The regD signal uses the PJM data on Jul 13th, 2016. According to the statistical analysis on the regD signal in 2016 [26], we setω mile = 2.7. According to the simulation results based on historical data, we conservatively assignω score = 0.92. Four types of GESs are considered in this paper, i.e., EES, EV, FFA and IVA, whose parameters are shown in Table.1, where U(a,b) indicates a uniform distribution between [a,b] and t res denotes the response cycle of GES. The profile of the outdoor temperature in simulation cases is shown in Fig.A1. To evaluate the control effect, the LA has to estimate the baseline power of the GES cluster when all GESs are uncontrolled, which will be denoted by P base hereafter. Many papers have studied the estimation method of the baseline load, e.g., a statistical based method proposed in [27]. In this paper, since the aggregate model of the GES cluster is available to the LA, it can estimate P base in each hour by solving the following optimization problem: min P base ,S agg 24 ∑ k=n S 2 agg,k s.t.P base,k = M 1,k S agg,k+1 + M 2,k S agg,k + M 3,k , ∀k − 1 ≤ S agg,k ≤ 1, ∀k P agg,min,k ≤ P base,k ≤ P agg,max,k , ∀k.
(35) Solution of Eq.(35) is referred to as the baseline case in this paper.

Case 1: only participate in the energy market
To evaluate the effectiveness of the aggregate dynamic model and the DoS-equality control, an optimization problem simplified from Eq. (31) is solved which only considers the energy market: s.t.P sch,k = M 1,k S agg,k+1 + M 2,k S agg,k + M 3,k , ∀k − 1 ≤ S agg,k ≤ 1, ∀k P agg,min,k ≤ P sch,k ≤ P agg,max,k , ∀k.
LA coordinates GESs to make the aggregate power P agg track the scheduled power P sch . As illustrated in Fig.8, the proposed method has high tracking accuracy. In addition, the aggregate power P agg can vary up and down around the baseline power P base as the electricity price changes to reduce energy cost. Therefore, the GES cluster can be scheduled as virtual energy storage, since it can be charged by making P agg higher than P base and discharged by making P agg lower than P base . Figure 8. Aggregate power and electricity price in case 1. Fig.9 demonstrates the performance of the DoS-equality control. As shown in Fig.9(a), DoS of all CP-GESs, i.e., IVAs and EESs, can track λ * quite well as they are able to adjust their power continuously. For DP-GESs (including FFAs and EVs), note again that the average DoS of a cluster rather than an individual's DoS can follow λ * . It can be seen in Fig.9(b) that S FFA avg of FFAs tracks λ * well. However, for the EVs, S EV avg slightly fluctuates around λ * . This is because the number of EVs is small, making the statistical characteristics inconspicuous and the distribution of DoS not well aligned with the analysis in section 3.3. Therefore, a larger-scale DP-GES cluster yields better DoS-equality control effect.   10 shows how well the aggregate dynamic model defined in Eq.(29) fits the cluster with heterogeneous GESs. As can be seen, S avg of all GESs except EVs are very close to the aggregate state variable S agg at the end of every optimization cycle. Since the number of EVs is small, especially during 6:00-9:00 and 18:00-20:00 when some EVs are off-grid, S EV avg fluctuates around S agg with relatively large errors.
In addition to the small number of DP-GESs, some other factors may also lead to errors in the aggregate dynamic model, e.g., the charge/discharge efficiency is not considered in EES's dynamic model. In order to prevent the error from being accumulated, this paper adopts the rolling optimization method to mitigate impacts of these factors continually. As can be observed from Fig.10, when combined with the rolling optimization, the simple aggregate model in (29) could be a useful tool for an LA to capture the aggregate dynamic feature of large-scale GESs.

Case 2: participate in both energy and regulation markets
In case 2, the GES cluster participates in both the energy and regulation markets. By solving the optimization problem in Eq.(31), the scheduled power P sch and regulation capacity C reg can be obtained, as illustrated in Fig.11. For comparison purposes, the scheduled power in case 1 is also plotted in the figure, which is denoted by P sch1 . The scheduled power profile in case 2 has significant difference from that in case 1 because the LA should allocate the cluster's flexibility to two markets according to both electricity price (µ ele ) and regulation price (µ cap and µ mile ). Since a symmetric regulation signal is used in this paper, C reg can reach its maximum value when P sch ≈ (P agg,min + P agg,max )/2. It can be seen that when µ cap is relatively high, e.g., during 14:00-15:00, 19:00-21:00, the LA tends to maximize C reg to gain higher payments from regulation market. The target power P tar in this case is calculated by Eq.(30), and the tracking performance is shown in Fig.12. The hourly value of ω score and ω mile are shown in Fig.A2. According to the result, ω score can basically reach 0.95 under the proposed control framework. In comparison, when responding to regD signal, ω score of a hydroelectric generator can be 0.7∼0.8, while that of an electric energy storage can be higher than 0.9 [28]. Therefore, the simulation results demonstrate that the GESs discussed in this paper are very promising alternatives to provide fast and accurate frequency regulation services. Trajectory of DoS and clearing price λ * are illustrated in Fig.13. Compared with Fig.10, since the GES cluster also needs to respond to the rapidly changing regD signal, some fluctuations and sudden changes can be observed in λ * , which makes DoS of GESs unable to exactly follow λ * . However, DoS always tends to approach λ * , and thus the DoS-equality control is basically realized. In addition, at the end of some optimization cycles, e.g., at 10:00, 16:00, 24:00, the difference between S avg and S agg is a little larger than that in Fig.10. It is because regD signal is not exactly zero-mean, which affects the actual energy consumption in each optimization cycle and enlarges the difference. To analyse the economic benefits of the proposed method, the energy bill, regulation payments and total cost between the above cases (case 1 and case 2) and the baseline case are compared. Hourly calculated cost is shown in Fig.A3, and Table.2 lists the daily results. As can be observed, case 1 which only considers energy market can significantly reduce the energy bill compared with the baseline case. By optimally allocating flexility in two markets, case 2 further reduce the total cost compared to case 1, even though it receives a higher energy bill in the energy market. By providing fast regulation service, an LA obtains high payments from the regulation market, leading to a significant reduction in the total cost. Therefore, our method could achieve great economic benefits.

Response Performance of Individual GESs
The proposed control strategy ensures DoS-equality among GESs. While for different types of GESs, their response behaviour may be different due to their distinguishing features. To demonstrate this, pick one GES randomly from each type of GESs, and plot their response power in Fig.14. For CP-GESs, this paper allows different response cycles. For example, an EES adjusts its response power each 10s (t res = 10s), while an IVA adjusts its power every 1min (t res = 60s) considering its response ability, thus an IVA's response power is stair-shaped as shown in Fig.14(a).
DP-GESs adjust their response power by regulating the duty cycle. Among them, an EV's duty cycle at different time is basically similar, because its operation is irrelevant to external conditions and mainly depends on the user's charging pattern. In contrast, the required power of an FFA varies over time as it is significantly affected by environment conditions, e.g., outdoor temperature. For example, it can be observed from Fig.14(b) that the duty cycle of an FFA during 13:00-16:00 is higher than that in the rest of the day as more cooling energy is required during these time slots.
In addition, an EV should ensure that the electric energy reaches its target value E tar at departure time. The change of electric energy is shown in Fig.15, where E + = E exp + r%C EV , E − = E exp − r%C EV . As we can see, the hysteretic model and control strategy used in this paper can guarantee that the difference between the energy at departure time and the target energy E tar would not exceed ±r%C EV .

Conclusions
In this paper, a unified coordination method is developed for large-scale heterogeneous GESs to participate in both energy and regulation market.
A generalized state variable referred to as DoS is first defined for GESs. The dynamic models with a unified form are then developed for different GESs. In real-time control, a market-based coordination framework is adopted, and a DoS-equality control method is then developed by construction of generalized demand curves for both GESs operating at continuous power and GESs with discrete states. Based on the unified dynamic models and the DoS-equality control feature, a low-dimensional aggregate dynamic model for a GES cluster is derived. At last, an optimization model aiming to allocate the flexibility of a GES cluster into both the energy market and the regulation market is developed, which uses the aggregate model to significantly reduce the mathematical complexity of the optimization problem.
The control framework has unified uplink/downlink information interfaces and supports a tree-like structure in both real-time coordination stage and optimization control stage, which makes it flexible, scalable and suitable for the control of large-scale GESs. Simulation results demonstrate that the aggregate model well describes the dynamic behaviour of a GES cluster. Additionally, the real-time control method can track the target power accurately while satisfying diversified requirements of different GESs and ensuring control fairness. It is also shown that an LA could gain considerable energy bill savings and high payments by participating in both energy and regulation markets.
However, simulations in this paper are based on ideal communication system and perfect model parameter identification. Future work would study the robustness of our method under communication problems and model errors.

Conflicts of Interest:
The authors declare no conflict of interest

Abbreviations
The following abbreviations are used in this manuscript: Appendix A Derivation of the function g IVA P (T tar , t p ) The analytical solution of Eq.9 can be described as: where T a,i (0) and T o,i (0) denote the current indoor air temperature and outdoor temperature of IVA i; T a,i (t) denotes the indoor air temperature at time t.
To derive the electric power required to make T a,i change from T a,i (0) to T a,i (t) over a certain period of time t p , we let t = t d , T a,i (t p ) = T tar , and the required heat rateQ IVA i can be calculated bŷ According to Eq.(11), the required electric powerP IVA i can be derived bŷ