Optimal and resilient coordination of virtual batteries in distribution feeders

This paper presents a novel hierarchical framework for real-time, network-admissible coordination of responsive grid resources aggregated into virtual batteries (VBs). In this context, a VB represents a local aggregation of directly controlled loads, such as smart inverters, electric water heaters, and air-conditioners. The coordination is achieved by solving an optimization problem to disaggregate a feeder's desired reference trajectory into constraint-aware set-points for the VBs. Specifically, a novel, provably-tight, convex relaxation of the AC optimal power flow (OPF) problem is presented to optimally dispatch the VBs to track the feeder's desired power trajectory. In addition to the optimal VB dispatch scheme, a real-time, corrective control scheme is designed, which is based on optimal proportional-integral (PI) control with anti-windup, to reject intra-feeder and inter-feeder disturbances that arise during operation of the power system. Simulation results conducted on a modified IEEE test system demonstrate the effectiveness of the proposed multi-layer VB coordination framework.


I. INTRODUCTION A. Background and Motivation
Coordinated control of demand-side, distributed energy resources (DERs), such as grid-tied PV inverters, distributed battery storage, and thermostatically controlled loads (TCLs; e.g., water heaters and air conditioners) is part of the solution that supports a renewable energy future [1]- [4]. Much of the recent literature on the coordination of DERs has focused on distributed control methodologies to turn large-scale aggregations of DERs into dispatchable grid assets (similar to [5]- [7]). Since the aggregation of DERs is dispatched as a single entity by a centralized coordinator and is subject to power and energy limits, the DER fleet is often referred to as a "virtual battery" (VB) [8]- [10]. To control and dispatch the VBs, much of the literature has focused on optimizing the VB operation over ISO/TSO market signals, but this does not directly consider the underlying AC network and the distribution system operator (DSO) constraints (e.g., voltage or power limits), see for example [11], [12].
To avoid violating operational limits of the grid and to ensure system reliability with DERs at scale, coordination between a DSO and DER owners and aggregators will become critical. This has spurred a multitude of concepts and models for how DSOs can interact with DERs, aggregators, and whole-sale (transmission) markets [13]- [15]. In this manuscript, we focus on the so-called "Market DSO" model, e.g., see [13], where the DSO performs coordination and aggregation of DERs to deliver grid services. While such a setup could preclude independent DER aggregators (i.e., increases regulatory complexity), the model simplifies the interaction between wholesale market signals and the DSO and the ideas herein can be adapted further to enable independent Fig. 1. The multi-layer VB coordination framework. The blocks with the inter-and intra-feeder controllers and optimal VB dispatcher form the focus of this paper. The DSO runs an AC OPF about every minute to dispatch VBs with optimal set-points. The inter-feeder and the intra-feeder control action, which take place in real-time (about every 5 s and 100 ms respectively), are also executed by the DSO. Specifically, the intra-feeder control computation occurs at the feeder's service transformers monitored by the DSO (this requires broadcasting from the substation to the service transformers).
"grid-aware" DER aggregators [16]. For other market-based DER coordination schemes, "transactive energy" can engender holistic TSO-DSO-Aggregator participation of DERs [17]. Some of these schemes focus on broadcasting prices directly to devices. However, with large-scale participation of DERs, transactive energy is susceptible to harmful load synchronization effects, power oscillations, and volatile prices, as shown in [18]. Other grid-aware approaches include optimization-based methods to account for AC network constraints [19], where VB control is achieved by solving an optimization problem based on AC network models and tracking a Karush-Kuhn-Tucker (KKT) point that satisfies the KKT optimality conditions. However, for non-convex AC OPF, the KKT conditions may not be sufficient to guarantee global optimality. Other optimization schemes can provide market services with VBs without exact grid models nor real-time measurements [20]. However, these methods do not directly incorporate multi-period energy constraints and the KKT point can be sensitive to exogenous disturbances.
Formulations based on convex relaxations can provide guarantees on feasibility and optimality of solutions when the cost function is monotonic [21]. However, optimal tracking of a power reference signal has a non-monotonic cost function, which means that one cannot guarantee that the predicted solutions are network-admissible. To address these challenges concerning real-time, optimal, and network-admissible coordination of VBs, this paper presents a hierarchical (multi-layer) framework for coordinating demand-side flexibility in the form of VBs (please see Fig. 1). The hierarchical coordination consists of a novel, convex OPF relaxation, which is provably tight at optimality under realistic conditions and generates grid-aware, feasible set-points for the VBs. The optimization is represented by the box "Optimal VB set-point dispatcher" in Fig. 1.
To ensure real-time tracking of VB set-points despite forecast and modeling errors and to ensure resilience, the hierarchy is augmented with a real-time controller to correct the feeder economic set-points. The real-time controller borrows concepts from wide-area control (WAC) [22], including (local) droop [23] and (regional) automatic generation control, and adapts them to dynamically managing VB power in distribution feeders. This results in a VB coordination scheme that enables a feeder head-node to optimally track a power reference while correcting, in real-time, for unexpected small (local, intra-feeder) and large (regional, inter-feeder) disturbances.
This scheme thus effectively deals with stochasticity, both in the input disturbances to these VBs (e.g., solar variations) -through the PI control mechanism, which rejects these disturbances, as well as that in the VB parameters themselves -by re-solving the OPF problem about every minute, and re-tuning the corrective controller gains about every 5 minutes with updated parameters. Unlike most of the existing literature, the presented framework explicitly and effectively accounts for AC (radial) networks and energy-constrained VB resources and is robust and resilient to grid disturbances. Similar to [9], [24], we consider VBs at the scale of 100-200 DERs per VB that are managed locally (e.g., from the same neighborhood) via load control (with full state information) and subject to lags from device dispatch and communications.
Prior work on hierarchical control of DERs in microgrids (e.g. [25]) has mainly considered using frequency and voltage droop characteristics to generate active and reactive power set-points for DERs using local measurements of frequency and voltage and compensating for the deviations, but in this paper, we compensate for the deviation in the head node power of the feeder from the economic set-point, thus taking into account the economic trajectory. Moreover, in this paper, the design of the proportional and PI gains is done based on the grid topology and device constraints, unlike local droop-based control strategies that are generally not cognizant of the network. While [26] develops a local (proportional) controller that incorporates network conditions into gains to minimize voltage deviations with active power injections, it does not consider system-wide power tracking objectives such as an economic trajectory that satisfies voltage limits across the feeder. Reference [27] suggests a Distributed Averaging PI (DAPI) control strategy to ensure proportional power-sharing and economic optimality, but since it requires extensive communication between the DERs, it may not be feasible on a large scale. Moreover, while the droop coefficients in [27] are chosen proportionally according to DER power limits, state-of-charge limits are not considered, and the coefficients are not optimized to minimize head node power deviations from the economic trajectory. In this paper, we overcome the above limitations through a novel hierarchical control of DERs.
To be more specific, the "optimal VB set-point dispatcher" in Fig. 1 considers a convex relaxation of the (balanced) AC network model. Furthermore, the convex OPF formulation is proven tight at optimality, which guarantees that the prediction of future physical operating states of the grid and the VBs are accurate. This is achieved by decomposing the feeder head-node (i.e., substation) economic reference into the aggregate VB dispatch, net-demand, and approximated total feeder line losses. The key feature here is the use of a first-order prediction model for the line losses to simplify the problem. It is shown through analysis and with simulations that this VB-plus-losses reference can achieve a tight optimal solution under practical conditions. To provide the formulation with updated estimates of feeder losses and network topology, it is assumed that a distribution system state estimator (DSSE) is available [28] (not considered in this work).
The VB set-points from the above optimization problem are further augmented with a real-time controller, which enables tracking under modeling and net-load variations due to, for example, solar PV injections. The controller is divided into intra-feeder and inter-feeder operations. The intra-feeder, droop-like proportional controller is implemented to correct the optimal VB set-points in real-time to track the feeder's head-node power reference signal. The inter-feeder proportional-integral (PI) controller with anti-windup allows for real-time disturbance rejection in the event of larger disturbance across multiple feeders, such as malfunctions in the VB's cyber network or the network structure, thus, ensuring resilience. The only measurement necessary for this purpose is the head-node power from the feeders.

B. Original Contributions
In summary, the main contributions of this paper are as follows: • A hierarchical framework adapted from wide-area control for coordinated control and optimization of VBs in radial, balanced distribution feeders, that considers multiple spatial and temporal scales. • A novel second-order conic AC OPF formulation of a multi-period optimization of VBs designed for tracking the desired power reference at the head-node of the feeder. In this formulation, the feeder structure is taken into account and the second-order cone relaxation is proven to be exact even under significant reverse power flow and with a non-monotonic cost function, which is an improvement over the present state-of-science in radial network OPF presented in [29]. • A practically relevant, real-time, corrective VB control scheme for intra-feeder and inter-feeder fluctuations in net-demand and localized communication outages, that requires only the measurement of the head node power, and is dynamically self-adjusting based on estimates of VB power and energy. • Systematic, simulation-based analysis on IEEE distribution feeders is performed, illustrating the validity of the approach.

C. Practical implementation and Data management
In this subsection, we explain how the presented framework (in Fig. 1) can be feasibly implemented in a practical scheme. In the proposed DSO-centric methodology, we assume that the Market DSO acts as the coordinator and aggregator of DERs and manages the entire scheme. The DSO may also use technology services to manage the DERs via a software interface (e.g., the VB-DER interface in Fig. 1), but this is to help the DSO reduce costs and manage constraints. The DSO runs an OPF every minute or so for each feeder based on the economic head node power reference. In this regard, it is reasonable to assume that the DSO has access to SCADA data and is aware of the grid topology and receives VB state of charge (SoC) estimates from the VB-DER interface to run the AC OPF. For the VB interface, the only information required is the corrected VB power set-point. Apart from this, the VB interface does not require any other information and it is not involved in solving an OPF.
The re-tuning of controller gains takes place by executing a nonlinear optimization problem by the DSO about every 5 minutes and only if the model parameters/grid structure change. For that purpose, the requirements are the updated grid topology, VB models, the operating point from the optimal dispatcher, and the VB power and SoC every minute, but only the head node power measurement is required in real-time (on the order of 100 ms). The VB power and SoC estimates are provided by the VB-DER interface. While every other network communication can be accomplished through SCADA, broadcasting head-node power measurements in real-time is not possible through SCADA, which operates on the order of 2-5 seconds. However, real-time automation controllers (e.g., SEL RTAC [30]) are sufficient and can process distribution substation power measurements and broadcast the tracking error to the feeder's VBs via the DSO's Gigabit Ethernet or fiber-optic networks on a timescale of about 100 ms [31].

A. Mathematical Notation
Consider a radial, balanced distribution network as a graph G = {N , E}, where N is the set of nodes and E is the set of branches, such that (i,j) ∈ E, if nodes i,j ∈ N are connected, and |E| = n, |N | = n+1. Node 0 is assumed to be the head-node (i.e., substation) node with a fixed voltage V 0 and define N + :=N \{0}. Let V i be the voltage phasor at node i and I ij the current phasor in branch (i,j) ∈ E. Then, we define v i := |V i | 2 and l ij := |I ij | 2 . Let S ij =P ij +jQ ij denote the sending end power flow from bus i to bus j where P ij and Q ij denote the active and reactive power flows respectively and let s i = p i +jq i denote the power injection into bus i where p i and q i denote the active and reactive power injections, respectively. Next, denote r ij and x ij as the resistance and reactance of the branch (i,j) ∈ E, respectively, which gives complex branch impedance z ij =r ij +jx ij .
Then, based on the DistFlow model for radial networks [32], the variables (s,S,v,l,s 0 ) at any time-step k are described by the following equations: Apart from the nonlinear relation (4) of l to S and v, (1)-(3) represent a linear relationship between the nodal power injections s, the branch power flows S, and the nodal voltages v. Thus, in an AC OPF optimization formulation, (4) would be a non-convex equality constraint, which begets a non-convex formulation. To overcome the non-convex formulation, we present an approach in Section II-C that accounts for the nonlinearity introduced by (4) via a secondorder cone relaxation [33], [34]. However, first, we introduce the VB concept and its corresponding dynamic model, which augments the AC network constraints in the OPF formulation.

B. Virtual battery model
In this paper, a VB is considered an energy-based model of a dispatchable aggregation of a relatively small number of controllable DERs (e.g., 100-200 distributed loads like ACs, which thus have the capacity to both consume and produce energy). This representation of DERs as an aggregate VB is adequate, based on several works in the literature (e.g., [8], [9], [35]). Moreover, while doing this aggregation, it is important to consider the human-in-the-loop in these flexible demand resources, which manifests itself in the form of State-of-Charge (SoC) of VBs [36], based on Quality of Service (QoS) constraints (e.g., temperature/comfort limits). It is assumed that the DERs reside in the low-voltage secondary network, while the local VB coordination (computation and control) takes place at a nearby primary feeder node. Owing to the small scale of aggregation, the DERs that make up a VB are controlled directly via utility Gigabit ethernet (e.g., IEC 61850) or wireless cellular/Wi-Fi connection. From literature, it has been shown that a VB endowed with a DER control policy, such as a priority-based switching of DERs [37], can be well-described by a first-order dynamic model of the VB's energy state, (e.g., see [8], [9], [38]): where at node i ∈ N + , B i (t) is the VB's state of charge (SoC), p b,i (t) its active power output, and p in,i (t) its desired total active power. The upper (lower) bound of the SoC is given by B i (B i ) while the VB's upper (lower) power limits are denoted p b,i (p b,i ). Note that (5) captures the SoC dynamics with α b,i as the energy dissipation rate. The coupling between the VB's ability to change output power and the VB's control and communication of DERs and their response is generalized with (6) as a lag-and-delay model. In that model, T d,i is the time delay associated with communicating with the ith VB and τ i is the time constant of the first-order model (similar to [39] for example). Specifically, (6) was formed by taking note of the following facts: i) The DERs that compose a VB turn on/off (possibly) sequentially, and there are power electronic components present inside each VB, both of which contribute to a net lag τ i ; ii) There are communication delays (generally of the order of 200 ms) between the head node of the feeder and each VB [40], [41], and delays associated with disaggregating the control signal into device-level signals [37]. The delays we consider in the VB model are in fact both these types of delays lumped together. The effect of time delays on stability is investigated in Section III-D. For this work, the battery charge and discharge efficiencies are assumed to be unity. Inclusion of nonunity battery efficiency requires binary variables to avoid simultaneous charging and discharging. A detailed description of the battery model with non-unity efficiency and analysis on avoiding simultaneous charging and discharging is provided in [42], [43]. This full VB model described in (5)- (8) is employed in simulations involving the real-time, corrective control scheme detailed in Section III. However, for the optimal dispatch to be presented next, since we optimize the VB set-points on a minutely timescale and α b,i ≈0 in realistic settings, it is reasonable to consider a simplified, discretized VB model for set-point optimization. This reduced, predictive model is below: valid for any k and where ∆t is the discretization timestep. Next, we augment the discrete-time VB model with the AC network model to formulate the feeder head-node reference power tracking OPF problem.

C. Conventional convex formulation
Convex relaxation techniques have gained popularity recently due to the existence of global optimality guarantees for AC OPF problems [29], [44]. In this subsection, we present a traditional secondorder cone programming (SOCP) relaxation of the AC power flow equations to formulate a convex, multi-period reference tracking OPF problem, (P1). Decoupling the cost function in the form: , the optimization problem can be expressed as: for discrete time-step k ∈ T over a prediction horizon T := {0,...,T − 1} and where the power injection s i at a node i ∈ N + is constrained to be in a pre-specified, compact, convex set S i ∈C.
The prediction horizon length T is chosen subject to availability of forecasts and depending upon the dynamics of VBs [8]. In our formulation, we choose T = 10. The set S i ∈ C depends upon the VB constraints, e.g., in case of inverters this set is given by: where S i is the apparent power limit of the inverter. In (11), the cost function (11a) minimizes the deviation of the head-node power p 0 [k] ∈ R from the economic head-node reference trajectory p econ ∈R is the active power demand at node i ∈ N + , and P S,i [k] ∈ R + is the solar PV generation at node i ∈ N + . We are optimizing over the VB dispatch, p b,i [k] ∈ R ∀i ∈ N + , which appears in power balance constraint (11c) while reactive power demand, Q L,i [k]∈R, is used in (11d). The second-order cone relaxation of the nonlinear equation (4) is given in (11e); (11f) and (11g) provide constraints on power injection and voltage magnitudes.
Several works in literature such as [29] provide conditions under which the second-order cone relaxation is exact for distribution networks. If an optimal solution of (P1) w * = (s * ,S * ,v * ,l * ,s * 0 ) is feasible for OPF, i.e., w * satisfies (4), then w * is global optimum of OPF and (P1) is said to be exact. Theorem 1 in [29] provides conditions for the SOCP problem in (P1) to be exact, however, it requires the part of cost function f 0 (p 0 )=(p 0 −p econ 0 ) 2 to be strictly increasing, which is not the case when tracking a reference power signal. Thus, even under the conditions provided in [29], (P1) may not be exact. To overcome these shortcomings of (P1), we propose a novel method for convexifying the AC OPF while ensuring an exact solution at optimality. Specifically, we utilize a linearized approximation of line losses and through this obtain a cost function that is strictly increasing in p 0 in order to satisfy the conditions in [29].

D. Reformulated convex formulation
To reformulate (P1), we consider each piece of the twopart composition of feeder's predicted head-node power, . Specifically, we employ a first-order approximation of total predicted line losses, L 1 [k], in the cost function (via p i -to-loss sensitivity factors, which we denote by ζ i ), and prove that with approximated losses, the solution is tight at optimality. Thus, the predicted grid response to an optimized VB dispatch is AC feasible. To achieve the above, we use p econ VB and p econ 0 mentioned earlier. This leads to a multi-objective reference tracking problem, similar to the form of a linear quadratic regulator (LQR) from optimal control [45]: subject to: (1)-(3), (9), (10), and (11c)-(11g), for discrete time-step k ∈ T over a prediction horizon T := {0, ... , T − 1}. The parameters α, ∈ R + are chosen appropriately with 1; L 1 [k] ∈ R is the first-order estimate of the total feeder line losses at time-step k and L 0,k ∈R + is the loss estimated for the operating point at time k. The term p 0 results in a tight relaxation as will be shown next in Theorem II.1, whereas the ) represents the change in network loss due to change in active power injection at node i ∈ N + , with p i [0] being the nominal injection. The factors ζ i ∈ R provide the first-order change in feeder losses due to changes in VB power injections. Similar power transfer distribution factors are often used in transmission system analysis but have recently been adapted for distribution networks [46].
Remark. The formulation in (P2) can easily be extended to account for solar curtailment as a control variable resulting in a more general formulation. If P econ C ∈R + represents the curtailment reference trajectory and P C,i ∈R + represents the solar curtailment at node i ∈ N + , then i∈N + (P S,i −P C,i ) is the net solar output and i∈N + P C,i − P econ C is the error in tracking the curtailment trajectory.
In the next section, Theorem II.1 proves that under practical conditions, the (P2) has a zero duality gap.

E. Exactness of formulation (P2)
In order to explain the notation in Theorem II.1, consider L:={l ∈N | ∃k ∈N such that k →l}, which denotes the collection of leaf nodes in the network. For a leaf node l ∈L, let n l +1 denote the number of nodes on path P l , and suppose P l ={l n l →l n l−1 →...l 1 →l 0 } with l n l =l and l 0 =0.
Also, define a + :=max{a,0} for a∈R and let I 2 denote the 2×2 identity matrix, and define vectors u i :=col{r ij ,x ij } and matrices are upper bounds on P ij and Q ij and are chosen so that A i only depends on the SOCP parameters (r,x,p,q,v).
Furthermore, let (Ŝ,v) denote the solution of the Linear DistFlow model, then where P i denotes the unique path from node i to node 0. Since the network is radial, the path P i exists and is unique. Physically, S ij (s) denotes the sum of power injections s h towards node 0 that go through line (i,j). Note that (Ŝ(s),v(s)) is affine in s, and equals (S,v) if and only if line loss z ij l ij is 0 for (i,j) ∈ E. Then based on the DistFlow model define: which denotes the power injection region wherev(s) is upper bounded by v. Since v(s) ≤v(s) (Lemma 1 in [29]), the set S volt is a power injection region where voltage upper bounds do not bind. Then based on this notation, Theorem II.1 below proves the exactness of (P2).
Theorem II.1. The SOCP problem (P2) is exact if the C1 and C2 conditions given in Theorem 1 of [29] are satisfied: C1: A ls A ls+1 ...A lt−1 u lt > 0 for any l ∈ L and any s,t such that 1≤s≤t≤n l ; C2: Every optimal solution w * = (s * , S * , v * , l * , s * 0 ) satisfies s * ∈S volt Proof. The cost function of the optimization problem (P2) can be expressed as: As f 0 in the cost function in (12a) is strictly increasing, the SOCP formulation satisfies all the conditions provided in Theorem 1 of [29] and hence the proof is a direct application of Theorem 1 in [29] under conditions C1 and C2. This concludes the proof.
The term p 0 is added to satisfy the additional condition in Theorem 1 of [29], where the cost function must be increasing with respect to p 0 . Note that the inclusion of this term in the cost function affects the optimal solution. However, > 0 can now be made arbitrarily small (per the proof of Theorem II.1), which ensures that the impact on the optimal solution is negligible.
C1 can be checked apriori and efficiently since A and u are simple function of (r,x,p,q,v) that can be computed in O(n) time and there are no more than n(n+1)/2 inequalities in C1. For practical parameters ranges of (r, x, p, q, v), line resistance and reactance r ij ,x ij << 1 per unit for (i,j) ∈ E, line flowP ij (p),Q ij (q) are on the order of 1 per unit for (i,j)∈E and voltage lower bound v i ≈1 per unit for i∈N + . Hence, A i is close to I for i∈N + , and therefore C1 is likely to hold. As has been shown in [29], C1 holds for several test networks, including those with high penetration of renewables.
To show the practical restriction of condition C2, Fig. 2 shows the increase inv with increase in reverse power flow due to increased VB injections. From the figure, it can be seen that the condition is valid for VB injections up to more than 400% of demand, compared to the base case of 20%. It can also be seen from Fig. 2 thatv matches the actual voltage v very closely due to the low impedance of the IEEE-37 node system resulting in small loss term z ij l ij . However, with solar PV penetration and an increase in impedance values, the maximum VB injection limit will reduce.
The next section presents simulation results that illustrate the effectiveness of (P2).

F. Optimal VB dispatch and head node tracking simulation
Simulation tests for the optimal reference tracking with VBs were conducted on the balanced version of the IEEE-37 node test feeder [47] to compare the conventional (P1) and the proposed (P2) formulations against the actual AC load flow from Matpower [48]. Simulation results in Fig. 3a show the reference tracking results in Matpower achieved through the flexibility of VBs using (P1) and (P2). The figure shows that (P2) can track the reference trajectory whereas (P1) cannot. For (P2), the error in tracking at each step change in the reference trajectory is due to the first-order loss approximation used in (P2). Since the loss approximation is updated every time-step, the effect on tracking error is small and corrected quickly as shown in Fig. 3a, which means that (P2) represents a reference-tracking OPF formulation that can effectively dispatch VBs while guaranteeing network admissibility. On the other hand, for the convex (P1), the non-zero duality gap creates a mismatch between the predicted power flow values and the actual AC power flow, which results in the sub-par tracking illustrated in Fig. 3a. Specifically, (P1) predicts perfect tracking, but the realized AC head node power does not match the grid reference, which results in suboptimal use of VB resources. The comparison of the aggregate State of Charge (SoC) obtained through (P1) and (P2) is shown in Fig. 3b. Clearly, (P1) predicts a different SoC trajectory than (P2) due to the non-physical solution of (P1).
Remark (High voltage conditions). The voltage condition in Theorem II.1 (C2) is not restrictive for practical distribution networks as can be seen from Fig. 2. Importantly, condition (C2) can still be satisfied under large reverse power flows, which occurs in feeders with significant penetrations of batteries or solar PV generation. To illustrate the effect of reverse power flows, we present simulation results in Fig. 4. For example, in Figs. 4a and 4b, it is shown that the predicted active head-node power matches the actual power while satisfying the voltage constraints. Of course, reverse power flows and high voltage conditions are related, which means that we are assuming that appropriate DER hosting capacity studies have been conducted to inform operations and avoid high voltage conditions. Nonetheless, theoretically, there are reverse power flows for which condition (C2) is violated at optimality, which means that the convex relaxation in (P2) may not be tight. Figs. 4c and 4d illustrate the effects of a non-tight solution and show that the mismatch between the predicted and actual head-node power can lead to voltage violations due to predicted (relaxed, fictitious) losses that ensure a feasible solution in (P2), but are not realized in the physical feeder and, thus, reduce the head-node power further. To guarantee an exact solution that is always physically meaningful, we could include additional constraints to (P2) that capture condition (C2) implicitly (e.g., augment (P2) with the LinDist formulation's voltage variables,v i , andv i 's upper voltage bound), which ensures that v <v i , if (P2) is feasible [29]. Alternatively, we could just project (P2)'s optimal solution onto the AC feasible set and accept the loss of optimality.
Note that the solve time for the above optimization problem is typically less than a second for a feeder with 37 nodes. Our prior work on large scale three-phase systems has shown that the OPF scales well and can be solved in under a minute [43]. The next section develops the real-time control of VBs to achieve tracking of grid reference by countering variations in net-demand and also to provide support under contingencies. In Section IV, the OPF is integrated with the real-time corrective control action to show the effectiveness of the combined approach through simulation results.

III. RESILIENT & CORRECTIVE CONTROL OF VIRTUAL BATTERIES
The optimal VB set-point dispatcher depicted in Fig. 1 solves (P2) about every minute and provides optimal active power set-points to all VBs in the feeder, such that the economic power trajectory for each feeder is tracked optimally. The aggregate of DERs that make up these VBs is then expected to provide and maintain those power set-points until the next set-point update. However, due to the nature of distribution feeders with solar PV, there will inevitably be short-term fluctuations in net-demand, which act as disturbances within a feeder. The flexibility available from the VBs can be used to mitigate these intra-feeder disturbances, by correcting the set-points provided by the optimal VB set-point dispatcher. Furthermore, for a utility with multiple feeders connected to a substation, one feeder may suffer from larger disturbances, e.g., forecast errors not accounted for in (P2), cyberattacks on the VBs' or DERs' communication channels, and changes to network topology from local outages. In these cases, it becomes important for the system to be resilient [49] and maintain the economic set-point provided by the market despite such inter-feeder disturbances. In this section, we hence provide a resilient and corrective control mechanism for mitigating these intra-feeder and inter-feeder disturbances by leveraging the flexibility of VBs. This ensures that feeders with high penetration of solar PV can be effectively dispatched to provide energy market services. This "real-time" corrective action will be executed in the order of a few hundred milliseconds.
The ideas proposed here are similar to existing ideas in wide-area control, including primary (droop) control and automatic generation control (AGC) of frequency [22], but they are adapted for distribution system operations and regulating power. For example, the intra-feeder control mechanism is needed to respond to small disturbances within a feeder without invoking VBs in other feeders and hence needs to operate on a fast time scale (similar to primary frequency control). The inter-feeder control, on the other hand, needs to operate only when there is a large contingency that necessitates a response from VBs present in all the feeders, and thus its response is expected to be slower than intra-feeder control, which is similar to automatic generation control (AGC). Unlike primary frequency control and AGC, we also need to take into account saturation in the energy and power of the VBs while designing for and implementing the real-time controllers.

A. Overview of Intra-Feeder Control of Virtual Batteries
The purpose of the intra-feeder control scheme is to reject shortterm disturbances (like solar PV/demand fluctuations) that enter the nodes inside a feeder and maintain the economic market setpoint at the feeder's head-node (i.e., substation). Our proposed intrafeeder control scheme essentially consists of a bank of proportional controllers, one to control each VB in the feeder. For this purpose, the only measurement required is the active power at the head-node of the feeder, which is readily available at the substation. The block diagram in Fig. 5 represents one feeder with n VBs and disturbances in the form of uncontrollable, unscheduled nodal injections. The VBs are modeled as mentioned in Section II-B where the power and energy bounds of the VBs are represented by the saturation blocks. P uf denotes the corrected economic reference head-node power for this feeder, as calculated by the inter-feeder controller described in the next section, and needs to be tracked by this intra-feeder control scheme. P h denotes the head node power of the feeder. The corrected VB set-point for the rth VB, p in,r , is then obtained as follows: where P set,r refers to the set-point provided by the optimal set-point dispatcher about every minute, K r is a proportional controller that is further adjusted by a factor K adj,r (as defined in (18), see below). The net VB power, p b,r , (computed from (6)) is injected into the "Feeder" block (which represents the feeder with nodal active power injections as input and its head node active power as output) at the respective location. The "Feeder" block is linearized when designing the gains, as will be described later. The plant that is controlled here thus consists of the VBs, together with the grid. Estimates of the SoC signal, B r , computed by the DSO based on DER models, and the net VB power signals, p b,r , obtained from the VB interfaces, are used to dynamically modify the adjustment factor for the proportional gain (described below), K adj,r about every minute (this adaptive and dependent behavior of K adj,r is indicated by slanted arrows in Fig. 5). This adjustment factor ensures that an unreasonable value of power that might deplete or fully charge the VB is not demanded from the VB. The computation of the control input is performed at the nodal service transformers that are monitored by the DSO. To obtain the response from the set of DERs in practice, however, we are assuming that a device level algorithm (like the priority-based scheme in [30]) is in place to ensure that the DERs can provide the response that is desired from them. This device-level dispatch can be performed every few hundreds of milliseconds.
To ensure a good disturbance rejection capability, we select the proportional gains, K r , optimally, using an approach similar to standard LQR. Specifically, we assume that the system (i.e., the feeder) is affected by nodal injections that are stochastic. This is reasonable since solar PV and demand fluctuations are typically random. The objective is then to reduce the effect of these injections on deviating Fig. 5. Intra-feeder Control. The adjustment factors K adj,r are adaptive, and depend on the VB power, p b,r and the SoC estimate, Br, (obtained through the VB interface) as given by (18). Note that the signal P uf is updated about every 5 minutes and Pset,r are updated about every minute, thus operating at a slower time-scale compared to the other signals.
the head-node power from the economic reference. Alternatively, we can treat the reference signal itself as stochastic and reduce the standard deviation of the tracking error. To do that, first, the system is linearized by 1) ignoring the saturation blocks, and 2) linearizing the AC power flow equations about an operating point set by (P2). The latter is done by finding the sensitivity of the head node active power of the feeder to the active power injection at the particular nodes where the VBs are situated. The time-delay blocks are also ignored, but to account for time-delays, the controller design is done in such a way that the gain-crossover frequency of the open-loop system, or the bandwidth of the closed-loop system (without delays), is less than 1/5T d,max , where T d,max is the maximum time delay in the system. Effectively, this ensures that the delay-free system and the delayed system behave similarly. Then, assuming that the system is excited by zero-mean wide-sense stationary Gaussian inputs, we minimize the sum of the variance of the error, P e , denoted by σ 2 Pe , and a weighted sum of the variances of the control inputs to each VB, P u1 , ..., P un , denoted by σ 2 Pu1 , ..., σ 2 Pun respectively: Here, σ Pe (K) and σ Pur (K) are related to the standard deviation of the reference via the H 2 norm of the transfer function, and to the gains K r through the Lyapunov equation, AΣ + ΣA + BB = 0, where A, B are system state matrices and Σ is the state covariance matrix. See [50] for details. ρ r >0 is a constant penalty parameter that we design to be inversely proportional to the power capacity of the rth VB. This penalizes power extraction from VBs that have a lower capacity to provide power output, thus resulting in a constraint-aware controller. The above nonlinear optimization problem can be solved efficiently (typically within six seconds for 30 VBs or less) using a path following algorithm [51]. Note that there is a trade-off to be considered while choosing these constants. If ρ 1, the control effort is less penalized, and hence, the resulting gains would be large (allowing larger control effort), and thus, VBs would saturate before they can assist in recovering the total head node power during a major disturbance.
On the other hand, if ρ 1, the control effort is penalized, and hence the resulting small gains may not be sufficient for better disturbance rejection capability. Thus, these constants should be chosen by taking into account this trade-off.
As a VB's energy state approaches its full capacity, the charging rate should be proportionately reduced, and when the energy state becomes empty, the discharging rate should be proportionately reduced. This helps to avoid a sudden step-change in power to zero when the VB saturates (either empty or full capacity). To achieve this, the proportional VB controller gains, K r , are then passed through an adjustment factor K adj,r (similar to standard gain scheduling) that is updated about every minute as follows: where B r and B r represent the lower and higher bounds on the state of charge of the rth VB respectively, B r represents the state of charge of the rth VB, p b,r is the power output of the rth VB, w is a constant greater than 1, and B rl and B rl are some chosen constants such that B r <B rl <(B r +B r )/2<B rl <B r .

B. Overview of Inter-Feeder Control
The intra-feeder control action described above is useful for mitigating small nodal disturbances that arise inside one feeder. However, to mitigate larger disturbances, such as if many of the VBs in one feeder saturate, the VBs inside one feeder may not be sufficient, and the flexibility of VBs in all the feeders will be required to reject the larger disturbance until (P2) updates the set-points. The block diagram of the real-time control system for mitigating these interfeeder disturbances is shown in Fig. 6. It is essentially a PI control scheme that corrects the economic set-point references to the m intra-feeder control systems, the structure of each of which is shown in Fig. 5. To ensure that the PI control is effected only when there is an appreciable change in the total head node power, a dead zone is also implemented before the PI controller. Moreover, a standard antiwindup mechanism is implemented inside this control scheme. This ensures that if all VBs inside all the feeders saturate in power, the PI control action is cut-off to avoid any integral windup. This control scheme requires measurements from the head node powers from all the feeders. Since there is communication overhead involved in receiving the head node powers from all the feeders, the inter-feeder control action is proposed to be effected only about every 5 s.
The working of the inter-feeder control system is described as follows. The sum of measured head node active powers from all feeders, denoted by P h,net , is compared with the total economic market set-point for all feeders, P econ,net . Then, for each time step k (executed about every 5 s), the control input to the ith intra-feeder control system, P uf,i [k], is computed as follows: where P u,net is the output from the saturation block: and P e,net is the tracking error, P ed is the deadzone limit, P e,net is the output of the deadzone, K w is the anti-windup gain, K p is the proportional controller gain and K i the integral gain in the PI controller, P u,net is the unsaturated output of the PI controller, P u,net is the upper limit of the sum of control inputs to the intra-feeder control systems, P u,net is the lower limit, K fi is a scaling factor, and P econ,i is the economic reference for the ith feeder. P u,net and P u,net are computed by assuming that the VBs are at their power limits.
To penalize the extraction of power from feeders with lower capacity to provide power, the inter-feeder gains K fi are chosen in proportion to the total capacity of the feeder with which they are associated. If the feeder i has a higher capacity, K fi is assigned a higher value, and if it has a lower capacity, K fi is assigned a lower value. Specifically, K fi = P i /P , where P i is the power capacity of the ith feeder, and P is the total power capacity of all feeders.
The PI gains are chosen to satisfy the desired requirement of phase margin and settling time. A high phase margin ensures that the system remains robust to VB failure, while a low settling time of less than a minute ensures that any large disturbance is quickly rejected before the optimal dispatcher provides new set-points to VBs every minute. To do this, first, the system is linearized in the same manner mentioned in Section III-A. Then, the phase margin of the open-loop transfer function from P e,net to ∆P h,net , as well as the settling time of the entire closed-loop system due to the response to a step input applied at ∆P econ,net , are computed for different values of K p and K i , and the PI gains for which the phase margin is high and settling time is low is selected. Fig. 7a shows an example of this process of selection of PI gains on a two-feeder system with 2 VBs in each feeder. The values of K p and K i were varied, keeping their ratio fixed, to study their effects on the phase margin and settling time. It can be seen that when K p =K i (blue curve), there are some values of K p for which the system has a low settling time (< 30 s, as indicated by the green dotted vertical line) but the phase margin varies significantly from 20 • to 80 • . In that region, there is a value of K p and K i for which the settling time is the lowest, and phase margin is above 60 • , indicating the optimal value for the PI gains. Other values of K p and K i , as those obtained from the K i = 0.1K p and K i = 10K p curves, lead to sub-optimal performance.

C. Order and Frequency of Selection of Gains
It is important to select/adjust the gains for both the intra-feeder and inter-feeder control mechanisms described above without hampering the real-time operation of the power system. If all the gains in both the inter-feeder and intra-feeder control mechanisms were chosen simultaneously, then the resulting optimization problem (if formed) would become intractable, hampering the real-time operation of the power system. Since it is important for VBs inside a feeder to respond to local disturbances quickly, without relying on the other feeders, we decouple the computation of the gains and design them in two stages: the proportional controllers inside each feeder are first chosen, considering each feeder to be separate, and then the PI controller to control all feeders is selected. This ensures the proper time-scale separation between intra-feeder and inter-feeder operation that was discussed before. Another benefit of this approach is that the optimization becomes more manageable as well. Moreover, to ensure that the gains accurately reflect the current feeder structure, we propose that the gains be updated at the end of every 5 min. This renders our control scheme adaptive, of self-tuning regulator type. Table I shows illustrative times (run on a laptop with 2.2 GHz Quad-Core Intel Core i7 processor and 16 GB RAM) required for the PI controller tuning every 5 min, as a function of the number of feeders. It can be seen that the time scales linearly with the number of feeders but is well less than 60 s. Real-time validation of this scheme on cyber-enabled devices through an OPAL-RT grid simulator is an ongoing work and will be the topic of a future publication.

D. Effect of Intra-Feeder Controller Gains and Delays on Stability
In this subsection, an analysis of the effect of time-delays and controller gains on the stability of the control system is presented. Consider the intra-feeder control diagram shown in Fig. 5, but with two controllable VBs (i.e., n = 2). For analysis, let this system be linearized as mentioned in Section III. The time-constants of the VBs are chosen to be τ 1 = τ 2 = 1 s. The region of stability was plotted (Fig. 7b) using the system transfer function from P econ to P h , replacing the delay with a third-order Padé approximation (to extract poles). Fig. 7b shows the effect of a delay in the response of one of the VBs on the region of stability. The delay in the second VB, T d2 , was varied from 0 to 200 ms, 500 ms, and 1 s (these values are chosen for the sake of illustration -in reality, delays may be smaller/larger), keeping no delay (T d1 =0) in the first VB. The regions of stability are shown (in Fig. 7b) by different shades of green, with the lightest corresponding to no delay and the darkest corresponding to a delay of 1 s. Note that the regions depicted by the lighter shades of green include those depicted by the darker.
It can be seen that when there are no time-delays in either VB, the region of stability is a half-space. Indeed, it can be shown by analyzing the poles of the transfer function that the region of stability is given by A p1 K 1 + A p2 K 2 > −1 (this provides an analytic bound of permissible controller gains). It also indicates that there is a competition between the VBs: if one of the controller gains increases (and is hence more responsive to changes in head node power), the other must decrease to ensure stability. With time-delays, however, the region of stability is no longer a half-space but reduces to a triangular region (in the case of 2 VBs). Furthermore, this reduction is larger for VBs with greater communication delays.
Thus, the above analysis indicates that for stability, allowable controller gains are limited by 1) the presence of other VBs, which introduce competitive behavior, 2) communication delays, and that only a more restrictive set of gains is permissible for VBs with more communication delays. For an arbitrary grid, once the structure is known, one can do a similar analysis and determine the range of allowable controller gains that maintain system stability, which would be important for the DSO.

IV. SIMULATION RESULTS
To test the real-time VB corrective control mechanism along with the optimal tracking, 10 IEEE-37 node feeders (single-phase equivalents) were simulated (Fig. 8), with VBs at certain randomly picked locations inside each feeder, and active and reactive power noise sources at certain randomly picked locations, as shown in Table II. Both the OPF presented in Section II and the real-time control mentioned in Section III are simulated together in this section, with the real-time control using the full VB model (5)-(8) at every time-step and the OPF using the reduced VB model (9)-(10) every minute. The parameters τ and T d of the VBs are also listed in Table II. Each feeder was assumed to have the same nominal active and reactive power demand that is present in the first phase of the standard IEEE-37 node feeder (which is about 727 kW in active power demand and 357 kVAR in reactive power demand). A unity power factor is assumed for the VBs. The sample time for the simulations was assumed to be 500 ms. The anti-windup gain was assumed to be 1, and the dead zone limit in the tracking error to be 72.7kW, which is 10% of the total base demand in each feeder. The active power limits of the VB was set to ±24.2 kW when there were 6 VBs in the feeder, and at ±29 kW when there were 5 VBs in the feeder, such that the total capacity of the VBs in each feeder represents 20% of the nominal demand. The energy bound of each VB was assumed to be 97 kWh (where there were 6 VBs in the feeder) and 116 kWh (where there were 5 VBs), thus, representing a maximum 4-hour continuous operation of each VB at its power limit. This is a reasonable scenario for many DERs, including water heaters and residential batteries. The controllers were designed according to the process mentioned in Section III. We first illustrate the rejection of intra-feeder disturbances such as fluctuations in solar PV, which represent clouds obstructing solar PV panels. Since (P2) operates every minute, such disturbances would have to be mitigated within a minute. Also, since such a disturbance is local, other feeders should not be actuated to reject such a disturbance. The results for a 1-minute simulation are shown in Fig. 9. After small modeling errors were corrected for by two consecutive runs of (P2) at t=60 s and t=120 s, a step disturbance was introduced at around 125 s into all the nodes of Feeder 1 at the active and reactive power noise locations shown in Table II. This step disturbance can represent, for example, a big cloud cover obstructing the sun. From 145 s to 180 s, random noise was injected into the same nodes instead of a step disturbance. This can represent small clouds intermittently covering the sun. The blue curve denotes the total actual head node power from all feeders. The orange dashed line depicts the economic market set-point desired to be met by all the feeders. The total head node power predicted by employing VB set-points resulting from solving (P2), every minute, is shown by the red dotted line. The optimization horizon for (P2) is 10 minutes. The purple dash-dotted line shows the head node power if the VBs were optimally dispatched, but no real-time control mechanism, as described in Section III, was in place. The results show that in the presence of a step disturbance, the proportional gains inside the feeder try to bring the power back towards the economic set-point. However, due to the nature of proportional control, there is a steady-state error. Moreover, if the noise is random, the intra-feeder control mechanism reduces the standard deviation of the random noise, as in Fig. 9, where the standard deviation reduces from 10.8 kW to 8.5 kW, a reduction by 26%. Note that inter-feeder PI control is not activated in the above simulation since the disturbances are all within the dead zone of 72.7 kW.

B. Inter-feeder disturbance rejection
In this subsection, we illustrate the effect of rejecting an interfeeder disturbance and thus, the resilience of the VB coordination framework. It is assumed that due to an adversarial cyberattack on a feeder, all the VBs in that feeder are set to their power limits. Fig. 10a shows the results of a 5-minute simulation with optimal PI controllers. At t = 150 s, one of the ten feeders is attacked in an above-mentioned manner. As a result of this attack, the power delivered by VBs suddenly change from 36.1 kW to 145.4 kW (an increase by 302.7%). It is seen that PI control brings the head node power back to the economic set-point. Because other feeders also participate in this major disturbance, the head node power from them also changes after the disturbance, as illustrated in Fig. 10b.
Remark (Voltage values). In both the intra-feeder and inter-feeder disturbance rejection, the control mechanism successfully keeps the voltages within acceptable limits. In the examples shown above, voltages in the nodes of the feeders are all above 0.97 pu, with the maximum voltage being 1.0 pu (which is the fixed head node voltage that assumes reactive power support from tap changing transformers is separately available).

C. Integration with Multiple runs of OPF
In this subsection, a 2-min simulation is performed (Fig.  11) to show the full hierarchy of the OPF and the real-time corrective control action in the case of VBs undergoing cyberattack. Specifically, two of the VBs in both the 9th and 10th feeder are assumed to be saturated (via attack) after 15 s, resulting in the PI control mechanism being activated, and obtained power (blue line) being restored to the desired value (yellow line). Note that without real-time control (purple dash-dotted line), restoration would not have been achieved. At the next run of the OPF at 1 min, however, even with no real-time control, the optimal set-point dispatcher adequately re-dispatches remaining VBs to result in power close to the desired power. Note that the spike in the blue line at 1 min exists because the VB set-points change, and so does the net real-time control action (which is the sum of set-point provided by the optimal set-point dispatcher and the corrective control).

V. CONCLUSIONS AND FUTURE WORK
In this paper, analysis and simulation results have been presented in support of a novel framework for large-scale coordination of DERs to support deep penetration of renewable energy. The explicit consideration of (temporal) energy and (spatial) grid constraints and the economic and reference-tracking (techno-economic) objectives have been achieved via a spatio-temporal decomposition approach that leverages information on demand-side flexibility to disaggregate grid economic trajectory into reference control signals for virtual batteries in distribution feeders. A convex optimal power flow (OPF) formulation has been presented that ensures a provably tight optimal dispatch of virtual batteries (VBs) to track an economic power trajectory. Resilient and real-time control techniques to adapt and recover from intra-feeder and inter-feeder disturbances have been analyzed and designed. To show the effectiveness of the decomposition approach and the real-time control mechanism, simulation results have been conducted on a modified IEEE-37 test system. Future work will incorporate reactive power control of VBs and extend the regulation of voltage with inverters. Further, the market economic problem will be considered explicitly [52] to study the coupling between the market layer economic problem and the feeder constraint aware dispatch. Future work will also try to extend this approach to the three-phase unbalanced system operation [43]. Finally, we are interested in extending the multi-period (P2) to a robust formulation to trade off conservativeness of dispatch and the probability of voltage and VB violations [53].