Bulk Electric Power System Risks From Coordinated Edge Devices

As smart load adoption grows on the electric power system, potential for losing load diversity increases, possibly in ways that impact system stability. Cloud computing resources are able to coordinate large amounts of behind-the-meter loads and resources. Inadvertent or malicious actions could potentially result in gigawatts of load, distributed across large regions, acting nearly simultaneously. We study the resulting impacts of such a perturbation, which were previously recognized, with improved fidelity and granularity using a physically-based power system and demand model. The ResStock tool was used to calculate residential air conditioning load at more than 3,000 locations across the Western Interconnection, corresponding in time to heavy summer and light spring loading. Under an assumption that one cloud platform managed smart thermostats controlling 10%, 15%, or 20% of residential air-conditioning, calculated load steps could be injected into Positive Sequence Load Flow dynamic simulations. These load-driven effects were coupled with two classes of distributed generation ride-through to evaluate the potential for cascading outages. We found frequency deviations in the spring case far exceed the credible contingency event, leading to widespread distributed generation loss, while voltage depressions during the summer loading lead to widespread distributed generation loss and system separation.


I. INTRODUCTION
T HE RAPID adoption of connected end-use loads presents an opportunity to reach previously unachievable levels of energy efficiency and dynamic demand integration with grid operations. As markets adapt to Federal Energy Regulatory Commission rulings 745 [1], 841 [2], and 2222 [3], grid-interactive loads will become more valuable, further accelerating their adoption.
Accompanying the benefits, such devices could have potential vulnerabilities which are exploited to impact electric power system operations. When central cloud computing resources provide command and control access to large numbers of energy-consuming devices, there are numerous new risks to consider. Inadvertent algorithmic glitches or platform errors could cause simultaneous actions across all coordinated devices [4]. Human error, either unintentional or intentional [5] could result in undesired or unexpected responses.
As we rely increasingly on central control platforms and algorithms, legal responsibility for such malfunctions is not well-established [6].
The cloud interface also presents an attack surface for a malicious actor. The compromise of individual edge devices has been widely studied and published, including the penetration of individual devices within the internet of things [7]. An edge device botnet have been used to launch a large-scale denial of service attack on the bulk internet [8]; compromised edge devices were leveraged to pivot into the Ukranian power grid and force a significant outage [9]; and we are aware that U.S. utilities are increasingly concerned about operational disruptions occurring through their digital networks and are bolstering their cybersecurity postures [10].
Regardless of how it occurs, there is potential risk that simultaneous digital operation of large quantities of grid-edge devices could result in negative impacts on the physical world. This paper presents scenarios over which today's utilities and grid operators have no control-disruption of the power system by modifying large amounts of aggregate demand.
This work presents analysis of coordinated demand assets-residential air conditioners-to expose a currently unmitigated risk to the U.S. power grid, which may be further exacerbated by distributed generation ride through adherence. Our work extends the state-of-the-art [11] by using three best-in-class models: the Western Interconnection, the U.S. residential building stock, and residential air conditioners.
The contributions of this work are summarized as follows: • The realistic potential of the coordinated control of cloud connected loads to yield a significant, instantaneous change of system net load is studied.
• The spatial and granular dispersion of residential air conditioning capacity across the western US is established, as is the dispatch at two specific times.
• Dynamic simulations are performed to study a realistic emergent system disturbance: an instantaneous residential demand step comprised of individually small while collectively significant loads.
• The cascading impacts of a collective residential air conditioning load step upon solar photovoltaic distributed generation is investigated depending on the adherent ride through criteria, and a follow on vulnerability is highlighted. A recent report [11] found a significant risk from coordinated residential air conditioners (RACs) when applied to a test distribution feeder and an approximate model of Poland's grid. Other researchers identified a method of attack on an urban power system by manipulating building occupants with false demand response messages [12]. Some recent work has investigated how bulk grid cybersecurity requirements might be applied to aggregated end-use loads that provide market services [13].

LIST OF ACRONYMS
In this paper, we look specifically at loads that may be small individually, but that are operated by one or more cloud computing platforms, including those managed by a third-party aggregator. The example used is cloud-connected RACs in the western United States; as the internet of things grows, other loads-e.g. water heaters, electric vehicles, batteries, connected solar inverters-are likely to comprise similar risk to power grids. There are more than 17 million connected thermostats deployed in North America today (and this is expected to grow to at least 60 million by 2027 [14]). Each of these operates a major load, the RAC, which is the category that dominates utility peak load planning. Of the 118 million households in the United States [15], 90% have air conditioning [16], which means roughly 15% of households (and presumably a significantly higher fraction of RAC demand, roughly 25%) is currently served by connected thermostats.

B. ResStock MODEL
To accurately determine the energy consumption of RACs and their maximum potential energy consumption, we used ResStock [17]-a tool for simulating the U.S. single-family housing stock. It provides sampled probability distributions for each feature of U.S. homes along with scaling multipliers, which together form a representative simulation model of a region's residential buildings. A large number of samples (typically 350,000) are used to ensure that the results are representative of the considerable variety in the housing stock. Annual simulation of the housing stock using Ener-gyPlus [18] was performed using actual 2012 weather data for 216 sites across the contiguous United States. This diversified approach allows for a geographically-accurate representation of, among other things, the RAC energy consumption and the distribution of RAC types and sizes. ResStock outputs include the RAC energy consumption each hour.

C. WESTERN INTERCONNECTION MODEL
The U.S. western interconnection (herein annotated WECC) is the synchronous power system serving 14 states in the western United States, as partitioned by the dashed line in Figure 1. The WECC summer peak demand as of 2019 is approximately 150 gigawatts (GW), which consists of around 80 million electricity customers connected through nearly 120,000 miles of transmission lines [19].
The Western Wind and Solar Integration Study (WWSIS) Phase III [20] investigated the transient stability of the WECC with a large share of renewable energy sources. Ten-year planning cases from 2012 were populated with high shares wind, utility and distributed solar photovoltaics (PV), and concentrating solar power. A light spring (LSP) case-April 4, 2012, at 7:50 a.m. Pacific Standard Time (PST)-and a heavy spring (HS)-August 25, 2012 at 1:20 p.m. PSTfrom WWSIS were used in this study. The dynamic model used for the WWSIS study, which was built within the Positive Sequence Load Flow (PSLF) simulation environment [21] on top of the current WECC dynamic model of the time, was used for this work. The authors have extensive experience with this PSLF WECC model [22]- [24]. Table 1 provides pertinent generation and load data-including the composite load with generation (CMPLDWG) dynamic model representation-for these two cases.

A. CONTROLLABLE RESIDENTIAL AIR CONDITIONING LOAD ESTIMATION
The initial and maximum controllable RAC loads were calculated as follows. ResStock was used to determine operating (online) RAC load accounting for diversity. The max RAC loads are calculated using the RAC-rated equipment cooling capacities from the ResStock sampling routines, the known outdoor conditions for that hour from the actual meteorological year weather file, and the energy input ratio and capacity coefficients [25], which are the same coefficients used in the ResStock model.
To increase the geographic granularity, RAC load was grouped by weather file location, building vintage bin, building size bin, heating fuel, and cooling type, and then disaggregated into census tracts by merging these grouped results using information from [26]. The results were then downselected to include only the WECC territory census tracts.
Census tract loads were assigned to the nearest bus based on latitude and longitude, splitting the load equally among collocated buses. The majority of the loads were allocated to existing RAC motor load in the HS and LSP dynamic models used. However, when changes had to be made to the dynamic model loading to fully represent allocated RAC load, the bus's real-time RAC load capacity was capped at 40% of the total load at the bus, which maintains an existing peak RAC penetration in the initial dynamic model. Once that capacity was reached, loads were assigned to the next-closest bus that still had available capacity. For the HS case, RAC motors were initialized at some CMPLDWG to accommodate ResStock loading.
This load-to-bus assignment was used to assign the online and max RAC loads for the HS and LSP cases, with a summary of a few of the state totals and systemwide totals  provided in Table 2. Figure 1 shows the allocated RAC by bus across the WECC system for the heavy summer case, as well as the fraction of RAC online at the load flow time (e.g., 1:20 p.m. PST). During the LSP case, only 10 megawatts (MW) of RAC are online. The associated WECC system image of LSP RAC online is not presented because of a lack of RAC online. The LSP maximums are similar, but less because of thermodynamics, than those for the HS case shown in Figure 1. These data are summarized and broken up by a few of the dominant RAC load states in Table 2.

B. CONTROLLABLE RESIDENTIAL AIR-CONDITIONING IMPLEMENTATION
More than 75% of the load in the HS and LSP PSLF dynamic models is modeled with the CMPLDWG dynamic model; the components of this model are shown in Figure 2. The three motors are defined via generic model parameters to capture distinct, aggregate types of motor load. The static load is voltage-and frequency-dependent, whereas the electronic load is a constant real/reactive power load. Photovoltaic distributed generation (DG) is also included, which reacts according to the PSLF generic solar inverter model (PVD1) model defined in [27]; the voltage and frequency ride-through VOLUME 8, 2021 characteristics are shown in Figure 6. The DG operates at a unity power factor, with constant power set point. A fourth motor is in the CMPLDWG with the purpose to model the characteristics of single-phase RAC motor load. As shown in Figure 2, for this study, the air conditioning load is partitioned into three types, nonresidential air conditioning, noncontrollable RAC, and controllable RAC. All are still maintained as a single motor in the dynamic model.
To implement the steps in the controllable RAC, PSLF motor loading factor of motor D (LFmD) is adjusted to reflect the commanded step in the controllable AC. This results in a step in the real power consumption of the AC motor D, proportional to the increase from the initialized LFmD (i.e., if LFmD is initialized at 1.0, then changing this to 1.2 during a dynamic simulation will increase the real power consumption by 20%). The loading factor is intended to capture the operating efficiency of the motor at the initialized loading, with changes impacting other values besides the real power consumption. Centering these changes at a unity loading factor minimizes these secondary impacts. This method does not model the in-rush currents that would be exhibited by an air conditioner upon start up [28]. For this reason, these simulations may be considered conservative because the in-rush currents amount to reactive power, leading to larger voltage dips than those caused by a step in real power alone. Figure 3 shows real power changes for a step off followed by a step to max at the San Mateo CMPLDWG model for three different controllable fractions of the RAC. The implemented fractions in this study are 10%, 15%, and 20%. In Figure 3, we see a consistent quantity of nonresidential AC and noncontrollable RAC that remains constant depending on the controllable fraction. The controllable RAC shows a step to zero real power, followed by a step to the max of the controllable fraction. It is evident that even for only 20% controllable, the real power consumption of the controllable RACs after stepping to max is more than one-third of the total RAC load at the San Mateo bus. Figure 4 shows a similar string of controllable RAC steps at the San Mateo bus for the LSP case, but because there is no online RAC at the time of the load flow (7:50 a.m.), there is no RAC to turn off. This results in a rough doubling of the total air conditioner load if only 10% controllable of the RAC units are stepped to max.

C. SIMULATION METHODOLOGY
The aggregate controllable RAC load for each selected fraction of controllable RACs-10%, 15%, and 20%-is shown in Table 3. Note that the total load for the LSP case is 117 GW, 38 VOLUME 8, 2021  whereas the total load for the HS case is 192 GW. The result is that the frequency nadir is expected to be lower for the LSP case for the same size load change as compared to the HS case. For a baseline comparison, the system frequency responses for varying controllable fractions are compared against the frequency response following the WECC credible contingency of a double Palo Verde unit trip equating to a loss of 2,750 MW of generation.
Load shedding driven by either under-or over-frequency power system controllers is not modeled in these dynamic simulations in order to preserve a relative impact, between the controllable fractions. Traces representing the frequency load shed values are included for comparison in the time series plots when the system frequency approaches or crosses these limits. These traces are constructed using the Standard PRC-006 criteria as outlined in [29] and [30].
For all simulations an averaged system frequency is presented as a time series. The system frequency is calculated as a rated mega-volt-ampere (MVA) weighted average of all synchronous generator frequencies across the system. Additionally, the volt-sec (VS) and volt-sec-distributed generation (VSDG) metrics as introduced in [31], and exercised in [22]- [24], are calculated for each of the simulations. The VS metric integrates the deviation from a set voltage range of all buses in the system. The VSDG metric is similar, except it weights the deviation by the quantity of DG on the For this study, the voltage deviation boundaries applied to the VS/VSDG metrics are for voltages below 0.95 per unit (p.u.), and above 1.05 p.u. Each dynamic simulation is started and run in steady state for 1 second before the one of the following control schemes are implemented: • HS: Controllable RAC step off at 1 second • HS: Controllable RAC step off at 1 second, step to max at 40 seconds • HS and LSP: Controllable RAC step to max at 1 second • HS and LSP: Controllable RAC step to max at 1 second, step off at 40 seconds Note that the capability of a RAC unit to turn on and off, or vice versa, in rapid succession, was not investigated explicitly for this work.

IV. CASCADING IMPACTS ON DISTRIBUTED GENERATION
More than 9 GW of DG is online throughout the system during the HS case, and more than 7 GW for the LSP case. Figure 5 shows the location of DG for the HS case throughout the WECC, with each circle representing a CMPLDWG in the dynamic model, the circle size its capacity in megawatts, and shading the penetration as a proportion of the simultaneous native load. Although these dispatches were based on projections, California alone has an estimated 8 GW of behindthe-meter DG in the form of PV as of 2019 [32]. Thus, the magnitude of DG located in the power system models used is a reasonable representation of the contemporary generation portfolio, although the locational aspect is less easily validated.   During frequency or voltage disturbances, the output is scaled according to the implemented ride-through criteria; this process is shown graphically in Figure 6. The defining characteristics are Vt0, Vt1, Vt2, Vt3, and Vrrecov for voltage, and Ft0, Ft1, Ft2, Ft3, and Frrecov for frequency. The four graph points define linear relationships between the fraction online and measured voltage or frequency. From these measurements, the fraction of online DG is derived from the two curves in Figure 6. The recovery fraction defines what fraction is permitted to recover after the frequency/voltage recover. While this model generates a dynamic response, it heavily obscures the temporal ride-through aspects, such as momentary cessation and recovery time, and the individual responses of DG located on a voltage diverse feeder, with a variety of control schemes originating from different manufacturers, which can heavily impact the collective system response [22], [31].
Two sets of ride-through parameters are used for these simulations, which are shown in Table 4 and correspond to Figure 6. The first set was taken from the WWSIS study, where the parameters were chosen to ignore any potential DG loss due to abnormal frequencies or voltages. Recent work by Pacific Gas and Electric (PG&E) [33] used the second set of parameters. For this work, two simulations, one for each parameter set, were performed for each defined RAC controllable fraction. The WWSIS parameters permit an examination of the frequency nadir (lowest frequency following generation loss) without incorporating the dynamics lost DG, while the PG&E parameters show how the cascading impacts related to DG loss influence the system response.
Note that for the PG&E parameters, the DG lost due to under-or over-frequency is unable to recover. Here, lost indicates DG that goes offline during the simulation, and recovered indicates DG that returns online before the end of the simulation.

V. RESULTS
The following are the results for the control schemes implemented in the LSP and HS cases. Each control scheme was performed at 10%, 15%, and 20% controllable RAC fractions, as well as each of the ride-through criteria, WWSIS and PG&E, as outlined in Table 4. A total of 38 (6 control schemes x 3 controllable fractions x 2 ride-through criteria + 2 Palo Verde contingencies) PSLF simulations were run.A time series of the system frequency is presented as the primary result, with all six relevant traces (or seven with the Palo Verde frequency for the step-to-max schemes) on each plot. The WWSIS ride-through results are shown with solid traces, and the PG&E ride-through with dashed traces. The trace color (green, blue, and magenta) corresponds with the controllable fraction (10%, 15%, and 20%, respectively). Where applicable, the frequency load shed limit is shown as a dotted red trace. Additional result analysis, i.e. VS/VSDG scores, is provided when insightful.

A. LIGHT SPRING
The system frequency traces for the LSP step-to-max simulations are shown in Figure 7. At 10%, no DG is impacted and therefore the two traces for the two different ride-through criteria are identical. The nadirs for the 15%-and 20%controllable fractions with WWSIS ride-through are below the double Palo Verde unit trip contingency. For the 20% case, the frequency deviates below 59.5 and exhibits oscillations upon recovery. Two frequency relays were triggered, removing a small amount of generation; the frequency oscillations indicate potential system instability. Table 5 provides the DG loss statistics for these three simulations, showing that there was an impact only for the 20% controllable case (a slight loss of 55 MW for the WWSIS ride-through), with a full recovery before the simulation ends.
For the case of 15% and 20% with the PG&E ride-through, all of the DG (7,310 MW) is lost, with none recovering. This loss in generation is nearly three times greater than the double-unit Palo Verde contingency. With a load increase of 3.5 and 4.7 GW in the 15% and 20% cases, respectively, followed by a loss of 7.3 GW of generation, very large frequency deviations are expected. This loss is evident in the frequency traces in Figure 7, as both traces show a precipitous drop in frequency occurring a few seconds after the RAC step to max. Thus, this is a cascaded impact with the RAC step to max influencing the system DG. The 20% case dips into under frequency load shedding territory. Both simulations diverge.  Light spring step-to-max-to-off results are similar to the stepto-max case, due to no RAC initially online.

B. HEAVY SUMMER 1) STEP OFF
For the step-off control scheme in the HS case, the controllable fraction of RAC is simply turned off and a frequency rise results, with the peak frequency growing with the magnitude of load reduction. The frequency traces are identical to the first 40 seconds of step off to max simulation in Fig. 10, and are not reproduced here. The frequency deviations are well within any frequency load shedding criteria. No DG is impacted during these simulations, and therefore the PG&E ride-through results perfectly match the WWSIS ride-through results (i.e., 0 MW lost and 0 MW recovered for all six simulations).

2) STEP TO MAX
The HS step-to-max simulations present a few unexpected frequency responses, although none deviate significantly from the nominal. For the WWSIS ride-through, the expected decrease in the nadir results with growing controllable fractions. At 20%, the frequency response is nearly the same as that for the double Palo Verde unit trip contingency, which is reasonable considering that the load increase of 3,400 MW is near the generation loss amount of 2,750 MW. It is likely that the spatial dispersion of the RAC is partly responsible for why the 20% nadir is not lower (even though it exceeds the  Palo Verde contingency by 650 MW), as was discovered in the WWSIS study itself when comparing a dispersed DG trip versus a double Palo Verde unit trip at the same power loss [20]. The WWSIS ride-through criteria results in a maximum loss of 105 MW of DG with 95 MW recovering, due to voltage deviations.
The PG&E ride-through yields a frequency nadir that is actually higher for the 10% simulation, and in fact the frequency responses are inverted for the 15% and 20% simulations. Here, the VS metric can help to understand what is occurring. The low-voltage (below 0.95 p.u.) VS and VSDG scores are presented in Table 7. Because these six simulations are the same time length, the same controllable fraction VS and VSDG scores are comparable.
The PG&E score is larger, indicating that either a larger number of buses experienced similar voltage deviations, or a similar number of buses experienced larger voltage deviations. The larger VSDG scores in Table 7 correlate with larger losses in DG as shown in Table 6. The growing scores indicate much larger and/or more time-persistent voltage depressions as compared to the WWSIS simulations. Figure 9 is a distribution of voltages for buses with voltages below 0.95 p.u., for the PG&E 20% controllable case. The black trace represents the quantity of buses with voltages below 0.95 p.u. captured by the distribution, whereas the red scale traces provide percentiles of the distribution. It is evident that more than 600 buses (out of 21,000 total) have voltages below 0.75 p.u. for more than 40 seconds. We see that as the simulation continues, the voltages improve VOLUME 8, 2021  significantly. Further investigation shows that the total static load in the system, which is voltage dependent, decreased significantly following the large loss in DG, which occurred just after the RAC step. Additionally, a large quantity of motor load was temporarily tripped offline because of low voltages. This explains the frequency rise, because although the RAC stepped up considerably, and generation was lost, the total decrease in load (RAC -(static + motor load)) exceeded this DG loss, resulting in the frequency rise. Thus, the RAC step to max causes an initial voltage decrease, causing a loss in DG, which causes a further voltage decrease that eventually results in a net decrease in load because of the voltage-dependent static load and motor tripping. The WWSIS study noted that a loss in real power DG comes with a fourfold increase in reactive power load [20]. Thus, such a significant voltage depression is not surprising following the loss of 1.7 GW of DG.

3) STEP TO OFF TO MAX
If all controllable RAC units are off, and then stepped on, the resultant load increase is roughly double that of turning all controllable RAC units on with half online at the start of the simulation. The system frequencies are shown in Figure 10. The first 40 seconds are identical to the step-off results, whereas the frequency decline following the step to max is larger than just the step-to-max simulations. For the PG&E simulations, the DG losses are similar between the step-tomax, and step-off-to-max simulations for the 10% and 15% cases, although far less recovers in the 15% case. For the 20% case, all DG is lost at 9,680 MW and the system separates,  as evident by the sharp (ambiguous) frequency increase in Figure 10.

C. RESULTS SUMMARY
In the light spring conditions, when the total loading on the system is far less than the heavy summer case and the online RAC is nearly zero, large frequency swings occurred following controllable RAC steps. With the PG&E ride-through criteria, all DG was lost following frequency deviations beyond the ride-through criteria in the 15% and 20% controllable cases, which resulted in larger frequency deviations into frequency load shed territory followed by simulation divergence. For the heavy summer conditions, the frequency swings were smaller as compared to the spring conditions for two reasons: roughly half of the RAC was online at the start of the simulations, resulting in a smaller amount of a load swing, and the total load on the system was much higher. Regardless, it was found that large-voltage depressions occurred during steps to max, which caused large DG losses that exacerbated the system response. In a few of the cases, the voltage depressions were significant enough to cause an overall decrease in load because of the decrease in voltage-dependent load and some low-voltage motor tripping. For the step-to-off-to-max 20% controllable PG&E ride-through criteria simulation, all 9 GW of DG were lost and the system separated.

VI. CONCLUSION
Within the assumptions and constraints of the methods used, our results indicate that centrally-coordinated load can pose a catastrophic risk to the bulk power grid, and multiple bulk power system risks can result if this load is synchronized in a step-wise manner. In particular, it was found that even modest fractions of controllable residential air conditioner load, during times when this load is not in use such as the Spring, could create a load/generation imbalance far exceeding the WECC credible contingency, resulting in system instability. These risks may be compounded depending on the distributed generation ride through response, which might lead to widespread disconnection due to the local voltage and frequency response to the synchronized load perturbations.

VII. FUTURE WORK
There are numerous opportunities to expand and enhance this work. Motor inrush from air conditioners was not considered, but will have an additional deleterious impact under some conditions mentioned. The dynamic model should be expanded to include additional bulk power grid dynamics. Our work should be expanded to include cloud vulnerabilities associated with other grid-edge devices-such as electric vehicles, rooftop solar photovoltaics, and batteries. It is also worth considering potential telecommunication limitations on triggering this many internet-connected devices simultaneously, and how network latency might impact these results-if at all.
Lastly, this work only explored power grid impacts that could result from a centrally-controlling cloud. A more comprehensive risk assessment could define the relative benefits of mitigating vulnerabilities at these type of cloud coordinators, versus other potential threat vectors such as utility supervisory control and data acquisition systems. This type of assessment is left for future work.