Identification and Validation of Virtual Battery Model for Heterogeneous Devices

TL;DR

This paper develops a first-order virtual battery model to accurately represent the aggregate flexibility of heterogeneous distributed energy resources, enabling improved demand side control for grid services.

Contribution

It introduces a novel identification framework for virtual battery parameters in heterogeneous DER ensembles, including time-varying power limits calculation.

Findings

Effective virtual battery modeling for heterogeneous DERs.

Improved demand response performance through optimized virtual battery parameters.

Validation on ensembles of air conditioners and water heaters.

Abstract

The potential of distributed energy resources in providing grid services can be maximized with the recent advancements in demand side control. Effective utilization of this control strategy requires the knowledge of aggregate flexibility of the distributed energy resources (DERs). Recent works have shown that the aggregate flexibility of DERs can be modeled as a virtual battery (VB) whose state evolution is governed by a first order system including self-dissipation. The VB parameters (self-dissipation rate, energy capacity) are obtained by solving an optimization problem which minimizes the tracking performance of the ensemble and the proposed first-order model. For the identified first order model, time-varying power limits are calculated using binary search algorithms. Finally, this proposed framework is demonstrated for different homogeneous and heterogeneous ensembles consisting of air conditioners (ACs) and electric water heaters (EWHs).

Figures

Tables

Table 1. Table 1: VB parameters for various ensembles

Type	# ACs	# EWHs	a (1/h)	C (KWh)	$x_0$ (KWh)
Homogeneous	100	0	0.0328	56.181	0
Homogeneous	100	0	0.033	61.515	43.682
Homogeneous	0	120	0.0855	71.954	0
Homogeneous	0	120	0.262	80.913	71.508
Heterogeneous	100	120	0.0987	126.869	117.807

Equations

$$\dot{x}(t)$$

$$C_1$$

$$P^{-}(t)$$

$$\text{mbox}{minimize}_s$$

$$T^{set}-\delta T/2\le T(t+1)\le T^{set}+\delta T/2$$

$$T_w^{set}-\delta T_w/2\le T_w(t+1)\le T_w^{set}+\delta T_w/2$$

$$|u_i(t)-s^{\top }P|\le \epsilon$$

$$s\in \{0,1\}$$

$$\text{mbox}{minimize}_s{W}_1||T(t+1)-T^{set}|{|}_2^2+{\sum }_i{s}_i(1-s_i)+W_2||T_w(t+1)-T_w^{set}|{|}_2^2$$

$$\text{mbox}subjectto0\le s\le 1,\text{ and}\text{mbox}Eqs.\text{eqref}eq:constrain{t}_{A}C-\text{eqref}eq:constrain{t}_{p}ower$$

$$x(t)=\sum _i\frac{T_i(t)-(T_i^{set}-\delta T/2)}{{\eta }_i/C_i}$$

$$x_w(t)=\sum _i\frac{T_{w_i}^{set}+\delta T_w/2-T_{w_i(t)}}{1/C_{w_i}}$$

$${\text{minimize}}_{a,C_1,C_2}||\log |B(a,C_1,C_2,x_0)-F||{|}_2$$

$$\text{subject to}B(a,C_1,C_2,x_0)\le F$$

$$B=[b(u_1(t),\varphi ,\bar{\tau })\dots b(u_n(t),\varphi ,\bar{\tau }){]}^{\top },$$

$$F=[f(u_1(t),\bar{\tau })\dots f(u_n(t),\bar{\tau }){]}^{\top }$$

$$P_{reg}=P_{base}+{\tilde{P}}_{reg}$$

Full text

1 Introduction

The last few years have seen a significant increase in integration of renewable energy into the electricity grid, and the intermittent nature of renewable energy causes uncertainty in power generation. Moreover, on the distribution side, there is an increased visibility of thermostatically controllable loads (TCLs) such as ACs and EWHs due to advancements in power electronics, communication capabilities that enable remote monitoring/controlling of TCLs. With increased renewable penetration, these advancements allow TCLs to provide several grid services such as demand response, frequency regulation, frequency response, and tracking regulation [1, 2, 3, 4, 5, 6, 7, 8]. In this work, we focus on characterizing the aggregate flexibility of the homogeneous and heterogeneous ensemble of ACs and EWHs using a first-order virtual battery (VB) model. This work is a preliminary step toward characterizing heterogeneous energy resources such as ACs, EWHs, Solar PV, electric vehicles and real battery as a single VB model.

The allowable power flexibility for any device is characterized by a set whose values are allowable deviations from the nominal power. Furthermore, the aggregate flexibility of the ensemble of heterogeneous devices is obtained by taking the Minkowski sum of the flexibility sets corresponding to each device [9]. Due to the complexity involved in identifying the exact Minkowski sum of the ensemble, the authors in [9] have come up with a method that gives an outer estimate of the aggregate flexibility using geometric methods. Further, authors in [10] have extended the geometric method by incorporating the uncertainties into characterizing the stochastic aggregate flexibility of the ensemble. Another way of characterizing the aggregate flexibility of the ensemble is by representing it with an equivalent possibly a first-order linear system, VB [4, 11, 12]. The state of charge (soc) of the VB denotes the aggregate flexibility of the ensemble and is captured simply in terms of battery parameters such as power limits, initial soc, self dissipation rate and capacity.

There exists several characterizations of aggregate flexibility using VB models in the literature [13, 4, 12, 11]. These works include characterizing a VB model for a wide range of systems from small residential TCLs to complex systems such as commercial building heating, ventilating, and air-conditioning (HVAC) loads. The thermal mass in the residential and commercial buildings enables them to store thermal energy and be used to provide grid services such as frequency regulation by varying the device power consumption. Traditionally, such services are provided by varying the generation at the generators. Similar to a real battery, the VB also has self dissipation rate, capacity, and power limits as parameters. The objective is to identify these VB parameters corresponding to the collection of homogeneous or heterogeneous TCL devices.

Two of the recent and most interesting works in this area are by Hughes, et al. [14, 12]. The authors compare the characteristics of a linear VB model with the actual nonlinear ensemble to identify the VB parameters by solving an optimization problem. This is interesting in the sense that they presented a realistic case study to illustrate the potency of the proposed approach. The VB model discussed here is inspired by [12] and follows it closely to characterize a VB model for an ensemble of homogeneous and heterogeneous devices.

The main contributions of this work can be summarized as follows. We formulated an optimization problem for identifying the device state (ON/OFF) in any given ensemble, such that the ensemble can successfully track the selected regulation signal, while satisfying the device level temperature constraints. Utilizing device level information from the aforementioned optimization problem, the VB parameters are identified by minimizing error in performance between the actual device ensemble and the equivalent VB model. In contrast to the existing works, we propose analytic expressions for the initial estimate of the VB state for different type of ensembles. Binary search algorithms are used to identify time-varying power limits of the identified VB model. Finally, applying the proposed framework, VB system is identified and validated for different homogeneous and heterogeneous ensemble of ACs and EWHs.

2 Virtual Battery Characterization

We begin this section with the description of the dynamics of the first order VB system.

[TABLE]

where $x(t)\in\mathbb{R}$ denotes the state of charge (soc) of the VB at time $t$ with the initial soc $x_{0}$, $a$ denotes the self dissipation rate, and the lower and upper energy limits of the VB are denoted by $C_{1}$ and $C_{2}$, respectively. The regulation signal, $u(t)$ acts as an input to the VB and must always lie within the time-varying power limits $P^{-}(t)$ and $P^{+}(t)$. This first-order VB model can be applied to characterize the aggregated flexibility of distributed energy resources (DERs) and many building loads [4, 11, 12], where the unknown VB parameters are denoted by vector $\phi=[a,C_{1},C_{2},x_{0},P^{-}(t),P^{+}(t)]$.

Homogeneous and heterogeneous ensemble of DERs can be represented using this first order VB model. For example, if an electric vehicle or a real battery is considered, then they appear parallel to the VB and an equivalent battery model can be derived. However identifying VB parameters when an electric vehicle is in the ensemble, requires further investigation, due to presence of additional constraints (such as charging constraints), coming from the electric vehicle. The inclusion of solar PV results only in additional power limits as this is a free instantaneous source of energy, and hence does not result in an increased battery capacity. Although this scenario changes if solar PV is connected to a battery storage. The scope of this work extends to ACs, EWHs, and a heterogeneous ensemble of ACs and EWHs. We refer the readers to [15, 8] for the first order hybrid models of ACs and EWHs. In what follows, we describe the method to identify the VB parameters, $\phi$.

3 Identification of VB Parameters

The VB parameters are computed in two steps: the first step involves finding the time at which the ensemble fails to track a regulation signal and the second step minimizes the difference in performance between the actual DER ensemble and the equivalent VB model [12]. We begin the discussion on control design problem to track a regulation signal.

3.1 Tracking regulation signals

For a given ensemble of ACs/EWHs, the solution to the following optimization problem gives a control law indicating which devices have to be on ON/OFF state such that the regulation signal is tracked at all times and user defined set-points are satisfied.

[TABLE]

where $\varepsilon>0$, $u_{i}(t)$ is the $i^{th}$ regulation signal at time $t$ and $T(t),T_{w}(t)$ denote the vector of temperatures corresponding to AC and EWH devices at time $t$ respectively. The temperature set points and dead band temperature limits for AC (EWH) devices are given by $T^{set}$ ($T^{set}_{w}$) and $\delta T(\delta T_{w})$.

The constraints, Eqs. (2) and (3) refer to the temperature bounds on AC and EWH devices. Notice that this optimization problem is nonconvex due to the binary constraints on the optimization variable $s_{i}$ (ON/OFF status of each device). We make the following two assumptions on the optimization problem to make it computationally tractable.

• •

relax the binary constraints on $s_{i}$ to a convex constraint, $0\leq s_{i}\leq 1$ for all $i$, and

• •

adding the additional cost on $s_{i}$, i.e., $\sum_{i}s_{i}(1-s_{i})$.

The binary relaxation of the optimization variable $s_{i}$ ensures that the optimization problem is convex. This results in a value of $s_{i}$ lying in the set $[0,1]$ which essentially implies that the power applied to either AC or EWH device is continuous which is not correct (refer Eqs. (8) and (9) in [15] or Eqs. (11) and (14) in [8]). Hence, we introduce the additional (nonconvex) cost on the optimization variable, $s_{i}$ by forcing it to take either [math] or $1$ value which results in the following control design problem.

[TABLE]

where $W_{1}$ and $W_{2}$ are some positive weights on the respective costs.

For each regulation signal, the control design problem, Eq. (6) is solved to identify whether the ensemble can be able to track the regulation signal at every time step. If at any time point, the solution does not exist then it implies that the ensemble failed to track the regulation signal. Although this process results in increased computation time, it ensures that the VB model thus obtained is accurate. PJM regulation signals [16] are considered in this work and sample normalized regulation signals corresponding to a two-hour time window can be seen in Fig. 1.

Any feasible regulation signal corresponding to an ensemble is computed by adding the normalized regulation signals to the baseline power consumption of the ensemble. The baseline power consumption is the actual power consumption of the ensemble without any regulation signal being applied. More discussion on this computation is reserved until Section 4.

The criterion for finding the VB model parameters is a comparison of the violation times of the ensemble and the linear first-order VB model. We recall the notion of violation times [12, 14] in the following subsection.

3.2 Violation times

The violation time of the ensemble is defined as the minimum time instant at which the control design problem, Eq. (6) fails to find a solution [12]. Similarly, the violation times for the first-order VB model can be computed when constraints given in Eq. (1c) are violated. Corresponding to the regulation signal, $u_{i}(t)$, the violation time function for the ensemble and VB are respectively denoted by $f(u_{i}(t),\bar{\tau})$ and $b(u_{i}(t),\phi,\bar{\tau})$.

In practical computation of these violation times, a time horizon is chosen over which the regulation signal is considered. Therefore when the ensemble succeeds in tracking the regulation signal, the violation time is considered to be the maximum time in which the regulation signal exists. We refer the readers to [12, 14] for a detailed discussion on computation of violation times. Next, the discussion on VB parameter identification follows.

3.3 VB parameter identification

The time varying power limits of the VB are identified by extending the binary search algorithm proposed in [14]. The algorithm to compute $P^{+}(t)$ is shown below.

[TABLE]

Similarly, the algorithm to find $P^{-}(t)$ can be formulated by changing $\beta$ to $-\beta$ and $\gamma$ to $-\gamma$ in the above algorithm. The initial condition of the VB system corresponding to the AC ensemble is denoted by $x_{0}=x(t)|_{t=0}$ where

[TABLE]

Similarly, the initial condition to the first-order VB model corresponding to an EWH ensemble is given by $x_{w0}=x_{w}(t)|_{t=0}$ where

[TABLE]

Finally in the case of VB for a heterogeneous ensemble (ACs and EWHs), the initial condition for the VB can be obtained as a sum of initial conditions defined in Eqs. (7) and (8).

The self dissipation rate and battery capacity bounds of the VB are obtained as a solution to the following optimization problem [12].

[TABLE]

where the vectors $B$ and $F$ are given by

[TABLE]

[TABLE]

and $n$ indicates the number of regulation signals. The cost in this optimization problem always tries to minimize the difference between the violation times of the linear first-order VB system and the actual ensemble. The motivation behind the optimization problem, (9), is that if the first-order VB model does not cause any violations, then the ensemble also does not cause any violations.

Note that, in [12], the initial condition $x_{0}$ of the VB system is identified as a solution to the optimization problem given in (9). By making the initial soc of the VB an optimization variable, the VB parameters thus obtained are conservative. In this work, we overcome this conservatism by giving an analytical expression for computing the initial soc for ACs and EWHs. In the following section, we consider an ensemble of AC and EWH devices and identify their corresponding VB parameters by applying the proposed framework.

4 Simulation Study

In this section, we present the VB parameters for different ensembles and compare the values thus obtained with existing VB parameter computation methods. The first step in computing the VB parameters of an ensemble is to compute the actual violation times by solving the control design optimization problem (shown in Eq. (6)) with respect to a set of regulation signals. The regulation signals are computed by combining the baseline power consumption and normalized PJM regulation signals. Discussion on the normalized PJM regulation signals is given in Section 3.

The baseline power consumption for the ensemble of AC devices is computed by allowing the ensemble to evolve in the absence of a regulation signal. In this work, a time window of 2 hours is considered and hence the baseline power is computed for 2 hour time horizon. The mathematical representation for the computation of regulation signal is given by

[TABLE]

where $\tilde{P}_{reg}=\gamma P_{base}^{mean}\hat{P}_{reg}^{PJM}$; $\gamma$ is the regulation factor; $P_{base}^{mean}$ is the average of the baseline power; and $\hat{P}_{reg}^{PJM}$ is the normalized PJM regulation signal. This computation ensures that the regulation signal obtained in this way is a feasible input for the ensemble. The regulation signal and baseline power corresponding to an ensemble of 100 AC devices can be seen in Fig. 2. The 2-hour time window is arbitrary and the resolution of the regulation signal is chosen to be $1$ second.

4.1 VB Model: Homogeneous Devices

Here we discuss the computation of VB parameters for two types of homogeneous devices, ACs and EWHs. The AC ensemble consists of $100$ devices and EWH ensemble consists of $120$ devices. The temperature evolution of AC and EWH devices is given in Eqs. (8) and (9) [15] respectively. The parameters such as set-point temperature, outside air temperature, thermal resistance and thermal capacitance for the $100$ AC devices is chosen such that no two devices have the same parameters. Similarly, the device parameters for all the EWHs are chosen. A medium water flow profile is considered for all the EWHs. The baseline power corresponding to $100$ AC devices is shown in Fig. 2.

We first identify the time-varying power limits for each ensemble. The upper and lower power limits ($P^{+}(t)$ and $P^{-}(t)$) for AC and EWH ensembles are computed by applying the binary search algorithm discussed in Section 3. Each regulation signal must satisfy the power limits before asking the ensemble to track it.

Here, we consider $200$ regulation signals for each ensemble and the ones violating the power limits are discarded. The remaining set of regulation signals are then used to identify the violation times. Time-varying power limits corresponding to the AC ensemble and some regulation signals satisfying the power limits can be seen in Fig. 3. Next, the control design optimization problem given in Eq. (6) is solved in Julia by using the Ipopt solver to obtain the violation time wrt each regulation signal. The solution to the optimization problem (Eq. (6)) identifies which AC device has to be in ON or OFF state such that the regulation signal is tracked by the ensemble. As a result, the violation times are computed for AC as well as EWH ensembles.

The violation times for AC and EWH ensemble are now used to solve for the VB parameters such as self dissipation rate and capacity (see optimization problem, Eq. (9)). It can be seen from the optimization problem given in Eq. (9) that the optimization variable, $s_{i}$ which decides the ON or OFF status of any device is allowed to take values between 0 and 1. However, the additional cost on the optimization variable, $s_{i}$ ensures that the resultant values are close to either [math] or $1$ and thereby trying to reduce the error due to convex relaxation.

Therefore, to identify the (binary) status of each device, the solution to the optimization problem Eq. (6), $s_{i}$ is rounded off to its closest integer. Thus the total power of the ensemble can be computed and is shown along with the regulation signal in Fig. 4(a).

Furthermore the distribution of percentage relative error between the regulation signal and the ensemble power is shown in Fig. 4(b).

From Fig. 4(b), it can be seen that the maximum relative error is very small and it shows the ensemble has tracked the regulation signal at every time step (even after relaxing the constraint set to convex and rounding off later to find the device status). The distribution of relative error percentage corresponding to every regulation signal with and without additional cost on the decision variable $s$ is shown in Fig. 5.

We next discuss the difference between the existing VB parameter identification approach such as [4] and draw comparisons with this work. The initial state of VB in [4] is considered to be zero and this assumption leads to a conservative estimate of the VB corresponding to the ensemble considered. Further in [12], the authors make the initial soc an optimization variable which also leads to conservative VB parameters. To see the conservatism and dependency of VB parameters on initial condition, we consider two different scenarios for computation of VB parameters. The initial condition in the first scenario is considered to be zero and in the later one, the initial state of VB is computed analytically as shown in Eq. (7) for AC devices and Eq. (8) for EWH devices.

The initial condition from the analytic expression is zero only under a special case where the initial temperature of every device is equal to their lower (upper) point of the temperature hysteresis band of AC (EWH) devices. It can be seen that as expected, the zero initial condition for ACs and EWHs led to a smaller VB capacity as shown in Table 1. The VB parameters shown in Table 1 are obtained by solving the optimization problem, Eq. (9) in Julia using Ipopt solver. Further, for computational convenience, it is assumed that $C_{2}=-C_{1}=C$.

We now present results that validate the VB parameters and we first consider the AC ensemble for this purpose. A regulation signal is considered and the VB state, energy corresponding to the AC ensemble is computed through the proposed framework. Also, the VB state is computed by extending Eq. 7 for all time points and we refer to this as the analytical VB state. Figure 6 shows the comparison of the VB states computed using analytical and optimization based methods. Figure 6 also show the regulation signal for which the VB state is computed and RMS error between the VB states obtained analytically and through the proposed approach.

Next, we show a comparison of violation times between the actual ensemble following nonlinear dynamics and the first order linear VB system to further validate the VB models identified. Figure 7 shows the comparison of violation times for EWHs ensemble and also for the case of heterogeneous ensemble (ACs and EWHs). It can be seen that the violation times for the VB tend to be very small due to the fact that there are only two optimization variables ($a$ and $C$) for any given set of regulation signals and for every regulation signal, the violation time constraint (constraint in Eq. (9)) has to be satisfied.

4.2 VB Model: Heterogeneous Devices

This section consists of VB parameter identification for heterogeneous devices. In this work, we restrict ourselves to a heterogeneous ensemble with AC and EWH devices. However, with less difficulty, it can be extended to a case with real battery, electric vehicle and solar PV generation along with ACs and EWHs.

The heterogeneous ensemble consists of $100$ AC devices and $120$ EWH devices. The time-varying power limits are computed by applying the one-sided binary search algorithm discussed in Section 3. The regulation signals which lie inside the power limits are chosen and violation times of the ensemble corresponding to those regulation signals are identified. The violation times are computed by solving the control design optimization problem, Eq. (6).

The violation times are then used to compute the VB parameters such as self dissipation rate and capacity. The initial state of the VB is identified by adding the initial state corresponding to the AC and EWH ensemble (see Eqs. (7) and (8)). Figure 8 show that the VB parameters such as self dissipation and capacity saturate when the number of regulation signals are increased.

5 Concluding Remarks and Future Scope

This work consists of VB parameters computation for a heterogeneous ensemble of devices such as ACs and EWHs. A control design optimization problem is developed which identifies the device status (ON/OFF) in the ensemble while tracking the regulation signal and simultaneously respecting the device level temperature constraints. The initial state of the VB is identified analytically and time-varying power limits are computed using the proposed binary search algorithm. The self dissipation rate and capacity of the VB are obtained as a solution to the optimization problem which minimizes the tracking performance of the ensemble and first order VB system. Finally, by applying the proposed framework, VB parameters for homogeneous/heterogeneous ensembles are identified and an extensive validation study shows the robustness and effectiveness of the proposed method. Further, it overcomes the conservatism with existing VB parameter identification methods.

In the subsequent publications, we would like to extend the VB computation with ACs and EWHs to include real batteries, electric vehicles and solar PV inverters. Furthermore, the capacity bounds for the VB can be made more tighter by considering the lower and upper bounds of the VB capacity as decision variables. Moreover, finding the device parameters with some accuracy is always a challenge and remains to be explored.