Reference-Tracking Control Policies for Packetized Coordination of Heterogeneous DER Populations

This manuscript presents design and analysis of a set of reference-tracking control policies for large-scale coordination of distributed energy resources (DERs) and quantifies tracking errors that arise due to heterogeneity in the power ratings for a fleet of DERs. In particular, the relay-based, reference-tracking control strategy that underpins much of packetized energy management (PEM) is augmented to uniquely leverage PEM’s energy packet request mechanism to optimize the number of accepted requests and to explicitly consider the quality of service (QoS). In addition, tracking errors from modeling a heterogeneous fleet of packetized DERs with a group of homogeneous macromodels are analytically derived for relevant PEM information scenarios. Finally, simulation-based analysis validates the results and shows that PEM is suitable for providing load balancing and ramping services for the grid.


I. INTRODUCTION
B ALANCING demand and supply is a fundamental problem in power system operations. For decades, the operating paradigm has been supply follows demand, which has been implemented via the hierarchical primary, secondary, and tertiary frequency regulation schemes. However, with increasing penetrations of variable renewable generation, such as solar PV and wind, and internet-enabled, connected, and controllable appliances (e.g., distributed energy resources or DERs), it is now technically and economically feasible for flexible demand to follow a variable supply (i.e., provide dynamic grid services) [1]- [3]. However, to control large fleets of DERs to provide grid services requires a distributed controller implementation, which has been known since the late 1970s [4], [5]. More recently, the notions of demand dispatch [6] and controllable loads [7] have reinvigorated the field of controlling DERs. In fact, it has become evident that large-scale coordination of DERs will require randomization Manuscript  to avoid harmful effects of synchronization [8]- [11]. Thus, this article focuses on so-called distributed, randomized control policies for demand dispatch that embed control logic both at the DER (device) and at the coordinator, which is also known as aggregator or virtual power plant (VPP). Information is then exchanged between the DERs and the coordinator, which informs the device to probabilistically transition its operational state (e.g., consume, standby, and/or supply) and permits the coordinator to compute control inputs regulating the controllable net demand (total consumption minus variable supply). However, coordinating large fleets of DERs is challenging and requires careful distributed controller design that is cognizant of one or more of the following: 1) the devices' finite energy limits (e.g., temperature bounds for thermostatically controlled loads or TCLs); 2) ON/OFF cycling rates; 3) communication costs; 4) fleet heterogeneity within and between device classes; and 5) privacy concerns around data sharing. It is within this context that a set of reference-tracking controllers are presented here for a new demand dispatch method called packetized energy management (PEM). PEM can coordinate diverse DER populations within a single ("coupled") controller to track an aggregate power regulation signal, while guaranteeing local DER-specific QoS constraints, and is described as follows.
1) A DER estimates its local need for energy, e.g., energy state of charge or SOC. 2) If the SOC is within a predefined range of comfort, the DER probabilistically requests from the coordinator to consume energy at a fixed rate (e.g., 4 kW) for a pre-specified epoch (e.g., 5 min), which begets an energy packet (e.g., 0.33 kWh). If the SOC falls below a pre-defined lower threshold the DER automatically and temporarily opts out of PEM to guarantee QoS and reverts to a conventional control mode (e.g., charges) until the SOC is returned within limits and after which it returns to PEM operation. Similarly, the DER also opts out if the SOC exceeds a pre-specified higher threshold.
3) The asynchronous coordinator or VPP either accepts or denies the DER's packet request on an arrival basis in relation to the current reference tracking error. If the request is denied, go to i ); else, consume the energy packet and then go to i ). The packetizing concepts that inspired PEM focused on coordinated EV charging under a fixed transformer constraint [12], [13]. Recently, PEM has been adapted for coordination of water heaters [14] together with batteries [15] and a This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ macro-level (aggregate) model has been characterized qualitatively and quantitatively in [16]- [18]. Most recently, Duffaut Espinosa and Almassalkhi [19] provides the most complete description of the macromodel and QoS requirements. In this article, the focus is entirely on design of control policies for the PEM aggregator that dynamically optimizes the number of accepted requests and the QoS under different information scenarios. However, PEM is not the first article to explore DER control and coordination at scale.
Methods based on mean-field models and control provided an early approach for coordinating large-scale population of DERs, including TCLs [20], [21]. The seminal work in [10] designed a feedback linearization scheme around a state-bin transition model where the aggregate power tracking error input together with estimated bin states were used to compute an updated state-bin transition probability vector that was broadcast to all devices. The feedback linearization overcomes a state-input bi-linearity in the controller, which requires that an accurate estimate of the states is available to compute the control input. The work herein also includes bi-linear input-state term, however, the control signal for PEM is a scalar signal that represents the proportion of available devices to transition (i.e., communicated only to a small sub-set of available devices). Mathieu's work has been advanced to consider higher-order TCL models [22], multiple market objectives with a virtual energy storage model [23], and communication delays [24].
The approach in [11] also considers a randomized control policy for a homogeneous population of DERs, but avoids a state-bin transition model by instead using a mean-field model based on a linearization around the average device's pre-determined baseline transition matrix. The linearization enables the use of a single scalar that can be broadcast to update the shape of the local transition probabilities. Recently, the group has been active in extending the work to consider estimation of the QoS (e.g., a proxy for state of charge and on-off device cycling rates) with and without local opt-out control [25], [26], intra-class DER heterogeneity via filtering [27], and virtual energy storage models [28]. Related work on stochastic control and modeling of heterogeneous TCLs can be found in [29] where the control strategy combines stochastic transitions with temperature dead-band (hysteresis) control to provide primary and secondary frequency regulation services for the grid. While the aggregate (macro) model in our work is a state-bin transition model similar to that in Mathieu, PEM's approach replaces the controllable probabilityto-transition curve at the flexible load with a fixed probabilityto-request curve that engenders a stochastic request process input to the PEM coordinator. The control strategy then is to determine the (scalar) proportion of requests to accept, which can be implemented as a probabilistic response to requesting devices or via recent Internet-of-Things (IoT) methods as a binary (accept/deny) response to each device. PEM has used the latter binary approach in earlier work on PEVs and TCLs [13], [14].
PEM's device-driven request-response mechanism was inspired by early TCP/IP and ALOHA network protocols. Similar approaches can be found in [30], which presents a decentralized packetized coordination scheme that is analyzed via queueing theory to determine QoS under a fixed, known resource budget. The hierarchical controller proposed in [31] is similar in spirit to PEM, but considers a non-randomized, priority-based scheme wherein a device regularly transmits a transition fitness score (e.g., in [0, 1]) instead of a request, which the coordinator uses to rank and prioritize the individual devices for control. In addition (unlike PEM), the controller uses the ranking to form on-to-off and off-to-on transition priority queues from which they can construct aggregate up/down flexibility curves that are used to transmit and update device-level set-points to ensure a rapid autonomous frequency response. This device-driven scheme is well suited for smaller populations due to the device-by-device prioritization and regular information exchanges. Another distributed DER control scheme replaces the randomized device-level transition with a deterministic, cyclical decision process for set-point control of TCLs [32]. Other related works include [33] where the effects of dynamic constraints and modeling errors are quantified.
Thus, the article herein first proposes control architectures for diverse fleets of PEM enabled DERs by leveraging the author's prior work on a PEM macromodel [16]- [19]. These architectures allow the coordinator to either generate a common control signal for different DER-types with the objective to track a market reference signal, or to obtain a distinct control signal for each type of DER that simultaneously tracks the aforementioned reference as well as satisfies a secondary objective such as QoS guarantees, minimize number of accepted requests etc. Then, the synchronous nature of the macromodel is analytically characterized and compared with the asynchronous request-response mechanism in PEM. Finally, a group-based approach to handle DER intra-class heterogeneity is provided as well as bounds on tracking error. The specific contributions of this article are as follows: 1) Two different DER coordination architectures are proposed and validated for coordinating diverse fleets of DERs. In the first architecture, denoted the coupled architecture, each packet request identifies as either a charge request or a discharge request only, that is the VPP is unaware of the DER-type the request is coming from. In this case, the VPP accepts the same proportion of requests across all DER classes. The second architecture is the decoupled architecture in which the DER-type is included in the packet request. The VPP then optimizes charge and discharge acceptance rates for each DER class separately. 2) Novel control policies are designed for tracking reference power signals that are applicable to both coupled and decoupled architectures. These policies have as secondary objectives the maximization or minimization of the number of packet requests to be accepted (input optimization) and guaranteeing QoS close to a specific energy state (state optimization). One can further co-optimize requests (input) and QoS (state), which bears resemblance to optimal controller problems. With respect to several classes of DERs, the control policies are validated for both architectures and performance metrics are provided. 3) Analytical characterization of the relationship between the controllers implemented in PEM's asynchronous request-response mechanism and the corresponding synchronous representation in the aggregated bin transition macromodel that describes the mean-field behavior. This type of characterization is useful when implementing the control schemes herein under coupled and decoupled architectures in cloud-based serverless implementations of PEM. 4) Analysis of parametric heterogeneity unique to PEM's packet request-response mechanism is conducted by studying uncertainty in the estimated rated power of requests. The tracking performance of a large group of DERs subject to parametric heterogeneity is compared against g smaller groups of DERs with the same rated power of request and modeled by the macromodel described in Section II-B. Analytical bounds on tracking errors are then calculated for the g groups. This article is organized as follows. Section II summarizes PEM and the macro-model. In Section III, reference-tracking controllers that maximize/minimize the acceptance rate of packet requests is presented along with a controller that maintains QoS. In Section IV, parametric heterogeneity with respect to DERs' rated power is discussed and the worst-case tracking error under steady-state conditions is presented. Conclusions and future research directions are given in the final section.

II. MODEL FOR PACKETIZED ENERGY MANAGEMENT
This section provides context for the main results of this article and summarizes the work in [14], [17], and [19].

A. Packetized Energy Management
For each DER under PEM, quality of service (QoS) is ensured by utilizing a probabilistic request rule based on the local dynamic state. In addition, if the QoS falls below a certain user-defined threshold, then the DER is capable of exiting PEM temporarily. Synchronization is mitigated by introducing randomization in the DERs' request rates based on the local state variables. The typical closed loop block diagram of a PEM system is shown in Fig. 1.
Let the energy state of the nth DER in a fleet be denoted as z n . The equation for this DER is given by the following discrete-time dynamic model where f n is a 1-D mapping (usually linear or bilinear). Energy transfer rates of the nth DER are given by P rate c,n and P rate d,n when charging (c) and discharging (d), respectively. The hybrid state φ n corresponds to the set of modes {c, sb, d} associated with {charge, standby, discharge}, respectively [14], [17]. In this article, the focus is on coordinating electric water heaters (EWHs), bidirectional energy storage systems (ESSs), and electric vehicles (EVs). The EWHs and EVs are seen by aggregators as having charging dynamics only whereas ESSs corresponds to both charging and discharging dynamics. The parameter w n ∈ R maps the end-use consumption to the energy state and is explained next for each DER class considered in this article.
The EWH end-use consumption w n corresponds to the hot water extraction from the tank. Since water extraction always leads to loss of heat energy, therefore, w n ≤ 0. In this article, Poisson random pulses are used to model the end-user process as detailed in [17], [19], and [35].Unless otherwise stated, the EWH parameters used in this article are the same as in [19,Tab. 1]. The ESS units are assumed to be installed for the purpose of backup power by the end-user and are managed by the aggregator (e.g., Tesla or VPP herein) to provide grid services. This setup is based on the state of Vermont's largest utility, Green Mountain Power, and their residential Tesla Powerwall battery program [35]. In this program a customer can pay $15/month per ESS (5 kW/13 kWh). This program has been popular as many rural customers opt for multi-day backup capability with two ESSs and, today, Vermont has more than 2000 PowerWalls under control. Therefore, it is assumed in all the simulations presented in this work, the ESS are rated at 5 kW/13.5 kWh, charge and discharge efficiency of around 95% (roundtrip of 92%) along with the set-point of 70% [19] and the background usage is absent. Finally, the driving patterns of the EV driver are modeled using a two-state Markov chain consisting of two modes; driving and parking [19]. Note that in [19], standby was used instead of parking. This driving pattern then determines the end-user consumption of EVs considered in this work. The average drive time of the EV driver is chosen to be approximately 30 min and the EVs are assumed to have an electric driving range of 150 miles with the driving efficiency of 7 milesper-kWh [19]. The request probability in PEM is presented next.
The probability that the nth DER with dynamic state z n [k] ∈ [z n , z n ] and desired set-point z set n ∈ (z n , z n ) over time k (for discretization time-step t) makes a request is given by a cumulative distribution function over the range of admissible dynamics states. For this purpose, an exponential distribution encoding three desirable conditions has been chosen; i ) no request (i.e., Pr(z n ) = 0) from device n for z n ≥z n ; ii) guaranteed request (i.e., Pr(z n ) = 1) from device n for z n ≤ z n ; and iii) Pr(z n ) = p R when z n ≡ z set n for a designed meantime-to-request value, 1/ p R , at the device's desired set-point, z set n , during interval t. This design is aimed at attracting a device's state towards the set-point. The probability of request is then given by where μ(z n [k]) > 0 is a variable rate parameter dependent on the local dynamic state. For charging energy packet requests where m R > 0 [Hz] is a design parameter that defines the mean time-to-request (MTTR) for z n = z set n which is fixed at the mid-point of (z n , z n ). A similar expression follows for μ(z n [k]) in the case of discharging packets. Fig. 2 shows the request probability curves for charging and discharging packets together with the probability of remaining in standby during coordination. While (2) is one possible choice for probability of request, however, this design can be changed or incorporated in the control (as done in [10] and [11]) and is the topic of ongoing work.

B. PEM Macro-Model
The state bin transition macro-model for a large homogeneous population of DERs of the same load type is presented next. A macro-model for a diverse population is comprised of a finite number of different type DER homogeneous populations that work together under the same VPP for demand dispatch. Hence, consider a population of DERs described by (1) with common underlying state space. To create a finite-state abstraction of the entire population evolution, the state space is discretized in a manner that the main features of the system are preserved. Furthermore, the size of the population is chosen so that the law of large numbers holds and the system behavior is approximated by a single averaged effect. This means that the macromodel is robust against the effects of individual DERs as the number of devices increases [36], [37]. The PEM macromodel is introduced next.
The transition probabilities between bins are determined from the dynamical system equations of the DERs comprising the population with respect to the type of end-user events for the particular class of DERs [19]. LetX = {x 1 , . . . , x N }, where each element is called a state and constitute a consecutively ordered partitioning of the continuous state space Z in which the DERs evolve. Denote by q j the probability of being in state x j and q := (q 1 , . . . , q N ) T the probability distribution overX . For the particular case of DERs having hybrid one dimensional dynamics as in (1), an interval [z min , z max ] within Z is divided into N consecutive bins each corresponding to a bin state inX , where x i ∈X corresponds to the interval [z i−1 , z i ) ⊂ [z min , z max ]. Since (1) includes three types of dynamics (charge/standby/discharge), the state space for the system consists of the union of the three identical copies of X given by X = X c ∪ X sb ∪ X d . At time k, the probability mass function of the system is q = (q c , q sb , q d ) with q c = (q 1 c , . . . , q N c ) and q sb and q d defined similarly. Note that q contains the percentage of the population in each state of X so that, 1 3N q = 1, where, 1 3N ∈ R 3N is a vector of all ones. Therefore, if N e is the total number of DERs and N i e,c is the number of devices in state x i c , then N i e,c = q i c N e . Similarly, the percentage of N e that is charging and discharging, and the total power demand of the system due to charging and discharging packets are In this article, symmetric charging and discharging is assumed, that is, P rate c and P rate d are equal. In case these are not the same, their average is taken instead for the sake of simplicity. Henceforth, the subscript is dropped and P rate is used for both charge and discharge rated powers. Furthermore, let T i req,c (T i req,d ) be the probabilities associated with making a charge (discharge) request using (2) in state x i sb . Then the probability of requesting a charging packet and not a discharging packet is given by . Furthermore, it should be noted that simultaneous charging and discharging requests can occur with probability T i req,c,d = T i req,c T i req,d . However, no request is submitted in that case. That is, making a charge or discharge request are mutually exclusive events. The total number of charging and discharging requests are In Section IV, uncertainty in P rate due to heterogeneity in the DER population's rated power is introduced and analyzed.
The discrete-time equation modeling the population dynamics of a PEM system is given by and . Furthermore, the 2-D control signal β := (β c , β d ) defines the proportion of charge/discharge packet requests out of the total charge/discharge requests denoted as n c r /n d r , that are accepted (i.e., defines how many DERs transition from standby (sb) to charge (c) and discharge (d)). It should be noted here that β c ∈ [0, 1] and β d ∈ [0, 1] since the total number of accepted charging and discharging requests cannot exceed n c r and n d r , respectively. In addition, is the proportion of DERs in (c) and (d), respectively, whose energy packets end now and return to (sb). The matrices in (4) (omitting time dependence for all βs) have the form d} is a multi-diagonal matrix containing the probabilities of staying, going to higher energy states and going to lower energy states, M c,sb and M sb,c are responsible for transferring any DER that exceeds z max from c to sb and any DER that falls behind z min from sb to c, respectively, and similarly M d,sb and M sb,d provide the transition probabilities from d to sb and from sb to d including the probabilities of uncontrollable events, T req,c, The computation of the transition probabilities has been addressed in [17] and [19].
Modeling the evolution of the number of active DER charging and discharging packets in the system requires introducing two sets of timer states (c and d). That is, given packet epoch δ, the sampling time step t, and two timer states vectors x p,h ∈ R n p with n p = δ/t and h = {c, d}, the timer dynamics are given by where c,d and B p,h ∈ R n p ×N is responsible for allocating the new charge/discharge population into their corresponding charge/discharge timer states. The number of charging/discharging packet requests received by the VPP is then n h r := 1 Nq h sb with h = {c, d} and 1 N = (1, . . . , 1) ∈ R N . The timer provides the formula for the percentage of DERs whose packet expires. That is,

C. Opt-Out Dynamics
Finally, the opt-out control mechanism that ensures QoS for each DER under PEM is introduced. This is achieved by augmenting q in (4) with opt-out states q ⊕ . That is, q is redefined as q := (q ⊕ , q ) and whereM is trivially augmented with a diagonal block identity matrix and zeros everywhere else (leaving the states q opt unaffected by β and β sb ) of M that has all rows and columns corresponding to states higher than the prespecified PEM re-entry bound removed, M pem provides the transition probabilities of leaving PEM and M ⊕ pem provide the transition probabilities of reentry PEM. Fig. 3 depicts a complete diagram of a DER population under PEM. In this diagram, controlled transitions among charging, standby, and discharging states are regulated by the VPP whereas uncontrolled transitions 1 are also permitted. For a detailed discussion of QoS guarantees, see [19].
Including opt-out dynamics completes the PEM macromodel that captures the aggregate dynamics of a large enough population of DERs. Theoretically, there exist modeling errors as a result of the finite-state abstraction [33]. However, the intended use of the macromodel in this work is to design reference tracking control policies for specific objectives such as satisfying QoS. In this regard, empirical evidence is provided here that shows that the chosen population size is sufficient for the Law of Large Numbers to hold. For that purpose and without loss of generality, consider a population of EWHs when all requests are accepted, that is, β c = 1. In this situation, all DERs get what they need without imposing restrictions from the coordinator. The error in the total power demand between the PEM macromodel and simulated EWH fleet in steady state, called modeling error in power in this work, is compared for different population sizes. The mean of the modeling error in power was observed to be close to zero whereas the standard deviation decreases with the size of the population as expected. Fig. 4 shows the standard deviation normalized with respect to the size of the population and the rated power. It can be seen that the standard deviation decreases proportionally to 1/ √ N , where N is the size of the population. The normalized standard deviation of the modeling error for 250 DERs is approximately 2.5% that corresponds  to ±6.25 EWHs and for 5000 DERs is approximately 1.2% amounting to ±60 EWHs, which is an acceptable tolerance given the fact that the coordinator has no information about the DERs state. Therefore, in the simulations presented in this article, the population consists of a minimum of 250 DERs of the same type. Moreover, the effect of modeling errors on the design of control policies presented in this work will be explored in future publications.

D. Control Architectures for Diverse DERs
While coordinating different classes of DERs within a fleet, the type of DER (e.g. TCL, ESS etc.) information can be included in the packet request along with the size of the packet (e.g. 4.5 kW). This gives rise to two VPP architectures, namely, coupled and decoupled. If the DER type is included in the request, then the VPP generates a separate control input for each class within the fleet. Whereas, in the absence of DER type information, the VPP generates a single control input for the entire fleet. The coupled implementation is illustrated in Fig. 5 for a fleet consisting of TCLs, ESSs, and EVs.

III. REFERENCE-TRACKING CONTROL POLICIES FOR PEM
This section presents novel reference-tracking control policies for PEM that explicitly considers the aggregate packet request rates and QoS. Specifically, control policies are implemented that optimize the number of accepted requests and preserves the QoS. This represents an extension from the request-maximizing policy presented in [17] and [18]. However, in order to track a reference while satisfying QoS for a population, a power baseline or nominal response for PEM needs to be established first.

A. Nominal Response of DER Fleet Under PEM
The flexibility available to a VPP for tracking a grid or market signal is subject to QoS constraints on the population of DERs under PEM. In particular, there exists a constant reference signal (or nominal demand), which maintains the average SOC above the desired value (to preserve QoS), and above/below which the average SOC increases/decreases. That is, the nominal demand represents a simple sustainable trajectory and produces the VPP's nominal response. A definition of the nominal response of a PEM system is provided and discussed next.
Definition 1: The nominal response of a fleet of DERs under PEM, given by (6), is the minimum constant power signal for which QoS is sufficiently satisfied for the average DER.
This nominal response provides a reference about which flexibility can be measured and establishes a relationship between aggregate VPP demand, QoS, and average SOC, which is the first step in the development of a dispatchable virtual energy storage model [23], [38]. This nominal response is characterized in PEM by the nominal control β * = (β * c , β * d ) that is the solution of the following constrained optimization problem where C dem,i q * i = P i is the total packetized nominal demand from the i th class of DERs at steady state q * i . The vector x i v ∈ X i contains the parameters associated with the discretized bin values of the i th DER class (e.g., temperatures, SOCs, etc.) and z i set is the desired set point for the dynamic state of the i th DER class.
1) For TCLs and EVs, the optimal solution to (7) is attained for β c ∈ [0, 1] due to n i=1 C dem,i q * i being a monotonically increasing function of β c (i.e., less packets accepted is less power consumed) and (q * i ) x i v in (7c) is also monotone in β c (i.e., accepting less requests leads to a lower average SOC, z i set ). However, including the ESS class introduces β d ∈ [0, 1] to regulate the acceptance rate of discharging packet requests from batteries that seek to discharge. This additional degree of freedom means that multiple ESS solutions can achieve the same net-power (charge minus discharge) without affecting the ESS average SOC. Hence, with the ESS class, the optimal solution is not unique. Work is ongoing to formally prove and characterize the optimal solution under various implementations of PEM and the scope herein is meant to introduce properties of the nominal response to motivate the reference tracking controller developed in the next section. 2) When a PEM system is in equilibrium (power and SOC) due to a fixed control input β, the internal signal β sb becomes time invariant, and packet interruptions are then negligible since the average SOC is near the set-point as shown in [18]. Specifically, it was also shown that β sb = (β − c , β − d ) ≈ (1/n p , 1/n p ) when tracking admissible trajectories.
3) In the case of coupled implementation (7) generates β * for the entire fleet (TCLs, ESSs and EVs), that is, a single nominal β * is obtained for all DER types. 4) The optimization problem presented in (7) is clearly non-convex, however, ongoing efforts consider convex re-formulations. However, those results are outside the scope of this manuscript. 5) The nominal control input β * also depends upon the parameter m R in the probability of request. Recall that m R is a design parameter that defines the meantime to request. Consider first the nominal charging acceptance rate β * c . Fig. 6 shows that increasing m R results in an increase in β * c in both decoupled and coupled implementations. This makes sense because increasing m R leads to DERs requesting less often. Therefore, the VPP has to accept more requests to maintain the fleet at discharge acceptance rate β * d . Nevertheless, further work is required for fully characterizing the behavior of m R and its effect on control effort that will be the focus of future publications. All simulations in this manuscript uses m R equal to the packet length.  = (0.302, 0). Note that the components of β for ESSs are close to zero since the background noise process is ignored (i.e., batteries are used only to provide flexibility). For EVs introduced in [19], an average driving time of 30 min and a driving occupancy of 20% was used. The nominal control for the coupled implementation yields β = (0.302, 0.208) that amounts to average TCL temperature, ESSs SOC, and EV SOC of 52.4 • C, 70% and 80%, respectively. Note that EVs in this manuscript can only request charging packets and do not participate in vehicleto-grid (V2G) discharging events. This means that the average SOC for a population EVs can only decrease during periods when the population drives more than it charges (on average). The main takeaway from Fig. 7 is that when the fleet starts in a nominal steady state, then the diverse DER fleet maintains its SOC under nominal acceptance rate, β * .

B. Reference-Tracking Controller With Diverse DERs
While the previous section focused on nominal control of a constant trajectory, this section provides explicit control policies for real-time reference-tracking by VPP. In fact, herein a diverse VPP is illustrated, i.e., a single coordinator for homogeneous groups of DERs. The control law that achieves the desired reference tracking determines the proportion of charging and discharging requests that are accepted (i.e., β = (β c , β d )) and is straightforward: 1) measure the tracking error and 2) accept requests so as to minimize the error. This logic is similar in effect to ACCEPT/DENY relay control for the individual packet requests. However, at the macro-level with an incoming stochastic stream of requests, the control effectively becomes continuous in [0, 1] where zero (0) means that all requests are denied and one (1) implies all requests are accepted.
Note that since each packet request can also include its corresponding charge/discharge energy transfer rate, e.g., P rate h,n , which provides the macromodel with an accurate estimate of average packet height P rate h := 1/β h n h r β h n h n=1 P rate h,n for h = {c, d}. Let P ref be the reference power signal provided by some grid or market condition and let P dem be the actual VPP net-load (charging minus discharging packets). In a coupled implementation of PEM, the two scalar control inputs β c and β d are generated by the VPP so that these minimize the tracking error. However, unlike other DER coordination schemes, PEM's unique packet request mechanism allows the VPP to consider maximizing or minimizing the number of accepted requests without affecting the aggregate VPP netdemand.
The scheme is shown in Fig. 5 for a diverse DER fleet consisting of TCLs, ESSs, and EVs. Despite the fact that there is only one VPP for all the DER populations, the VPP can generate either one control signal β = (β c , β d ) for all the DER populations (coupled case) or multiple control signals designed for the specific DER population (decoupled case). The latter will obviously require more information (e.g., a request could include pinging the VPP the DER type). In the coupled case, for the purpose of tracking, the VPP solves the following optimization problem at time-step k: where is the tracking error of the VPP that includes the effect of expiring packets in terms of n sb,c and n sb,d which is the number of DERs that move from charge and discharge modes to standby modes respectively. Thus, after solving (8) the VPP controls are given by Observe that (F (χ)) 2 is only positive semi-definite so the optimal solution is not unique as formulated above. Thus, one can add a secondary objective to maximize/minimize the total number of accepted packet requests (i.e., χ * c + χ * d ). The solution to (8) is provided by the theorem below.
Theorem 1: Define the available packet budget as n error = /P rate c when ≥ 0 and n error = ||/P rate d otherwise. The solution (χ * c , χ * d ) to (8) that maximizes and minimizes the number of accepted requests is given in Tables I and II,  respectively. Proof: Due to the simplicity of the proof, the focus will be in the case where there are enough requests to track the signal and for maximizing number of requests and > 0. The other cases follow directly. When there are enough requests to track the reference signal, it follows that F (χ) = 0. Then, . This solution corresponds to a family of solutions for χ c ≥ /P rate c since discharging > 0 when this number is less than n c r , otherwise χ c = n c r . In both cases, the total number of accepted requests is maximized. This completes the case shown in Table I for > 0 and α r := n c r > n error . The case for n c r ≤ n error is obtained by obeying the constraints on the number of accepted requests in the problem.
The mechanics for the decoupled case are similar to the coupled one, except that n h r for h = {c, d} is no longer the total number of requests across all DER classes, but the number of requests coming from each DER class separately. This decoupled case implies an implicit disaggregation of the reference signal at each time instant. The work in [39] also deals with the dissagregation of the reference signal in the frequency domain so that DER types are pre-selected to track the corresponding pieces. This reduces tracking to considering individual DER types tracking their own reference independently of each other. The decoupled case here relies on (8) and a secondary objective (e.g., maximizing or minimizing requests in Theorem 1, which dissagregates the reference at each time step. In this case, if suddenly some DER types become unavailable (e.g., TCL peak time of the day or EVs start driving en masse due to the end of the working day) then the other DER populations automatically begin to track the part of the reference that was left unattended by the offline DERs. In other words, the decoupled case performs a dynamic dissagregation of the reference. However, when dealing with DER classes with different rated powers, decoupled implementation can potentially prioritize one class of DERs over the others. Consider for example two DER classes with rated powers 4 and 5 kW, tracking a common reference signal while maximizing the total number of accepted requests. The decoupled VPP is inclined towards accepting more requests from 4 kW DER class than the class with rated power 5 kW since the objective is to accept as many requests as possible. It is worth pointing out here that although the grid access of the 5-kW class is reduced, VPP's main objective of tracking the reference signal is not compromised. Therefore, to avoid prioritization of 4-kW DER-class over the other, the objective function in (9) can be regularized by penalizing accepting requests from the 4 kW DER class. Another approach involves a QoS-optimizing policy that prioritizes access to the grid based on the energy needs of the DER class.
A QoS-optimizing control policy for accepting packet requests seeks to track a reference signal while limiting the deviation of the average SOC from the desired setpoints. Compared with the two previously described Minimize/Maximize requests policies, the QoS-optimizing control policy consists of QoS (state) and number of accepted requests (input) that resembles an optimal controller problem [28]. The decoupled QoS-optimizing control policy aims to find the optimal solution, χ QoS c,d , such that where the objective function is given by  ) . The latter is obtained from the nominal problem of (7) for the decoupled VPP case. Finally, ν nom := (ν nom TCL,c , ν nom ESS,c , ν nom ESS,d , ν nom EV,c ) and the scalars ν represent weights on the different objectives. For the coupled VPP case, the box constraint in (9) reduces to 0 ≤ χ c ≤ n c r and 0 ≤ χ d ≤ n d r , where n c r and n d r are, respectively, the total charging and discharging requests of all the DER classes combined. Naturally, one requires to know the current state distribution q of the respective populations. In [18], an extended Kalman filter (EKF) was constructed for the state-bin transition model providing the dynamics of DERs populations. The same EKF formulation is employed to estimate QoS TCL [k], QoS ESS [k] and QoS EV [k] at time k in a manner that the distance between the average dynamic state at k and a pre-established set point is minimized by χ QoS c,d . It is important to highlight here that the scheme is not manipulating individual DER's SOC but choosing a control input β such that the aggregated SOC behaves as desired.
Figs. 8-10 illustrate the decoupled control policies for maximizing, minimizing, and QoS-optimizing the number of accepted packet requests for a diverse VPP with 1000  TCLs, 1000 ESSs, and 250 EVs. The reference signal in the simulations represents a de-trended (with respect to the nominal response) and scaled AGC signal obtained from [40] and shows that the VPP is able to provide approximately ±0.5 MW of flexibility. Furthermore, it should be noted that in Fig. 9 the average power consumption of ESSs is greater than zero where the average SOC remains close to the setpoint. This is because of the nonunity battery efficiency [19] These optimization problems were solved using fmincon in Matlab on a MacBook Pro with a 2.2-GHz processor and 16-GB memory and the solution was obtained on average within 45 ms.
1) Reference Tracking Performance: The root mean square (RMS) tracking error between the reference signal and the total power output is reported in Table III for the coupled and decoupled implementations of PEM. The one-step-ahead up (down) flexibility for a single DER population with n c r charging requests and n d r discharging requests specifies the maximum increase (decrease) in power demand P dem due to that DER population. Specifically, it is defined as follows: )P rate,⊕ where n sb,h , for h = c, d is defined in Section III-B, n ⊕ is the number of DERs opting out, n sb ⊕ the number of DERs leaving opt-out mode to sb mode and P rate,⊕ is the average power rate of DERs in opt-out mode. In the decoupled implementation, the total upward (downward) flexibility is obtained as the sum of all individual DER classes upward (downward) flexibilities within a fleet. The coupled implementation, however, considers the aggregated charge and discharge requests, aggregated charging and discharging expiring packets as well as the aggregated number of DERs leaving and entering the opt-out mode to compute the up (down) flexibility using the average P rate c and P rate d rates within the DER classes in the fleet. The system flexibility is a function of β in that f lex [k] corresponds to β c = 1 and β d = 0 and f lex [k] corresponds to β c = 0 and β d = 1. In addition to RMS tracking error, the average of the one-step-ahead up/down flexibility is provided in Table III. The RMS tracking error and one-step-ahead flexibility suggest that minimizing the number of accepted requests improves flexibility -the QoS policy is similar too and was implemented with (ν TCL ν ESS , ν EV ) = (19, 1, 1) and ν nom = (45, 1, 10, 10). The main difference between the coupled and decoupled implementations of the three VPP control policies is in the available downward/negative flexibility. That is, ESSs can be used to control down flexibility in the decoupled case while EWHs and EV can only control flexibility upwards. The coupled case does not allow for such behavior since the rate of acceptance of charging packets is the same for all DER types, which limits the down flexibility provided by ESSs.    Table IV that the average dynamic state of each DER class tends to be around the desired set point, which ensures that the average DER satisfies QoS requirements while the population is tracking the reference. Specifically, Table IV shows the RMS tracking error values of the average QoS with respect to the desired set point.
However, the average dynamic state could be misleading if the DER populations show too much variation over the dead-band or beyond it. In this regard, Table V presents the average standard deviation for each DER population under the control laws described in Theorem 1. One can infer from these values that while TCLs and ESSs populations are relatively close to their average dynamic state (≈ 1 • C and around 9% for ESSs), the EV population have an almost constant standard deviation of around 14.5%, which is a result of EVs uncontrollably transitioning to driving (i.e., discharging) events. Note that future work will focus on control mechanisms for allowing PEM to not only regulate with respect to the average population QoS set-point but to actively reduce the variance in the population (beyond opt-out control), while also tracking the reference. Table V illustrates that managing a fleet relative to the QoS set-points does not guarantee an overall reduction in variance.
Note that when the reference signal is centered around the nominal power level of the system, then QoS set-points can be attained. Furthermore, as mentioned in Section III-B, due to the presence of multiple DER types, several combinations of charging and discharging requests can be accepted so that the tracking error is zero. However, the QoS objective of (10) determines the number of requests to accept that ensures QoS as well as the control input β is close to the nominal obtained from (7). In the absence of χ c,d [k] − χ nom c,d term in (10), the VPP has no information about the nominal control input and hence determines χ c,d that minimizes QoS deviation from the set-point at the next time-step depending only on the current SOC of each population. This results in prioritizing one class of DERs whose SOC is farthest from the set-point at the current time step. Moreover, the χ c,d obtained without the nominal information is not guaranteed to minimize QoS deviation over longer time periods. Recall that in PEM, the DERs whose packet requests have been accepted, charge or discharge for a time equal to at least the packet length. As a result, oscillations are produced in the power consumption of individual DER population while the aggregate still tracks the provided reference with minimal tracking error. One can potentially avoid such oscillations by optimizing QoS over a prediction horizon in a model predictive control (MPC) formulation, similar to the one in [28]. However, this is not trivial due to the bin based population model, studied in this work, being nonconvex and having a large number of states. Further work is required to develop a low-order model that admits an MPC type formulation and is the topic of ongoing research. Nevertheless, including the nominal information along with QoS ensures that the VPP obtains χ c,d that satisfies QoS while the control input is close to the desired nominal values.
Finally, Fig. 11 shows the tracking of a signal that is well below the nominal power and, therefore, the system starts losing energy as shown in the SOC plots. Clearly, DERs start opting out as designed by the PEM scheme around minute 350. Different DER classes opt-out at different rates. The figure also shows that as the overall SOC of the system saturates, the system is unable to keep tracking the signal. Again, to be clear, the SOC of the EV population is decreasing because insufficient charging requests are accepted by the VPP compared with the discharging from end-use driving events.

IV. PARAMETRIC HETEROGENEITY OF DERS IN PEM
Thus far, a diverse VPP has been defined as one that coordinates different classes of DERs (e.g., EWH, ESS, and PEV) via coupled or decoupled macromodels as shown in Fig. 5. However, a single macromodel, as described in Section II-B, represents a population of homogeneous DERs that are all of the same type (e.g., EWHs with the same parameters). Furthermore, it is valuable to be able to employ the macromodel to accurately represent intra-DER-class heterogeneity. The goal of this section is to analyze and quantify the tracking error when a parametrically heterogeneous population of one specific class of DERs (with respect to their rated power parameter) is approximated by a group of homogeneous macromodels. The reason to focus on the rated power parameter is due to PEM's unique request-response mechanism for coordination that uses an estimate of rated power to compute the control action, e.g., the computation of β c and β d in Section III requires knowledge of P rate c and P rate d . This introduces uncertainty in the coordination of an intra-DER-class heterogeneous population. Hence, the tracking error is due to the fact that control inputs are computed with respect to a finite discrete group rather than the continuous underlying distribution that defines the heterogeneous parameter variations.
Consider a population of residential EWHs from different manufacturers whose temperatures evolve as where z n is the average tank temperature [ • C] corresponding to DER n, c p is the specific heat constant, ρ is the density of water when close to 50 • C and L n is liters in the container tank. Furthermore, z a,n ≡ z a is the ambient temperature, τ n ≡ τ is the standing loss time constant to ambient temperature, η n ≡ η for all EWHs and w n is the end-use process as described in Section II-A. The parameters are the same as in [19, Table 1]. The mechanical input u φ n ,n = P rate c,n [kW], if logic state φ n = c (charging) and u φ n ,n = 0, if φ n = s (standby). Note that (11) is a particular case of (1). In addition, the model in (11) is heterogeneous only in P rate c,n . Each EWH is simulated individually according to the request-response mechanism of PEM described in Section II and is referred to as agent-based simulations. The rated power heterogeneity considered in this section is introduced next.
The rated power for the nth EWH is assumed to be uniformly distributed over some fixed power rating interval. This forms a rated-power-heterogeneous population of packetized EWHs. Each EWHs is simulated in an agent-based simulation environment and the population is divided into groups based on their rated power operating under the same VPP and denoted as group-based simulations. The VPP receives measurements of the total power demand of the population and determines the number of requests to be accepted. This decision depends upon the information available to the VPP regarding the groups. The objective, therefore, is to compare the behavior of a controlled rated-power-heterogeneous population of packetized DERs divided into g groups, with that of an aggregation of g homogeneous macromodels. The rated power associated with the i th macromodel then corresponds to the average rated power of that group. The set of g intra-DER-class groups under PEM are in an arrangement similar to that in Fig. 5.
To illustrate the idea, let a rated-power-heterogeneous population over the power interval [4,12] kW be used to track a variable regulation signal. An increasing finite set of evenly separated values P := {P 1 , . . . ,P g } withP j ∈ [4,12] kW is used to group the devices. That is, an EWH with rated power P i belongs to the j th group ifP j makes P i −P j be the smallest with respect to all other values in P. Fig. 12 shows the tracking results for a rated-power-heterogeneous population of EWHs compared against 1, 2, 4 groups having rated power evenly distributed over [4,12] kW while tracking the same regulation signal. Note here that g = 1 indicates a group characterized by rated power of 8kW while g = 4 implies groups with rated powers {5, 7, 9, 11} kW. The key takeaway is that the VPP makes decisions with respect to the rated power valueP j of the groups whereas, the rated power of EWHs is uniformly distributed over the mentioned interval. It is also important to highlight that for EWHs there is no difference between the control strategies in Section III since there are no discharging packets available.
As expected, the tracking error is lower when more groups are used. In Fig. 12, the VPP assumes that all requests from the i th group has rated powerP i instead of the exact value which is uniformly distributed. Clearly, if an accurate day-ahead or multi-hour prediction of flexibility from a VPP was critical for feed-forward control (e.g., optimizing VPP's response to a peak demand event), then g > 2 groups are helpful. However, if an observer is implemented to provide an updated estimate of the states, say, every minute, then it is foreseeable that one or two groups could suffice for control. Estimation and prediction of QoS is an important topic and represents ongoing work that is outside the scope of this manuscript. For the sake of simplicity, the analysis hereafter is limited to EWHs, therefore, only charging requests are considered. Furthermore, in Fig. 12, the VPP responds sequentially to packet requests whereas in the macromodel, the VPP gathers all requests within a time interval and determines the total number of requests to accept as discussed next.

A. Comparing Synchronous Model With Asynchronous Reality
By design, PEM is asynchronous and, therefore, responds to each energy packet request in a sequential manner. This means that the controller keeps accepting requests while demand power is less than the reference. This is similar to event-driven, real-time, cloud-based but serverless implementations of large-scale IoT data processing applications. A request in that scenario is an event that triggers a response from the load coordinator amounting to accept or reject the packet request [41]. On the other hand, the macromodel from Section II-B is implemented synchronously and, therefore, gathers all requests over a time step interval, say [kt, (k + 1)t], and then computes the proportion of packet requests Fig. 12. (Top) Reference-tracking comparison between a VPP agent-based (micro-)model with 10 000 heterogeneous EWHs with rated power uniformly distributed over [4,12]kW and g groups each with 10 000/g EWHs evenly divided into groups of constant rated power. (Bottom) The tracking error for different groupings of the 10 000 EWHs.
to accept, β, using the guidelines in Section III. This section establishes a relationship between the sequential and gathered implementations of PEM in the presence of parametric heterogeneity as modeled by g homogeneous groups.
First, consider a homogeneous population of DERs and compare sequential and gathered implementations of PEM. In the synchronous (gathered) case, the PEM coordinator receives requests as a queue. More precisely, let n t be the number of requests received over a time step interval and n a ≤ n t be the number of accepted requests. The sequence of requests is denoted by s r := (r i ) i≤n t , r i is the i th request received and P i is the associated rated power of such request. Denoting by S n t the group of all bijections from the set of n t elements on itself, a permutation is an element of S n t from a given set of elements. An element ∈ S n a is written as If the permutation is realized in a random fashion, then it is called a random permutation. It then follows that another realization of the request process during a time interval produces a random permutation (s r ). If, in addition, the population is homogeneous with the same rated power for each request, then (s r ) is indistinguishable of s r because the output power they contribute is the same. Now, given a sequence s r , asynchronous (sequential) PEM accepts the first n a requests based on order of arrivals. On the other hand, gathered PEM accepts request based on the entire (unordered) sequence s r by computing the percentage of n t that balances demand and the reference. Note in this case that there are n t n a sequences that can be chosen by the coordinator from the n t received requests. In fact, all these sequences form an equivalence class with respect to their output power. Moreover, the output power from this equivalence class is exactly the same as that of the sequential PEM acceptance methodology. This fact indicates that the macromodel for a homogeneous population of DERs can accurately represent PEM.
The case is more complex for heterogeneous populations approximated by g groups of homogeneous DERs. Under sequential PEM, incoming requests may have different rated powers, one can think of g homogeneous sequences s r,i for i = 1, . . . , g. A request received at any particular time can come from any of the g groups depending upon their probability of request. To model such a string of request, let n t,i be the total number of request in group i and r i j be the j th request of the i th group at time τ i j ∈ R. The set τ i := {τ i 1 , . . . , τ i n t,i } is naturally ordered since τ i j < τ i j +1 . For a set A = {a 1 , a 2 , . . . , a n } with a i ∈ R, the ordering operation ord < (.) ∈ S n permutes A into ord < (A) = {a ord < (1) , a ord < (2) , . . . , a ord < (n) } such that a ord < (i) < a ord < (i+1) . Definingτ = ord < (τ 1 ∪ · · · ∪ τ g ), the sequence s r = (r i ) i≤n t such that r i 2 is the request associated toτ i ∈τ constitutes the string of request received sequentially by the coordinator. The coordinator now accepts the first incoming n a requests from s r in similar fashion to the homogeneous case. For the gathered mechanism when g > 1, the acceptance rate is computed as in the coupled case of diverse populations in Section IV. That is, one computes the overall accepting rate as β = n a n t (12) where n a is as before the total number of accepted requests gathered from the g groups without any specific order of arrival. Then β is applied equally to all the groups. The objective is now to relate the sequential and gathered mechanisms when g > 1. To achieve that, first observe that the power output of the accepted request now depends specifically on how the requests from different groups arrive and from the sequence of accepted requests of length n a . For example, the first n a requests could come from the population with minimum rated power, P 1 , or from the population with maximum rated power, P g . Clearly, n a P 1 = n a P g . Nevertheless, the distribution characterizing all the possible sequences of length n a drawn from the n t requests from the g groups obey a multivariate hypergeometric distribution (MHD). Considering that requests from the same group are indistinguishable, the MHD gives the likelihood of a particular sequence to appear in a realization of the requesting process. Let n a,i be the number of accepted requests in group i . The acceptance rate for the i th group is β avg,i = n a,i /n i , where n i denotes the total number of requests received from group i (with mild abuse of notation). From the properties of the MHD [42], the expected number of accepted requests per group is E[n a,i ] = n a n i /n t . The average of the expected accepting rates is then n a n t = n a n t which is the same as (12). In fact, for every i , one has that E[β avg,i ] = E[n a,i ]/n i = n a /n t , which is constant and independent of i . Hence, the expected packet acceptance rate over the set of all possible sequences of length n a obtained from the received request sequence s r with g type of requests is equivalent to the even application of (12) to all the groups. This result will allow the computation of bounds on tracking error using the PEM macromodel and the control schemes presented in Section III.
2 r i can come from any of the s r i sequences forming s r .

B. Group-Based Tracking Error and Information Scenarios
The objective now is to quantify the VPP tracking error by considering each group being modeled by a macromodel that has an estimate of the group's DERs' rated powers. Depending on the distribution of the rated powers (P i ), the macromodel groups can differ in their population count (N i ). As discussed in Section III, the VPP accepts/denies energy packet requests based on the tracking error (or gap) where P dem represents, as before, the measured power demand of the entire DER population, n i sb,c is the total number of packets expiring from i -th group, and g i=1 P i n i sb,c is the power that will be lost due to n t sb,c := g i=1 n i sb,c packets expiring since only EWHs are considered and there are no discharging packets. The information regarding the total number of expiring packets is always available to the VPP from the decision making process and from the timer states, as described in Section II-B, which act as ledger for the expiration of packets.
Recalling that the current total number of packet requests from the i th group is n i and that the total number of requests from g groups is n t = g i=1 n i . In a similar manner, the total number of EWHs in the entire population is where N i is the number of DERs in the i th group. Without loss of generality, assume P 1 < P 2 < · · · < P g , where P i is the rated power used by the i th group. Under this setting, the VPP coordinator receives a total n t number of requests from the packetized EWHs and the total power demand P dem . However, the VPP does not know the power associated with each of these requests. Therefore the rated power of the requesting population P is estimated by the VPP in order to compute a scalar control input β (i.e., the proportion of requests to accept) The focus hereafter is when 0 < P gap ≤ Pn t . From Section III, the VPP's control objective tracks the reference to minimize P gap . Thus, the VPP accepts βn t requests and packetized demand increases by Note in (14) that the coordinator's estimate of the requesting power rating P is used explicitly to reduce P gap (i.e., track the reference) since individual rated powers are not communicated to the VPP. Thus, from (14) and (15), the values N i , P i and n i affect the accuracy in which devices supply power in order to make P gap as small as possible. Defining tracking error P err = P gap − P acc , then it follows that It should be noted here that the VPP assumes the packet size of all requests to be P which is an estimate depending upon the information available to the VPP. In this regard, three information scenarios are considered as illustrated in Fig. 13. Furthermore, analytical bounds on tracking error in each of these scenarios are explicitly derived and validated next. 1) Full Information: It is available at all times: i.e., P i , N i , and n i are known for all the groups. Hence, the VPP's estimate of P is computed as the weighted average P * = g i=1 P i n i /n t . Since only packets of fixed power are allowed, the main source of mismatch between accepted packets and P gap is roundoff error. The roundoff errors is considered coming from two cases: sequential and gathered implementations of PEM. Since, the VPP makes decisions with full knowledge of the power of each request, the error incurred P * err in sequential PEM is bounded by Relation (17) is due to the fact that the last accepted request brings P gap within half a packet. However, in the gathered PEM case, the proportion of accepted requests is applied equally to all the groups, therefore Equation (18) tells that each group incur into a roundoff error of half the size of its corresponding P i . Hence, gathered PEM with groups is not implemented in the field since it amplifies the roundoff error. 2) Static Information: It is available offline: only N i and P i are known for all the groups. The resulting estimate of rated power for the entire population is then In this scenario, the coordinator is unaware of the rated power of individual requests, therefore, it estimates P gap from the total number of requests n t and the total number of expiring packets n t sb,c as where P n t sb,c and P n t are the coordinator's best estimate of the demand reduction due to packets expiring and total power associated with packet requests respectively, in the time-step k. The control β is then computed from P gap and P as β = P gap P n t (21) where β = 0 if P gap < 0 and β = 1 when P gap > P n t . The increase in demand after accepting β n t requests is P acc := P β n t = P gap .
Note that (22)  is an roundoff error due to PEM as well as the error introduced due to difference realizations of the request sequence presented in Section IV-A. Thus, using (13), (17) and (20), the VPP tracking error is where, P err, is the error due to different realizations of request sequence and is explicitly derived in the next paragraph. 3) Least Information: It is available offline: only P i is known for every group. Specifically, the VPP has no information about requests from each group (n i ) as well as the composition of each group (N i ). Then a pragmatic estimate of the rated power for the entire population is the simple average, i.e., P † = 1/g g i=1 P i . The VPP's estimate of P † gap is obtained by replacing P with P † in (20). The corresponding tracking error P † err is similar to (23) but with P replaced by P † and the corresponding P † err, instead of P err, which is derived next. In addition, one can provide conservative bounds for tracking errors P * err , P err , and P † err . With the availability of Full information, the only tracking error is due to roundoff error. The error bounds follow directly from (17) as The bounds on tracking error in Static and Least information cases are obtained by quantifying P err, and P † err, , respectively, in (23). Consider first the Static information case where P and P err are the inferred rated power by the coordinator and the corresponding error caused by such a choice respectively. Least information case will then follow from the same procedure. The error due to the assumption that expiring packets are of the same rated power P is characterized by (23). However, as described in Section IV-A, different realizations of the request sequence provide different output powers to balance P gap . The statistics of the MHD is used to bound the error due to realizations of the request sequence s r , denoted as P err, in (23). The difference is now that the number of requests per group n i are unknown. But one can still rely on the fact that N i is known in the Static information case. That is, under the assumption that the g groups have similar probabilities of request (see Section II-A), the estimated value of n i isn i = n t N i /N t . The expectation value of n a,i is E[n a,i ] = n ani n t = n a N i N t and its variance is Var[n a,i ] =n i n t 1 −n i n t n a n t − n a n t − 1 n a n t − n a n t − 1 (25) where (25) follows directly from the MHD [42]. Furthermore, observe in (25) that all quantities are known except for n a . A straightforward maximization of (25) gives n a,max := n t 2 .
Note that n a,max is independent of i , therefore one has that (25) is upper bounded by σ 2 := Var[n a,i ]| n a =n a,max . Furthermore, defining β σ = g √ σ /n t , it follows that the power accepted is redefined as where P err, := P g √ σ and contains the error due to the variations in request acceptance from each group. Using (21) and (26) one can obtain (23) after including roundoff error.
As an example, consider the Static information case with g = 2, N t = 2500, N 1 = 2000, N 2 = 500, and n t = 20. Then, n a,max = 10, σ = 0.84 and β σ = 0.0918. The latter additional acceptance rate amounts to P err, = 1.8 P extra error in power accepted. In general, one can include Kβ σ with K ∈ N in (26) instead of just β σ depending on how much of the distribution tail of n a,i needs to be captured. It is shown by the simulations in the next section that K = 1 suffice in capturing the tracking errors described previously.
The tracking error bounds for the Static information case are now obtained as follows: P err in (23) achieves its minimum (P err,min ) when the total number of expiring packets n i sb,c belong to group with the smallest rated power (P 1 ). This occurs when n i sb,c = 0 for i = 1 (i.e., n t sb,c = n 1 sb,c ). Similarly, P err is maximized (P err,max ) when n i sb,c = 0 for i = g (i.e., n t sb,c = n g sb,c ). The maximum and minimum bounds are P err,max := P * err,max + P g − P n t sb,c + P err, , P err,min := P * err,min + P 1 − P n t sb,c − P err, .
For the case of Least information, the bounds on tracking error are obtained by replacing P with P † in (26) through (28) and P † err, := P † g √ σ . Particularly for the least information case, one has assumed that N i = N t /g since there is no other offline information available. The next section validates the tracking error bounds for the three information scenarios.

C. Simulations
Consider a heterogeneous population consisting of 4000 EWHs with rated power 6 kW and 2000 EWHs with rated power 10 kW. Fig. 14 shows the results of the agent based simulations of the heterogeneous population tracking a scaled and shifted reference signal. In Fig. 14, it can be seen that the full information case (with P * ), results in the tracking error |P * err | ≤ P 2 = 5 kW according to (24) and the RMS tracking error is 1.96 kW. The Static information case produces tracking error bounds according to P err,max into (27) and (28). The tracking error is always within P err,max and P err,min and the RMS tracking error is 12.3 kW, whereas the maximum and minimum tracking error observed is 36 and −45.5 kW respectively. These bounds are a function of the number of expiring packets and hence a degree of variation is observed in Fig. 14. The RMS values of the maximum and minimum bounds (P err, max andP err, min , respectively) are also plotted with dashed lines in Fig. 14. The RMS tracking error for the Least information case is 12.9 kW and the maximum and minimum values of the tracking error are 39.6 and −34.7 kW respectively. The mean error due to the Least information case's estimate is larger than the one for the Static information case, but they are not expected to disagree by much since P = 8 kW and P † = 7.3 kW. Furthermore, the dashed lines in Fig. 14  The second simulation aims to illustrate the effect on tracking when populations vary significantly in rated powers. Consider two groups (P 1 , P 2 ) = (3, 12) kW with (N 1 , N 2 ) = (2000, 500) EWHs. The results are shown in Fig. 15 and the corresponding errors are tabulated in Table VII. The RMS tracking error in the presence of Full information is 1.64 kW which is within ±6 kW according to (24). For the Static information case, the coordinator's inferred power rating of the population is P = 4.8 kW which is biased towards the group with rated power P 1 . The RMS tracking error is 17.3 kW and the maximum and minimum observed errors are 64.5 and −55.5 kW, respectively. Similarly, for the Least information case, the RMS tracking error is 19.5 kW along with the maximum and minimum errors of 65.2 and −56.8 kW, respectively. Furthermore, it should be noted that the error bounds in the Static information case are tighter compared to the Least information case. The reason is that the coordinator with Least information infers power rating P † = 7.5 kW and ignores the size of each group. Furthermore, it should be noted here that although P † is less accurate as compared to the Static information case, however, the worst case bounds account for such inaccuracies by considering all expiring packets either belong to the group with rated power 3 kW (P † err,min ) or to the group consisting of 12-kW EWHs (P † err,max ). Finally, as in Fig. 14, the dashed lines in Fig. 15 are the RMS values of the maximum and minimum bounds on tracking error.
As expected, the Full or Static information cases outperform the Least information case, which highlights the performance gains possible in PEM by including static and dynamic information in each request sent to the VPP. Of course, more information per packet request requires larger data exchanges and leads to larger communication overhead for PEM, which represents a possible trade-off between communication and control. This topic represents future work.

V. CONCLUSION
Reference-tracking control policies were analyzed and presented for aggregated, diverse DER populations operating under PEM. This allowed for real-time coordination of VPP resources in a power balancing scenario. Specific to PEM, a QoS-preserving nominal response was defined and used to compute an effective power reference signal that enables effective and sustained tracking. Several scenarios were presented to validate the control schemes. Among these, were different cases of information sharing among DER groups of packetized DERs. In addition, control policies were presented that optimize (maximize/minimize) the number of accepted packet request and prioritize QoS with respect to a desired set point for the DERs. The latter case highlights interesting, undesirable effects arising from the DER population dynamics constrained to maintain specific QoS levels and requires prioritization of flexible DER classes. Performance metrics developed include RMS tracking error, up/down power flexibility, and average standard deviation of the DER populations. Additionally, a grouping strategy for handling intra-class DER heterogeneity under PEM was developed. First, sequential and gathered accepting mechanisms were studied and related to each other showing that gathered acceptance constitutes the averaging scenario for all possibilities of sequential accepting. Finally, bounds on reference tracking performance were quantified analytically and validated with simulations to describe the effects of approximating intra-class heterogeneous populations of DERs (relative to their rated power) with groups of homogeneous macromodels.
Future work will consider including m R and packet size δ in the optimization problem (7) to obtain the nominal control β * (m * R , δ * ). This will allow the VPP to select the best probability of request curve corresponding to each DER class. Furthermore, the analysis of the error due to the aggregation processes as in [33] for the case of EWHs will improve the quantification of the flexibility provided by PEM. One can then expand the heterogeneous assumption to include thermal parameters (e.g., tank capacity and insulation). Finally, state-estimation and predictive control policies for PEM based on the average state of charge in the form of observer design and a low-order virtual battery model [38] will also be investigated. The state estimator will build on preliminary results in [18] where a Kalman filter performed well for the macro-model under some simplifying assumptions.

DISCLOSURE
The author M. Almassalkhi is a Co-Founder of startup Packetized Energy, which seeks to bring to market a commercially viable version of Packetized Energy Management.