Resilience of the Electric Grid through Trustable IoT-Coordinated Assets

a]Vineet J. Nair
b]Priyank Srivastava
c]Venkatesh Venkataramanan
d]Partha S. Sarker
d]Anurag Srivastava
e]Laurentiu D. Marinovici
f]Jun Zha
g]Christopher Irwin
h]Prateek Mittal
i]John Williams
h,1]H. Vincent Poor

a]Anuradha M. Annaswamy

[a]Department of Mechanical Engineering, MIT, Cambridge, MA 02139
[b]Department of Electrical Engineering, Indian Institute of Technology Delhi, India
[c]National Renewable Energy Laboratory, Golden, CO
[d]West Virginia University, Morgantown, WV
[e]Pacific Northwest National Laboratory, Richland, WA
[f]Larsen & Turburo Digital Energy Services, San Francisco, CA
[g]Office of Electricity, Department of Energy, Washington, DC
[i]Department of Civil and Environmental Engineering, MIT, Cambridge, MA 02139
[h]Princeton University

This paper provides a framework for achieving power grid resilience against cyber-physical attacks through the coordination of resources at the grid edge that are trustable and resilient. It is proposed that such coordination can be enabled through a suite of market operators suitably located in a distribution grid. We use this coordination to validate the mitigation of attacks of different levels of severity, with attack magnitudes that range from 5 to 40% of the total peak load. Both grid-connected and islanded cases are studied. In all cases, we show that grid resilience can be obtained through a combination of locally available flexible assets and reconfiguration of the grid topology.

A.M.A, V.N., P.S., V.V., and H.V.P. designed research; A.M.A., V.N., P.S., and V.V. performed research; V.N., P.S., V.V., and P.S.S. developed the analytical tools, V.N., L.M., J.Z., and V.V. performed validation;
A.M.A., V.N., P.S., V.V., and A.S. wrote the paper; C.I., P.M., and J.W. provided critical inputs and valuable guidance.
The authors declare no conflicts of interest.

1To whom correspondence should be addressed. E-mail<EMAIL_ADDRESS>
The electricity grid has evolved from a physical system to a cyber-physical system with digital devices that perform measurement, control, communication, computation, and actuation. The increased penetration of distributed energy resources (DERs) that include renewable generation, flexible loads, and storage provides extraordinary opportunities for improvements in efficiency and sustainability. However, they can introduce new vulnerabilities in the form of cyberattacks, which can cause significant challenges in ensuring grid resilience.
We propose a framework in this paper for achieving grid resilience through suitably coordinated assets including a network of Internet of Things (IoT) devices. A local electricity market is proposed to identify trustable assets and carry out this coordination.
Situational Awareness (SA) of locally available DERs with the ability to inject power or reduce consumption is enabled by the market, together with a monitoring procedure for their trustability and commitment. With this SA, we show that a variety of cyberattacks can be mitigated using local trustable resources without stressing the bulk grid. The demonstrations are carried out using a variety of platforms with a high-fidelity co-simulation platform, real-time hardware-in-the-loop validation, and a utility-friendly simulator.

This manuscript was compiled on August 27, 2024

The electricity grid is going through a rapid transformation in an effort toward deep decarbonization. Large synchronous generators powered by fossil fuels such as oil, natural gas, and coal are being phased out in favor of solar and wind-based generation. While the latter enables the necessary move towards a reduced carbon footprint, it brings two major challenges in ensuring reliable and resilient delivery of electricity to the end-user. The first of these is the temporal signature of these renewables  – the amount of generation varies with time, both in terms of intermittency and uncertainty. The second is that these are distributed and large in number. A strong enabler of the scale of the DERs is IoT, which denotes a network of physical devices such as water heaters (WHs), air-conditioners, and electric vehicles (EVs), as they enable automated and fast operation of various loads. And their pervasiveness brings in complexities of heterogeneity, decentralization, and scale. In order to ensure the reliability of the grid despite these challenges, a precise coordination of these DERs, both in space and time, has to be carried out. In particular, power balance of generation and consumption has to be ensured at all locations and at each instant. These challenges are being overcome using a pervasive cyber layer that senses, communicates, coordinates, and enables the requisite power injection and consumption throughout the grid.

In addition to reliability, an essential property of the electricity grid is its resilience [1]. This central property, which denotes the ability of the grid to withstand and recover quickly to supply critical loads following a major disruption, such as an outage, a natural calamity, a cyberattack, or a cascading failure, is paramount, even with increased penetration of DERs. In this context of ensuring resilience, the very transformations that enable deep decarbonization, including the development of cyber-grid infrastructure, adoption of IoT devices, use of dynamic renewable energy sources, and increased electrification of transportation, could also introduce new vulnerabilities. Cyberattacks can disclose, deceive, or disrupt crucial information, thereby causing significant damage, ranging from small outages to brownouts and blackouts.  Recent reports [2, 3, 4, 5]
indicate the ubiquity, ease, and scale of cyberattacks on sensitive industrial environments including supervisory control and data Acquisition (SCADA), operational technology (OT), and industrial control systems (ICS), underscoring the importance of ensuring resilience to such adversaries.

By and large, most of the information for power grid operations flows through utility-controlled communication networks which are more reliable and resilient than commercial networks, and utilize commercial telecommunications services for other informational needs such as accessing the internet and communicating with customers. Such a tight separation is challenged by the increased information flow which becomes necessary with a stronger presence of a cyber-layer, which in turn is necessitated due to increased coordination and automation at the grid edge. What have remained as tight closed systems thus far, may have to relax their boundaries, introducing complexities in the underlying communication. While air gaps and protections will always be important and included, imperfect protections are inevitable as complexity increases. With increased penetration of instrumentation and automation, motors and generators may be manipulated by adversaries to open and close at will. Another point to be noted is that with increased complexities due to intermittent and uncertain generation and consumption, utilities alone cannot cater to all needs; public and private partnerships may be necessary. It is therefore extremely important to design an appropriate cyberinfrastructure that ensures that the lights stay on, despite increased communication, which may be between disparate stakeholders. The focus of this paper is on such a distributed decision-making framework.

Given the size and complexity of the problem of cyberattacks, providing a complete resilience framework for the entire power grid is a tremendously difficult task. In this paper, we propose a first step, of providing SA to the grid operators in a distribution grid, with SA corresponding to the knowledge of local DERs in terms of their location and the amount of power generation that they are able to provide, as well as a resilience score (RS) that the operators can make use of to provide resilience. We argue that this first step, of providing SA, is enabled through a local electricity market (LEM) structure that consists of operators at different voltage levels in a distribution grid. This market structure is proposed to be local, across the distribution grid, electrically co-located with primary and secondary circuits, with operators scheduling all DERs at the corresponding nodes in a given region. The market will also include IoT-coordinated assets (ICAs), with the assumption that the ICA will have computing capability and the ability to exchange information. The overall framework, EUREICA (Efficient, Ultra-REsilient, IoT-Coordinated Assets), is the innovation in the proposed cyberinfrastructure, and will be shown to lead to SA made available to operators placed hierarchically at various locations, thereby providing an important first step in ensuring resilience.

LEMs have been addressed in several studies including [6, 7, 8, 9, 10, 11], with real-field implementations beginning to be reported [12, 13], all of which show the feasibility of a local market structure, and its advantages compared to alternate solutions that are designed to encourage full participation of DERs [14, 5]. The LEM structure that we propose in this paper builds on that in [6]. The resilience of the electricity grid to cyberattacks has been explored in a very large number of studies (see [3, 15, 16, 17] and references therein), with new results appearing continuously. Broadly, these approaches can be categorized into detection and isolation of the attack [18], prevention of the attack, and resilience in the presence of attacks. For large-scale attacks such as those described in [4, 19, 20], these methods are inadequate; it may be near-impossible to identify the attacker but rather that an attack has occurred. Prevention of the attack can be enabled through varying levels of access and authorization [21] and monitoring, isolation, and protection at the component level [4]. However, as the scale, location, and number of IoT devices in particular, and DERs in general grow, it becomes exceedingly difficult to completely prevent attacks. Ensuring resilience, especially in the face of large-scale attacks, for a large-scale system such as the electricity grid, is exceedingly difficult; current literature has either focused on systems at a small scale or with low levels of renewables. The EUREICA framework that we propose will provide SA that detects that an attack has occurred, and with this SA, deploys trustable ICAs in order to mitigate the impact of the attack, and ensure grid resilience through a distributed decision-making strategy.

The distributed decision-making in EUREICA is enabled through an LEM. The same market structure [6], which has been shown to lead to grid reliability [22] and provide grid services such as voltage support [23] in addition to overall power balance, is demonstrated in this paper to ensure grid resilience against cyberattacks using local trustable DERs. In particular, the results will show that local resilience is attainable through SA of locally available ICAs that have the ability to inject power or reduce consumption as well as a procedure for monitoring their trustability and commitment. The demonstrations are carried out using a variety of platforms such as (i) GridLAB-D which enables the simulation of distribution grids with high fidelity, (i) the advanced research on integrated energy systems (ARIES) platform that includes a real-time digital simulator (RTDS) and enables hardware-in-the-loop (HIL) validation, and (iii) General Electric's advanced distribution management system (ADMS) [24], distribution operations training simulator (DOTS), and DER integration middleware (DERIM).

The problem statement is laid out in <ref>, and the local market structure together with the SA and Resilience score (RS) it will facilitate, is described in <ref>. Various use cases that constitute the various attack surfaces on the grid are described in <ref>. The main results are presented in <ref>, followed by a summary in <ref>.

§ PROBLEM STATEMENT

A typical path of electricity delivery to end users traverses generation, transmission, and distribution. Distribution substations connect to the transmission system and gradually step down the voltage from 44kV or higher to 33kV (denoted as a primary network), then down to 11.2/4.6kV (denoted as a secondary network), and further down to 110V or 220V, depending on the specific region in the world. While the 20th century witnessed distribution systems operating as simple distribution lines as vehicles for sharing the electricity from transmission networks, today’s distribution systems are increasingly becoming heavily integrated with distributed energy resources, that correspond to resources that are located closer to the load, including renewable generation, some of which may be behind the meter [5], batteries, and flexible consumption units. This in turn is causing distribution systems to become more independent, and to be required to take on increased responsibilities of services such as grid reliability and grid resilience.  Other examples of DERs are distributed photovoltaics (DPVs) like rooftop solar systems, combined heat and power plants, electric vehicles, and diesel generators. DERs vary in size, from DPV systems that range between 1-1000kW in size to larger ground-mounted solar farms that range to several MW. With technological advances in power electronics and associated smart inverters as well as protection systems, fewer restrictions are being placed on the size and locations of the DERs, providing an opportunity for them to play stronger and more central roles in grid reliability and resilience [25].

Over the past years, DERs have been shown to be increasingly useful in providing key grid services such as volt-var control [26]. The central idea in these explorations is that key information is exchanged, in a distributed manner, between suitable individual components in the primary and secondary networks, coordinated both in space and time, thereby allowing local control over power injection and reduction of load at key locations and instants. Such a correct operation of the complete distribution network is predicated on this key information reaching the recipients in a secure manner. This sets the stage for malicious attacks that can disconnect and disrupt the overall grid by impairing key components.

Several attacks on power systems have been recently reported [27, 28, 29, 30, 31, 32, 19, 4] on the central control systems, key nodes in the distribution grid, or at the devices at the end-user level. Those at the device end, denoted as MadIoT (Manipulation of Demand via IoT) attacks, correspond to a botnet at a secondary network node that causes the corresponding load to change abruptly. If this node corresponds to a high-wattage device, and the attack is coordinated through malware that simultaneously corrupts a large number of these devices, an argument can be made that it can cause frequency instabilities, line failures, and subsequently a severe disruption on the overall power grid. Building on the results in [31, 32], the results in [19] show that even with realistic load profiles, a strategically coordinated attack can show a better success rate than in [31, 32] requiring fewer compromised IoT devices without triggering well-established protection systems. The well-known attack studied in [4] on the other hand is at the central control system level, which was a well-planned strategic attack that led to a power outage affecting 250,000 customers over a significant period of time. The question we address in this paper is: How can we use a cyberinfrastructure with IoT-Coordinated Assets (ICA) to support grid resilience against cyber-attacks?

The specific approach that we propose to circumvent the anomalous scenario consists of two steps: (1) Enable improved visibility over the grid and net power injections available at various nodes through a hierarchical market structure with operators at the primary network and secondary network nodes; (2) Enable the market operators to determine an RS computed through monitoring of various features of the communication network. Steps (1) and (2) together provide SA to the grid operators (as shown in <ref>). Our central thesis is that through this SA, operators can determine that an attack has occurred and take appropriate steps to mitigate the impact of the attack in a timely manner. The system operators and resilience managers are suitably co-located with the electrical assets so as to respond quickly through a distributed decision-making framework. The framework therefore avoids the computational pitfalls of a centralized architecture while still underpinned by a substrate of communication, sensing, and actuation. The overall solution is also well-placed to integrate with the existing grid operational and market structures, helping accelerate its adoption in the field.

Summary of attack scenarios and use-cases, LA = load alteration attack, DG = distributed generator attack.

Attack Number     Attack type     Attack surface     Grid connection     Power flow model     Grid model      Scale of attack [kW]

1a                       LA                     PMA                   Grid-connected               Current injection             Unbalanced, 3-phase      36

1b                        DG                       PMA                   Grid-connected               Current injection             Unbalanced, 3-phase     45

1c                        DG                       SMA                         Grid-connected               Current injection             Unbalanced, 3-phase      157

2a                        DG                       PMA                   Grid-connected               Branch flow                   Balanced, single-phase     261

2b                        DG                       PMA                   Grid-connected               Branch flow                   Balanced, single-phase     650

3                      DG                       PMA                        Islanded                     Current injection             Unbalanced, 3-phase     2500

§ OUR APPROACH: LOCAL ELECTRICITY MARKETS

A Hierarchical LEM for a Distribution Grid. The resilience infrastructure utilizes the dual market layer consisting of PM-SM.

LEM co-located with distribution grid. This shows a primary and secondary feeder distribution network based on the modified IEEE-123 node test case.

In order to provide visibility into a distribution grid, we propose an LEM that is hierarchical (see <ref>) in nature and electrically collocated with a radial network. The starting point for the overall LEM is a distribution system operator (DSO) that oversees several substations in the distribution grid with multiple primary and secondary markets downstream and acts as their representative in its transactions with the wholesale electricity market (WEM). The substation connects the distribution grid to the high voltage transmission grid at the point of common coupling (PCC) (node 150 in <ref>). The dual-layer downstream of the substation, consisting of a primary market (PM) and a secondary market (SM), is the core of the resilience infrastructure of this paper. The PM consists of Primary market operators (PMOs) and Primary market agents (PMAs). The PMA at each of the primary nodes either own a DER at a primary feeder node or are aggregators representing DERs at the secondary feeder level and below.  In the latter case, the PMA plays a second role as an SM operator (SMO) and coordinates with SM agents (SMAs). The PMO, PMAs/SMOs, and the SMAs are located at the coupling between the substation and the primary feeder, primary feeder nodes, and secondary feeder nodes, respectively (see <ref>). The PM and SM operate at medium- and low-voltage levels, respectively. The DSO supervises the entire distribution grid - our proposal could be viewed as an expansion of the current responsibilities of a DSO, which comprise grid maintenance and grid reliability, to include market oversight and regulation as well. In this sense, the role of the DSO would be analogous to that of existing independent system operators for transmission grids [33]. SI <ref> illustrates a possible implementation of the dual-layer LEM for the city of Boston, MA, USA.

Formally, we define SA at an operator $x$ as the tuple

\begin{align}\label{eq:sa}
\text{SA}_x=\{\text{ICA}_x,\text{RS}_x\},
\end{align}

where ICA$_x$ stands for the IoT-coordinated assets and denotes the generator and/or consumption flexibilities of DERs under the purview of agent $x \in \{\text{SMA}, \text{SMO}\}$, and RS$_x$ denotes their resilience scores, to be defined in <ref>. We will show, in what follows, that $\text{RS}_x$ can be determined based on the asset's market performance and security against possible attacks.

We now show how the LEM made of the PM-SM layers will allow the computation of SA. The operation of a distribution grid is challenging due to its scale, complex topology, and presence of various active  DER assets and fixed load nodes. We separate this complex task by having the PM focus on grid-specific costs and constraints while the SM focuses on consumer-centric costs and constraints. We assume that PM and SM clear once every 5 minutes and 1 minute respectively. The main reason for this separation of timescales is that the SM typically needs to monitor fewer assets than the PM, and is closer to DER devices (such as rooftop solar and batteries) and therefore may need to operate at a faster timescale than a PM. The starting point for both markets is the submission of bids by the corresponding agents. Bids for the SM are submitted by the SMAs exogenously, whereas bids for the PM are computed by the SMOs via the SM. Accordingly, we describe the operation of the SM before going into the details of the PM.

§.§ Secondary market

The operation of the SM consists of three sequential stages: bidding, clearing, and monitoring. We denote $\mathcal{N}$ to be the set of all SMOs in the network and $\mathcal{N}_i$ to be the set of all SMAs under a given SMO $i \in \mathcal{N}$.

§.§.§ SM bidding

During the bidding phase, each SMA $j \in \mathcal{N}_i$ submits  a bid $\mathcal{B}_j^{iS}$ defined as

\mathcal{B}_j^{iS} =\{P_j^{i0},Q_j^{i0},\underline{P}_j^i,\underline{Q}_j^i,\overline{P}_j^i,\overline{Q}_j^i,\beta_j^{iP},\beta_j^{iQ}\}. \label{eq:sm_bid}

$P_j^{i0}$ and $Q_j^{i0}$ denote the baseline active and reactive injections of SMA $j$, along with the upward ($\overline{P}^i_j,\overline{Q}^i_j$) and downward flexibility ($\underline{P}^i_j,\underline{Q}_j^i$). $\beta_j^{iP}$ and $\beta_j^{iQ}$ denote the disutility parameters associated with providing active and reactive power flexibility, respectively. It should be noted that Bid ${B}_j^{iS}$ requires SMA $j$ to have a realistic estimate of its energy profile for the next 1 minute. Since it is not always trivial to predict future power availability, agents deploy a decentralized federated learning (FL)-based framework [34] to determine their bids. Using FL helps ensure that the privacy of the participating agents is preserved and the computational aspects of the prediction algorithm scale well as the number of agents increases. Further details on the FL implementation can be found in SI <ref>. SI <ref> summarizes details of the overall LEM.

§.§.§ SM clearing

Once the SMO $i$ has received bids from the participating SMAs, it clears the market with active and reactive power injection setpoints ($P^{i*}_j,Q_j^{i*}$) and the corresponding retail tariffs ($\mu_j^{iP*},\mu_j^{iQ*}$).
In addition, the SMO also solves for the optimal flexibility ranges ($\delta P_j^{i*},\delta Q_j^{i*}$) for $j \in \mathcal{N}_i$.
The SMO clears the markets with the following objectives: (O1)  maximization of aggregate resilience $f_i^1$, (O2) minimization of the net cost to the SMO, $ f_i^2$, (O3) maximization of total flexibility $f_i^3$ that the SMO can extract from all its SMAs and (O4) minimization of the disutility of the SMAs $f_i^4$, arising from flexibility provision (see SI <ref> for all details). This gives rise to a multiobjective constrained optimization problem:

\begin{align}
& \min_{\textbf{y}^S_i} f_i^S=\{f_i^1, f_i^2, f_i^3, f_i^4\}^\top \\
& \text{s.t.} \underline{P}^i_j + \delta P^i_j \leq P^i_j \leq \overline{P}^i_j - \delta P^i_j \; \forall j \in \mathcal{N}_i, \; \forall \; \text{constraints} \\
& \underline{Q}^i_j + \delta Q^i_j\leq Q^i_j \leq \overline{Q}^i_j - \delta Q^i_j \\
& \delta P^i_j, \; \delta Q^i_j \geq 0, 0 \leq \mu_j^{iP} \leq \overline{\mu}^{iP}, 0 \leq \mu_j^{iP} \leq \overline{\mu}^{iQ}  \;  \\
& \sum_{t_p}\sum_{t_s} \sum_{j \in \mathcal{N}_i} \mu_j^{iP}(t) P_j^i(t) \Delta t_s \leq \sum_{t_p} \mu^{iP^*}(\hat{t}_p) P_i^*(\hat{t}_p) \Delta t_p  \label{eq:budgetP} \\
& \sum_{t_p}\sum_{t_s} \sum_{j \in \mathcal{N}_i} \mu_j^{iQ}(t) Q_j^i(t) \Delta t_s \leq \sum_{t_p} \mu^{iQ^*}(\hat{t}_p) Q_i^*(\hat{t}_p) \Delta t_p \label{eq:budgetQ} \\
& \sum_{j \in \mathcal{N}_i} P_j^i(t_s) = P^{i*}(\hat{t}_p), \quad \sum_{j \in \mathcal{N}_i} Q_j^i(t_s) = Q^{i*}(\hat{t}_p) \label{eq:PQbalance}
\end{align}

The constraints include capacity limits and operational bounds on SMA injections (including flexibilities), budget balance constraints, price ceilings, and lossless power balance. Note that, here we do not account for all the power physics, these will be considered in the PM. The decision variables consist of the P and Q injection setpoints as well as retail tariffs for each SMA i.e. $\textbf{y}^S_i = \{\textbf{y}^{iS}_j\} \; \forall j \in \mathcal{N}_i$ where $\textbf{y}^{iS}_j = [P_j^i, Q_j^i, \delta P_j^i, \delta Q_j^i, \mu^{iP}_j, \mu^{iQ}_j]$. We note from the choice of $f^1$ that the solution of <ref> requires the resilience scores $\text{RS}^i_j$. This is assumed to be communicated by the secondary resilience manager (SRM) to the SMA, the details of the SRM are addressed in the next section.

In general, the optimization problem in <ref> has multiple solutions known as Pareto points, with each solution prioritizing different objectives.
However, since the objective functions have different units, instead of finding the Pareto solutions, we use a hierarchical approach proposed in [23] where the SMO optimizes one objective at a time in descending order of importance. While optimizing the subsequent objective functions, additional constraints on the degradation of prior objectives are added to the optimization problem (see [6] for details). The cleared market schedules $\textbf{y}^{S*}_i$ are sent by the SMO to their corresponding SMAs, as well as to their SRM.

§.§.§ SM monitoring and resilience scores

The final stage in the SM is monitoring. During the market operation, the responses of each SMA $j$ to the market schedules, in terms of its actual DER injections $\hat P_j^i$ and $\hat Q_j^i$ are suitably monitored by its corresponding SRM.

Sequence of communication steps and events leading to SA with an LEM. The red arrows indicate the entities and communication links that would be affected by an attack. A more detailed diagram can be found in SI <ref>.

In addition to the market operators, we also propose the addition of two new entities, which we denote as the primary resilience manager (PRM) and the SRM, both of which provide grid functionalities, with the PRM located at the primary circuit level and the SRM at the secondary level, as shown in <ref>. With the market clearing providing the first step of awareness in the form of power available at each of the nodes at the secondary and primary level, the PRM and SRM monitor the actual injections, determine corresponding scores of commitment, trustability, and resilience (to be defined below), and communicate them using protected channels to the PRM. Not only do these entities enable a separation between grid-specific decision-making from market-specific decisions, but they also provide a pathway for mitigating the impact of any attacks that can occur through the addition of local resources, as will be shown in the following sections.

In this monitoring stage, the SRM assigns each SMA an RS that is updated constantly based on its performance in the market and susceptibility to being compromised. The RS is a weighted combination of its commitment score (CS) and trustability score (TS). Formally, for an agent $j$

\begin{align*}
\text{RS}_j= \alpha \text{CS}_j + (1-\alpha) \text{TS}_j,
\end{align*}

where $\alpha \in (0,1)$ is a parameter chosen by the SRM.
The CS and TS are defined below.

* Commitment score (CS).
The CS of an agent measures its reliability in executing its cleared schedules and is updated at every SM clearing instance.
The first step in updating $CS_j$ for each agent $j \in \mathcal{N}_i$ is the computation of any relative deviation between the cleared schedule and its executed value over the past market period. A moving average is then computed to account for the past performance. Finally, a min-max normalization across all the SMAs is performed to keep $CS_j \in [0,1]$ for all $j$ (see SI <ref> for further details).

* IoT Trustability score (TS).
The TS captures the possibility of the agents (or the devices underneath them) being compromised. TSs are computed using an FL-based anomaly detector and like the CS, past values are used again to compute a weighted moving average. However, unlike the CS, which solely depends on power injections, the TS is a cyber-power metric [35] that also takes into account the associated cyber information, e.g., packet length, arrival time, and communication protocols, etc. (see SI <ref> for further details).

In summary, the overall operation of the SM allows the generation of the schedules $\textbf{y}^{S*}_i =[P_j^{i*},Q_j^{i*},\delta P_j^{i*},\delta Q_j^{i*},\mu_j^{iP*},\mu_j^{iQ^*}] \; \text{and} \; \text{RS}^i_j \; \forall \;\text{SMAs} \; j$, all of which provide $\text{SA}^i_j$ for the SRMs corresponding to all SMAs $j$ at the primary node $i$. Similar measures of the resilience of large-scale networks to attacks can be found in [36].

§.§ Primary market

PM transactions happen between the PMO and the PMAs. Similar to the SM, the operation of the PM also consists of bidding, clearing, and monitoring.

§.§.§ PM bidding

Here, we describe the link between the PM and SM. As noted previously, the PM is cleared every 5 minutes while the SM operates more frequently at 1-minute intervals. Before each PM clearing, the SMOs (or PMAs) aggregate the schedules and cleared flexibilities of all their SMAs resulting from the most recent prior SM clearing (at the lower level) to submit their flexibility bid to the upper-level PM. All market bidding and clearing for both the SM and PM are based on forecasts (assuming perfect foresight) and for the very next period. The complete bid submitted by each SMO $j \in \mathcal{N}$ into the PM $\mathcal{B}_i^{P}$ defined as:

\begin{align}
\mathcal{B}_i^{P} =\{P_i^{0},Q_i^{0},\underline{P}_i,\underline{Q}_i,\overline{P}_i,\overline{Q}_i,\alpha_i^P, \alpha_i^Q, \beta_i^P,\beta_i^Q\}. \label{eq:pm_bid}
\end{align}

This includes the nominal power available, the lower and upper bounds, net generation cost coefficients, and flexibility disutility parameters (see details in SI <ref>). We note that standalone PMAs such as a large industrial facility, a community solar farm, or an EV charging station, may also be present. In this case, the PMA would directly bid into the PM on its own instead of aggregating over SMAs.

§.§.§ PM clearing

At each PM clearing instance, an optimal power flow (OPF) problem is solved to optimize the PMO's objective while satisfying all grid physics and network power flow constraints. For simplicity, in this work, we consider the cost functions of all the PMAs (or SMOs) to be quadratic. The objective function
utilized is a weighted linear combination of (i) maximization of social welfare, (ii) minimization of total generation costs, and (iii) minimization of electrical line losses (see SI <ref> for details of these functions). The total cost includes paying the locational marginal price (LMP) $\lambda$ for importing power from the transmission grid at the PCC, as well as the payments to local generator PMAs that provide net positive injections into the PM. We divide by suitable base values to convert all quantities to per unit (between 0 and 1 p.u.). Thus, it is reasonable to combine all the terms into a single objective function using a simple weighted sum.

With the objective function thus defined, the constraints are determined by the choice of the power flow model used to describe the system. Since the original alternating current OPF (ACOPF) is inherently nonconvex and NP-hard, we need to convexify the problem to make it more tractable. In this study, we considered two different approaches for this convexification. The first is a branch flow (BF) model or nonlinear DistFlow [37] based on a second-order conic program (SOCP) convex relaxation - this is a simpler implementation that is valid for radial and balanced networks. The second is a linear current injection (CI) model [38] based on a McCormick envelope convex relaxation that is more generally applicable to unbalanced and meshed grids common in distribution systems (in addition to radial, balanced), although this adds some overhead due to certain pre-processing steps needed.

We deployed both of these models for different use cases considered in this paper, as shown in <ref>. Further details of the BF and CI approaches are provided in SI <ref> and SI <ref>, respectively. The exact set of decision variables $\textbf{y}^P_i$ for each PMA $i$ differs slightly depending on the OPF model used. Both models solve for the nodal power injections and voltages. However, the BF model only considers branch currents while the CI model also considers nodal current injections. BF also models all variables as only having a single phase while the CI models these as three-phase, complex phasor quantities. For simplicity, we have only included the single-phase formulations in the paper thus far. However, these can easily be extended to the complex three-phase representation by simply modifying all variables to 3-dimensional complex vectors instead of scalars. A three-phase extension of the SM optimization is also given in SI <ref>. Thus, the decision vectors for the BF model are given by $\textbf{y}^{P,BF}_i = [P_i^G, Q_i^G, P_i^L, Q_i^L, v_i, I_{ik}] \; \forall i \in \mathcal{N}, \; (ik) \in \mathcal{E}$ and for the CI model, $\textbf{y}^P_i =[P_i^{\phi},Q_i^{\phi},V_{i}^{\phi,R},V_{i}^{\phi,I},I_{i}^{\phi,R},I_{i}^{\phi,I},I_{ik}^{\phi,R},I_{ik}^{\phi,I}]$, where phases $\phi \in \{a,b,c\}$ are the phases and $\mathcal{E}$ is the set of all network edges or branches.

§.§.§ Distributed optimization for PM clearing

Since the number of nodes (and hence the number of PMAs) in a primary feeder could be arbitrarily large, rather than using traditional centralized optimization solvers, we employ a distributed proximal atomic coordination (PAC) algorithm [39] to solve the OPF using peer-to-peer communication between the agents. This also helps preserve data privacy since each PMA only needs to exchange limited information with its immediate neighbors. A distributed approach also enables the PMAs to clear the market independently of the PMO, alleviates the communication burden, and reduces latencies since PMAs do not need to send all their data to a centralized entity, thus allowing for scalability. This is achieved by a process called atomization wherein the overall global optimization problem is decomposed into several local optimization problems called atoms for each PMA.
The constraints can also similarly be decoupled. However, certain network constraints also depend on other PMAs' variables. To deal with this we include additional coupling or consensus constraints to ensure consistency. We also used an enhanced variant of PAC known as NST-PAC that employs Nesterov (NST) acceleration and has enhanced privacy features by further masking the variables exchanged between atoms (i.e., the PMAs) [40]. After a sufficient number of iterations, both the PAC and NST-PAC algorithms provably converge to globally optimal and feasible solutions $\textbf{y}^{P*}= \{\textbf{y}^{P*}_i\}$ for each of the PMAs. Further details on PAC and NST-PAC are provided in SI <ref>. These cleared market schedules are communicated by the PMAs to their respective SRMs as well as to the PMO.

§.§.§ PM monitoring and resilience scores

During the actual market operation, the injections $\hat P_i$ and $\hat Q_i$ from the DERs at PMA $j$ are monitored by their SRM. These could be either from standalone PMAs or aggregated information from all the SMAs at a given PMA. The SRM also assembles resilience scores $\text{RS}_i$ for each PMA $i$. This is done through aggregation (via a weighted average) of $\text{RS}^i_j \; \forall j \in \mathcal{N}_i$. The RSs for standalone PMAs can also be directly computed at the SRM using their monitored injections. $\textbf{y}^P_i$ and $\text{RS}_i$ thus provide complete SA at each PMA node $i$. All SRMs send this information to the PRM so that the PRM has complete SA of all PMAs. This SA can then be used to redispatch the ICAs in both the PM and SM to mitigate the impact of various attacks. Further details on the mitigation strategy can be found in SI <ref>.

§.§ Reconfiguration paths

The final tool that we use in our proposed EUREICA framework is the determination of reconfiguration in the wake of islanding which can occur if an attack or natural disaster causes an entire section of the grid to be disconnected from the main grid. Reconfiguration can also In such cases, an algorithm that determines a self-sustaining operation of the islanded system, which is enabled by reconfiguration paths with suitable switch settings, is essential. We propose a reconfiguration algorithm (see SI <ref> for details) that considers power flow feasibility, available distributed generators (DGs), critical load, as well as RS information, to determine switching actions to restore specific sections of the distribution feeder. The reconfiguration paths will be determined based on the available amount of generation and the amount of critical load to be supplied, which is obtained through the SA provided by the EUREICA framework. In addition, the TSs are used at the secondary feeder level to intelligently disconnect non-critical loads, thus enabling the maximum restoration of critical loads. Once the feasible paths are determined for the optimal selection of loads, the RSs for all feasible paths are computed, and the most resilient path is implemented in the system.

§ USE CASES

In this section, we present use cases that illustrate how SA can be leveraged to ensure grid resilience in a distribution grid with a high penetration of DERs. We consider four different attack scenarios, all of which are motivated by the two large-scale attacks in  [4, 31] on power grids. Disruptive attacks are assumed to occur in the form of (a) a sudden loss of generation, and/or (b) a sudden increase in load, at multiple vulnerable locations. All use cases are simulated using an IEEE 123-node test feeder (see SI <ref>); extensions to more realistic and larger networks [41] can be implemented similarly.

§.§ Attack 1

In this attack, it assumed that a small percentage of generation or load resources at either the primary or secondary feeder level are compromised. In particular, it is assumed that these units are offline due to either an outage, natural calamity, or malicious cyber-attacker using elevated privileges to disconnect the units. In addition to the generation shortfall, it is assumed that the communication link between the market operators (PMO/SMO) and the resilience managers (PRM/SRM) is also affected by a denial of service (DoS) attack, which compromises the availability of a resource (see [42] for an attack which occurred on an sPower installation in Utah). Attack 1 draws inspiration from [4], where a malicious attacker used (i) elevated and unauthorized access to disconnect several resources, and (ii) severed communication links, to hamper operator visibility and response. While these attacks occurred at the transmission level, it is feasible that a similar impact can be had by targeting distribution grid entities, especially with the larger attack surface provided by grid-edge devices. Independently, it is possible that IoT load devices such as heating, ventilation, and air conditioning (HVAC) devices, WHs, EV chargers, or refrigerators may be attacked as well as noted in [19]. Elements of both of these types of attacks are explored here in two different cases, 1a and 1b.

§.§.§ Case 1a

In this case, the grid is assumed to be subjected to a sudden increase in load at the primary feeder level level (SMO or PMA) due to malicious agents. There are several large loads connected to the primary feeder such as commercial buildings or industries, and a malicious agent can manipulate the loads in these entities to affect the grid. Typically, the grid would rely on the margin provided by grid inertia to mitigate the effect of a sudden load increase. However, in a case where the grid’s resources are stretched, such as a cold snap or similar natural hazards, it is imperative that the grid-edge IoT resources be tapped to mitigate this condition. Examples of this scenario are already seen in operations, such as requests from grid operators in Alaska, Texas, and others in response to cold snaps. The operators requested customers to reduce their power consumption to support large critical loads such as chillers in hospitals. Furthermore, increased DER penetration will also lead to a loss of inertia, currently provided largely by large coal and gas plants. Case 1(b) details the performance of the proposed framework from this generation shortfall, even when the PRM does not have complete observability in the system.

§.§.§ Case 1b

Here, several generating resources are assumed to be unavailable at the primary feeder level (i.e. SMO or PMA). There are several scenarios that motivate this case – for example, in the case of several cloudy days in a row (affecting wind power), or unforeseen maintenance on generating units, the grid operates at a lower margin than under normal conditions. There is also the case of a malicious actor disconnecting generation resources. Here, the grid experiences a generation shortfall, and in combination with the DoS attack, the system operator (PRM) loses observability.

§.§.§ Case 1c

In this case, the grid is subjected to a sudden increase in load and/or corruption of distributed generation from the IoT devices, in a coordinated fashion directly at the secondary feeder. DER IoT devices will soon be operated via cloud-based service mechanisms that allow them to be controlled remotely. Thus, a sufficiently motivated malicious actor could gain control of a large number of these to suddenly reduce generation or increase load in a coordinated fashion. We simulated a case in which a large number of DGs (such as solar PV smart inverters) are attacked at the SMA level.

§.§ Attack 2

A larger-scale attack is assumed to occur at the distribution grid level in the form of several DGs being corrupted, causing them to go offline. The scale of this attack is assumed to be such that the impact is felt even in the transmission grid. We will explore how SA by the PRM and SRM helps mitigate this impact. Similar to attack 1.0, this use case combines elements of both [19] and [4]. The similarity to the latter is that the corruption is inserted in the form of outages of large DGs, while that to the former is that it introduces oscillations at the transmission level. For this purpose, we will utilize the well-known Kundur 2-area test system used to understand the transient and dynamic transmission-level impacts [43]. In particular, we will assume that there is an outage in one of the two areas (Area 2) that is load-rich, which introduces additional stress on the tie-line connecting the 2 areas (see SI <ref> for a diagram of the 2-area system).

§.§ Attack 3

The substation transformer is located at node 150, which is connected to the main transmission grid under normal operating conditions. However, under this attack, the distribution grid is islanded from the main grid at node 150. This could be due to a multitude of factors – such as wildlife tripping the transmission line from the substation to the distribution system, or a cyber-attack (i.e., integrity or disruption attack) that trips the circuit breaker from the main grid. With the increased SA introduced through our framework, we will demonstrate that the distribution system loads can be picked up in a coordinated fashion.

§ RESULTS AND DISCUSSION

In this section, we focus on a few specific attack scenarios 1a, 2b, and 3 described in <ref>, along with some results for attacks 1b and 1c. The remaining scenarios are elaborated in SI <ref> and SI <ref>. With the numerical simulation setup described in <ref>, we present details of how each of these attacks is mitigated using the proposed EUREICA framework. Note that for all attacks except attack 2.0, we use the mitigation strategy described in SI <ref>. For attacks 2a and 2b, we use the algorithm in SI <ref> instead. In addition to market simulations, we also validated our results using high-fidelity software at the Pacific Northwest National Lab (PNNL), LTDES, and the National Renewable Energy Lab (NREL). Technical details for each validation platform can be found in SI <ref>. For each attack, we focus on a specific subset of results to highlight the most pertinent findings and insights. Complete details including the remaining market simulation results for all attack scenarios can be found in SI <ref>, along with the complete set of validation results in SI <ref>.

§.§ Mitigation of attack 1a

We note that in attack 1a, loads are compromised leading to an increase in the power import from the bulk grid. It is also assumed that the communication from all SRM to the PRM is disrupted, while the communication from the PRM to theSRM remains intact. That is, the PRM loses observability but is still able to communicate the redispatch of the new coefficients to the SRM. We do not consider the case when such observability is not lost, a discussion of which is beyond the scope of this paper. With the redispatch, the PM-SM framework identifies all of the new trustable PMAs (through the SA computations described in SI <ref>), which will provide the injections needed to fully mitigate the attack, and the overall power balance is thus met at all points in the distribution grid.

The steps in mitigation are as follows: 10 SMO nodes are attacked, resulting in a total increase in load (generation shortfall) of 36 kW for the entire feeder as seen in <ref>. A large number of flexible load nodes across the entire feeder help with mitigation by curtailment and shifting as in <ref>. Flexible load curtailment at individual SMO nodes ranges from a minimum of 0.55 kW to a maximum of 7.8 kW reduction per primary feeder node - using a combination of resources like HVAC, WHs, batteries, and EVs to reduce the net load. There is a 123 kW decrease in power import after mitigation as seen in <ref>. The new SMO setpoints from the PM redispatch are then disaggregated amongst their SMAs during the following SM redispatch, with an example for SMO 77 shown in <ref>.

Net load at attacked nodes.

Feeder power import from main grid.

Effect of attack 1a and mitigation.

Curtailment of flexible loads for attack 1a mitigation.

Dis-aggregation of setpoint changes (from the PM) for SMO at node 77 across its 3 SMAs (in the SM) on phase B, after attack 1a mitigation.

The outputs from the PM-SM market framework were sent to the DERIM interface using which the effect on the total net load at the substation feeder head could be determined with the DERIM-ADMS-DOTS software platform (see SI <ref> for an overview of the validation process). It is clear from <ref> that without the intervention of EUREICA, the impact of the attack is a 37 kW jump in the feeder demand; in contrast with EUREICA, the feeder demand is cut by 94 kW. Moving further ahead from the attack timestep, the feeder net load eventually approaches back to the same value as if there hadn't been an attack. See SI <ref> for further details on this validation. See SI <ref> and SI <ref> for the other validation results using the HELICS and ARIES platforms, respectively.

LTDES validation of attack 1a in the DERIM-ADMS platform, showing total power import at the substation around the attack time at 13:00.

§.§ Attack 1b validation based on resilience

Here, we briefly highlight the effects of resilience scores on the mitigation of attack 1b, where a number of DGs are attacked. Full details on this attack can be found in SI <ref>. Here, we focus only on how the RSs of SMOs and SMAs influence which resources are utilized to mitigate the attack. The RSs of the flexible SMOs are plotted against their absolute and relative levels of net load curtailment in <ref> and <ref>, respectively. We see that in relative terms, the curtailment is generally distributed evenly to ensure that no single SMO is disproportionately affected. However, if the PMO does need to utilize more flexibility from certain SMOs, it generally calls upon more reliable ones with higher RSs. The absolute amounts of curtailment vary for each SMO based on their baseline load. This also holds while dis-aggregating SMO setpoints at the SM level, where the SMO allocates greater flexibility to SMAs with higher RSs, as seen in SI <ref>. See SI <ref> and SI <ref> for validation results of this attack using DERIM-ADMS-DOTS and HELICS, respectively.

Absolute curtailment and RS.

Relative curtailment and RS.

Distribution of absolute and relative amounts of load curtailment across the flexible net load SMOs, along with their corresponding resilience scores.

§.§ Contributions of SM and PM to attack 1c mitigation

Attack 1c is a more distributed attack where individual SMAs are attacked directly. Here, we show how both the SM and PM flexibility are needed to fully mitigate the attack. <ref> shows the contributions of the SM and PM toward attack mitigation. We see that for most of the SMO nodes, both the SM and PM flexibility play a significant role in reducing the net load compared to the post-attack case. At the SM level, we utilize the available upward flexibility of any SMAs with remaining online DGs and use the downward flexibility of all net load SMAs. At the PM level, we utilize the downward load flexibility of the SMOs (which are all net loads after the attack). Further details on this attack can be found in SI <ref>.

Contributions of SM and PM flexibility for attack 1c mitigation.

§.§ Mitigation of attack 2

Here, we describe the mitigation of two attacks at the primary feeder level that are relatively broader in scope, one is a medium-scale and the second is a large-scale attack. Both are disruption attacks where the attacker shuts down one or more of the large DGs in the network. We only consider a single primary market time step to study the effects of an instantaneous attack. Mitigation can use P dispatch from batteries, P and Q curtailment from flexible loads, limited P dispatch from PV, Q support from smart inverters (connected to PV and batteries), as well as conventional dispatchable fossil fuel sources like diesel generators. We only present results for the large-scale attack here, details on the medium-scale attack can be found in SI <ref>.

§.§.§ Large-scale attack 2b

Here, we adopted a top-down approach in emulating an attack and started with a Kundur 2-area transmission model, with the attack occurring in Area 2 (see SI <ref>) which consists of a load of 1767 MW. Noting that Area 2 can be broken down into 552 IEEE-123 feeders, each with approximately 3.2 MW, we assume that an attack that compromises about 650 kW of generation, occurs in each of these 552 feeders. This in turn corresponds to an overall shortfall of 359 MW at the transmission level. We introduced this 650 kW shortfall in the form of a generation loss at four nodes 25, 40, 81, and 94, in each of the 552 primary feeders. The only remaining SMO with significant generation capability is at node 67. With the same procedure as outlined in the previous scenarios, the use of our proposed EUREICA framework leads to the results in <ref>. In order to mitigate the attack, we need to leverage the upward generation flexibility of the remaining SMO 67 to increase its output injection after attack mitigation, while the net injections for all the other four attacked SMOs drop to zero as seen in the left plot. The right plot shows the new SMA schedules resulting from the revised SM clearing.

Mitigation of large-scale attack 2b.

However, due to the larger scale of the attack, re-dispatching the generator SMOs is no longer sufficient to fully meet the shortfall. Furthermore, as seen in <ref>, we are not able to utilize all the upward flexibility of the remaining online SMO 67 since its dispatch is limited by power flow constraints, on nodal voltages and line currents in particular. Thus, we also need to perform some shifting and curtailment of high-wattage flexible loads. These could include EVs and thermostatically controlled loads like HVAC and WHs. In addition, it could also involve some discharging of battery storage systems to reduce the net load. The distribution of net load reductions across the remaining SMOs is shown in <ref>, with a total decrease of around $14 \%$ as seen in <ref>. From <ref> we also see that the attack would have potentially increased the power import from the transmission grid by over $37 \%$, but the combination of increased local generation and load curtailment helps keep the imported amount almost the same as before. Further details on this attack can be found in SI <ref>.

Flexible load curtailment for large attack mitigation.

§.§.§ Effects at the transmission level

The overall impact of the generation shortfall and mitigation using EUREICA is simulated in the RTDS using a proxy where the individual feeders are not modeled, but the aggregated effect is studied at the transmission scale. A combined shortfall of 359 MW, corresponding to a simultaneous compromise and outage of 650 kW in all 552 primary feeders in Area 2 triggers a frequency event (see <ref>). Left unchecked, this can potentially lead to drastic load shedding or parts of the system being blacked out. To mitigate this situation, the power flow from Area 1 to Area 2 needs to be increased, which was observed in the RTDS, through the action of the governor system, which responds in the timescale of seconds, by increasing generation from the other generators present in the system proportionately based on a droop value. This increases the power flow from the generation-rich Area 1 to Area 2. However, changing the tie-line power flow creates a frequency imbalance, resulting in the system frequency oscillating, and settling at a lower/higher frequency, as shown in <ref>. With the EUREICA framework, the frequency mismatch is mitigated by suitably leveraging the flexibility of the remaining generation as well as demand response (DR) mechanisms from flexible loads at both the SMO and SMA levels (see <ref>). Once the governor response is completed and the system settles at a sub-optimal frequency, a combination of intelligent DR and generation redispatch in Area 2 facilitated by the EUREICA framework allows the system frequency to be restored to normal, ensuring grid resilience, avoiding system stress and increased operational costs.

Response without EUREICA; system settles at sub-optimal frequency

Frequency response with EUREICA; system settles at 60 Hz following demand response and load shedding enabled by the EUREICA framework

§.§.§ Key system metrics, economic, and distributional impacts

In our simulations, we find that attack mitigation comes at the expense of increased operational costs for the PMO since it needs to dispatch more expensive local resources to a greater extent, rather than importing cheaper power from the main grid (at the LMP rate). The PMOs and SMOs also need to adequately compensate agents for the critical flexibility they provide. As shown in <ref> for attack 2b, the attack increases the system operating costs by around $7\%$, and the mitigation steps raise the cost by over $31\%$, both relative to the pre-attack case. However, the PMO could recoup this through other revenue streams and cost savings. For example, the transmission system operators may compensate PMOs for locally containing attacks. Being able to leverage local DER flexibility through markets could also reduce the amount of auxiliary backup generation that the PMO needs to maintain, and lower the reserves it may have to otherwise procure from capacity or ancillary service markets. The PMO in turn could also redistribute some of these benefits among the SMOs and SMAs.

Summary of metrics for large-scale attack 2b scenario.

Pre-attack     Post-attack     Attack mitigation
Power import from main grid [kW]     1,325                   1,821 (+37.4%)          1,328

Total cost [$]                            10,752                  11,500 (+7%)            14,156 (+31.7%)

Total load [kW]                            2,064                   2,023 (-0.02%)          1,775 (-14%)

§.§ Mitigation of attack 3

EUREICA IEEE 123-node feeder for reconfiguration module validation.

We now consider the attack scenario where the distribution grid is isolated from the transmission system. In such a case, the distribution grid is fed through an alternate circuit such as from node 350 (see <ref>). A typical response in such a case is that the distribution grid breaks into several “zones” - creating smaller islands where only a portion of the load is fed through any DERs that may be present. We show below that with the increased awareness provided by the EUREICA framework, a much higher percentage of consumers remain unaffected, by suitably leveraging the DERs at node 48, the microgrid system connected at node 65 (marked by the red circles in <ref>, and DR methodologies. In order to ensure feasibility and supply-demand balance with islanding, we also introduce two large diesel generators located at nodes 48 and 65 which may only be called upon when the feeder is islanded. Three cases are presented.

§.§.§ Critical loads distributed across the feeder

In this case, through the proposed resilience-based IoT load restoration with DR optimization strategy (see SI <ref>, a feasible reconfiguration path is computed to open or close tie switches and completely or partially shed non-critical grid edge loads using reconfiguration to allow the available generation resources to cover approximately $30\%$ of total load in the system. As seen in <ref>, with almost $70\%$ of the load shed (second graph from the top) between 13:00 and 14:00, and batteries only allowed to discharge, if possible, to supply extra energy (third graph from the top), the burden on the diesel generators is significantly alleviated as they only need to ramp up to about $230$ KW. These results were validated using the HELICS co-simulation platform at PNNL (see SI <ref> for details). Additional validation results using HELICS and AIRES are included in SI <ref> and SI <ref>, respectively).

Demand and DG with resilience-based reconfiguration during attack 3.

§.§.§ Critical loads aggregated in a single zone

In this case, the SA from EUREICA helps the reconfiguration algorithm to disconnect or open the switches 18-135 and 151-300 to island zone 3 and pick up only the critical loads in this zone using the DG at node 48, which is a total of 430 kW. The results from this case are shown in <ref>. These results were validated using the DERIM and ADMS-DOTS software at LTDES (see SI <ref> for details). Additional results are included in SI <ref>.

Primary node load change between 12:59 (before) and 13:00 (after attack).

§.§.§ Mitigation with a military microgrid

We assume that there is a military microgrid at node 66 in the primary circuit, which serves as a backup directly in the distribution system. Under current regulations, defense critical systems have to be disconnected and isolated in the event of contingencies. Since EUREICA has the ability to identify trusted resources, our thesis is that there is confidence in the security of this resource as well as in meeting the power flow requirements, making it feasible to use this additional resource for attack 3 mitigation. First, the fault is isolated using reconfiguration based on the algorithm described in SI <ref>. The reconfiguration algorithm returns the most resilient path for implementation, in this case, only one feasible path is present so that is chosen. This islands the feeder by opening the switch between nodes 150 and 149 and connecting the switches to the DG and microgrid at nodes 48 and 66, respectively. Then, a combination of $\approx$ 1.7 MW from the microgrid at node 66, 560 kW from the DG at node 48, and customer-side DR is utilized to pick up approximately 80% of the total load of the feeder. Some results from this case are shown in <ref>, validated using the ARIES platform at NREL (see SI <ref> for details). The complete set of results can be found in SI <ref>.

Microgrid response after reconfiguration.

DG response after reconfiguration.

System response after reconfiguration with microgrid.

§ SUMMARY

We have proposed a framework, EUREICA, for achieving grid resilience through the coordination of IoT-Coordinated Assets that are trustable. A local electricity market that has been previously shown to lead to grid reliability and provide services such as voltage support and overall power balance, is leveraged in this framework to ensure grid resilience. The local market accomplishes this through SA to co-located operators. This SA consists of information about DERs and their power injections, as well as their levels of trustability, commitment, and resilience. With this SA, we have shown that a range of cyberattacks can be mitigated using local trustable resources without stressing the bulk grid. The demonstrations have been carried out using a variety of platforms with high fidelity, hardware-in-the-loop, and utility-friendly validation software.

This work was supported by the US Department of Energy under Award
DOE-OE0000920 and the MIT Energy Initiative. We gratefully acknowledge several useful discussions with Karan Kalsi at PNNL, Rob Hovsapian at NREL, and Jayant Kumar at LTDES. V.N. would like to acknowledge the support of a summer internship at NREL.

[1]
NASEM, Enhancing the Resilience of the Nation's Electricity System.
(The National Academies Press, Washington, DC), (2017).

[2]
CHERNOVITE's PIPEDREAM Malware Targeting Industrial Control Systems (ICS) (2022) available at <https://www.dragos.com/blog/industry-news/chernovite-pipedream-malware-targeting-industrial-control-systems/> (accessed 10 April 2024).

[3]
T Nguyen, et al., Electric power grid resilience to cyber adversaries: State of the art.
IEEE Access 8, 87592–87608 (2020).

[4]
DE Whitehead, K Owens, D Gammel, J Smith, Ukraine cyber-induced power outage: Analysis and practical mitigation strategies in Proceedings of the 2017 70th Annual Conference for Protective Relay Engineers (CPRE).
(IEEE), pp. 1–8 (2017).

[5]
NASEM, The Role of Net Metering in the Evolving Electricity System.
(The National Academies Press, Washington, DC), (2023).

[6]
VJ Nair, V Venkataramanan, R Haider, AM Annaswamy, A hierarchical local electricity market for a der-rich grid edge.
IEEE Transactions on Smart Grid 14, 1353–1366 (2022).

[7]
L Kristov, P De Martini, JD Taft, A tale of two visions: Designing a decentralized transactive electric system.
IEEE Power and Energy Magazine 14, 63–69 (2016).

[8]
T Chen, Q Alsafasfeh, H Pourbabak, W Su, The next-generation us retail electricity market with customers and prosumers—a bibliographical survey.
Energies 11, 8 (2018).

[9]
S Bjarghov, et al., Developments and challenges in local electricity markets: A comprehensive review.
IEEE Access 9, 58910–58943 (2021).

[10]
T Sousa, et al., Peer-to-peer and community-based markets: A comprehensive review.
Renewable and Sustainable Energy Reviews 104, 367–378 (2019).

[11]
T Pinto, Z Vale, S Widergren, eds., Local Electricity Markets.
(Academic Press), pp. i–iii (2021).

[12]
Utility Dive, Coned virtual power plant shows how new york's rev is reforming utility practices; available at <https://www.utilitydive.com/news/coned-virtual-power-plant-shows-how-new-yorks-rev-is-reforming-utility-pra/421053/> (accessed 28 march 2024) (2016).

[13]
A Lüth, J Weibezahn, JM Zepter, On distributional effects in local electricity market designs—evidence from a german case study.
Energies 13, 1993 (2020).

[14]
R Haider, et al., Reinventing the utility for distributed energy resources: A proposal for retail electricity markets.
Advances in Applied Energy 2, 100026 (2021).

[15]
SM Dibaji, et al., A systems and control perspective of CPS security.
Annual Reviews in Control 47, 394–411 (2019).

[16]
CD Brummitt, RM D’Souza, EA Leicht, Suppressing cascades of load in interdependent networks.
Proceedings of the National Academy of Sciences 109, E680–E689 (2012).

[17]
Z Li, M Shahidehpour, F Aminifar, A Alabdulwahab, Y Al-Turki, Networked microgrids for enhancing the power system resilience.
Proceedings of the IEEE 105, 1289–1310 (2017).

[18]
Y Li, X Wei, Y Li, Z Dong, M Shahidehpour, Detection of false data injection attacks in smart grid: A secure federated deep learning approach.
IEEE Transactions on Smart Grid 13, 4862–4872 (2022).

[19]
T Shekari, AA Cardenas, R Beyah, MaDIoT 2.0: Modern High-Wattage IoT botnet attacks and defenses in Proceesdings of the 31st USENIX Security Symposium (USENIX Security 22).
(USENIX Association, Boston, MA), pp. 3539–3556 (2022).

[20]
P Eder-Neuhauser, T Zseby, J Fabini, G Vormayr, Cyber attack models for smart grid environments.
Sustainable Energy, Grids and Networks 12, 10–29 (2017).

[21]
H He, J Yan, Cyber-physical attacks and defences in the smart grid: a survey.
IET Cyber-Physical Systems: Theory & Applications 1, 13–27 (2016).

[22]
VJ Nair, P Srivastava, A Annaswamy, Enhancing power grid resilience to cyber-physical attacks using distributed retail electricity markets in Proceedings of the 2024 IEEE/ACM International Conference on International Conference on Cyber-Physical Systems (ICCPS).
(IEEE), (2024).

[23]
VJ Nair, A Annaswamy, Local retail electricity markets for distribution grid services in Proceedings of the 2023 IEEE Conference on Control Technology and Applications (CCTA).
(IEEE), pp. 32–39 (2023).

[24]
SS Amiri, M Rahmani, JD McDonald, An updated review on distribution management systems within a smart grid structure in 2021 11th Smart Grid Conference (SGC).
(IEEE), pp. 1–5 (2021).

[25]
KA Horowitz, et al., An overview of distributed energy resource (der) interconnection: Current practices and emerging solutions, Technical report (2019).

[26]
P Srivastava, et al., Voltage regulation in distribution grids: A survey.
Annual Reviews in Control 55, 165–181 (2023).

[27]
K Zetter, Inside the cunning, unprecedented hack of ukraine’s power grid (2018) available at <https://www.wired.com/2016/03/inside-cunningunprecedented-hack-ukraines-power-grid/> (accessed 11 April 2024).

[28]
J Kennedy, Dragonfly: Western energy sector targeted by sophisticated attack group (2017) available at <https://www.symantec.com/blogs/threatintelligence/dragonfly-energy-sector-cyber-attacks> (accessed 10 April 2024).

[29]
RM Lee, MJ Assante, T Conway, ICS defense use case: Analysis of the cyber attack on the Ukrainian power grid, (Electricity Information Sharing and Analysis Center, SANS ICS), Technical report (2016).

[30]
M Zeller, Myth or reality – does the aurora vulnerability pose a risk to my generator? in 64th Ann. Conf. for Protective Relay Engineers.
pp. 130–136 (2011).

[31]
S Soltan, P Mittal, HV Poor, BlackIoT: IoT botnet of high wattage devices can disrupt the power grid in Proceedings of the 27th USENIX Security Symposium (USENIX Security 18).
(USENIX Association, Baltimore, MD), pp. 15–32 (2018).

[32]
B Huang, AA Cardenas, R Baldick, Not everything is dark and gloomy: Power grid protections against IoT demand attacks in 28th USENIX Security Symposium (USENIX Security 19).
(USENIX Association, Santa Clara, CA), pp. 1115–1132 (2019).

[33]
WW Hogan, Independent system operator: Pricing and flexibility in a competitive electricity market.
Center for Business and Government, JF Kennedy School of Government, Harvard University, MA (1998).

[34]
V Venkataramanan, S Kaza, AM Annaswamy, Der forecast using privacy-preserving federated learning.
IEEE Internet of Things Journal 10, 2046–2055 (2022).

[35]
PS Sarker, SK Sadanandan, AK Srivastava, Resiliency Metrics for Monitoring and Analysis of Cyber-Power Distribution System With IoTs.
IEEE Internet of Things Journal 10, 7469–7479 (2023).

[36]
CM Schneider, AA Moreira, JS Andrade Jr, S Havlin, HJ Herrmann, Mitigation of malicious attacks on networks.
Proceedings of the National Academy of Sciences 108, 3838–3841 (2011).

[37]
DK Molzahn, IA Hiskens, A Survey of Relaxations and Approximations of the Power Flow Equations.
Foundations and Trends R in Electric Energy Systems 4, 1–221 (2019).

[38]
G Ferro, M Robba, D D'Achiardi, R Haider, AM Annaswamy, A distributed approach to the Optimal Power Flow problem for unbalanced and mesh networks.
IFAC-PapersOnLine 53, 13287–13292 (2020).

[39]
JJ Romvary, G Ferro, R Haider, AM Annaswamy, A Proximal Atomic Coordination Algorithm for Distributed Optimization.
IEEE Transactions on Automatic Control 67, 646–661 (2022).

[40]
G Ferro, M Robba, R Haider, AM Annaswamy, A Distributed-Optimization-Based Architecture for Management of Interconnected Energy Hubs.
IEEE Transactions on Control of Network Systems 9, 1704–1716 (2022).

[41]
R Meyur, et al., Ensembles of realistic power distribution networks.
Proceedings of the National Academy of Sciences 119, e2205772119 (2022).

[42]
First cyberattack on solar, wind assets revealed widespread grid weaknesses (2019) available at <https://www.utilitydive.com/news/first-cyber-attack-on-solar-wind-assets-revealed-widespread-grid-weaknesse/> (accessed 14 April 2024).

[43]
P Kundur, Power system stability.
(CRC Press New York) Vol. 10, (2007).

ACOPFAlternating Current OPF
ADMSAdvanced Distributed Management System
ADPAnomalous Data Point
ARAnomaly Ratio
ARIESAdvanced Research on Integrated Energy Systems
CARCumulative Anomaly Ratio
CICurrent Injection
CSCommitment Score
[CVSS]Common Vulnerability Scoring SystemCVSS
CMOConsumer Market Operator
DCVSDevice and Communication Vulnerabilities present at the Secondary
DERDistributed Energy Resources
DERIMDER Integration Middleware
DGDistributed Generator
DOTSDistribution Operations Training Simulator
DoSDenial of Service
DPVsDistributed Photovoltaics
DRTSDigital Real-Time Simulation
DSODistribution System Operator
DSRDistribution System Resiliency
EUREICAEfficient, Ultra-REsilient, IoT-Coordinated Assets
EVElectric Vehicle
FLFederated Learning
AHPAnalytic Hierarchy Process
GEGeneral Electric
HELICSHierarchical Engine for Large-scale Infrastructure Co-simulation
HVACHeating, Ventilation, and Air Conditioning
ICAsIoT-Coordinated Assets
ICSIndustrial Control Systems
IEEEInstitute of Electrical and Electronics Engineers
IoTInternet of Things
LEMLocal Electricity Market
LMPLocational Marginal Price
MadIoTManipulation of Demand via IoT
MCDMMultiple-Criteria Decision-Making
NRELNational Renewable Energy Lab
NVDNational Vulnerability Database
OPFOptimal Power Flow
OTOperational Technology
PACProximal Atomic Coordination
PCCPoint of Common Coupling
PMPrimary Market
PMAPrimary Market Agent
PMOPrimary Market Operator
PNNLPacific Northwest National Lab
PNRPrimary Node Resiliency
PRMPrimary Resilience Manager
RSResilience Score
RTDSReal-Time Digital Simulator
SASituational Awareness
SMSecondary Market
SMASecondary Market Agent
SMOSecondary Market Operator
SOCPSecond-Order Conic Program
SRMSecondary Resilience Manager
STNRSecondary Transformer Node Resiliency
TSTrustability Score
WHWater Heater

§ IEEE 123-NODE DISTRIBUTION NETWORK

<ref> shows the IEEE 123-node test feeder kersting2001radial, which is the distribution network model used to simulate and validate all our attack scenarios.

IEEE 123-node test feeder network.

§ EXAMPLE INSTANCE OF LEM

Here, we sketch out a possible (hypothetical) instantiation of our LEM for the city of Boston, MA, which is located in the New England (NE) region in the US. <ref> shows the IEEE 39-bus transmission system athay1979practical which is a synthetic representation of the entire NE region, with a peak load of 6254 MW and a total of around 7 million homes. This corresponds to $\approx$ 162 MW and 180,000 homes per bus. Given that the IEEE 123-node distribution feeder has a peak load of roughly 3.6 MW, we can estimate that there will be 44 such primary feeders per transmission bus and 4100 homes per feeder. Thus, the city of Boston with a total of 300,000 homes bostonhouses, will be served by 73 primary feeders across 2 transmission buses.

IEEE 39-bus transmission system.

<ref> shows a breakdown of different entities to form a hierarchical LEM for Boston. Note that the main market operators and agents that are relevant for this work are marked in green.

Example of hypothetical LEM for the city of Boston, MA.

§ DETAILS OF THE LOCAL ELECTRICITY MARKET

<ref> shows the inputs and outputs for different levels of the hierarchical LEM. For both the SM and the PM, the inputs consist of the baseline power injections and flexibility bids, while the outputs are the market schedules (setpoints for power injections) and their associated flexibility ranges, along with the corresponding electricity prices of tariffs.

Overall inputs and outputs in the LEM.

§.§ Objective functions for optimization in the Secondary Market

The four objective functions considered in the SM clearing are defined as:

O1. Maximization of aggregate resilience, $f_i^1$, given by the following, where $RS_j^i$ denotes the resilience score of SMA $j$ under SMO $i$

$$f_i^1 = - \sum_{j=1}^n RS^i_j ((P_j^i - P_j^{i0})^2 +  (Q_j^i - Q_j^{i0})^2)$$

O2. Minimization of net cost, $f_i^2$ to the SMO for running the SM

$$ f_i^2  = \sum_{j=1}^n \mu_j^{iP} P_j^i + \mu_j^{iQ} Q_j^i$$

O3. Maximization of total flexibility, $f_i^3$ that the SMO can extract from all its SMAs

$$f_i^3 = - \sum_{j=1}^n (\delta P_j^i + \delta Q_j^i)$$

O4. Minimization of disutility of the SMAs, $f_i^4$ arising from flexibility provision
$$f_i^4 = \sum_{j=1}^n \beta_j^{iP}(P_j^i - P_j^{i0})^2 + \beta_j^{iQ}(Q_j^i - Q_j^{i0})^2.$$

§.§ Three-phase SM optimization problem

\begin{align}
& \min \sum_{j \in \mathcal{N}_{J,i}} \{f_{j,1}^i,f_{j,2}^i,f_{j,3}^i,f_{j,4}^i\} \label{eq:cost} \\
& f_{1,j}^i \succ f_{2,j}^i \succ f_{3,j}^i \succ f_{4,j}^i \label{eq:ranking}, \; \; \Phi = \{a, b, c\} \\
& f_{j,1}^i = -C^i_j \left(\sum_{\phi \in \Phi}(P_j^{i,\phi} - P_j^{i0,\phi})^2 + (Q_j^{i,\phi} - Q_j^{i0,\phi})^2 \right)\nonumber \\
& f_{j,2}^i  = \sum_{\phi \in \Phi} \mu_j^{iP,\phi} P_j^{i,\phi} + \mu_j^{iQ,\phi} Q_j^{i,\phi} \nonumber \\
& f_{j,3}^i = -\sum_{\phi \in \Phi}(\delta P_j^{i,\phi} + \delta Q_j^{i,\phi}) \nonumber\\
& f_{j,4}^i = \beta_j^{iP}\sum_{\phi \in \Phi}\left(P_j^{i} - P_j^{i0}\right)^2 + \beta_j^{iQ}\sum_{\phi \in \Phi}\left(Q_j^{i} - Q_j^{i0}\right)^2 \nonumber \\
& \text{subject to:} \nonumber\\
& P_j^{i,\phi} - \delta P_j^{i,\phi} \geq \underline{P}_j^{i,\phi} \; Q_j^{i,\phi} - \delta Q_j^{i,\phi} \geq \underline{Q}_j^{i,\phi}  \label{eq:PQ_lowerlim} \\
& P_j^{i,\phi} + \delta P_j^{i,\phi} \leq \overline{P}_j^{i,\phi}, \; Q_j^{i,\phi} + \delta Q_j^{i,\phi} \leq \overline{Q}_j^{i,\phi}  \label{eq:PQ_upperlim} \\
& \delta P_j^{i,\phi}, \; \delta Q_j^{i,\phi} \geq 0, \; 0 \leq \mu_j^{iP} \leq \overline{\mu}^{iP}, 0 \leq \mu_j^{iP} \leq \overline{\mu}^{iQ}  \; \label{eq:flex_fairness} \\
& \sum_{t_s}^{t_s + \Delta t_p} \sum_{j \in \mathcal{N}_{J,i}} \sum_{\phi \in \Phi} \left(\mu_j^{iP,\phi}(t) P_j^{i,\phi}(t) + \mu_j^{iQ,\phi}(t) Q_j^{i,\phi}(t)\right) \Delta t_s  \nonumber \\
& \leq \sum_{\phi \in \Phi} \left(\mu_i^{P^*,\phi}(\hat{t}_p) P_i^{\phi^*}(\hat{t}_p) + \mu_i^{Q^*,\phi}(\hat{t}_p) Q_i^{\phi^*}(\hat{t}_p)\right) \Delta t_p \label{eq:budgetPQ} \\
& \sum_{j \in \mathcal{N}_{J,i}} P_j^{i,\phi}(t_s) = P_i^{\phi^*}(\hat{t}_p), \quad \sum_{j \in \mathcal{N}_{J,i}} Q_j^{i,\phi}(t_s) = Q_i^{\phi^*}(\hat{t}_p) \label{eq:PQbalance1}
\end{align}

§.§ Bidding in the Primary Market

The bid submitted by each SMO $j \in \mathcal{N}$ into the PM $\mathcal{B}_i^{P}$ defined as:

\begin{align}
\mathcal{B}_i^{S} =\{P_i^{0},Q_i^{0},\underline{P}_i,\underline{Q}_i,\overline{P}_i,\overline{Q}_i,\alpha_i^P, \alpha_i^Q, \beta_i^P,\beta_i^Q\}. \label{eq:pm_bid1}
\end{align}

\begin{align}
& P_i^0(t_p) = \sum_{j \in \mathcal{N}_i}  P_j^{i*}(t_p),\; Q_i^0(t_p) = \sum_{j \in \mathcal{N}_i}  Q_j^{i*}(t_p) \label{eq:aggregation}  \\
& \underline{P}_i = \sum_{j\in \mathcal{N}_i} P^i_j - \delta P_j^{i*}, \overline{P}_i =\sum_{j\in \mathcal{N}_i} P_j^{i*} + \delta P_j^{i*} \nonumber \\
& \underline{Q}_i = \sum_{j\in \mathcal{N}_i} Q_j^{i*} - \delta Q_j^{i*}, \overline{Q}_i = \sum_{j\in \mathcal{N}_i} Q_j^{i*} + \delta Q_j^{i*}\nonumber
\end{align}

In the above, (i) $P_i^0, Q_i^0$ denote the baseline active and reactive power injection bids of the SMOs, (ii) ($\underline{P}_i, \underline{Q}_i$) and ($\overline{P}_i, \overline{Q}_i$) denote the downward and upward flexibilities (around their nominal values) in active and reactive power, respectively, (iii) $\alpha_i^P, \alpha_i^Q$ are the local net generation costs, and (iv) $\beta_i^P, \beta_i^Q$ are the flexibility disutility parameters of SMO $i$ for P and Q, respectively. The SMO computes (iii) and (iv) as weighted averages of all their SMA retail tariffs, and SMA disutility parameters, respectively, as follows:

\alpha_i^P = \frac{\sum_{j \in \mathcal{N}_i \mu_j^{iP*} |P_j^{i*}|}}{\sum_{j \in \mathcal{N}_i |P_j^{i*}|}}, \beta_i^P = \frac{\sum_{j \in \mathcal{N}_i \beta_j^{iP*} |P_j^{i*}|}}{\sum_{j \in \mathcal{N}_i |P_j^{i*}|}}

§.§ Objective functions for optimization in the Primary Market

We define the following functions

\begin{align}
f^P(\textbf{y}^P) &=\sum_{i \in \mathcal{N}} f^P_i(\textbf{y}^P_i) = \sum_{i \in \mathcal{N}} \Big[f_{i}^{\text{Load-Disutil}}(\textbf{y}^P_i) \nonumber \\
& + f_{i}^{\text{Gen-Cost}}(\textbf{y}^P_i)\Big] + \xi\Big[\sum_{(ki) \in \mathcal{E}} f_{ki}^{\text{Loss}}(\textbf{y}^P_i)\Big] \label{eq:cost:1} \\
f_{i}^{\text{Load-Disutil}}(\textbf{y}^P_i)&= \beta_{i}^{P}(P_{i}^{L}-P_{i}^{L0})^2+\beta_{i}^{Q}(Q_{i}^L-Q_{i}^{L0})^2 \label{eq:cost:2} \\
f_{i}^{\text{Gen-Cost}}(\textbf{y}^P_i)&= \begin{cases}\alpha_{i}^{P} (P_{i}^{G})^2 +\alpha_{i}^{Q} (Q_{i}^{G})^2, & \\
\lambda_{i}^{P} P_{i}^{G}+\lambda_{i}^{Q} Q_{i}^{G},\text{if } i & \text{ is PCC}
\end{cases} \label{eq:cost:3}\\
f_{ki}^{\text{Loss}}(\textbf{y}^P_i)&=R_{ki}|I_{ki}|^2 \label{eq:cost:4}
\end{align}

The objective function used in <ref> used is a weighted linear combination of (i) maximizing social welfare in <ref>, (ii) minimizing total generation costs in <ref> and (iii) minimizing electrical line losses in <ref>. The total cost includes paying the locational marginal price (LMP) $\lambda$ for importing power from the transmission grid at the point of common coupling (PCC), as well as the payments to local generator PMAs that provide net positive injections into the PM. We divide by suitable base values to convert all quantities to per unit (between 0 and 1 p.u.). Thus, it is reasonable to combine all the terms into a single objective function using a simple weighted sum. The hyperparameter $\xi$ controls the tradeoff between penalizing line losses versus optimizing for other objectives. The coefficients $\alpha_i, \beta_i$, are communicated by each PMA $i$ as part of their bids, while $\xi$ is a global hyperparameter common to all PMAs, and determined by the PMO. Here, $R$ denotes the network resistance matrix and $\mathcal{E}$ denotes the set of all edges in the network.

§.§ Computation of commitment scores

We describe here the details of computing the commitment reliability score, mentioned in <ref>. From the SM clearing, the SMAs $j$ are directed by their SMO $i$ to keep their net injections within the intervals $[P_j^{i*} - \delta P^{i*}_j, P_j^{i*} + \delta P_j^{i*}]$. We first compute the deviations (if any) in their actual responses $\hat{P}^i_j$ from this range, where $\llbracket \cdot \rrbracket$ denotes the indicator function:

\begin{align}
e^{iP}_j(t_s) & = \llbracket\hat{P}^i_j > \overline{P}_j^{i*}\rrbracket (\hat{P}_j^i - \overline{P}_j^{i*}) + \llbracket\hat{P}^i_j < \underline{P}_j^{i*}\rrbracket(\underline{P}_j^{i*} - \hat{P}^i_j) \nonumber \\ & + \llbracket\underline{P}_j^{i*} \leq \hat{P}^i_j \leq \overline{P}_j^{i*}\rrbracket\max (\hat{P}^i_j - \overline{P}_j^{i*},\underline{P}_j^{i*} - \hat{P}^i_j) \label{eq:commit1}
\end{align}

We then obtain relative deviations by comparing these with the magnitudes of their corresponding baseline setpoints:

\begin{equation}
\overline{e^{iP}_j} (t_s) = \frac{e^{iP}_j(t_s)}{|P_j^{i*}(t_s)|}, \quad \overline{e^{iQ}_j} (t_s) = \frac{e^{iQ}_j(t_s)}{|Q_j^{i*}(t_s)|} \label{eq:commit2}
\end{equation}

These are then normalized to unit vectors to compare the deviations among all SMAs overseen by the SMO. This allows the SMO to assess their relative performance across all its SMAs.

\begin{equation}
\widetilde{\mathbf{e^{iP}}}(t_s) = \frac{\mathbf{\overline{e^{iP}}}(t_s)}{{\|\mathbf{\overline{e^{iP}}}}(t_s)\|}, \; \widetilde{\mathbf{e^{iQ}}}(t_s) = \frac{\mathbf{\overline{e^{iQ}}}(t_s)}{\|\mathbf{\overline{e^{iQ}}}(t_s)\|} \label{eq:commit3}
\end{equation}

The scores are then updated, with the score being increased when the SMAs follow their contracts and decreased otherwise:

\begin{align}
C^i_j(t_s) & = \begin{cases} 1 &\mbox{if } t_s = 0 \\
C^i_j(t_s - 1) - \frac{\widetilde{e^{iP}_j} (t_s) + \widetilde{e^{iQ}_j} (t_s)}{2} & \mbox{if } t_s > 0 \end{cases} \label{eq:commit4}
\end{align}

Finally, we perform min-max normalization across all the SMAs' scores to ensure that $0 \leq C_j \leq 1 \; \forall \; \text{SMAs } j$.

\overline{C}^i_j = \frac{C_j^i - \max_j C^i_j}{\max_j C^i_j - \min_j C^i_j}

§ TRUSTABILITY SCORES AND RESILIENCY METRICS

§.§ Computation of IoT trustability scores

The IoT trustability score (TS) is computed utilizing the federated self-learning concept Sarker2023ResiliencyIoTs. Anomalies in IoT data are the key factor in the formulation of the IoT TS. Another contributing factor is the IoT device's market commitment history. More details of the features are shown in <ref>.

Features considered for each type of data.

Data Source     Features
c]@c@IoTs network
packet      c]@c@Source/Destination IP,  Source/Destination port,
Packet length, Protocols, Intra-packet arrival time
HVAC     c]@c@Timestamp, Load, Indoor temperature,
outdoor temperature, Temperature setpoint,
Indoor area, Building thermal insulation
PV     c]@c@Timestamp, Power generation,
Rating, Solar irradiance
Battery     c]@c@Timestamp, Charging/Discharging rate,
SoC, KW capacity
EV     c]@c@Timestamp, Charging rate,

To detect anomalies, we learn the IoT data's expected behavior and prediction for the short time steps. Predicted data is compared with measured data for anomaly detection. For prediction, we are using an autoencoder neural network for federated unsupervised learning. There will be one autoencoder model for each IoT device to train on its physical data and one more autoencoder model to train only on IoT network packet data. For each type of data, there is a tolerance value $T_{err}$ for the relative error (RE). If any data point ($DP$) crosses $T_{err}$, then that is flagged as an anomalous data point ($ADP$). So, for any reporting time period $\Delta t$, non-anomaly ratio ($NAR$) is calculated using,

\begin{equation}
NAR=1- \frac{\text{Total ADP number over }{\Delta t}^{}}{\text{Total DP number over }{\Delta t}}
\end{equation}

Next, the cumulative non-anomaly ratio ($CNAR$) is computed, where $T$ is the fixed total time period and is always divisible by $\Delta t$.

\begin{equation}\label{cnar}
CNAR_t=\sum_{j=1}^{\frac{T}{\Delta t}}\frac{T}{j\Delta t} NAR_{t-j\Delta t}
\end{equation}

IoT Trustability Score ($TS$) for time $t$ and building/house $i$ is calculated:

\begin{equation}\label{its}
TS_{t,i}=w_t\times NAR_t+ w_{t-}\times\frac{CNAR_t}{CNAR_{max}}
\end{equation}

\begin{equation} \label{wtits}
(I) \; w_t \geq w_{t-} \quad \quad
(II) \; w_t + w_{t-}=1
\end{equation}

Here, $CNAR_{max}$ is calculated using (<ref>) with the maximum $NAR$ being $NAR=1$ for the whole time period $T$. Finally, to get the overall $TS_t$ of any observation node with IoTs at time $t$, we average the $TS_{t,i}$ of all the clients $i$ of that observation node to calculate $TS_t$:

\begin{equation}
TS= \frac{\sum_{i=1}^{M} TS_{t,i}}{M}
\end{equation}

where $M$ is the total number of clients or buildings/houses at that observation node.

§.§ Secondary Transformer and Primary Node Resiliency Metric (STNR and PNR)

Secondary transformer node resiliency (STNR) is computed using multiple resiliency factors and TS.

\begin{equation}
STNR_j = \prod_{i=1}^{n_c} F_i^{W_i}
\end{equation}

where $n_c$ is the total number of factors for the category of the secondary level node, $F_i$ is the value for each factor, and $W_i$ is the normalized weight for each factor. These factors $f$ influencing resiliency are determined and assigned weights to aggregate into the PNR score. Factors that can be determined directly from the secondary level configuration are described in <ref>. All the device and communication vulnerabilities present at the secondary (DCVS) level of a primary node are identified using the national vulnerability database (NVD) booth2013national. Then DCVS factor is calculated as,

\begin{equation}
DCVS=\frac{1}{\sum_{i=1}^{N_s} CVSS_i}
\end{equation}

where $N_s$ is the number of total vulnerabilities present at the secondary level. Here, the common vulnerability scoring system (CVSS) is one of several methods to measure the impact of vulnerabilities in devices known as Common Vulnerabilities and Exposures (CVE). It is an open set of standards used to assess the vulnerability of software and assign severity along a scale of 0-10. The National Institute of Standards and Technology (NIST) analyzes all identified vulnerabilities and enlists these in the NVD. In case of the absence of any vulnerability, DCVS will be equal to $1$.

Weight assignment and aggregation are managed by fuzzy multiple-criteria decision-making (MCDM), specifically the fuzzy analytic hierarchy process (Fuzzy AHP). A weighted average of the STNR results in the primary node resiliency (PNR):

\begin{equation}\label{eq:pnr}
PNR_k =  \frac{\sum_{j=1}^n (STNR_j \times W_j )}{\sum_{j=1}^n W_j }
\end{equation}

where $W_i$ is the weighted coefficient for the $i^{th}$ secondary feeder node.

§.§ Distribution System Resiliency (DSR)

Let $F=(f_{ij})\in R_+^{m \times n}$ be the factors value matrix, where $f_{ij}$ is value of factor $i$ of primary node $j$. The higher the value of $f_{ij}$, the more the node will contribute to the resiliency metric in regard to that factor. Now, following the data envelopment analysis (DEA) method, each node $p$ can choose a set of weights $w^p=(w^p_1,...w^p_m)$, where, $ \sum_{i=1}^{m}w^p_i=1$. Now the relative contribution (RC) of the node $p$ to the total contribution of all the nodes towards DSR, as measured by node $p$’s weight selection can be evaluated as,

\begin{equation}\label{eq:ratio}
RC^p=\frac{\sum_{i=1}^{m}w^p_i f_{ip}}{\sum_{i=1}^{m}w^p_i \sum_{j=1}^{n}\left (f_{ij}\right)}
\end{equation}

Now, each node wants to maximize this ratio in <ref> to have the best set of weights so that they can contribute to the maximum possible value in DSR. Using the weight vector for each node, a combination of multiplicative and additive methods are used to get the DSR.

\begin{equation} \label{eq:dsr}
DSR= \sum_{j=1}^{n} \left (\prod_{i=1}^{m} \left (f_{ij}  \right )^{w^j_i} \right )
\end{equation}

Overview of the developed resilience score for the distribution system with IoTs

Details related to the computation of DSR are shown in <ref>.

§ POWER SYSTEM MODELS

§.§ Branch flow model

The branch flow OPF problem is formally stated as follows, where $R$ and $X$ denote the network resistance and reactance matrices respectively, $v$ and $I$ denote the nodal voltage magnitudes and branch currents respectively, and $\mathcal{E}$ denotes the set of all edges in the network. The primal decision variables here for each PMA $i$ are $\textbf{y}^P_i = [P_i^G, Q_i^G, P_i^L, Q_i^L, v_i, I_{ki}] \; \forall i \in \mathbf{N}, \; (ik) \in \mathcal{E}$.

\begin{align}
& \underset{\textbf{y}^P}{\min} \;\; f^{S-W}(\textbf{y}^P) \label{eq:pmo} \\
& \text{subject to:} \nonumber\\
& v_{i} - v_{k} = \left(R_{ki}^{2} + X_{ki}^{2}\right) |I_{ki}|^2 - 2 \left(R_{ki} P_{ki} + X_{ki} Q_{ki}\right)\nonumber\\
& P_{i}^{G} - P_{i}^{L} = -P_{ki} + R_{ki} |I_{ki}|^2 + \sum_{k:\left(i,k\right) \in \mathcal{E}} P_{ik}\nonumber\\
& Q_{i}^{G} - Q_{i}^{L} = -Q_{ki} + X_{ki} |I_{ki}|^2 + \sum_{k:\left(ik\right) \in \mathcal{E}} Q_{ik}\nonumber\\
& P_{ki}^{2} + Q_{ki}^{2} \leq \overline{S}_{ki}^{2}, \; P_{ki}^{2} + Q_{ki}^{2} \leq v_{i} |I_{ki}|^2, \;\underline{v}_i \leq v_i \leq \overline{v}_i  \nonumber\\
& \underline{P}_i^{G} \leq P_i^{G} \leq \overline{P}_i^{G}, \; \underline{P}_i^{L} \leq P_i^{L} \leq \overline{P}_i^{L} \nonumber\\
& \underline{Q}_i^{G} \leq Q_i^{G} \leq \overline{Q}_i^{G}, \; \underline{Q}_i^{L} \leq Q_i^{L} \leq \overline{Q}_i^{L} \nonumber\\
\end{align}

§.§ Current injection model

The primal decision variables for each SMO $i$ obtained by solving the optimization problem $\textbf{y}^P_i =[P_i^{\phi},Q_i^{\phi},V_{i}^{\phi,R},V_{i}^{\phi,I},I_{i}^{\phi,R},I_{i}^{\phi,I},I_{ik}^{\phi,R},I_{ik}^{\phi,I}]$ consists of (i) active ($P_i^{\phi*}$) and reactive ($Q_i^{\phi*}$) power setpoints (ii) real and imaginary components of nodal voltages ($V_{i}^{\phi,R*}, V_{i}^{\phi,I*}$) and current injections ($I_{i}^{\phi,R*}, I_{i}^{\phi,I*}$). Note that these are solved for each non-zero phase $\phi \in \mathcal{P} = \{a,b,c\}$. The CI-OPF problem formulation is given by:

\begin{align}
&\min _x f^{obj}(x)\\
&I^R=\operatorname{Re}(\mathrm{Y} V), \; I^I=\operatorname{Im}(\mathrm{YV}) \label{eq:ohmslaw} \\
&P_i^\phi=V_i^{\phi, R} I_i^{\phi, R}+V_i^{\phi, I} I_i^{\phi, I} \quad \forall i \in \mathcal{N}, \phi \in \mathcal{P} \label{eq:bilinearP} \\
&Q_i^\phi=-V_i^{\phi, R} I_i^{\phi, I}+V_i^{\phi, I} I_i^{\phi, R} \quad \forall i \in \mathcal{N}, \phi \in \mathcal{P} \label{eq:bilinearQ}\\
& (I_{ik}^{\phi,R})^2+(I_{ik}^{\phi,I})^2\le \overline{I_{ik}^{\phi}}^2 \quad \forall i\in \mathcal{N}, \phi\in \mathcal{P}, (ik) \in \mathcal{E} \label{eq:Ibranchlims} \\
& \underline{{V_{i}^\phi}}^2\le({V_{i}^{\phi,R}})^2 + ({V_{i}^{\phi,I}})^2 \le \overline{V_{i}^{\phi}}^2 \quad \forall i\in \mathcal{N}, \phi\in \mathcal{P} \label{eq:Vmaglims} \\
& \underline{P_i^\phi} \leq P_i^\phi \leq \overline{P_i^\phi}, \;  \underline{Q_i^\phi} \leq Q_i^\phi \leq \overline{Q_i^\phi}
\end{align}

where $Y$ is the 3-phase bus admittance matrix for the network, and $V$ and $I$ are matrices of nodal voltages and currents respectively. Problem <ref> is nonconvex due to bilinear constraints <ref>, and the ring constraint <ref> on voltage magnitudes. We obtain a convex relaxation by using McCormick envelopes (MCE), which represent the convex hull of a bilinear product $w=xy$ by using upper and lower limits on $x, \; y$. Thus, we replace the bilinear equality with a series of linear inequalities, denoted as $\text{MCE}(w) = \{w=xy: x\in [ \underline{x}, \overline{x}],  y\in [ \underline{y}, \overline{y}]\}$:

\begin{equation} \label{eq:mce}
MCE(w,\underline{x},\overline{x},\underline{y},\overline{y}) =
\begin{cases}
w\ge \underline{x}y+x\underline{y}-\underline{x}\underline{y} \\
w\ge \overline{x}y+x\overline{y}-\overline{x}\overline{y} \\
w\le \underline{x}y+x\overline{y}-\overline{x}\underline{y} \\
w\le \overline{x}y+x\underline{y}-\underline{x}\overline{y}
\end{cases}
\end{equation}

We introduce auxiliary variables for each of the four bilinear terms $\{a_i^{\phi},b_i^{\phi},c_i^{\phi},d_i^{\phi}\} = \{V_i^{\phi, R} I_i^{\phi, R},V_i^{\phi, I} I_i^{\phi, I},V_i^{\phi, R} I_i^{\phi, I},V_i^{\phi, I} I_i^{\phi, R}\}$ allowing us to convert constraints <ref> to linear constraints with MCE constraints on each of the auxiliary variables. We also need additional constraints on the nodal current injections and nodal voltages in order to define the MCE constraints. These voltage and current bounds can be determined by applying a suitable preprocessing method using the nodal $P$ and $Q$ limits from the SMO bids Ferro2020ANetworks. The resulting bounds will also implicitly satisfy constraints <ref> and <ref>. Thus, we can replace constraints <ref> with the following set of constraints in order to obtain the relaxed CI-OPF problem, which reduces to a linear program that can be solved easily. However, we do incur the overhead of computing the tightest possible $V$ and $I$ bounds to obtain a good convex relaxation, which in turn ensures that the relaxed solutions are feasible for the original problem.

\begin{align}
& P_i^\phi=a_i^{\phi} + b_i^{\phi}, \; \; Q_i^\phi=-c_i^{\phi} + d_i^{\phi} \quad \forall i \in \mathcal{N}, \phi \in \mathcal{P} \label{eq:linearPQ}\\
& \underline{I_{i}^{\phi,R}} \le I_{i}^{\phi,R} \le \overline{I_{i}^{\phi,R}}, \; \underline{I_{i}^{\phi,I}} \le I_{i}^{\phi,I} \le \overline{_{i}^{\phi,I}}  \label{eq:Ilims} \\
& \underline{V_{i}^{\phi,R}} \le V_{i}^{\phi,R}\le \overline{V_{i}^{\phi,R}}, \; \underline{V_{i}^{\phi,I}} \le V_{i}^{\phi,I}\le \overline{V_{i}^{\phi,I}} \label{eq:Vlims} \\
& a_i^{\phi} \in MCE(V_{i}^{\phi,R}I_{i}^{\phi,R},\underline{V_{i}^{\phi,R}},\overline{V_{i}^{\phi,R}},\underline{I_{i}^{\phi,R}},\overline{I_{i}^{\phi,R}}) \\
& b_{i}^{\phi} \in MCE(V_{i}^{\phi,I}I_{i}^{\phi,I},\underline{V_{i}^{\phi,I}},\overline{V_{i}^{\phi,I}},\underline{I_{i}^{\phi,I}},\overline{I_{i}^{\phi,I}}) \\
& c_{i}^{\phi} \in MCE(V_{i}^{\phi,R}I_{i}^{\phi,I},\underline{V_{i}^{\phi,R}},\overline{V_{i}^{\phi,R}},\underline{I_{i}^{\phi,I}},\overline{I_{i}^{\phi,I}}) \\
& d_{i}^{\phi} \in MCE(V_{i}^{\phi,I}I_{i}^{\phi,R},\underline{V_{i}^{\phi,I}},\overline{V_{i}^{\phi,I}},\underline{I_{i}^{\phi,R}},\overline{I_{i}^{\phi,R}}) \label{eq:mce_ineqs}
\end{align}

§ DISTRIBUTED OPTIMIZATION

For a given global optimization (primal) problem with equality and inequality constraints for $K$ number of nodes (or agents):

\begin{equation}
\min _x \sum_{i=1}^K f_i(x) \text { s.t. } G x=b, \quad H x \leq d
\end{equation}

We can decompose this into $\mathcal{S} = \{S_1, S_2,\dots S_K\}$ coupled optimization problems, known as atoms (representing each SMO $i$). We separate the vector of all decision variables $x$ into two sets: $\mathcal{L} = \{L_i, \forall i \in [K]\}$ and $\mathcal{O} = \{O_i, \forall i \in [K]\}$ which is a partition of decision variables into those that are owned and copied by atom $i$, respectively. We can similarly also decompose the constraints into sets owned by each atom $\mathcal{C} = \{C_i, \forall i \in [K]\}$. These variable copies across multiple atoms can then be used to satisfy coupled constraints and global objectives. Note that for a number $K$, $[K] =\{1, 2, \dots K\}.$

The decomposed (or atomized) optimization problem is shown in <ref>, where $a_j$ and $f_j(a_j)$ are the primal decision variables (both owned and copies) and individual objective functions corresponding to each SMO atom, respectively. $G_j$ and $H_j$ are the atomic constraint submatrices of $G$ and $H$, while $b_j$ and $d_j$ are subvectors of $b$ and $d$ of the right hand side constraint vectors $b$ and $d$, respectively. $B$ is the directed graph incidence matrix defining the owned and copied atomic variables. This incidence matrix allows us to fully parallelize the distributed optimization by defining coordination or consensus constraints, which enforce that all the copied variables for each atom $j$ must equal the values of their corresponding owned values in every other atom $i \neq j$. $B_j$ and $B^j$ denote the incoming and outgoing edges for atom $j$ respectively. Here

\begin{align}
& \min _{a_j} \sum_{j \in K} f_j\left(a_j\right) \label{eq:opt_atom} \\
& \text { s.t. } G_j a_j=b_j, \; H_j a_j \leq d_j, \; B_j a=0 \; \forall j \in [K]  \nonumber\\
& B_{im} \triangleq \begin{cases}
-1, & \text { if } i \text { is 'owned" and } m \text { a related "copy" } \\
1, & \text { if } m \text { is "owned" and } i \text { a related "copy" } \\
0, & \text { otherwise } \end{cases} \nonumber
\end{align}

The augmented Lagrangian is first atomized or decomposed for each node or SMO, introducing dual variables $\eta$ and $\nu$ corresponding to primal equality and coordination constraints respectively. Note that the inequality constraints are handled directly during the primal minimization step by appropriately defining the feasible set.

\begin{align}
\mathcal{L}(a, \eta, \nu) & =\sum_{j \in K}\left[f_j\left(a_j\right)+\eta_j^T\left(G_j a_j-b_j\right)+\nu_j^T B_j a\right] \nonumber\\
& =\sum_{j \in K}\left[f_j\left(a_j\right)+\eta_j^T\left(G_j a_j-b_j\right)+\nu^T B^j a_j\right] \nonumber\\
& \triangleq \sum_{j \in K} \mathcal{L}_j\left(a_j, \eta_j, \nu\right) \label{eq:aug_lang}
\end{align}

§.§ PAC algorithm

We can then apply the prox-linear approach of aljanaby2005survey to <ref> and obtain the proximal atomic coordination (PAC) algorithm Romvary2022AOptimization,Haider2020TowardGrids:

\begin{aligned}
& a_j[\tau+1]=\underset{a_j \in \mathbb{R}^{\left|T_j\right|}}{\operatorname{argmin}}\left\{\begin{array}{c}
\mathcal{L}_j\left(a_j, \bar{\mu}_j[\tau], \bar{\nu}[\tau]\right) \\
+\frac{1}{2 \rho}\left\|a_j-a_j[\tau]\right\|_2^2
\end{array}\right\} \\
& \mu_j[\tau+1]=\mu_j[\tau]+\rho \gamma_j \tilde{G}_j a_j[\tau+1] \\
& \bar{\mu}_j[\tau+1]=\mu_j[\tau+1]+\rho \hat{\gamma}_j[\tau+1] \tilde{G}_j a_j[\tau+1]
\end{aligned}

Communicate $a_j$ for all $j \in[K]$ with neighbors

\begin{aligned}
& \nu_j[\tau+1]=\nu_j[\tau]+\rho \gamma_j[B]^{O_j} a[\tau+1] \\
& \bar{\nu}_j[\tau+1]=\nu_j[\tau+1]+\rho \hat{\gamma}_j[\tau+1][B]^{O_j} a[\tau+1]
\end{aligned}

Communicate $\bar{\nu}_j$ for all $j \in[K]$ with neighbors.

The primal and dual variables are initialized as follows, $\forall j \in [K]$:

\begin{aligned}
& a_j[0] \in \mathbb{R}^{|| T_j \mid} \\
& \mu_j[0]=\rho \gamma_j \tilde{G}_j a_j[0] \\
& \bar{\mu}_j[0]=\mu_j[0]+\rho \hat{\gamma}_j[0] \tilde{G}_j a_j[0] \\
& \nu_j[0]=\rho \gamma_j[B]^{O_j} a[0] \\
& \bar{\nu}_j[0]=\nu_j[0]+\rho \hat{\gamma}_j[0][B]^{O_j} a[0]
\end{aligned}

§.§ NST-PAC algorithm

This work employs an enhanced, accelerated version called NST-PAC developed in Ferro2022AHubs. It is a primal-dual method incorporating both $L2$ and proximal regularization terms. The convergence speed is increased by using time-varying gains and Nesterov-accelerated gradient updates for both the primal and dual variables. The iterative NST-PAC algorithm consists of the following steps at each iteration $\tau$:

\begin{align}
a_j[\tau+1] & = \underset{a_j}{\operatorname{argmin}}\big\{\mathcal{L}_j\left(a_j, \hat{\eta}_j[\tau], \hat{\nu}[\tau]\right)\\
& +\frac{\rho_j \gamma_j}{2}\left\|G_j a_j-b_j\right\|_2^2 + \frac{\rho_j \gamma_j}{2}\left\|B_j a_j\right\|_2^2 \nonumber\\
& +\frac{1}{2 \rho_j}\left\|a_j-a_j[\tau]\right\|_2^2\big\} \nonumber\\
\hat{a}_j[\tau+1] &= a_j[\tau+1]+\alpha_j[\tau+1]\left(a_j[\tau+1]-a_j[\tau]\right) \nonumber\\
\eta_j[\tau+1]& = \hat{\eta}_j[\tau]+\rho_j \gamma_j\left(G_j \hat{a}_j[\tau+1]-b_j\right) \nonumber\\
\hat{\eta}_j[\tau+1]&= \eta_j[\tau+1]+\phi_j[\tau+1]\left(\eta_j[\tau+1]-\eta_j[\tau]\right) \nonumber\\
\text {Communicate } & \hat{a}_j \text { for all } j \in[K] \text { with neighbors } \nonumber\\
\nu_j[\tau+1]&= \hat{\nu}_j[\tau]+\rho_j \gamma_j B_j \hat{a}_j[\tau+1] \nonumber\\
\hat{\nu}_j[\tau+1]&= \nu_j[\tau+1]+\theta_j[\tau+1]\left(\nu_j[\tau+1]-\nu_j[\tau]\right)\nonumber \\
\text {Communicate } & \hat{\nu}_j \text { for all } j \in[K] \text { with neighbors } \nonumber
\end{align}

The algorithm further protects privacy by masking both the primal and dual variables. Masking is implemented by using iteration-varying and atom-specific parameters $\alpha_j[\tau]$, $\phi_j[\tau]$ and $\theta_j[\tau]$. Masking the dual variables (or shadow prices), in particular, is desirable since these may reveal sensitive data related to costs, operating constraints, or other preferences of SMOs. Instead, masked variables $\hat{a}$ and $\hat{\nu}$ are exchanged between atoms. By iteratively solving the local, decomposed optimization problems across all SMOs, NST-PAC (and PAC) provably converge to the globally optimal ACOPF (relaxed) solutions for the whole primary feeder Romvary2022AOptimization,Ferro2022AHubs.

§ FEDERATED LEARNING

When it comes to DER forecasting, the challenges of privacy as well as the requirement of large training data sets, can be met using a distributed machine learning paradigm, federated learning (FL) yang2019federated,li2020federated,lim2020federated,buildFL2020. FL is a machine learning framework where each device participates in training a central model without sending actual data, but only exchanges gradient information in the training phase and sends prediction estimates during deployment. A general overview of the proposed DER prediction process is shown in <ref>. In the figure, various IoT devices including those at a house level such as smart thermostats and smart washers, and energy-producing devices such as PV and EVs are considered. The future smart grid will include a wider range of devices that will be capable of computation and communication. The house-level devices are grouped under a home energy manager $H_i$ to enable aggregation at a house level while energy-producing devices such as EV and PV (grouped under $E_i$) are assumed to directly participate in a transactive environment. Both $H_i$ and $E_i$ can be considered typical IoT-based DERs connected to a power grid, whose energy consumption needs to be predicted. Using a communication infrastructure, our goal is to exchange information between the DERs in the bottom local layer and the global decision makers at the top layer in a private and secure manner so as to lead to an accurate prediction of the DER consumption/generation. The DER prediction is then utilized to formulate a grid service, to mitigate the load swings and  that occur in the distribution grid due to a lack of situational awareness. The FL-based DER prediction is used to anticipate the load swings and mitigate them by proactively controlling the DERs venkataramanan2022forecast.

An overview of the DER prediction process using federated learning

A typical process of DER-forecast can occur in the following manner. Collect the input-output pair $[x_t, \hat{P}^t(T)]$ for a federate $F_i$ for several samples $n$. The features used in this work are , where  denotes the actual power consumption, and $m$ denotes the minutes prior to time $T$. The number of samples $n$=2880, was obtained by collecting data every 15 minutes over a period of 30 days. The overall training schematic of the FL-based neural network is shown in <ref>.

The schematic of neural network training using federated learning is shown here. The steps (1)-(6) are repeated until $L^t \leq \epsilon$

§ DETAILED COMMUNICATION STEPS AND TIMELINE

See <ref> for more details on the communication among different market agents, operators, and resilience managers.

Diagram showing a more detailed communication scheme and steps for information exchange between the various market operators, agents, and resilience managers at the secondary and primary levels.

§ MITIGATION USING THE MARKET OPERATORS AND RESILIENCE MANAGERS

Our market framework consisting of the SM and PM provides situational awareness (SA) in the form of available power injections at various nodes at the primary and secondary levels. Once the market is cleared, during execution, the actual injections from the SMA and PMA are monitored by the SRM and PRM, respectively (see <ref>). These injections are then utilized by these managers to compute commitment scores, trustability scores, and resilience scores (RS), as shown in <ref> and <ref>. In what follows, we will show how the SA from the market operators and the RS from the resilience managers can be utilized to mitigate all the attacks described in <ref>.

As a result of continuous monitoring, any unexpected deviation from the agents' nominal performance in the form of change in the net injection at the PCC, raises a flag. Any such flag makes the operators shift from the nominal operating mode to the resilience mode. Minimal visibility regarding actual injections from all PMA is assumed to be available. Rather, we assume that each SRM only locally observes the actual injection from the corresponding SMA, and each SRM communicates that information to the PRM. More importantly, the attack scenarios considered also assume that this important communication to the PRM from all SRMs is completely sabotaged (as was the case in the Ukraine 2015-16 attacks). Despite this loss of communication, the PRM is able to step in and mitigate the attack as the flag raised is independent of this communication loss and is due to a physical impact of the agents' deviation from nominal performance. Subsequently, the PMO redispatches trustworthy PMAs so as to bring the power import from the bulk grid down to pre-attack levels. The new setpoints for the PMAs/SMOs are in turn suitably disaggregated to compute new setpoints for the SMAs through a re-dispatch by the SM. Before proceeding to the results, we propose a specific mitigation strategy that leverages the SA provided by our approach.

Timeline of attack detection and mitigation.

§.§ Algorithm (A) for redispatch by the PMO in a balanced network

We first consider a balanced, equivalent single-phase network using the BF model. The starting point for the overall mitigation sequence is the awareness that an attack has occurred. This is realized by the PRM in the form of a change in the net load from $P_{PCC}$ to $\overline{P}_{PCC}$, which denotes the net load from the entire primary feeder at the substation before and after the attack, respectively. This can be detected by the PRM at the substation or point of common coupling (PCC) since this is the power it imports from the main transmission grid. We propose a redispatch algorithm that the corresponding SMOs can carry out based on the ratio between these two values. To describe this redispatch algorithm, we begin with the cost function in <ref> for the PM ACOPF problem. For ease of exposition, this can be rewritten in a simplified manner as:

\begin{gather}
\sum_{i=1}^n \left(\frac{1}{2} \; \alpha_i P_i^{G^2} + \beta_i \left(P_i^L - P^{L0}_i\right)^2\right) + \xi \cdot losses \label{eq:update_obj}\\
\overline{\alpha}_i = \Delta_{\alpha} \alpha_i, \; \overline{\beta}_i = \Delta_{\beta} \beta_i, \; \overline{\xi} = \Delta_{\xi} \xi; \; \alpha, \beta, \xi, \Delta > 0 \label{eq:update1}\\
\Delta_{\alpha} = \Delta_{\beta} = \frac{|P_{PCC}|}{|\overline{P}_{PCC}|}, \; \Delta_{\xi} = \frac{|\overline{P}_{PCC}|}{|P_{PCC}|} \label{eq:update2}
% \alpha, \beta, \xi, \Delta > 0 \nonumber
\end{gather}

We note that a change in the power import from the main grid causes $\Delta_\alpha, \Delta_\beta, \Delta_\xi$ to deviate from unity. Suppose that several distributed local generator SMOs are attacked, as in Attack 1(a). This would increase net feeder load, i.e.  $|\overline{P}_{PCC}| > |P_{PCC}|$ (note that both $P_{PCC}, \overline{P}_{PCC} < 0$ since net loads are negative injections), thus causing $\Delta_\alpha<1$. Applying this cost coefficient update would lower the cost coefficients from $\alpha_i$ to  $\alpha_i^\prime$. This results in dispatching more local generation from remaining online SMOs instead of importing power from the bulk grid. As the SMOs also have information about the flexibility in their corresponding SMAs in the form of $\delta P^*, \delta Q^*$ (see <ref>), the overall hierarchical PM-SM market structure automatically provides the solutions of the new dispatch. Similarly, a value of $\Delta_{\beta} < 1$  reduces the disutility coefficients to encourage more demand response via load shifting and/or curtailment, by utilizing the downward flexibility provided by the SMOs bidding into the PM, and subsequently also by the SMAs bidding into the SM. In contrast to these two values, when the net import from the main grid increases, then $\Delta_{\xi} > 1$ penalizes electrical line losses more heavily in the objective function. As a result, the redispatch discourages imports from the transmission grid in favor of dispatching more local DERs. This is because distribution grids are more lossy (have higher resistance to reactance ratios), and hence prioritizing the loss minimization makes it more efficient to utilize local generation closer to the loads being served.

After deriving the multiplicative coefficient update factors $\Delta_{\alpha}, \Delta_{\beta}, \Delta_{\xi}$, the PRM can broadcast these common values to all the SRMs simultaneously, who in turn send them to their corresponding SMOs. The SMOs update each of their objective function coefficients using these factors and then perform distributed optimization to redispatch the PM, resulting in new $P$ and $Q$ setpoints for SMOs, along with new nodal distribution LMPs (d-LMPs). This is followed by each SMO also re-dispatching their SM, in order to disaggregate the new setpoints among their SMAs. A timeline of the key events is shown in <ref>.

§.§ Algorithm (B) for redispatch in an unbalanced, 3-phase network

For the unbalanced 3-phase case, we use a modified algorithm for the coefficient update. The update rule here is more sophisticated since in this case, the variables are now 3-phase vectors rather than scalars.

\begin{gather}
\Delta=\mathbf{P}_{PCC}-\overline{\mathbf{P}}_{PCC} \\
Z_i\left(\delta_i\right)=1+\frac{RS_i \Delta^{\top} \delta_i}{\mu \sum_i RS_i} \Longrightarrow \gamma_{i \delta}=\frac{1}{Z_i\left(\delta_i\right)} \\
\overline{\boldsymbol{\alpha}}_i=\gamma_{i \alpha} \boldsymbol{\alpha}_i, \quad \overline{\boldsymbol{\beta}}_i=\gamma_{i \beta} \boldsymbol{\beta}_i, \quad \overline{\boldsymbol{\xi}}=\left(\frac{\sum_i \gamma_{i \alpha}+\gamma_{i \beta}}{2n}\right)^{-1} \boldsymbol{\xi} \label{eq:update_ci_model}
\end{gather}

Note here that $\textbf{P}_{PCC}, \overline{\textbf{P}}_{PCC}$ are the 3-phase power imports from the tie line before and after the attack. $\boldsymbol{\alpha}_i, \boldsymbol{\beta}_i$ are $3\times 1$ vectors representing cost and disutility coefficient for each phase at SMO node $i$, and $\boldsymbol{\xi}$ is a 3-phase hyperparameter that penalizes line losses in the objective function. A DG attack that increases net load would result in $\gamma_{i\alpha}, \gamma_{i\beta} < 1$ and $\overline{\boldsymbol{\xi}} < \boldsymbol{\xi}$. Thus, these coefficient updates work using a similar intuition to algorithm (A) in that it favors local DER generation and load flexibility over transmission imports. A key difference here is that the PRM also takes into account the RS of each SMO during the redispatch so that it relies more heavily on resilient SMOs for attack mitigation. The PRM updates the coefficients $\alpha_i, \beta_i$ and $\xi$ to $\alpha_i^\prime, \beta_i^\prime$ and $\xi^\prime$, and sends the new coefficient values to all SRMs. The SRMs send these new objective functions to the corresponding SMOs, and the rest of the mitigation procedure follows in the same manner as in the previous section.

§ RESILIENCE-DRIVEN RECONFIGURATION ALGORITHM

All possible shortest paths are computed between each generation source and critical load pairs present within the system using the graph network. If the generation is not enough to supply the total critical load, then the algorithm searches for the next available generation. This will continue until the critical load demand is met. As the generation sources are assigned to critical loads, if any source's capacity is more than the assigned load, the source's partial remaining capacity will be utilized for other loads. Once we have all the feasible paths for reconfiguration, we will compute the resiliency metric for each path (see <ref> for how to compute the resiliency metrics and $PNR$), which will support the operator in finding the most resilient path to restore. The reconfiguration paths will be determined based on the stress levels of the grid and the corresponding degree of the failed SMA node, the tolerance bands and flexibility of the ICA, and the security levels and privacy needs of the SMA. This process is outlined in <ref>.

Resilience-based reconfiguration algorithm.

§ VALIDATION PLATFORMS

§.§ PNNL

Since the proposed framework is intended to be deployed in a large-scale distribution system with a high penetration of renewables, we choose our first validation platform to be driven by HELICS (Hierarchical Engine for Large-scale Infrastructure Co-simulation), an open-source cyber-physical energy co-simulation framework for energy systems trevor2024Helics. As the core engine of the platform, HELICS provides time management and data exchanges between the simulators, also known as federates. Moreover, through standard procedures and application programming interfaces (APIs), data exchange between federates is performed either as values or messages. Within this platform, GridLAB-D is the distribution system simulator used to simulate all of the use cases we discussed above. GridLAB-D is versatile, with the capability of simulating large 3-phase unbalanced distribution systems, with agent-based and information-based modeling tools and extensive data collection tools for end-use technologies and interface APIs for co-simulation Chassin_4517260. Thus, it allows us to modify the standard IEEE 123-node system to include (1) residential (either single or multi-family) and commercial buildings with or without heating, ventilation, and air conditioning (HVAC) systems, (2) edge devices with IoT connectivity including HVAC type appliances, small electronics and lighting, (3) distributed generators (DGs) such as photovoltaic (PV) panels, battery energy storage systems (BESSs), and diesel generators. The distribution model is also assumed to have all its loads connected via smart meters with load-shedding capability to regulate energy demand according to a distribution system operator’s command. Specifically, the IoT-enhanced version of the IEEE 123-node system consists of the assets detailed in <ref>.

GridLAB-D IEEE 123-node test feeder features with IoT-enhanced model.

Number     Capacity

Standard IEEE 123-node test feeder     Spot loads     $85$     $3,985.7$ kVA
EUREICA IEEE 123-node test feeder     Houses - Demand response (HVACs in all, WHs in 348)     $1,008$     variable ($4$ KW avg/house) ($20\%$ to $30\%$ critical)
Distributed generators (DGs)     $380$     $1,745.8$ kVA ($\approx44\%$ system penetration)
PVs     $207$     $880.84$ kVA
BESSs     $173$     $865$ kVA
Community PVs     $12$     $96$ kVA

§.§ NREL

The objective of validating the EUREICA framework at NREL is to evaluate the feasibility of implementing the framework in real-time. Since electrons flow in the grid in real-time, it is critical that the operations proposed should also function in real-time, and be in compliance with operational requirements. The Advanced Research on Integrated Energy Systems (ARIES) at NREL is a cutting-edge virtual emulation platform that encompasses actual Distributed Energy Resource (DER) hardware systems, such as wind turbines, photovoltaic (PV) arrays with controllers, batteries, and storage systems kurtz2021aries. The digital real-time simulation (DRTS) cluster of ARIES is used for validating the performance of the EUREICA modules. The overall validation platform is shown in <ref>. It constitutes 5 components - (i) IoT device virtualization (using Raspberry Pis and Typhoon HIL to characterize IoT devices), (ii) Communication emulation (using analog and network connections to emulate communication at the speed of actual communication in the field) (iii) hardware-in-the-loop interface (to provide increased fidelity for components under study) (iv) digital real-time simulation (using RTDS and Typhoon HIL to emulate the power grid in real-time) (v) Time synchronization (to bring together the hardware components using a time server to ensure accuracy of simulation). This validation platform is used to determine the performance of the EUREICA framework for all the modules.

ARIES-DRTS Validation Platform at NREL

§.§ LTDES

A training simulator-based platform was also used to validate the EUREICA framework. In particular, the General Electric (GE) ADMS DOTS (Advanced Distribution Management System-Distribution Operations Training Simulator), used for training operators and dispatchers, was used as the validation platform. Rather than users waiting to experience challenging events on the job, dispatchers are able to familiarize themselves with advanced application functionality and gain an understanding of how they interact with other subsystems of the ADMS. Into this ADMS-DOTS we introduce DERIM, a Distributed Energy Resource Integration Middleware, an interface that allows integration of various DERs with the ability to communicate as dictated by the EUREICA framework (see <ref>). This integrated system uses the same software components, programmatic and user interfaces as the real-time ADMS, and creates an effective training and testing environment to operate with the actual network model, data, and functions in a controlled and safe environment. <ref> shows an example of what the validation process looks like for attack 1a.

DERIM interface with ADMS-DOTS in the LTDES validation platform.

Attack 1a validation process workflow in the LTDES validation platform with DERIM and ADMS-DOTS.

§ KUNDUR 2-AREA SYSTEM

<ref> is a diagram of the Kundur 2-area transmission system commonly used as a test case to study dynamic stability, power interchange, oscillation damping, etc. The system contains 11 buses, four generators, and two areas. The two areas are connected with weak tie lines kundur2007power.

Schematic of Kundur 2-area power system

§ NUMERICAL SIMULATION SETUP FOR MARKETS

All use cases considered are based on an IEEE 123 test feeder (see <ref>), which is radial, unbalanced, and multi-phase. The feeder was modeled in the GridLAB-D environment (see SI <ref> for more details) and augmented to have a high penetration of DERs. For all attacks except attack 3, we assumed that the switch settings are in their nominal positions such that we have one primary feeder, with 85 active nodes with SMOs/PMAs (out of the 123 in total). A PMO was assumed to be at the slack bus (substation), at either 115 or 69kV, with the SMOs at 4.16kV, and each SMA at 120-240V. The flexibility bids for the SMAs and the SMOs were randomly generated, allowing each to offer flexibilities of up to $\pm 30\%$ around their baseline power injections Olsen2014Demand2020_2. We used 5-minute real-time market LMPs from the California ISO and assumed the Q-LMP to be 10% of the P-LMP. Note that for all attack scenarios except attack 3, we used the CI model to represent the feeder as is.

For attack 2, however, we considered a modified version of the feeder and deployed the BF model instead. Here, we modified the original IEEE-123 feeder to consider a case where we have a few large distributed generators (PV, batteries, diesel generators) concentrated at just five primary feeder (SMO) nodes numbered 25, 40, 67, 81, and 94. This is in contrast to the other attacks where there were instead a larger number of smaller DERs distributed throughout the network. Another distinguishing factor of this scenario is that the originally unbalanced feeder was converted to an equivalent balanced 3-phase model by (i) assuming all switches to be at their normal positions, (ii) converting single-phase spot loads to 3-phase, (iii) assuming cables to be 3-phase transposed, (iv) converting configurations 1 thru 12 to symmetric matrices and (v) modeling shunt capacitors as 3-phase reactive power generators Haider2020TowardGrids. Each SMO was assumed to have between 3-5 SMAs with the number chosen uniformly at random. Since the injection data in the original IEEE-123 model was only available up to the primary feeder node level, we artificially randomly disaggregated the injections at each SMO amongst its SMAs, which could be net loads or generators.

We then performed a co-simulation of both the PM and SM for all attack scenarios. We refer the reader to nair2022hierarchical,nair2023local for the behavior of this market structure for a nominal scenario when there is no attack. In what follows, we only consider the three attack scenarios described above.  We also note that our flexibility bids were synthetically created, so the resulting flexible ranges in our simulations may be quite large at times and not realistic in some cases. However, our proposed framework can be generally applied to cases where there is less DER flexibility as well.

§ MARKET SIMULATIONS OF OTHER ATTACKS

§.§ Mitigation of Attack 1b

We note that in Attack 1b, there is a loss in net generation, and therefore the power imported from the bulk grid increases. It is also assumed that the communication from all  SRM to the PRM is disrupted, while the communication from PRM to SRM remains intact. That is, the PRM loses observability but is still able to communicate the redispatch of the new coefficients to the SRM. We do not consider the case when such observability is not lost, a discussion of which is beyond the scope of this paper. With the redispatch, the PM-SM framework identifies all of the new trustable PMAs through the SA computations described in Section 3 with the overall power balance met at all points in the distribution grid.

The steps in mitigating this attack are as follows. Due to the attack, 45 kW of net-generation compromised as shown in <ref>. The PMO alerts other trustable PMAs/SMOs to redispatch their generation assets in the PM. Trustable PMAs/SMOs will curtail flexible loads to respond and mitigate the attack as in <ref>. The redispatch is also influenced by the resilience scores of different SMOs over time shown in <ref>. SMOs redispatch the SM which provides correct setpoints to all their SMAs. An example of this in <ref> shows how the SMO at node 35 disaggregates its new setpoint amongst its 3 SMAs.
As a result of mitigation, the total import from the main grid stays at the same level as shown in <ref>.

Generation with and without attack 1b.

Feeder power import.

Effects of attack 1b on SMO net generation and power import.

Dis-aggregation of changes in the setpoints for SMO (from the PM) at node 35 across its 3 SMAs (in the SM), resulting from attack 1b mitigation, along with each SMA's RS. All 3 SMAs are on phase A.

Curtailment of flexible loads for attack 1b mitigation.

Locational-temporal trends of RS across all flexible SMO nodes and over the whole simulation period of 24h.

§.§ Mitigation of Attack 1c

Attack 1c is a more distributed attack where individual SMAs representing secondary feeders are attacked directly. We considered a case where a large number of DGs including solar PV and batteries are attacked. A total of 53 SMA nodes with DGs were compromised and taken offline, resulting in a total loss in generation capacity of 157 kW. This leads to a decrease in the net injections across all the SMOs as seen in <ref> - there are no longer any SMOs with net generation after the attack and the loss of local generation also increases the net load at the SMOs. This leads to an increase in power import from the main transmission grid as in <ref>.

SMO injections.

3-phase power imports from the main grid.

Effects of attack 1c on SMO injections and power import.

In the case of all other attacks, the mitigation strategy involves the PM redispatch occurring 1st followed by the SM redispatch. There, only the PM is directly involved in attack mitigation while the SM is only used to disaggregate the new SMO setpoints amongst their SMAs. However, in the case of attack 1c, the SM redispatch occurs first at the secondary feeder level and is then followed by the PM redispatch at the primary feeder level. Thus, both the SM and PM are actively involved in attack mitigation here. We see in <ref> that we can partially mitigate the attack by leveraging the flexibility of SMAs in the SM. However, SM mitigation alone is not sufficient. We need to also utilize the inter-SMO flexibility in the PM to fully mitigate and restore the feeder import back down to the pre-attack level. A summary of the attack metrics is shown in <ref>.

Attack 1c summary.

Power import [kW]     Total net load [kW]

Pre-attack     1412           1457

Post-attack        1722       1716

SM mitigation only        1553       1547

SM + PM mitigation        1422       1417

We also compare the flexibility bids of the SMOs before and after the attack in <ref>. As expected, the net load of the bids generally increases across all SMOs due to the loss of local DGs at their respective SMAs. However, by leveraging their SMA flexibilities, the SMOs are still able to offer some flexibility to the PM to help mitigate the attack.

Comparison of SMO flexibility bids into the PM before and after the attack. The dashed and dotted lines indicate the baseline values while the shaded regions are the flexibility bids around the baseline.

§.§ Mitigation of Attack 2a

This corresponds to a case where there are five large distributed generators in the modified IEEE-123 system, one of which (at SMO node 94) is taken offline. Here, we see that the remaining four SMO nodes (25, 40, 67, 81) have more than enough remaining generation capacity to meet the shortfall caused by the attack. Without mitigation, the attack would have resulted in an additional import of about 261 kW from the main grid. However, by utilizing the upward flexibility of remaining SMOs, we are able to fully resolve the attack and bring the total power imported back to pre-attack levels. The left figure in <ref> shows the results of the PM dispatch before the attack and after attack mitigation for the five key SMO nodes of interest. The plot also shows the SMO's bids into the PM, with the dashed blue line being the baseline injection bid and the blue-shaded region representing the upward/downward flexibility around it. The right figure shows the results of the SM re-dispatch after the attack mitigation and PM re-dispatch for SMO 67 as an example. It disaggregates the new setpoints among its three SMAs, with SMA 1 being a net load while SMAs 2 and 3 are net generators.

Mitigation of small-scale attack 2a.