Robust hierarchical dispatch for residential distribution network management considering home thermal ﬂexibility and model predictive control

In the transactive energy (TE) paradigm, the devices of participative consumers, or prosumers, may be aggregated and employed to drive operational objectives at the network level. Home heating, ventilation and air-conditioning (HVAC) systems in particular are well-suited to modulate their behaviours based on both home thermal ﬂexibility and requests from the utility grid. This paper develops a robust, hierarchical power dispatch scheme in the context of a residential distribution network. The formulation couples a unique, multiphase linear distribution optimal power ﬂow (OPF) at the upper level with model predictive control (MPC)-based HVAC ﬂeet controllers at the lower level. The proposed approach is tested on nearly 2000 homes with a three-phase distribution network in an intraday market setting, where two major applications are explored and analysed.

distribution networks (ADN), operators and controllers can use price signals to modulate aggregated device response [11].
Although HVAC systems present an opportunity to optimize distribution market operations, the challenge is in the management of a large number of devices. Sections of the network with a high proportion of responsive devices require smart and effective ways to facilitate their participation. While purely centralized approaches may be overburdened with the management of so many devices, hierarchical and more distributed approaches offer a better management strategy. The hierarchical control approach put forth in this work offers an effective way to manage the interactions between the distribution network and home devices. It does so by both enabling prosumers to participate in the distribution market and allowing the ADN to operate its network in a secure and stable manner. This mutually beneficial operational integration between the ADN and its customers is a key goal of TE development.
TE methods are increasingly being incorporated into distribution operations and planning [12] and, as such, they have been a popular and important research topic. In the next section, some relevant studies will be briefly noted.

Relevant literature
This work falls within an application of TE known as buildingto-grid (B2G) and, more specifically, home-to-grid (H2G). In such schemes, HVAC fleets use signals from the network level to optimize their power consumption throughout the day. This is like the optimal scheduling of a fleet of electric vehicles (EV), where much groundwork has already been laid, e.g. [13][14][15][16].
While an EV fleet can export power back out to the grid, however, an HVAC fleet may only shift its power consumption in time.
In general, B2G/H2G schemes attempt to couple a distribution network with HVAC fleets in an optimal dispatch routine. As noted in [17], many approaches are available and typically vary by a few major features, including temporal versus incentive-based response; centralized versus decentralized decision-making; type and performance of optimization approach; and others. Some promising approaches have been put forth in general coordination of demand response with network and household operational constraints [18][19][20][21][22][23]. For example, [21] focuses on day-ahead scheduling problem for one individual home energy management system (HEMS); [22] integrates DR aggregators with distribution network operators for a secure and efficient scheduling; and [23] investigates an aggregator coordination mechanism to support distribution network. A common assumption of these solutions is a pre-trained building thermal model, which may not match the dynamics of real building operation. Hence, it may be challenging to incorporate more realistic consideration of thermal load operations into the models in the literature.
Since some form of distribution optimal power flow (OPF) routine is required in a B2G/H2G scheme, a major distinguishing aspect in research revolves around the distribution network model employed. Distribution networks require careful modelling due to their unique characteristics: Higher line resistance-to-reactance (r/x) ratios, un-transposed lines, radial power flows and imbalanced loading, among others [24]. Some works, e.g. [25], use full AC models to capture such distribution network behaviour. These approaches, however accurate, suffer from the well-known convergence and/or scalability challenges. Linear developments such as in [26] are thus preferred to mitigate these issues, due to their computational efficiency, relative accuracy and ease of incorporation into an optimization program.
Further, most prior works assume balanced network operation and single-line formulations. This may be problematic in practical implementations that require the distribution network to consider how device response impacts load and voltage balance. Therefore, multiphase-based approaches are better suited to characterize lower voltage distribution networks. One exception includes the B2G scheme laid out in [27], which requires an iteration algorithm to solve the problem. The proposed approach in this work, on the other hand, uses a fully linear, multiphase model that both characterizes the distribution network and can be solved in one shot.
Finally, although robust/stochastic distribution OPF formulations for energy management have been documented, e.g. [28][29][30], these are typically not employed in current B2G/H2G literature. Existing hierarchical approaches such as [31][32] mainly focus on a broad DR application instead of deep diving into specific characteristics as well as its implementation in practical B2G/H2G cases. As distribution networks stand as key facilitators for renewable energy integration, uncertainty in at least net bus power should be considered. In this work, interval variables are used to characterize uncertain parameters, and a computationally efficient, robust optimization program is used to solve the distribution OPF.

Contributions and organization
The general contribution of this work is a unique and practical H2G integration scheme, which couples a multiphase ADN model with participating HVAC fleets. The proposed approach blends the strongest available features via the following key areas: • Distribution model -a unique and fully linear multiphase distribution system model is employed to characterize threephase power flows and voltages. • Robust distribution OPF -an enhanced robust optimization formulation is used to incorporate uncertainties as interval variables and solve for a robust solution. • Overall control algorithm -a hierarchical control scheme is used to separate network optimization from the scheduling of participating HVAC units, which decouples the overall computation. Specifically, the upper level handles the distribution OPF routine, while the lower level, model predictive control (MPC)-based HVAC fleet controllers handle the optimal operation of participating HVAC units. Further, the case study shows how the proposed approach can be employed for multiple applications, including strategic peak demand shaving and black-start enhancements, while keeping the network in balance.
The structure of this paper is organized as follows: Section 2 presents the upper level problem formulation, which lays out the distribution system model, its validity, and how the model is incorporated into a robust distribution OPF. Section 3 presents the lower level problem formulation, which describes the home thermal modelling approach and the MPC-based HVAC fleet controller. Section 4 describes the overall control hierarchy and details the interaction between the upper and lower levels. Section 5 presents the test results of this approach on 1920 participating homes on a modified IEEE 13-bus three-phase distribution network. Two applications are explored in this section: Peak demand shaving and a black start enhancement, both in an intraday market context. Finally, Section 6 provides the conclusion, including some notes on future work.

UPPER LEVEL FORMULATION
The approach is built upon a hierarchical control structure, consisting of upper and lower levels. The upper level refers to the domain of the ADN, while the lower level refers to HVAC fleet controllers that control participating homes. While the upper level passes down a reference power signal for HVAC fleets to follow, the lower level passes back up information about HVAC fleet power flexibility. Note that there is one controller (load aggregator) per flexible node. Figure 1 gives a basic visual overview.
At the upper level, the operator wishes to optimally dispatch available resources without violating voltage and power limitations. In other words, the upper level needs to compute a distribution OPF routine. In the following subsections, a multiphase distribution system model and robust optimization formulation for distribution OPF is developed.

Multiphase distribution system model
As mentioned previously, linearized distribution system models are preferred, especially for intraday market applications, where computational efficiency is key. One approach is the well-known 'DistFlow' model, first proposed in [33][34] and more rigorously derived in [35]. Although the full version includes line losses, when these marginal quantities are neglected and bus voltages are assumed close to 1.0 p.u., the simplified version or 'LinDistFlow' equations are formed as: where Ω B is the set of all buses; ( j ), ( j ) is the subset of all children and parent buses of bus j , respectively; r i j , x i j are line i j positive sequence resistances and reactances, respectively; P i j , Q i j are line i j real power and reactive power flows, respectively; p i , q i are bus j real and reactive power injections, respectively; v i , v j are busses i and j voltages, respectively; and v 0 is approximated bus voltages, typically 1.0 p.u. Under normal operating conditions, these assumptions are acceptable. Due to its fully linear formulation and relative simplicity, the LinDistFlow model finds wide use in distribution system optimization and TE applications [36][37][38]. Like most market applications, the balanced operation is usually assumed. The approach in this work employs an extension of this balanced, single-phase version to a multiphase version, which will be referred to as the 'multiphase LinDistFlow' model [39]. In this development, buses range from 1 to N, where N represents the total number of buses, and bus 1 represents the root (main) substation. Lines are then numbered from 2 to N, with the line index corresponding to its 'to' bus. Figure 2 depicts a simple three-phase radial distribution system, for reference.
In brief, a set of LinDistFlow equations is applied to each phase in the network. The result is a set of equations describing a radial network for each phase: where each expression is identical to (1a)-(1c), but now indexed by phase d and associated with the set of nodes Ω N . Note that, in this work, a 'node' refers to a bus-phase pair, while a 'branch' refers to a line-phase pair. For better readability, subsequent equations will omit the phase index d .
In this formulation, if one phase is not present from one bus to another, that equation is simply omitted. Each voltage equation (2c) still uses positive sequence impedances, which can be considered as an approximation of the self and mutual impedances given in branch impedance matrices. The numerical validity of this approach is discussed in the following section.

Numerical verification of multiphase LinDistFlow
To test the validity of the multiphase LinDistFlow approach, it was compared with the power flow results of the forwardbackward sweep (FBS) full power flow method [24]. Error is computed from nodal voltage calculations for each approach, with the FBS results taken as 'true.' The error was investigated along the lines of both loading condition and maximum voltage imbalance, the latter computed as [24]: where dV max i equals the maximum deviated voltage between each phase at node i and V avg i is the average per-phase voltage at the node i. The loading condition is defined as the total load applied across all phases as a proportion of the base case, normal loading condition in per unit. Figure 3 presents the results, which depicts error ranges as a function of imbalance and loading condition.
Numerous loading scenarios were applied in a simple test network to investigate how various network conditions affect the accuracy of the multiphase LinDistFlow model. Expectedly, at near-zero imbalance (i.e. a balanced network) there is minimal error. The error increases in a roughly linear fashion as the network becomes more imbalanced and has an increased range at higher loading conditions. Although this approach is erroneous at high levels of imbalance and heavy loading, for smaller levels the error is minor. For instance, the green dots in Figure 3 show that the error is between 0-2% when the load is up to 2.0 of the base case and the imbalance is up to 5%, and the blue dots show that the error is between 2-4% when the load is 2.0 and imbalance is 5-10%. Since high imbalance is unacceptable to distribution network performance [40] and generally mitigated, the multiphase LinDistFlow model should prove valid under normal operating conditions. The following section describes how this model fits into the overall upper level distribution OPF routine.

Distribution OPF based on multiphase LinDistFlow
An optimal power flow routine attempts to minimize the energy cost of a network while satisfying the operational constraints of the grid [41]. The upper level attempts to minimize a weighted cost function of power importation and flexible load. In this sense, it is multi-objective. Further, the distribution OPF is developed dynamically, where each term has an associated temporal component. Multiple time steps are used to optimally dispatch given future (forecasted) conditions. This 'receding horizon' approach benefits the operator by allowing for some anticipative ability. Thus, the objective function has the following form: where p g i,t is the real power generation at node i at time t ; p fl i,t is the flexible (dispatchable) load at node i at time t ;c g i,t cost of generation at node i at time t ;c fl i,t is the 'cost' of flexible load dispatch at node i at time t ; and i is the comfort weight factor at node i, which is included to capture comfort requirement variations at each flexible node. The T U hor index refers to the number of future time steps considered at the upper level, while the N g and N f indices refer to the number of generation and flexible nodes, respectively.
Of note is the negative sign of the flexible load term: When the cost to import power (the locational marginal price (LMP)) is less than the revenue of flexible load, it will drive the flexible nodes to absorb as much power as possible, and vice-versa when the cost to import is higher.
Nodal power injections are broken down into generation and load components, which are further broken down into dispatchable and non-dispatchable terms. Thus, the following equation is used to model real power node injection (only real power dispatch is considered in this work): where p g i,t is any dispatchable generation at node i at time t ; p re i,t is the renewable (non-dispatchable) injection at node i at time t ; p l i,t is the inflexible load demand at node i at time t ; and p fl i,t represents the flexible load demand term at node i at time t . The objective function is combined with the multiphase LinDistFlow model given in (2a)-(2c) to produce power balance and nodal voltage equality constraints for each phase. As noted previously, in this model each phase of the network is solved using its own set of equations. To link each phase, voltage balance inequalities have been added to the set of constraints: where v an i , v bn i and v cn i are the line-to-neutral voltages phases a, b and c at bus i, respectively; and K v is a maximum voltage imbalance constant that is set by the operator. As to be shown in Section 5, these inequalities help drive voltage (and power) balance throughout the network.
All variables are subject to typical upper and lower limits: where Ω G is the set of generation nodes; Ω F is the set of flexible nodes; p  Due to the inherent variability and randomness of load and renewable generation, net bus injection uncertainty must be considered in a distribution optimization routine. In addition, as flexible load capacities are estimated from the lower level and are also subject to variability, these parameters are considered uncertain. Thus, in this work, two optimization parameters are made uncertain: Nodal net injections, given by p re i,t − p l i,t , and flexible load capacities, given in the terms p f l min i and p f l max i . The following subsection describes how uncertain parameters are modelled as well as how the overall distribution OPF formulation is transformed into a robust, interval-based optimization program for the upper level dispatch routine.

Robust solution based on interval model for uncertain parameters
Uncertain parameters are modelled using interval variables, which are ordered pairs of real numbers that represent an upper and lower bound of a continuous range. Interval variables offer the comparative benefit that they do not require probability distribution functions -only their respective set of bounds. For each uncertain parameter, expected values are used to create a forecast trajectory. Interval terms are then created by computing a lower and upper bound, which correspond with the amount of error associated with each forecast. Thus, an uncertain parameter is transformed into two trajectories at each time step, one representing its lowest possible values and one representing its highest possible values. The uncertainty bounds on interval parameters increase over time; this is based on the practical assumption that uncertainty increases the farther out the forecast is made. In this work, the forecast error increases linearly over 11-time steps from 0 to 10% for both renewable power output as well as load. The 10% error is chosen for this 12-hour window based on the wind forecast error distribution in [59]. The quantities are then summed to make a net uncertain load term.
With the inclusion of these random variables, the base distribution OPF model developed in the previous section is transformed into a robust optimization program with uncertain parameters modelled as intervals, which has the following general form [42]: where vectors of interval terms are given in As mentioned in [43], the indeterminate constraint (8b) cannot be processed directly. Therefore, the original optimization problem needs to be converted to its equivalent deterministic form, first by splitting up the base problem into two subproblems, one representing the lower bound (best case) outcome and one representing the upper bound (worst case) outcome. Although the lower bound sub-problem results in a basic linear program, the upper bound (robust solution) results in a computationally challenging max-min program. While prevalent solution approaches employ traversing algorithms, the approach in this work uses an enhanced technique [43][44], which is briefly described as follows: • First, strong duality theory is used to transform the upper bound, max-min sub-problem into a single maximization problem. • Next, the Big-M method is used to convert combinatorial constraints derived from (8b) into associated inequalities with included binary variables. • Finally, the problem is reformulated into a single mixedinteger linear program, which can now be given to a standard solver.
These upper bound results are taken as the robust optimization solution and applied at each iteration of the upper level routine. Thus, due to the final linear structure of the upper bound sub-problem, the full routine is stable and computationally efficient. This upper level robust scheduling result will interact with the lower level problem, which is elaborated in the following.

LOWER LEVEL FORMULATION
The lower level describes the function of the HVAC fleet controller. Each controller computes the cycling (i.e. the on-off switching) of its participating HVAC units, at a defined sample rate, in line with a power reference signal sent by the upper level. Fleet controllers are applied at each node in the distribution network that has participating homes and are considered flexible nodes by the upper level. The number of homes at each flexible node may be arbitrary; as shown in Section 5, homes are drawn randomly from a distribution to populate the network. The lower level optimization routine is based on unique, dynamic home thermal models, described in the following.

System-identified home thermal modelling
A variety of models are available to characterize home thermal behaviour [45][46]. The following widely accepted approach is adopted in this study: where Ω H is the set of homes, temp in i,t is the indoor temperature of home i at time t , u hvac i,t is the on/off status of the HVAC unit of home i at time t , and temp out t is the outdoor temperature at time t . The parameters a, b and g uniquely describe the thermal dynamics of each home. In this work, a system identification approach is used to develop unique home thermal models for each home.
Although model parameters a, b and g can simply be assumed, in this work a system identification procedure is used to parameterize each home model from historical data. System identification is a well-known optimization procedure to fit a known model to a set of training data [47]. This approach allows for real-world data to be used in building each home's unique thermal model.
In this work, the highly detailed GridLAB-D [48] home objects were used to simulate home thermal behaviour, the data of which was used as the training set for the system identification procedure. Figure 5 gives an example of indoor air temperature data for a GridLAB-D home object (measured) with its system-identified model (predicted).
Once a nominal set of parameters has been calculated, an HVAC fleet can be created by adding several homes to a node and then slightly randomizing the parameters for each home.

HVAC fleet controller
The lower level controller is built upon an MPC formulation, which is a control strategy that utilizes dynamic plant models with an optimization procedure to drive the controlled parameter towards a reference signal [49]. In an MPC implementation, the controller optimizes the control input for the current timeslot considering predicted disturbances. Like the upper level, this is also a receding horizon approach. An MPC-based controller may be formulated in multiple ways to drive the system to a reference signal [50]. In the approach used in this work, the objective function has been modified from the classical formulation [51] to drive the HVAC fleet to follow a reference signal passed down from the upper level, which is a flexible load dispatch quantity. In this formulation, each controller attempts to minimize a penalty function that is proportional to the deviation in the sum of HVAC powers from this reference signal. The objective function is developed as: where P is a constant penalty factor and s k is an introduced slack variable for each lower level time k. The controller attempts to minimize this function from current time step t = k to some defined lower level time horizon t = k + T L hor . The slack variable is then used in the following tracking constraints: where p fl i,k is the flexible demand signal generated from the upper level transformed to lower level time k; u hvac j is the on/off status of home j at time k; p hvac rated is the rated real power draw of HVAC units; and N N i is the total number of participating homes at flexible node i.
Each fleet controller uses the objective function (10), participating home thermal models (9) and slack variable constraints (11a)-(11b) to compute optimal set points at each iteration. A synchronicity constraint is also included to ensure that a high proportion of homes do not simultaneously turn on to cause power spikes. This constraint is developed as: (12) where sync max is a fraction that limits simultaneous maximum HVAC on status. The output of each controller includes a sequence of HVAC set points at each participating home for certain future horizon, as well as updated limits on the amount of flexible load (p f l min i and p f l max i ) that can be dispatched by the operator over its next optimization horizon. The maximum capacity term is simply estimated from the maximum amount of power that can be absorbed by the aggregation of homes at each node. Likewise, the minimum amount of flexible load required is an estimation based on the minimum power required to maintain an acceptable level of comfort. As noted previously, these terms are converted into interval parameters for use in the upper level robust distribution OPF routine.

ALGORITHM AND HIERARCHICAL INTERACTION
The overall control is hierarchical with uncertain parameters considered at the upper level. Therefore, this approach is called robust hierarchical dispatch. There are two interacting control layers, each with its own respective iteration times. It is important to note here that, although this robust dispatch scheme is centralized, each level handles its own domain: The upper level dispatches generation and flexible load at each node and lower level fleet controllers determine how to dispatch their respective HVAC fleet to follow a requested power signal. In this way, the centralized controller will not be overburdened.
In brief, after the upper level completes the robust distribution OPF routine, it sends flexible load dispatch quantities to each flexible system node. Each lower level fleet controller uses these flexible dispatch signals as a power reference in setting HVAC on/off points over its control horizon. Each controller then updates the upper level with p f l min i and p f l max i terms before the next iteration. Figure 6 diagrams the algorithmic approach.
The upper level solves every 15 minutes using a three-hour time horizon. Each fleet controller then uses 15-minute flexible load dispatch signals over a one-hour time horizon to optimize three-minute HVAC set points over that period. In terms of MPC design parameters, each lower level routine uses a sampling rate of three minutes, a control period of five-time steps, and a prediction horizon of 20-time steps. For ease of reference, Table 1 summarizes each level's time parameters.
The algorithm can be applied and started at any point in a simulation study, so long as system data is valid and inappro- priate form. Thus, the approach can be considered a 'rolling window' optimization procedure.

CASE STUDY
The presented approach is implemented in the context of a residential distribution network, where the operator wishes to drive certain objectives, while meeting the comfort requirements of its customer base. The network may import power at the main substation, which is modelled as a generator, and can dispatch flexible load based on the capabilities of each node's HVAC fleet. This work assumes that a portion of the network's customer base wishes to participate in the electricity market, whose HVAC units will be controlled by their respective HVAC fleet controller. To compensate participants, the operator could pass through realized cost savings via a rebate or credit, for example. The simulation environment is built upon Matlab R2019b [52], along with the open-source, agent-based distribution net- work simulator, GridLAB-D. As noted previously, GridLAB-D home object data was used as input to the system identification process used to build home thermal models. CPLEX 12.8 [53] optimization functions were used to solve the upper level distribution OPF routine, while Gurobi [54] was used to solve each lower level HVAC control problem. After primary simulation results were collected, GridLAB-D was used to compute full AC power flows and calculate network imbalance and nodal voltages. The hardware environment is a laptop with Intel® Core™ i7-8650U 1.6 GHz CPU and 16.00 GB RAM.

System description
The residential distribution network is modelled as a modified IEEE 13-bus three-phase distribution network based on [55]. Homes were randomly drawn from a normal distribution at each node, totalling 1920 participating homes throughout the network. Table 2 shows the distribution of homes for each bus and phase, as well as totals across phases and buses. Homes are described by heterogeneous thermal models with unique temperature bands between 67-74 • F. Figure 7 gives a statistical representation of average temperature bands at each bus and phase. The diversity of comfort bands at each home allows a variety of thermal storage capacities available at each flexible node.
Apart from bus 1, the main substation, every bus has inflexible load demand, non-dispatchable renewable PV generation and flexible load determined by each node's number of participative homes. The main substation stands as the only dispatchable generator, which represents the network drawing in power from the utility. The cost to import power follows a typical summer day LMP curve [56], which has high and low of 43.0 $/MWh and 16.0 $/MWh, respectively. The peak pricing period in this case occurs from 2 to 6 PM. Additional major upper level input parameters are summarized in Table 3.  Additional major lower level input parameters are summarized in Table 4. Note that the nominal thermal model coefficients were derived from the system identification procedure described in Section 3.1, and randomized slightly for each home applied to the network (the a parameter was constrained to 1.0).
Two major applications are explored in the case study, which are presented in the following sections. In each case, the simulation is applied at a noted time step, and runs over the course of the day. The proposed dispatch scheme, called the 'flexible case,' is compared with the base case in multiple areas, including total power importation, average home temperatures at each node, and maximum power imbalance. The base case represents normal, non-demand response operation, where HVAC units simply follow their unique set points in a narrow temperature band. Note that base case homes follow a temperature setpoint of around 71 • F. Although two different applications are explored, the objective function and formulation are not modified for either. The first case represents an economic dispatch scenario under normal operation, where network constraints are not driving the outcome. In the second case, the operator must first satisfy network constraints, specifically voltage and power inequalities, before optimizing dispatch economically. The comprehensive formulation makes the approach applicable in multiple contexts.

Load shaping with power balancing
In the first application, the network operator wishes to shape total load demand throughout the day to reduce the network peak, while also balancing the network. The operator accomplishes this through the objective function described in the problem formulation: As the revenue of flexible load dispatch falls below the cost of power import, which occurs during the system peak load condition, the operator will minimize the dispatch of flexible load. Typical residential load curves were applied at each bus and varied slightly by phase to achieve some system diversity [57]. Some solar generation was also applied at each node to produce an uncertain node injection term. Each term is transformed into an interval parameter (in line with the description in Section 2.2) and incorporated into the upper level robust distribution OPF routine. Figure 8 shows total power importation at the main substation throughout the simulation period.
The simulation begins around 7 AM and lasts until the day end. As the base case peaks to near 7 MW, the flexible case reduces the peak to 6.5 MW, resulting in a peak reduction of 7.14%. The operator can achieve the reduction through the optimal dispatch of flexible resources through the HVAC fleet controller and home temperature bands. If this peak shaving occurs during a peak pricing period, it can be expected that significant cost savings may be achieved. The average home temperature for each node is shown in Figure 9.
It can be seen from Figure 9 that the homes straddle their upper temperature bound to reduce power during the peak pricing period. Further, the fleet of homes pre-cool during the event in an anticipative manner. This is the benefit of the receding horizon optimization approach.
Finally, Figure 10 depicts maximum power imbalance (i.e. the maximum percentage deviation of branch real power flows from the average that shares the same to and from buses) at any branch throughout the simulation period.
The flexible case maintains less power imbalance. This is partly due to the voltage balance constraints described previously, as well as each controller driving its respective fleet to closely follow a reference signal, which inherently keeps the system more balanced. Thus, the proposed approach can drive load shifting while simultaneously keeping the network balanced.  There are multiple benefits of more balanced operation in a distribution system, including more efficient utilization of tapchanging transformers and reduced system losses [58].

Black start enhancement with power balancing
In the second application, the network operator wishes to restore power after an outage. In this scenario, there has been a previous power outage during a summer day from 1 to 2 PM, in which all home temperatures have risen well above their upper temperature thresholds. The operator attempts to black start the network at 2 PM. With the absence of flexible operation, as power is restored, there is an initial power surge due to the synchronous actuation of HVAC units (along with inflexible load demand coming back online). This synchronous actuation in the base case could cause voltage and other operational violations at several nodes in the system, which could lead to additional outage time. The proposed dispatch scheme, on the other hand, mitigates this by engaging flexible resources to meet security constraints -in this case, upper and lower voltage limits. The results of this scenario are shown in Figures. 11-14, the first of which depicts the main substation power import.
As the base case results in a peak power draw of around 10.6 MW, the flexible operation does not allow a draw of over 8.1 MW, a peak surge reduction of around 23.6%. Figure 12 gives nodal voltages by each bus and phase.
As the flexible case mitigates the initial power draw and smooths out the power consumption over time, voltage limits are not violated, and normal operation is restored. Figure 13 shows average home temperatures at each node.
The base case operation, shown on the left, cools homes to the reference set point as quickly as possible by consuming maximum power. The flexible case, on the other hand, cools homes over a longer period to mitigate the power surge. Although some homes will require more time to arrive at their respective FIGURE 11 Total power importation at the main substation set points, the extra time required is marginal, and the customer base will benefit from avoiding another possible power outage. The flexible case also restores the network in a more balanced way, as seen in Figure 14. Again, the proposed approach can simultaneously modify flexible load and keep the network more balanced.
It is expected that, without flexibility, an additional outage could occur due to the initial power surge. Thus, the proposed robust dispatch scheme could aid reliability metrics. Where strict power quality standards drive utilities to prevent any network outages, incorporating flexible loads in a smart manner into a black start scenario not only enhances operation, it may be necessary.

CONCLUSIONS
A robust hierarchical dispatch model is proposed in this paper to solve a multiphase distribution OPF at the upper level with uncertainties considered. The underlying optimization model is developed from a unique multiphase LinDistFlow formulation. The formulation is a receding horizon approach such that optimal decisions consider forecast disturbances. At the lower level, an MPC-based HVAC fleet controller is developed to control each participating home at each node. These HVAC fleet controllers similarly use a horizon of data to anticipate both disturbances as well as the trajectory of the reference signal. The reference signal is passed down to each flexible node from the upper level, and the lower level sends back flexible load capacities.
A 2000-home case study was conducted on a modified IEEE 13-bus three-phase distribution network to test the proposed approach in an intraday market context. Two applications were considered in the case study: Load shaping for peak reduction with power balancing, and a black start enhancement with power balancing. In each case, the flexibility enabled by the  Further enhancements could be made to this approach. Although the focus of this work is concerned with home HVAC systems, other devices could be explored and likely solved with the hierarchical approach. For example, an ensemble of HVAC units, EVs and water heaters could be engaged to follow a reference signal from the upper level optimization procedure.
Further, the authors would like to explore and compare the various multiphase distribution models that have and are being developed. Some contexts may require more accuracy than is available from the multiphase LinDistFlow approach.

ACKNOWLEDGEMENT
This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the US Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).