# Operational and Economy-Wide Impacts of Compound Cyberattacks and Extreme
Weather Events on Electric Power Networks

Charalampos Avraam, Luis Ceferino, Yury Dvorkin

## Abstract

The growing frequencies of extreme weather events and cyberattacks give rise
to a novel threat where a malicious cyber actor aims to disrupt stressed
components of critical infrastructure systems immediately before, during, or
shortly after an extreme weather event. In this paper, we initiate the study
of Compound Cyber-Physical Threats and develop a two-stage framework for the
analysis of operational disruptions in electric power networks and economy-
wide impacts under three scenarios: a Heatwave, a Cyberattack, and a Compound
scenario when the Cyberattack is timed with the Heatwave. In the first stage,
we use a bilevel optimization problem to represent the adversarial rationale
of a cyberattacker in the upper level. In the lower level, we model
disruptions in the electric power network using an optimal power flow model.
In the second stage, we couple the disruption of electricity supply with a
Computable General Equilibrium model to elucidate the impacts on all economic
sectors. For the New York Independent System Operator, we find that a 9%
demand increase in a Heatwave may not lead to unserved load. The Cyberattack
can lead to 4% of unserved electric load in Long Island, while the Compound
scenario can increase unserved electric load in Long Island to 13% and affect
almost 600,000 customers. Our results show that the activity of state and
local government enterprises can decrease by 30% in the Compound scenario. We
conclude that the vulnerability of federal, state, and local government
enterprises to electricity disruptions can affect a broad range of
populations.

##### Keywords

heatwave, cyberattack, compound risks, CGE, bilevel optimization, cyber-
physical systems

##  1 Introduction

Information communication and control (ICC) technologies enhance the
efficiency and flexibility of infrastructure operations [18], but expose
critical infrastructure components to cyberattacks. In infrastructure systems
with strong coupling of ICC technologies and physical operations, malware
infiltration can result in large disruptions [47]. According to the Industrial
Control Systems Cyber Emergency Response Team, the annual number of data
breaches and exposed records in the U.S. grew sixfold and twofold,
respectively, between 2005 and 2020. [70]. Across all economic sectors,
cyberattacks increase in frequency and intensity [30], but energy
infrastructure is at particular high risk. Energy was the second most targeted
economic sector in 2015 [70].

Cyberattacks can be more likely and damaging when infrastructures are facing
stresses than in normal operating conditions. For example, the U.S. Federal
Bureau of Investigation noted in 2022 that cyberattackers may target U.S. food
infrastructure more aggressively during the critical planting and harvest
seasons, e.g., February-March and September-October for grains [26]. Since the
surge of COVID-19, data breaches against healthcare organizations have been
persistently increasing, compared to 2018, and affected a record-high of 45
million individuals in 2021 [15]. The more recent attack against oil terminals
in Germany and in the Amsterdam-Rotterdam-Antwerp area in late January and
early February 2022 [75] happened during Europe’s peak gas demand season [52].

Extreme weather events also damage critical infrastructures. For example, they
damage electric power substations, utility poles, and distribution lines [13,
14, 53]. Damaged electric power components can disrupt drinking water pumps,
sewage pump stations, telecommunications cellular tower sites, causing
substantial economy-wide financial losses. Hurricane Isaias left almost
800,000 utility customers without power in New York in 2020 [41] and Hurricane
Ida more than 1.2 million customers in Louisiana in 2021 [34], largely due to
damaged utility poles. A substation explosion caused by Hurricane Sandy left
most of Lower Manhattan without power for more than four consecutive days in
2012 [49], affecting also the transportation and telecommunications
infrastructure systems [69].

Disruptions in the electric power sector can also propagate to food
infrastructure systems, e.g., in the case of Winter Storm Jonas in 2016 [8].
In addition, disruptions in the electric power sector propagate across sectors
and cause economy-wide losses. Estimates of the economy-wide losses of the
2003 Northeastern Blackout range between $4-10 billion [68]. Nationally, the
10-year economic losses in the US due to extreme weather events have
quadrupled between 1980-2020, averaging $120 billion per year in 2016-2020
[65].

The growing frequency of extreme weather events and cyberattacks manifest a
novel threat where a malicious cyber actor disrupts stressed components of
critical infrastructure immediately before, during, or shortly after an
extreme weather event. Therefore, this paper aims to initiate the
investigation of compound cyberattacks and extreme weather events, hereinafter
referred to as Compound Cyber-Physical Threats, in critical infrastructure
systems. We focus on the electric power sector due to its importance for the
operation of other critical infrastructure systems and its vulnerability to
cyberattacks and extreme weather events. Recognizing the potential of Compound
Cyber-Physical Threats, the U.S. National Guard in coordination with the City
of Houston in 2018 and the State of Indiana in 2021 designed and conducted
three-day drills of a cyberattack striking during an extreme weather event
[31, 33]. However, there does not exist a framework to (a) understand the
compounding effect of a cyberattack timed with an extreme weather event
compared to a cyberattack under normal operating conditions, (b) identify the
spatial distribution of critically-vulnerable infrastructure components within
an electric power network, and (c) generate realistic, synthetic scenarios of
Compound Cyber-Physical Threats for service disruptions and their impacts on
the economy. Thus, this paper investigates the following research questions:

  * •

What is the disruption potential of Compound Cyber-Physical Threats for
electric power networks?

  * •

Which regional electric power network components are more vulnerable to a
cyberattack timed with an operational disruption caused by an extreme weather
event?

  * •

How do operational disruptions caused by a Compound Cyber-Physical Threat
against electric power network components affect the activity of other sectors
in the economy?

In this paper, we identify vulnerable electric power network components across
New York State and within New York City, and assess economy-wide operational
and price disruptions due to Compound Cyber-Physical Threats across sectors.
We focus on New York City because it is the most densely populated city in the
U.S. [67], hence a regional disruption of electricity supply can impact a
larger population and cause more severe cascading failures to interdependent
infrastructure systems compared to other cities. Moreover, the New York State
is the most exposed U.S. state to physical and cyber risks according to
Lloyd’s and the Cambridge Centre for Risk Studies [16] that measure the impact
of all risks on the Gross Domestic Product (GDP@Risk in Billion US$).

The assessment of the operational disruption of a Compound Cyber-Physical
Threat requires integrating an electric power network model capturing the
physical and economic principles of the electric power network operations in
New York State, with a cyberattacker model, capturing the adversarial
rationale of the cyberattacker. We use an Optimal Power Flow model (OPF)
model, which is representative of operational practices in the New York
Independent System Operator (NYISO) electric power network, to capture the
response of the system operator [22] to a heatwave. We focus on heatwaves
because they stress the grid by increasing electricity demand. Heatwaves will
be exacerbated by global warming [57], rendering them more costly [17] and
dangerous for vulnerable populations [39, 19] in future climates. A bilevel
optimization problem models the adversarial rationale of a cyberattacker
anticipating the system operator’s response [12, 36].

To further assess the economy-wide impacts, we couple the cyberattacker
bilevel optimization problem with a Computable General Equilibrium model
(CGE), which can capture flow of money and goods across all sectors in an
economy [60]. The assessment of sectoral disruptions caused by natural hazards
[58], terrorist attacks [59] and cyberattacks [20, 23] in existing literature
is exogenous and is based on a series of assumptions. For example, Rose et
al., [59], assume that completely disconnecting the transmission lines between
Los Angeles County and the rest of the electric power network results in 100%
unserved load for two weeks. The authors then constrain electricity supply in
a CGE to simulate the response of all other sectors. Instead, in our
framework, unserved load is the result of the bilevel optimization problem.
Following a failure, our framework considers the physical response of the
electric power network and the decisions of the system operator to derive
unserved load. For example, in an attack against transmission lines in our
framework, the system operator may choose to operate the electric power
network in islanded mode to support parts of the network using local
generation and avoid a complete blackout. Then, we inform the constraints on
regional electricity supply in a CGE in Rose et al., [59] based on the results
of the bilevel optimization problem. Hence, the coupling extends the existing
literature by relaxing the assumptions on the operational disruptions of
natural hazards and cyberattacks. We do not aim to quantify disaster impacts
in the economy, as in Markhvida et al., [46], but understand the response of
all economic sectors to operational disruptions in the electric power network.

The rest of the paper is organized as follows. Section 2 describes the
assumptions behind the Compound Cyber-Physical Threat scenarios and Section 3
introduces the proposed modeling framework and methodology. We identify
electric power network and inter-sectoral outcomes under combinations of
extreme weather events and cyberattacks in Section 4 and discuss resilience
implications and limitations in Section 5. We conclude in Section 6 with the
summary of the main findings and future research.

##  2 Scenario Design

We first describe the Baseline scenario, which will serve as a reference for
the Cyberattack, Heatwave, and Compound Cyber-Physical Threat scenarios. Among
different extreme weather events, New York City is particularly vulnerable to
Heatwaves [56]. First, although winter outages result in greater welfare
losses in the current climate and electrification levels, summer outages due
to Heatwaves can be more costly under warmer future climate [17]. Second,
higher temperatures increase power demand directly due to air conditioners
[7], which in turn strains the electric power network. The heat-driven strain
is particularly critical in power network which have their annual peak
consumption in summer, e.g., the NYISO network. Third, Heatwaves, or extreme
heat events, disproportionately affect the elderly [39], low-income [19] and
other vulnerable populations [21], Finally, densely populated cities,
including New York City, can become urban heat islands during summer, i.e.,
the cities are subject to increased temperatures compared to neighboring rural
areas [63]. In what follows, we describe the assumptions behind each scenario.

##### Baseline.

We obtain the 2019 installed capacity and electricity demand for each NYISO
zone (Figure 1) and an hourly summer load-duration curve for NYISO zones from
[37, 50]. We retrieve the capacity and impedances of transmission lines of the
NYISO electric power network from [37]. To the best of our knowledge, a 2019
estimate of the number of NYISO customers is not available. Therefore, in what
follows, we use a total number of customers equal to 19.8 million, reported by
the Federal Energy Regulatory Commission (FERC) on July 2020 [28].

Figure 1: Left: NYISO installed capacity in 2019 by generation technology.
Right: Maximum Demand and installed power generation capacity by NYISO Zone.

##### Heatwave.

We base the formulation of the Heatwave scenario on the 2022 analysis by
Rommitti and Wing [57]. The authors study the impact of increasing
temperatures on annual electricity demand by 2050 for 36 cities, including New
York City, using projections of future temperature from 21 global climate
models [35]. For New York City, they find that yearly mean peak electricity
demand can increase approximately by 5.6%, and the maximum increase can be
7.5%. However, actual and predicted peak electricity demand can differ. For
example, the short-term prediction of peak electricity demand in each node of
the NYISO deviated by 16% from the actual demand in 2020 [51]. Therefore, in
our Heatwave scenario, we assume that the increase in demand is 20% higher
than the maximum calculated increase in [57], i.e., a $7.5\cdot 1.2=9\%$
increase in hourly electricity demand for all NYISO zones.

##### Cyberattack.

The cyberattacker aims to maximize the operational disruption to the electric
power network by compromising power generation capacity, transmission
capacity, and limits on voltage angle differences. Since generation,
transmission, and end-use components use Supervisory Control And Data
Acquisition (SCADA) devices and Programmable Logic Controller (PLC) interfaces
[2], the cyberattacker can modify SCADA and PLC firmware [40, 72] and hardware
[73, 74] to execute unauthorized code. In our framework, the cyberattacker can
inject spoof measurements and modify control algorithms and messages exhanged
SCADA-SCADA, SCADA-PLC, PLC-PLC communications to constrain and manipulate
regional generation, flow and angle difference values and limits. The NYISO
oversees transmission components and participates in the cybersecurity
initiatives of the North American Electric Corporation [24]. Individual
generators do not necessarily participate in the North American Electric Power
Corporation. For that, we assume that generators are more vulnerable to
cyberattacks in the default Cyberattack scenario. Equivalently, it is
significantly more difficult, i.e., costlier, for a cyberattacker to
compromise a transmission component compared to generation capacity in the
default Cyberattack scenario where $\tilde{c}_{f}=5\cdot\tilde{c}^{g}$.
Because our knowledge is limited regarding the cyberattacker capabilities, we
conduct sensitivity analysis on the available cyberattacker resources
$(\tilde{b})$ and the relative difficulty of compromising transmission
components versus power generation capacity
$(\frac{\tilde{c}^{g}}{\tilde{c}^{f}})$.

##### Compound Cyber-Physical Threat.

In the Compound Cyber-Physical Threat scenario, the Cyberattack is timed with
the Heatwave. We assume that the cyberattacker has access to heatwave
advisories or predictions, which are possible up to 72 hours ahead of the
event [54] to prepare and launch the cyberattack. Moreover, we assume that the
cyberattack is perfectly timed with the heatwave to maximize unserved
electricity. In this scenario, electricity demand increases by 9% for all
hours and NYISO zones, similarly to the Heatwave scenario, and the
cyberattacker compromises transmission and generation components
simultaneously. Given the uncertainty regarding the cyberattacker
capabilities, we conduct sensitivity analysis on the available cyberattacker
resources and the relative difficulty of compromising transmission components
in this scenario as well.

##  3 Methods

In our framework, the cyberattacker attempts to maximize the disruption to the
electric power network by compromising generation capacity, power flow
capacity, or the limit on nodal voltage phase angle differences. Acharya et
al. [2] show how a cyberattacker can infer features of electricity
transmission and distribution systems using publicly available data. Thus, we
assume that the cyberattacker has complete knowledge of the electric power
network, including the location, capacity, moment of inertia, and type of
controller of generators; the impedance of transformers and transmission
lines; and the load damping constant for New York City. Under complete
knowledge, the cyberattacker can design their strategy by fully anticipating
the response of the NYISO and exacerbate the impact of the cyberattack. Hence,
the complete knowledge assumption is true also under a worst-case scenario and
typical in cyber vulnerability assessments [78, 44, 12]. We represent the
operation of the power system using a Direct Current Optimal Power Flow (DC-
OPF) model, commonly used in power network operations, e.g., for electricity
market dispatch [22], contingency analysis [11], storage siting [76], and
transmission planning [77]. We pose the DC-OPF as an optimization problem,
which regulates generation and power flow decisions within the power system to
minimize the system-wide operating cost. Since the cyberattacker aims to
maximize the electric power network disruption, the DC-OPF model is a helpful
modeling tool for the attacker to determine and analyze its actions.
Therefore, we formulate the cyberattacker problem as a bilevel optimization
problem where the upper level accounts for the objective and constraints of
the cyberattacker and the lower level includes the DC-OPF model (Figure 2).
The bilevel optimization framework allows us to capture the cyberattacker
rationale and the anticipated response of the system operator under the
cyberattack and increased electricity demand due to a Heatwave, which in turns
informs the design of the attack strategy against electric power network
components.

Disruptions in the electric power network decrease electricity supply to all
other sectors. We compute economy-wide losses for all sectors by coupling the
bilevel cyberattacker model with WiNDC, a CGE model [60]. Using WiNDC, we find
the response of all economic sectors to limited electricity supply, caused by
a Heatwave, Cyberattack, and Compound Cyber-Physical Threat. Coupling the
cyberattacker bilevel optimization problem with the CGE allows us to estimate
economy-wide losses, through changes in the Gross Domestic Product, and
identify the top-five impacted sectors in each scenario. In what follows we
detail the formulation of the bilevel problem, the solution algorithm, and the
economy-wide model. Figure 2 details the proposed methodology.

Figure 2: Proposed Framework. A Heatwave stresses the NYISO electric power
network, modeled as a DC-OPF. A cyberattacker chooses strategy $z$ which
maximizes unserved load $\tilde{C}(\cdot)$ while anticipating the response of
the system operator $x$. We input the percent change in regional electricity
supply into a Computable General Equilibrium Model (CGE) to find activity
adjustments of all other economic sectors under a disruption.

###  3.1 Lower-Level Decisions: DC-OPF Model

##### DC-OPF as a Linear Program.

We consider an electric power network with $G$ power generation technologies,
$N$ nodes and $E$ edges/transmission lines. The DC-OPF model computes power
generation $\left(g_{sh}\in\mathbb{R}^{G}\right)$, power flows
$\left(f_{sh}\in\mathbb{R}^{E}\right)$, and voltage angles
$\left(\theta_{sh}\in\mathbb{R}^{N}\right)$ within an electric power network
for all four seasons $s\in\mathcal{S}$ and hours $h\in\mathcal{H}$ of
representative days. Through the DC-OPF, we can also compute wholesale
electricity prices $\left(\pi^{d}_{sh}\in\mathbb{R}^{N}\right)$ and other
payments to power generators, e.g., rent
$\left(\bar{\rho}^{g}_{sh}\in\mathbb{R}^{G}\right)$ when a generator is at
capacity. Therefore, emergency response decisions within a DC-OPF model under
a disruption include change in power generation, change in power flows, change
in voltage angles $\left(\theta_{sh}\right)$, change in electricity prices
($/kWh), and unserved electricity demand
$\left(u_{sh}\in\mathbb{R}^{N}\right)$ when generation and transmission
capabilities can not cover electricity demand at a node in season $s$ and hour
$h$. Equation (1) describes the DC-OPF model, where variables
$\pi^{f},\underline{\rho}^{g},\bar{\rho}^{g},\underline{\rho}^{f},\bar{\rho}^{f},\underline{\rho}^{\theta},\bar{\rho}^{\theta},\underline{\rho}^{u}_{sh},\bar{\rho}^{u}_{sh},\delta$
are the dual variables of the respective constraints and are given in
parentheses. Table 3 in the Supplementary Material provides a detailed
description of all lower-level parameters and variables.

$\begin{array}[]{ll}\min\limits_{\begin{array}[]{l}g_{sh}\in\mathbb{R}^{G},f_{sh}\in\mathbb{R}^{E},\\\
u_{sh},\theta_{sh}\in\mathbb{R}^{N}\end{array}}&C\left(g,u\right)\end{array}$
(1d) s.t. (1e)
$\displaystyle\sum_{t\in\mathcal{T}}g_{sh}+u_{sh}-d_{sh}+A^{T}f_{sh}=0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ $\displaystyle\qquad\left(\pi^{d}_{sh}\right)$ (1f)
$\displaystyle f_{sh}+BA\theta_{sh}=0,$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$
$\displaystyle\qquad\left(\pi^{f}_{sh}\right)$ (1g)
$\displaystyle\underline{g}_{sh}\leq g_{sh}\,\,\,\leq\bar{g}_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$
$\displaystyle\qquad\left(\underline{\rho}^{g}_{sh},\bar{\rho}^{g}_{sh}\right)$
(1h) $\displaystyle\underline{f}_{sh}\leq f_{sh}\,\,\,\leq\bar{f}_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$
$\displaystyle\qquad\left(\underline{\rho}^{f}_{sh},\bar{\rho}^{f}_{sh}\right)$
(1i) $\displaystyle\underline{\theta}_{sh}\leq
A\theta_{sh}\leq\bar{\theta}_{sh},$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$
$\displaystyle\qquad\left(\underline{\rho}^{\theta}_{sh},\bar{\rho}^{\theta}_{sh}\right)$
(1j) $\displaystyle 0\,\,\,\leq u_{sh}\,\,\,\,\leq d_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$
$\displaystyle\qquad\left(\underline{\rho}^{u}_{sh},\bar{\rho}^{u}_{sh}\right)$
(1k) $\displaystyle\theta^{ref}_{sh}=0,$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$
$\displaystyle\qquad\left(\delta_{sh}\right).$ (1l)

The objective function (1e) is the cost observed by the system operator
$C(\cdot)$ for the 24-hour period, while equation (1f) enforces the power
balance constraint in all nodes and equation (1g) describes the relationship
between power flows and voltage angles. Matrix $A\in M_{E\times
N}\left(\mathbb{R}\right)$ in (1f) is the incidence matrix of the electric
power network and matrix $B\in M_{E\times E}\left(\mathbb{R}\right)$ in (1g)
is the diagonal matrix of the susceptances of the electric power network
edges. Equations (1h)-(1k) describe the upper and lower limits on power
generation, power flows, voltage angle difference in all edges, and unserved
demand respectively. Equation (1l) defines the reference node of the power
system. The system operator in our formulation minimizes the following
operating cost of the power system:

$\displaystyle C\left(g,u\right)=$
$\displaystyle\sum_{s\in\mathcal{S},h\in\mathcal{H}}\left(c^{g}_{sh}\right)^{T}g_{sh}+\sum_{s\in\mathcal{S},h\in\mathcal{H}}\left(c^{u}_{sh}\right)^{T}u_{sh}.$
(2)

##### DC-OPF as a Complementarity Problem.

The DC-OPF comprises a convex (linear) objective function and convex feasible
set, defined by a set of linear equality and inequality constraints, therefore
(1) is a convex optimization problem. Hence, a point that satisfies the Karush
Kuhn Tacker (KKT) conditions of (1) is also a minimizer of (1) [10]. Moreover,
we can recast the KKT conditions of a convex optimization problem as a Mixed
Complementarity Problem (MCP) [29]. The MCP reformulation allows us to
represent in a single problem the physical variables $(g,u,f,\theta)$ and
prices and rents $(\pi^{d},\pi^{f},\bar{\rho}^{g},\bar{\rho}^{f})$ of problem
(1). In addition, we will see in Section 3.2 that the MCP reformulation will
facilitate the formulation and solution of the cyberattacker bilevel problem.
The conditions (3) describe the MCP formulation of the KKT conditions of (1).

$\displaystyle 0=$ $\displaystyle\quad
c^{g,l}_{sh}-\underline{\rho}^{g}_{sh}+\bar{\rho}^{g}_{sh}-\pi^{d}_{sh}$
$\displaystyle\quad\bot\quad$ $\displaystyle g_{sh}\quad$
$\displaystyle\text{: free},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (3a) $\displaystyle 0=$
$\displaystyle\quad
A\pi^{d}_{sh}+\pi^{f}_{sh}-\underline{\rho}^{f}_{sh}+\bar{\rho}^{f}_{sh}$
$\displaystyle\quad\bot\quad$ $\displaystyle f_{sh}\quad$
$\displaystyle\text{: free},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (3b) $\displaystyle 0=$
$\displaystyle\quad
A^{T}\left(B\pi^{f}_{sh}-\underline{\rho}^{\theta}_{sh}+\bar{\rho}^{\theta}_{sh}\right)+e_{\text{ref}}\delta_{sh}$
$\displaystyle\quad\bot\quad$ $\displaystyle\theta_{sh}\quad$
$\displaystyle\text{: free},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (3c) $\displaystyle 0=$
$\displaystyle\quad
c^{u}_{sh}-\underline{\rho}^{u}_{sh}+\bar{\rho}^{u}_{sh}-\pi^{d}_{sh}$
$\displaystyle\quad\bot\quad$ $\displaystyle u_{sh}\quad$
$\displaystyle\text{: free},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (3d) $\displaystyle 0\leq$
$\displaystyle\quad g_{sh}-\underline{g}_{sh}$ $\displaystyle\quad\bot\quad$
$\displaystyle\underline{\rho}^{g}_{sh}$ $\displaystyle\geq 0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (3e) $\displaystyle 0\leq$
$\displaystyle\quad\bar{g}_{sh}-g_{sh}$ $\displaystyle\quad\bot\quad$
$\displaystyle\bar{\rho}^{g}_{sh}$ $\displaystyle\geq 0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (3f) $\displaystyle 0\leq$ $\displaystyle\quad
f_{sh}-\underline{f}_{sh}$ $\displaystyle\quad\bot\quad$
$\displaystyle\underline{\rho}^{f}_{sh}$ $\displaystyle\geq 0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (3g) $\displaystyle 0\leq$
$\displaystyle\quad\bar{f}_{sh}-f_{sh}$ $\displaystyle\quad\bot\quad$
$\displaystyle\bar{\rho}^{f}_{sh}$ $\displaystyle\geq 0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (3h) $\displaystyle 0\leq$ $\displaystyle\quad
A\theta_{sh}-\underline{\theta}_{sh}$ $\displaystyle\quad\bot\quad$
$\displaystyle\underline{\rho}^{\theta}_{sh}$ $\displaystyle\geq 0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (3i) $\displaystyle 0\leq$
$\displaystyle\quad\bar{\theta}_{sh}-A\theta_{sh}$
$\displaystyle\quad\bot\quad$ $\displaystyle\bar{\rho}^{\theta}_{sh}$
$\displaystyle\geq 0,$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (3j) $\displaystyle 0\leq$
$\displaystyle\quad u_{sh}$ $\displaystyle\quad\bot\quad$
$\displaystyle\underline{\rho}^{u}_{sh}$ $\displaystyle\geq 0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (3k) $\displaystyle 0\leq$
$\displaystyle\quad\left(d^{w}_{sh}+d^{o}_{sh}\right)-u_{sh}$
$\displaystyle\quad\bot\quad$ $\displaystyle\bar{\rho}^{u}_{sh}$
$\displaystyle\geq 0,$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (3l) $\displaystyle 0=$
$\displaystyle\quad\theta^{\text{ref}}_{sh}$ $\displaystyle\quad\bot\quad$
$\displaystyle\delta_{sh}$ $\displaystyle\text{: free},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}.$ (3m)

Following the literature on complementarity problems [45, 55, 64, 61], we use
“$\bot$” to describe the complementary relationship between the left-hand and
right-hand side conditions in each equation. The notation depicts that an
equilibrium satisfies the equations in both sides of all complementarity
conditions and at least one side is equal to zero in each condition. Following
the MCP notation, we also denote variables without lower bounds as “free”. For
ease of notation, we denote all optimal power flow variables by

$y=\left(g,f,u,\theta,\pi^{d},\pi^{f},\underline{\rho}^{g},\bar{\rho}^{g},\underline{\rho}^{f},\bar{\rho}^{f},\underline{\rho}^{\theta},\bar{\rho}^{\theta},\underline{\rho}^{u}_{sh},\bar{\rho}^{u}_{sh},\delta\right)\in\mathcal{Y}.$

We can write (3) compactly as

$\displaystyle 0=$ $\displaystyle\quad F_{1sh}(y_{1},y_{2})\quad$
$\displaystyle\bot$ $\displaystyle\quad y_{1sh}\quad$ $\displaystyle\text{:
free},$ $\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad
h\in\mathcal{H}$ (4a) $\displaystyle 0\leq$ $\displaystyle\quad
F_{2sh}(y_{1},y_{2})\quad$ $\displaystyle\bot$ $\displaystyle\quad
y_{2sh}\quad$ $\displaystyle\geq 0,$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad h\in\mathcal{H},$ (4b)

where

$y=\left(y_{1},y_{2}\right)\in\mathcal{Y}_{1}\times\mathcal{Y}_{2}=\prod_{s\in\mathcal{S},h\in\mathcal{H}}\mathcal{Y}_{1sh}\times\prod_{s\in\mathcal{S},h\in\mathcal{H}}\mathcal{Y}_{2sh},$

$y_{1sh}=\left(g_{sh},u_{sh},f_{sh},\theta_{sh},\delta_{sh}\right)\in\mathcal{Y}_{1sh}=\prod_{i=1}^{5}\mathcal{Y}_{1sh}^{(i)},$
$y_{2sh}=\left(\underline{\rho}^{g}_{sh},\bar{\rho}^{g}_{sh},\underline{\rho}^{f}_{sh},\bar{\rho}^{f}_{sh},\underline{\rho}^{\theta}_{sh},\bar{\rho}^{\theta}_{sh},\underline{\rho}^{u}_{sh},\bar{\rho}^{u}_{sh}\right)\in\mathcal{Y}_{2sh}=\prod_{i=1}^{8}\mathcal{Y}_{2sh}^{(i)}.$

Notice that $F_{1},F_{2}$ are linear functions. An equilibrium of (3) is a
point $y^{*}=\left(y_{1}^{*},y_{2}^{*}\right)$ that satisfies all
complementarity conditions (3a)-(3m). Since the original problem is convex,
then a point $\left(g^{*},f^{*},u^{*},\theta^{*}\right)$ that satisfies the
KKT conditions (4) is also a minimizer of (1) [10].

###  3.2 Upper-Level Decisions: Model of the Cyberattacker

The Cyberattack scenario designed in Section 2 assumes that the cyberattacker
maximizes the operational disruption in electric power infrastructure by
targeting SCADA devices and PLC interfaces. The cyberattacker strategy
includes compromising power generation capacity
$\left(z^{g}_{sh}\in\mathbb{R}^{G}\right)$, i.e., the upper limit,
transmission limits $(z^{f}_{sh}\in\mathbb{R}^{E})$, and limits on voltage
angle differences $\left(z^{\theta}_{sh}\in\mathbb{R}^{E}\right)$. Hence, the
physical power generation capacity, transmission, and limits on voltage angle
differences under a cyberattack are modeled as

$\displaystyle\underline{g}_{sh}\quad$ $\displaystyle\leq$ $\displaystyle\quad
g_{sh}$ $\displaystyle\leq$ $\displaystyle\quad\bar{g}_{sh}-z^{g}_{sh},\quad$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (5a) $\displaystyle\underline{f}_{sh}-z^{f}_{sh}\quad$
$\displaystyle\leq$ $\displaystyle\quad f_{sh}$ $\displaystyle\leq$
$\displaystyle\quad\bar{f}_{sh}-z^{f}_{sh},\quad$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (5b)
$\displaystyle\underline{\theta}_{sh}-z^{\theta}_{sh}\quad$
$\displaystyle\leq$ $\displaystyle\quad A\theta_{sh}\quad$ $\displaystyle\leq$
$\displaystyle\quad\bar{\theta}_{sh}-z^{\theta}_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}.$ (5c)

Given the ability of the cyberattacker to compromise components of electric
power infrastructure as modeled in (5a)-(5c), the bilevel diruption-maximizing
problem of the cyberattacker can be formulated as

$\begin{array}[]{l}\max\limits_{\begin{array}[]{l}z^{g}_{sh}\in\mathbb{R}^{G},z^{f}_{sh}\in\mathbb{R}^{E},\\\
z^{\theta}_{sh}\in\mathbb{R}^{E},y\in\mathcal{Y}\end{array}}\tilde{C}\left(y,z^{g},z^{f},z^{\theta}\right)\end{array}$
(6d) s.t. (6e) $\displaystyle
y\in\arg\min\limits_{y}\left\\{C(y;z^{g},z^{f},z^{\theta})\,\,|\,\,\eqref{eq:opf-
lp-mkt-clrng-cnstr}-\eqref{eq:opf-lp-pow-flow-cnstr},\eqref{eq:attacker-
compromised-capacity-gen}-\eqref{eq:attacker-compromised-capacity-
angles},\eqref{eq:opf-lp-unsd-lim-cnstrs}-\eqref{eq:opf-lp-ref-node-
cnstr}\right\\}$ (6f) $\displaystyle 0\leq z^{g}_{sh}\leq\bar{g}_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (6g) $\displaystyle 0\leq z^{f}_{sh}\leq\bar{f}_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (6h) $\displaystyle 0\leq
z^{\theta}_{sh}\leq\bar{\theta}_{sh},$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (6i)
$\displaystyle\sum_{s\in\mathcal{S}}\left(\tilde{c}^{g}_{sh}\right)^{T}z^{g}_{sh}+\sum_{s\in\mathcal{S}}\left(\tilde{c}^{f}_{sh}\right)^{T}z^{f}_{sh}+\sum_{s\in\mathcal{S}}\left(\tilde{c}^{\theta}_{sh}\right)^{T}z^{\theta}_{sh}\leq\tilde{b}_{s},$
$\displaystyle\quad\forall s\in\mathcal{S}.$ (6j)

The objective function of the cyberattacker in (6e) is a measure of the target
disruption $\tilde{C(\cdot)}$, while (6f) enforces that all lower-level
variables are solutions to the DC-OPF. Equations (6g)-(6i) capture that
compromised generation and transmission capacity, and limits on voltage angle
differences can not exceed the respective physical capacities and limits.
However, a cyberattacker can target intelligent electronic devices of
substations, e.g., transmission line protective relays [38], and generators,
e.g., manipulating the automatic generation control signals [3], to constrain
power flows and generation. Equation (6j) describes the capabilities of the
cyberattacker in each hour $h$ of each season $s$. Compromising generation and
transmission capacity, and voltage angle limits requires committing resources
relevant to the vulnerability of each component, e.g., work-hours, software
and hardware infrastructure, or cost of employing specialized personnel, with
corresponding costs
$\tilde{c}^{g}_{sh},\tilde{c}^{f}_{sh},\tilde{c}^{\theta}_{sh}$ with respect
to the total available resources $\tilde{b}_{s}>0$ in season $s$. In Section 2
we assume that in the default Cyberattack scenario, the cost of attacking
transmission $\left(\tilde{c}^{f}_{sh},\tilde{c}^{\theta}_{sh}\right)$ is
significantly higher than generation $\left(\tilde{c}^{g}_{sh}\right)$. Given
the difficulty of quantifying component vulnerabilities, we assign
$\tilde{c}^{f}_{she}=5c^{g}_{sht}$ for all $e\in\mathcal{E}$,
$t\in\mathcal{T}$ in the default Cyberattack scenario and conduct sensitivity
analysis in Section 4. In our formulation, the cyberattacker targets unserved
electricity demand, i.e.,

$\tilde{C}\left(y,z^{g},z^{f},z^{\theta}\right)=\sum_{s\in\mathcal{S},\in\mathcal{H}}\left(c^{u}_{sh}\right)^{T}u_{sh}.$

Notice that (6) is an optimization problem where (6f) yields an equilibrium of
the electricity market (1). Therefore, (6) is a Mathematical Program with
Equilibrium Constraints [45]. Finding a global optimum of an MPEC is a
combinatorial problem [45], thus computationally demanding. In the next
section we describe how to efficiently solve large-scale instances of (6).

###  3.3 Bilevel Optimization of the Cyberattacker: Solution Algorithm

In Section 3.1 we showed that finding a minimizer of the DC-OPF formulation
(1) is equivalent to finding a solution $y^{*}$ of the MCP (4). Therefore, we
can substitute (6f) with the respective complementarity conditions (4) and
problem (6) becomes a Mathematical Program with Complementarity Constraints
(MPCC), a special case of an MPEC [55], which allows us to model in a single
framework the physical and technological (primal), and economic (dual)
variables of the lower-level problem (6f). Problem (7) is the MPCC
reformulation of the bilevel optimizaiton of the cyberattacker in (6).

$\displaystyle\text{MPCC}:=$
$\begin{array}[]{ll}\max\limits_{\begin{array}[]{l}z^{g}_{sh}\in\mathbb{R}^{G},z^{f}_{sh}\in\mathbb{R}^{E},\\\
z^{\theta}_{sh}\in\mathbb{R}^{E},y_{1}\in\mathcal{Y}_{1},y_{2}\in\mathcal{Y}_{2}\end{array}}\sum\limits_{s\in\mathcal{S},\in\mathcal{H}}\left(c^{u}_{sh}\right)^{T}u_{sh}\end{array}$
(7d) s.t. (7e) $\displaystyle F_{1sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})=0$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (7f) $\displaystyle 0\leq y_{2sh}\quad\bot\quad
F_{2sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})\geq 0$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (7g)
$\displaystyle 0\leq z^{g}_{sh}\leq\bar{g}_{sh},$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (7h)
$\displaystyle 0\leq z^{f}_{sh}\leq\bar{f}_{sh},$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (7i)
$\displaystyle 0\leq z^{\theta}_{sh}\leq\bar{\theta}_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (7j)
$\displaystyle\sum_{s\in\mathcal{S}}\left(\tilde{c}^{g}_{sh}\right)^{T}z^{g}_{sh}+\sum_{s\in\mathcal{S}}\left(\tilde{c}^{f}_{sh}\right)^{T}z^{f}_{sh}+\sum_{s\in\mathcal{S}}\left(\tilde{c}^{\theta}_{sh}\right)^{T}z^{\theta}_{sh}\leq\tilde{b}_{s},$
$\displaystyle\quad\forall s\in\mathcal{S}$ (7k)

Although all functions in (7) are linear, the MPCC is a nonlinear optimization
problem since the complementarity conditions (7g) imply that either $y_{2,i}$
or $F_{2,i}(y_{1},y_{2})$ are zero for all $i\in|\mathcal{Y}_{2}|$. Finding a
global optimum of a nonlinear optimization problem can be computationally
challenging, and many commercial solvers compute a local optimum [62]. To
compute a global optimum of (7), we can recast the complementarity conditions
using binary variables, and reformulate the MPCC into a Mixed Integer Linear
Program (MILP) [45]. The MILP reformulation requires introducing as many
binary variables as complementarity constraints in (3), which in our problem
are in the order of tens of thousands111The complementarity conditions include
roughly $9\text{ (generation technologies)}\times 16\text{ (nodes)}\times
24\text{ (hours)}\times 3\text{ (variables
}g_{sh},\underline{\rho}^{g}_{sh},\bar{\rho}^{g}_{sh})\approx 10,000$
generation-related conditions and variables in (3a),(3e)-(3f); $17\text{
(edges)}\times 24\text{ (hours)}\times 3\text{ (variables
}f_{sh},\underline{\rho}^{f}_{sh},\bar{\rho}^{f}_{sh})\approx 1,000$ flow-
related conditions and variables in (3b),(3g)-(3h); and $16\text{
(nodes)}\times 24\text{ (hours)}\times 6\text{ (variables
}\theta_{sh},u_{sh},\underline{\rho}^{\theta}_{sh},\bar{\rho}^{\theta}_{sh},\underline{\rho}^{u}_{sh},\bar{\rho}^{u}_{sh},\delta_{sh})\approx
2,000$ node-related conditions and variables in (3c)-(3d), (3i)-(3m).. For the
MILP reformulation, we introduce the following binary variables and
parameters222Note that
$y_{2sh}=\left(\underline{\rho}^{g}_{sh},\bar{\rho}^{g}_{sh},\underline{\rho}^{f}_{sh},\bar{\rho}^{f}_{sh},\underline{\rho}^{\theta}_{sh},\bar{\rho}^{\theta}_{sh},\underline{\rho}^{u}_{sh},\bar{\rho}^{u}_{sh}\right)\in\prod_{i=1}^{8}\mathbb{R}^{\kappa_{i}}_{\geq
0}$.

$\gamma_{sh}\in\prod_{i=1}^{8}\left\\{0,1\right\\}^{\kappa_{i}},\quad
M_{sh}\in\prod_{i=1}^{8}\mathbb{R}_{\geq
0}^{\kappa_{i}},\quad\kappa_{i}=\left\lvert\mathcal{Y}^{(i)}_{2sh}\right\rvert,\qquad\forall
s\in\mathcal{S},\,\,h\in\mathcal{H}.$

The MILP reformulation of the MPCC in (7) is

MPCC-MILP:=
$\begin{array}[]{ll}\max\limits_{\begin{array}[]{l}z^{g}_{sh}\in\mathbb{R}^{G},z^{f}_{sh}\in\mathbb{R}^{E},\\\
z^{\theta}_{sh}\in\mathbb{R}^{E},y_{1}\in\mathcal{Y}_{1},y_{2}\in\mathcal{Y}_{2}\end{array}}\sum\limits_{s\in\mathcal{S},h\in\mathcal{H}}\left(c^{u}_{sh}\right)^{T}u_{sh}\end{array}$
(8d) s.t. (8e) $\displaystyle F_{1sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})=0,$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (8f) $\displaystyle 0\leq
F_{2sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})\leq\left(\mathbf{1}_{sh}-\gamma_{sh}\right)M_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (8g) $\displaystyle 0\leq y_{2sh}\leq\gamma_{sh}\cdot
M_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (8h) $\displaystyle 0\leq
z^{g}_{sh}\,\,\leq\bar{g}_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (8i) $\displaystyle 0\leq
z^{f}_{sh}\,\,\leq\bar{f}_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S}$
$\displaystyle\quad\forall h\in\mathcal{H}$ (8j) $\displaystyle 0\leq
z^{\theta}_{sh}\,\,\leq\bar{\theta}_{sh},$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (8k)
$\displaystyle\sum_{h\in\mathcal{H}}\left(\tilde{c}^{g}_{sh}\right)^{T}z^{g}_{sh}+\sum_{h\in\mathcal{H}}\left(\tilde{c}^{f}_{sh}\right)^{T}z^{f}_{sh}+\sum_{h\in\mathcal{H}}\left(\tilde{c}^{\theta}_{sh}\right)^{T}z^{\theta}_{sh}\leq\tilde{b}_{s},$
$\displaystyle\quad\forall s\in\mathcal{S}$ (8l)

where $\mathbf{1}_{sh}\in\mathcal{Y}_{2sh}$ is a vector of all ones and
“$\cdot$” in equations (8g)-(8h) is the element-wise product between the two
vectors of each bilinear term. Equations (8h), (8g) are equivalent to
(7g)333When $\gamma_{sh,i}=1$ then $F_{2sh,i}=0$ and $0\leq y_{2sh,i}\leq
M_{sh,i}$. Similarly, when $\gamma_{sh,i}=0$ then $y_{2sh,i}=0$ and $0\leq
F_{2sh,i}\leq M_{sh,i}$. Thus, equations (8h), (8g) enforce that
$y_{2sh,i},F_{2sh,i}(y_{1},y_{2})$ are non-negative and at least one is zero..
If $\left(z^{g*},z^{f*},z^{\theta*},y_{1}^{*},y_{2}^{*}\right)$ is a global
optimum of (7) and parameter
$M_{i}>\lvert\lvert\left(z^{g*},z^{f*},z^{\theta*},y_{1}^{*},y_{2}^{*}\right)\rvert\rvert_{\infty}$,
then a solution of the MILP in (8) is a global optimum of the original MPCC
(7) [27].

Finding a global optimum of the MPCC (7) is a combinatorial problem [45] and
requires significant computational power. In our formulation, the number of
binary variables of the MPCC-MILP (8) is in the order of tens of thousands,
which exacerbates the problem of available computational resources and long
computation time. The Nonlinear (NLP) reformulation of a MPCC addresses the
problem of computational power and is the most widely used alternative to the
MILP reformulation for large-scale MPCCs [42, 43]. Under the NLP
reformulation, we cast the complementarity conditions (7g) as equality and
inequality constraints (9g)-(9i) to derive the Nonlinear Programming (NLP)
reformulation of the MPCC (7), i.e.,

$\displaystyle\text{MPCC-NLP}:=$
$\begin{array}[]{ll}\max\limits_{\begin{array}[]{l}z^{g}_{sh}\in\mathbb{R}^{G},z^{f}_{sh}\in\mathbb{R}^{E},\\\
z^{\theta}_{sh}\in\mathbb{R}^{E},y_{1}\in\mathcal{Y}_{1},y_{2}\in\mathcal{Y}_{2}\end{array}}\sum\limits_{s\in\mathcal{S},\in\mathcal{H}}\left(c^{u}_{sh}\right)^{T}u_{sh}\end{array}$
(9d) s.t. (9e) $\displaystyle F_{1sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})=0$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (9f) $\displaystyle
F_{2sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})\cdot y_{2sh}=0$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (9g) $\displaystyle 0\leq
F_{2sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (9h)
$\displaystyle 0\leq y_{2sh}$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (9i) $\displaystyle 0\leq
z^{g}_{sh}\leq\bar{g}_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (9j) $\displaystyle 0\leq
z^{f}_{sh}\leq\bar{f}_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (9k) $\displaystyle 0\leq
z^{\theta}_{sh}\leq\bar{\theta}_{sh},$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (9l)
$\displaystyle\sum_{h\in\mathcal{H}}\left(\tilde{c}^{g}_{sh}\right)^{T}z^{g}_{sh}+\sum_{h\in\mathcal{H}}\left(\tilde{c}^{f}_{sh}\right)^{T}z^{f}_{sh}+\sum_{h\in\mathcal{H}}\left(\tilde{c}^{\theta}_{sh}\right)^{T}z^{\theta}_{sh}\leq\tilde{b}_{s},$
$\displaystyle\quad\forall s\in\mathcal{S}$ (9m)

Using a NLP solver to solve the MPCC-NLP is computationally more efficient for
many MPCC problems [64, 27]. However, the bilinear term (9g) implies that
conventional NLP solvers can converge to a stationary point or local optimum
of the MPCC-NLP (9) [43]. Moreover, MPCCs, like MPECs, can have disjoint
feasible sets [32], hence the solver optimum depends on the starting point. In
our MPCC, there exists a tradeoff between finding a global optimum using the
MILP reformulation and computational efficiency using the NLP reformulation.
The first is more computationally demanding while the second may converge to a
local optimum. To overcome these challenges, we will use a heuristic to
approximate the solution of the MILP. Consequently, we will initiate the MPCC-
NLP (9) at a point informed by the structure of our problem. Specifically, for
each season $s$ and hour $h$ in our formulation, functions $F_{1sh},F_{2sh}$
depend only on $y_{1sh},y_{2sh}$. Hence, only constraint (9m) couples
decisions across time periods, while all remaining constraints (9f)-(9l) and
all variables are decoupled in $s,h$. Moreover, the objective function is
additive. Therefore, we can approximate the MPCC-MILP reformulation in (8) by
maximizing a series of hourly problems $(c^{u}_{sh})^{T}u_{sh}$ for each $s$
and $h$ subject to an hourly budget

$\tilde{b}_{sh}=\frac{\tilde{b}_{s}}{H},\qquad\forall s\in\mathcal{S},\qquad
h\in\mathcal{H}.$

The decoupled MPCC-MILP(s,h) for each season $s$ and hour $h$ is

DMPCC-MILP(s,h):=
$\begin{array}[]{ll}\max\limits_{\begin{array}[]{l}z^{g}_{sh}\in\mathbb{R}^{G},z^{f}_{sh}\in\mathbb{R}^{E},\\\
z^{\theta}_{sh}\in\mathbb{R}^{E},y_{1}\in\mathcal{Y}_{1},y_{2}\in\mathcal{Y}_{2}\end{array}}\sum\limits_{s\in\mathcal{S},h\in\mathcal{H}}\left(c^{u}_{sh}\right)^{T}u_{sh}\end{array}$
(10d) s.t. (10e) $\displaystyle
F_{1sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})=0,$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (10f)
$\displaystyle 0\leq
F_{2sh}(y_{1},y_{2};z^{g},z^{f},z^{\theta})\leq\left(\mathbf{1}_{sh}-\gamma_{sh}\right)M_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (10g) $\displaystyle 0\leq y_{2sh}\leq\gamma_{sh}\cdot
M_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (10h) $\displaystyle 0\leq
z^{g}_{sh}\,\,\leq\bar{g}_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S},$
$\displaystyle\quad\forall h\in\mathcal{H}$ (10i) $\displaystyle 0\leq
z^{f}_{sh}\,\,\leq\bar{f}_{sh},$ $\displaystyle\quad\forall s\in\mathcal{S}$
$\displaystyle\quad\forall h\in\mathcal{H}$ (10j) $\displaystyle 0\leq
z^{\theta}_{sh}\,\,\leq\bar{\theta}_{sh},$ $\displaystyle\quad\forall
s\in\mathcal{S},$ $\displaystyle\quad\forall h\in\mathcal{H}$ (10k)
$\displaystyle\left(\tilde{c}^{g}_{sh}\right)^{T}z^{g}_{sh}+\left(\tilde{c}^{f}_{sh}\right)^{T}z^{f}_{sh}+\left(\tilde{c}^{\theta}_{sh}\right)^{T}z^{\theta}_{sh}\leq\tilde{b}_{sh},$
$\displaystyle\quad\forall s\in\mathcal{S},$ $\displaystyle\quad\forall
h\in\mathcal{H}$ (10l)

The solution
$\left(z^{g(1)}_{sh},z^{f(1)}_{sh},z^{\theta(1)}_{sh},y^{(1)}_{1sh},y^{(1)}_{2sh}\right)$
of the DMPCC-MILP (10) may not optimize the MPCC-MILP (8), since the decoupled
formulation assumes the attacker’s budget allocation to be homogeneous across
time periods. We then initiate the MPCC-NLP with
$\left(z^{g(1)}_{sh},z^{f(1)}_{sh},z^{\theta(1)}_{sh},y^{(1)}_{1sh},y^{(1)}_{2sh}\right)$
to compute the cyberattacker’s actions under a more disruptive budget
allocation. Table 1 describes the steps of the proposed solution algorithm.

Algorithm. |
---|---
Begin |  |  |
|  For $s\in\mathcal{S}$ |  |
|  | For | $h\in\mathcal{H}$
|  |  |  Find $\left(z^{g(1)}_{sh},z^{f(1)}_{sh},z^{\theta(1)}_{sh},y^{(1)}_{1sh},y^{(1)}_{2sh}\right)$ solution of:
|  |  | DMPCC-MILP(s,h):=$\max\limits_{\begin{array}[]{l}z^{g}_{sh},z^{f}_{sh},z^{\theta}_{sh},\\\ y_{1sh},y_{2sh}\end{array}}\left(c^{u}_{sh}\right)^{T}u_{sh}$
|  |  |  s.t. (8g)-(8k)
|  |  |  $\left(\tilde{c}^{g}_{sh}\right)^{T}z^{g}_{sh}+\left(\tilde{c}^{f}_{sh}\right)^{T}z^{f}_{sh}+\left(\tilde{c}^{\theta}_{sh}\right)^{T}z^{\theta}_{sh}\leq\tilde{b}_{sh}$
|  | End |
|  End |  |
|  Initialize MPCC-NLP with the decoupled DMPCC-MILP(s,h) solutions:
|  |  $\left(z^{g},z^{f},z^{\theta},y_{1},y_{2}\right)^{(0)}=\prod\limits_{s\in\mathcal{S},h\in\mathcal{H}}\left(z^{g(1)}_{sh},z^{f(1)}_{sh},z^{\theta(1)}_{sh},y_{1}^{(1)},y_{2}^{(1)}\right)$.
|  Find $\left(z^{g*},z^{f*},z^{\theta*},y^{*}_{1},y^{*}_{2}\right)$ solution
of:
|  | MPCC-NLP:=$\max\limits_{\begin{array}[]{l}z^{g}_{sh},z^{f}_{sh},z^{\theta}_{sh},\\\ y_{1sh},y_{2sh}\end{array}}\left(c^{u}_{sh}\right)^{T}u_{sh}$
|  |  |  s.t. (9f)-(9l).
|  |  |  $\sum\limits_{h\in\mathcal{H}}\left(\tilde{c}^{g}_{sh}\right)^{T}z^{g}_{sh}+\sum\limits_{h\in\mathcal{H}}\left(\tilde{c}^{f}_{sh}\right)^{T}z^{f}_{sh}+\sum\limits_{h\in\mathcal{H}}\left(\tilde{c}^{\theta}_{sh}\right)^{T}z^{\theta}_{sh}\leq\tilde{b}_{s}$
End |  |  |
Table 1:  Decomposition algorithm of MILP reformulation (8) of MPCC (7). Since
the objective function of the original problem is additive and all constraints
are disjoint in $s,h$ for fixed budget allocation, then we can retrieve an
approximate solution to the MPCC-MILP under a homogeneous budget allocation by
solving the DMPCC-MILP(s,h) for all $s\in\mathcal{S},h\in\mathcal{H}$. We then
feed the DMPCC-MILP solution as a starting point to the MPCC-NLP and allow for
flexible budget allocation between time periods to compute an optimal
disruption also across time periods.

###  3.4 Economy-Wide Model: Computable General Equilibrium

We evaluate the ripple effects of disruptions in electric power networks
across all economic sectors using the WiNDC model, developed and maintained by
the University of Wisconsin, Madison [60]. The economy-wide cost of blackouts
can be hard to derive. For example, in the 2003 Northeastern Blackout, the
U.S. Department of Energy estimated losses between $4 – $10 Billion in the
U.S. [68] following the event. However, finding the Value of Lost Load (VOLL)
under a power outage can be challenging and estimates vary. For example, VOLL
is different between industrial, commercial and residential consumers and the
estimation can be based on econometric modeling under multiple assumptions
[66]. System operators employ a variety of methods to compute VOLL, including
an ex-post macroeconomic analysis of the contribution of an event to the
reduction of the Gross Domestic Product [25]. Despite the shortcomings of
computing the VOLL using ex-post macroeconomic analysis, the Electric
Reliability Council of Texas, Inc. reports that the results can serve as a
“reference point”. Instead, we aim to provide an ex-ante estimation of how
transactions between sectors deviate from the Baseline using a Computable
General Equilibrium (CGE) model to improve the assessment of outages’ riple
effects on the economy.

The WiNDC model is a multi-regional, multi-sectoral CGE Model comprising 51
regions, including the 50 States and the District of Columbia denoted with
their two-letter code, and 71 economic sectors. In New York State, outside the
“Utilities” sector, the “Residential” sector consumes the most output of the
“Utilities” sector in New York State, followed by the “Construction”, “State
and Local Government Enterprises (SLGE)”, and “Pipeline Transportation”
sectors. The output of each sector requires intermediate inputs from other
sectors, labor and capital, and can cover regional demand and the demand for
exports to other regions and the rest of the world. Domestic prices drive the
optimal output and demand for production factors of each sector, while changes
in regional prices drive inter-regional and international trade in each
sector. Household consumption in WiNDC is driven by consumer preferences and
regional price changes allowing households to substitute between goods in
response to price changes. The model also allows for taxes and subsidies on
intermediate demand, and import tariffs. Since 2019, the WiNDC model has been
deployed to study economy-wide impacts to U.S. states of international food
trade policies, climate policy design, and the impact of environmental
regulation on regional labor markets. For climate policy design, Rutherford et
al., [60] expand the core WiNDC module to account for energy intensity of
sectoral supply chains to understand carbon leakage due to energy consumption
at the subnational level, across U.S. states, while Becker et al., [6] further
enhance WiNDC to understand the impact of climate policies on state-level
labor supply across demographic groups. In the food sector, Ballistreri et
al., [4], employ WiNDC to understand adverse impacts of The Phase One Trade
agreement between the United States and China across U.S. states [71].

Figure 2 describes the coupling of the bilevel optimization problem with the
CGE in our framework. First, the bilevel optimization problem assesses the
availability of electricity supply of the New York State electric power
network under a Heatwave, Cyberattack, or Compound Cyber-Physical Threat.
Then, WiNDC receives regional electricity supply as input and assesses the
response of all other sectors and potential economic losses to the New York
State economy. The implementation of the scenarios in Section 2 in WiNDC
required diverging from the core model in [60]. We adjust the base year of
WiNDC from 2014 to 2018 to ensure that trade across sectors in the CGE are
consistent with the operating year of the NYISO electric power network. For
that, we introduce an exogenous total factor productivity, $tfp_{r}$ for each
region $r$, which is set to the Gross Domestic Product (GDP) growth rate of
each state [5] in all quantities, e.g., $X_{r,g}:=tfp_{r}\cdot X_{r,g}$ in
equations (52)-(54) in [60] which become

$\displaystyle\sum_{s}Y_{r,s}\bar{ys}_{r,s,g}+\bar{yh}_{r,g}-tfp_{r}X_{r,g}\bar{s}_{r,g}=0$
(11) $\displaystyle
tfp_{r}X_{r,g}S^{D}_{r,g}-A_{r,g}D^{D}_{r,g}-\sum_{m}MS_{r,m}\bar{dm}_{r,g,m}=0$
(12)
$\displaystyle\sum_{r}\left(tfp_{r}X_{r,g}S^{N}_{r,g}-A_{r,g}D^{N}_{r,g}-\sum_{m}MS_{r,m}\bar{nm}_{r,g,m}\right)=0$
(13)

In WiNDC, set indices $r,s,g,m$ denote the WiNDC regions, sectors, goods
produced by sectors, and margin adjustments for trade and transport. Equation
(11) regulates state-level supply, where $\bar{yh}_{r,g}$ denotes fixed
household production, $\bar{ys}_{r,s,g}$ denotes baseline sectoral production
in region $r$, $Y_{r,s}$ is the change of sectoral production compared to
baseline, $\bar{s}_{r,g}$ is baseline regional supply, and $X_{r,g}$ is the
change of regional supply compared to baseline. Equation (12) regulates
allocation of supply to the state-level market, where $\bar{dm}_{g,m}$ denotes
baseline state-level margin supply, $MS_{r,m}$ denotes change of state-level
margin supply compared to baseline, $D^{D}_{r,g}$ denotes state-level goods
demand, $A_{r,g}$ is the regional allocation of goods, and $S^{D}_{r,g}$
denotes state-level supply. Equation (13) ensures the balance of quantities at
the national scale, where $S^{N}_{r,g}$ denotes national supply, $D^{N}_{r,g}$
denotes national goods demand, and $\bar{nm}_{r,g,m}$ denotes the baseline
national supply of margin sectors.

The New York State annual GDP in our model grows from $1.47 Trillion in 2014
to $1.57 Trillion in 2018, which translates to a daily GDP of $4.31 Billion.
The Compound Cyber-Physical Threat causes a shortage of regional electricity
supply, which leads to unserved demand. We implement the Compound Cyber-
Physical Threat scenarios in WiNDC by decreasing regional electricity supply
provided by the “Utilities” sector, similarly to the simulation of power
outages impacts using a CGE in Rose et al., [58]. The “Utilities” sector in
this version of WiNDC comprises electricity, natural gas, and water/wastewater
services. For that, we assume that a supply disruption in any of the three
subsectors translates to an equivalent supply disruption to the services of
the “Utilities” sector as a whole. The percentage decrease of regional
electricity supply (variable $X_{r,g}$ for state $r$ and sector $g$) is equal
to the percent of unserved load with respect to Baseline load in the NYISO
model. Figure 2 shows the coupling between the bilevel optimization problem
and the CGE which allows us to endogenously consider the operational
disruptions within the electric power network under each scenario when
assessing the response of all other economic sectors.

The computation of CGE parameters assumes yearly or 5-year steps. Therefore, a
CGE simulates the long-term substitution between intermediate inputs and
production factors for economic sectors, and between goods for households.
Calculating robust substitution elasticities under extreme events requires
often unavailable time-series data on the response of different sectors during
the extreme event. For that, substitution elasticities under infrastructure
disruptions are unavailable for most sectors, e.g., the food sector [48].
Nevertheless, conducting economy-wide analysis of a 24-hour disruption using a
CGE can provide a lower bound on ripple effects and economic losses. Since
sectors and consumers can substitute inputs and consumed goods easier in a CGE
framework, we expect that the economic losses will be greater under more rigid
substitution regimes and New York State would face at least losses reported
from the CGE.

##  4 Results

###  4.1 Vulnerability of Regional Infrastructure Components

In this section, we identify vulnerable components of the NYISO grid and
compare the magnitude and regional distribution of unserved load under the
Heatwave, Cyberattack, and Compound Cyber-Physical Threat scenarios as
described in Section 2. Figure 3 shows that the NYISO grid has the generation
and transmission capacity to withstand the Heatwave without shedding load,
while the Cyberattack under normal operating conditions results in 408 MW of
unserved load in a 24-hour period, which is 0.1% of total NYISO load or
approximately 19,800 customers444We estimate number of customers by
multiplying the percentage of unserved load with the total number of NYISO
customers, as reported in the Baseline scenario. However, timing the
Cyberattack with the Heatwave can compound the impact of the Cyberattack by
more than 3 times. Therefore, the Compound Cyber-Physical Threat results in
2401 MWh of unserved load in a 24-hour period, which is almost 0.4% of total
NYISO load and affects almost 79,200 customers. However, the impact varies
temporally, with the largest shedding at 17:00, that is during the peak load
hour. For comparison, under the Cyberattack, peak load shedding reaches 253
MW, while under the Compound Cyber-Physical Threat, peak load shedding reaches
722 MW.

Timing the Cyberattack excacerbates also the duration of load shedding. In the
Cyberattack scenario, load shedding happens for five hours, between
13:00-18:00, while in the Compound Cyber-Physical Threat scenario, load
shedding happens for 11 hours, between 10:00-21:00. Therefore, the Compound
Cyber-Physical Threat exacerbates both the impact and the duration of the
Cyberattack.

Figure 3: NYISO Unserved Load (MW). Under the Heatwave there is no unserved
load. During peak load at 17:00 and under normal operating conditions, the
Cyberattack compromises 253 MW (1%) of NYISO load. For the same time period,
the Compound Cyber-Physical Threat causes load curtailment of 722MW, which is
three times higher compared to the Cyberattack. Figure 4: Regional NYISO
Unserved Load (MW) at 17:00. Under the Heatwave there is no unserved load. The
Cyberattack under normal operating conditions compromises 4% of Long Island
load and 1% of NYISO load. The Compound Cyber-Physical Threat increases load
curtailment in Long Island to 13%, while NYISO curtailment is 2%.

Figure 4 shows the regional distribution of unserved load in each scenario
across all NYISO zones. We observe that the cyberattacker targets Long Island
power generation in the Cyberattack and Compound Cyber-Physical Threat
scenarios between 17:00-18:00PM, when the NYISO power network experiences the
highest load shedding in both scenarios. The reason is that Long Island power
generation is close to peak generation capacity at 18:00 on a summer day.
Moreover, Long Island connects with only one other NYISO zone and the
interconnection is close to capacity under the Baseline scenario. Therefore,
compromising power generation capacity has immediate impact on local unserved
load.

Figure 5: Cyberattacker strategy. Under the scenarios assumptions, the
cyberattacker compromises consistently more generation capacity than
transmission capacity. The cyberattacker allocates all of their budget to
compromise 3% of the total NYISO generation capacity during peak load
shedding.

Between compromising transmission and generation capacity in Long Island,
which are operated close to their limits in the Baseline scenario, the
cyberattacker chooses to compromise more generation capacity. Figure 5 shows
the strategy of the cyberattacker. Between 13:00-18:00 in the Cyberattack
scenario, and 10:00-21:00 in the Compound Cyber-Physical Threat scenario, when
load curtailment happens in both scenarios, the attacker allocates all their
resources to compromising generation capacity. Figure 5 reveals that the
available cyberattacker resources allow them to compromise 3% of total NYISO
generation capacity during load curtailment periods. The optimal cyberattacker
strategy arises from the distribution of cyberattacker costs across
components, where we have assumed that compromising generation capacity is
relatively cheaper compared to transmission capacity. However, our knowledge
on the cyberattacker costs distribution and budget availability is limited and
may change over time. For that, in the next section we test the sensitivity of
our results to varying cyberattacker costs and budget availability.

###  4.2 Infrastructure Vulnerabilities Sensitivity Analysis

####  4.2.1 Cyberattacker cost distribution

In Section 2, we describe how transmission components are often more protected
against cyberattacks compared to consumer and generation components. Hence, in
our formulation, we have assumed that the cost of compromising transmission
capacity is disproportionately greater than the cost of compromising
generation capacity for the same budget. In this Section, we explore the
operational impact of a cyberattack if the difficulty of compromising
transmission lines decreases in the future. For example, cybersecurity
protocols may not keep up with the enhanced skillset of cyberattackers. For
that, a cyberattacker may be able to exploit novel attack vectors which can
provide easier access to transmission components. In this section, we conduct
sensitivity analysis on the relative cost
$\left(\frac{\tilde{c}^{g}}{\tilde{c}^{f}}\right)$ compared to the default
relative cyberattacker cost
$\left(\frac{\tilde{c}^{g,0}}{\tilde{c}^{f,0}}\right)$, by decreasing the
cyberattacker transmission cost by 10% and increasing the cyberattacker
generation cost by 10% compared to their default values in each iteration,
i.e., for iteration $i\in\mathcal{I}=\left\\{1,2,\dots,6\right\\}$, the
relative cyberattacker cost $(\tilde{c}^{r})$ in each iteration is

$\tilde{c}^{r}=\frac{\tilde{c}^{g}}{\tilde{c}^{f}}=\gamma^{(i)}\left(\frac{\tilde{c}^{g,0}}{\tilde{c}^{f,0}}\right),\qquad\gamma^{(i)}=\frac{1-(i-1)\cdot
0.1}{1+(i-1)\cdot 0.1},$

where $\gamma^{(i)}$ is the relative cyberattacker cost ratio for iteration
$i\in\mathcal{I}$.

Figure 6 shows that increasing $\gamma$ increases the impact and duration of
load shedding. For $\gamma=1$, load shedding under the Compound Cyber-Physical
Threat happens between 11:00-18:00 and peaks at 0.7 GW at 17:00, while for
$\gamma=3$ load shedding spans all 24 hours and peaks at 6.6 GW. Unserved load
in all hours is larger under the Compound Cyber-Physical Threat than in the
Cyberattack. NYISO unserved electricity in a typical summer 24-hour period
grows exponentially in $\gamma$ for the Cyberattack and Compound Cyber-
Physical Threat scenarios, as shown in Figure 7. We also observe that greater
$\gamma$ decreases the deviation of unserved electricity between the
Cyberattack and the Compound Cyber-Physical Threat scenarios. A lower
transmission attacker cost allows for compromising more transmission capacity
with the same budget. For a small enough cyberattacker transmission cost, the
impact is extremely severe, up to 25% of unserved electricity for $\gamma=3$,
or 4.95 million customers, under both scenarios compared to no more than 1% of
unserved electricity for $\gamma=1$, or 198,000 customers. Therefore, for
$\gamma\geq 3$, the impact of the Compound Cyber-Physical Threat depends
vastly on the capabilities of the capabilities. Finally, the cyberattacker
strategy in Figure 8 provides more insights on the exponential growth of
unserved electricity for increasing $\gamma$. When $\gamma=3$, compromising
transmission capacity is cheap enough for the cyberattacker to compromise more
than 70% of NYISO transmission capacity and an average of 82% of NYISO
transmission capacity each hour under both the Cyberattack and Compound Cyber-
Physical Threat scenarios. By attacking transmission capacity, the Genesee
zone disconnects from neighboring zones. When Genesee is disconnected from
neighboring zones, compromised Genesee generation capacity has direct impact
on unserved load in Genesee. Therefore, the cyberattacker chooses to also
compromise 500 MW of Genesee generation capacity in all hours, or 1% of NYISO
installed capacity, leading to 46% of unserved load in Genesee.

Figure 6: Unserved load during a summer 24-hour period for varying $\gamma$
under the Cyberattack (left) and Compound Cyber-Physical Threat (right)
scenarios. Higher $\gamma$ increases the duration of unserved load and the
impact for all hours and scenarios. Figure 7: Unserved electricity during a
summer 24-hour period for varying relative cyberattacker cost ratio
$(\gamma)$. Unserved electricity grows exponentially in the Compound Cyber-
Physical Threat and Cyberattack scenarios. While unserved electricity is
higher under the Compound Cyber-Physical Threat compared to the Cyberattack,
the deviation between scenarios decreases for increasing $\gamma$. Figure 8:
Cyberattacker strategy for varying relative cyberattacker cost ratio
$(\gamma)$. For $\gamma=3$, the cyberattacker cost of compromising
transmission capacity is low enough to compromise more than 70% of
transmission capacity in all hours. The cyberattacker compromises 82% of NYISO
transmission capacity on average, while compromised generation capacity
remains almost unchaged.

Figure 9 shows the spatial distribution of load curtailment within NYISO
during a typical summer day for the time-period with longest load curtailment.
Electricity importing regions, i.e., regions whose average electricity demand
exceeds regional installed capacity, are the most affected. For example, the
Genesee (B), Mohawk (E), and Dunwoodie (I) regions face 56%, 48%, and 100%
load curtailment during the peak curtailment period at 17:00. Finally, for
$\gamma=3$, the Cyberattack timed with the Heatwave increases unserved load in
New York City by 10% compared to the Cyberattack under normal operating
conditions. In the Compound Cyber-Physical Threat scenario, the attacker
compromises transmission capacity between Capital (F) and Hudson Valley (G),
and Dunwoodie (I) and New York City (J), which isolates New York City from
upstate supply from Capital (F) and Hudson Valley (G). Unserved load increases
from 2 GW in the Cyberattack scenario to 2.2 GW in the Compound Cyber-Physical
Threat scenario.

Figure 9: Regional distribution of unserved load within the NYISO grid for a
summer day, at 17:00. The cyberattacker compromises less than 5% and more than
95% of transmission capacity for black and red lines respectively. tha For
$\gamma=3$, the attacker compromises two additional interconnections compared
to the Cyberattack scenario between Capital (F) and Hudson Valley (G), and
Dunwoodie (I) and New York City (J), leading to a 10% increase in New York
City unserved load during peak load shedding.

####  4.2.2  Availability of Resources to the Cyberattacker

In this section we explore the sensitivity of our results to the availability
of cyberattacker resources. We use the default cyberattacker costs defined in
Section 2, but increase the budget $(\tilde{b})$ by 20% compared to the
default budget $(\tilde{b}^{0})$ in each iteration $i\in\mathcal{I}$, i.e.,

$\tilde{b}=\beta^{(i)}\cdot b^{0},\qquad\beta^{(i)}=(1+(i-1)\cdot 0.2),$

where $\beta^{(i)}$ is the relative cyberattacker resources ratio for
iteration $i\in\mathcal{I}$.

Figure 10 shows that increasing the cyberattacker budget $\tilde{b}$ increases
the impact and duration of load shedding. For $\beta=1$, load shedding under
the Compound Cyber-Physical Threat happens between 11:00-18:00 and peaks at
0.7 GW at 17:00, while for $\beta=2$ load shedding spans all 24 hours and
peaks at 1.7 GW. Unserved load in all hours is greater under the Compound
Cyber-Physical Threat compared to the Cyberattack. NYISO unserved electricity
in a typical summer 24-hour period grows linearly in $\beta$ for the
Cyberattack and Compound Cyber-Physical Threat scenarios, as shown in Figure
11. Greater $\beta$ preserved the percentage deviation of unserved electricity
between the Cyberattack and the Compound Cyber-Physical Threat scenarios.
Figure 12 shows that the cyberattacker chooses a similar strategy as in
Section 4.1, but can cause greater damage due to the availability of
additional resources. The cyberattacker can cause similar damage by targeting
transmission and generation capacity for all time periods except for 17:00,
hence the choice to compromise transmission capacity. During the peak
electricity demand hour, close to 17:00, the cyberattacker can cause
additional unserved load by targeting capacity. Hence, during that hour, the
cyberattacker focuses all their resources on generation capacity. When the
cyberattacker doubles their budget, unserved load increases by 3%, or roughly
594,000 customers, in both scenarios.

Figure 10: Unserved load during a summer 24-hour period for varying $\beta$
under the Cyberattack (left) and Compound Cyber-Physical Threat (right)
scenarios. Higher $\beta$ increases the duration of unserved load and the
impact for all hours and scenarios. Figure 11: Unserved electricity during a
summer 24-hour period for varying relative cyberattacker resources ratio
$(\beta)$. Unserved electricity grows linearly in the Compound Cyber-Physical
Threat and Cyberattack scenarios. Unserved electricity is higher under the
Compound Cyber-Physical Threat compared to the Cyberattack and grows by 3% in
both scenarios when the cyberattacker resources double with respect to the
default budget ($\beta=2$). Figure 12: Cyberattacker strategy for varying
relative cyberattacker resources ratio $(\beta)$. For double resources
($\beta=2$), the cyberattacker chooses to allocate their budget between
compromising transmission and generation capacity in each hour similarly to
the default budget scenario ($\beta=1$), i.e., compromise generation capacity
during the peak load shedding period.

Figure 13 reveals that the spatial distribution of unserved load remains
insensitive to increasing the budget, but the magnitude of unserved demand is
magnified. Although the total curtailed electricity in the 24-hour period
increases to 3% when doubling the budget, Long Island unserved load increases
to 100% under the Cyberattack and 109% under the Compound Cyber-Physical
Threat compared to the Baseline. The Cyberattack scenario results show that
doubling the budget of the cyberattacker can compromise all of Long Island
generation. In the Compound Cyber-Physical Threat, the system operator sheds
109% of the Baseline load, i.e., the load under normal operating conditions,
and the additional 9% of the Heatwave-induced demand increase.

Figure 13: Regional distribution of unserved load within the NYISO grid for a
summer day, at 17:00. For $\beta=2$, the attacker has the capabilities to
compromise enough generation and transmission capacity in Long Island to and
cause a blackout in the whole region during normal operating conditions and
during the Heatwave scenarios.

###  4.3 Economy-Wide Impacts

In this section, we evaluate the ripple effects of the Compound Cyber-Physical
Threat across all economic sectors using the WiNDC model. In Figure 14 we can
observe the compounding nature of a timed Cyberattack and a Heatwave. Under
the Cyberattack and Heatwave scenarios there is no or little unserved load
respectively, similarly to the Baseline scenario, and economic losses are
negligible. However, under the Compound Cyber-Physical Threat, the New York
State economy incurs $2 Million losses. Increasing the cyberattacker budget
magnifies economic losses, up to $21 Million for a budget twice the default
cyberattacker budget. Under the Compound Cyber-Physical Threat, doubling the
cyberattacker budget causes a 9-fold increase of unserved load in the electric
power sector, from 2.4 GW to 21.8 GW, and a 10-fold increase of economy-wide
losses. Similarly, increasing the vulnerability of transmission components to
cyberattacks $(\gamma)$ increases economy-wide losses exponentially to $120
Million for $\gamma=3$, i.e., 2.3% of the daily GDP.

The economic impact can vary across sectors. From a total of 71 sectors, only
5 sectors report losses greater than 0.5% under the Compound Cyber-Physical
Threat with $\beta=2$. When $\gamma=3$, 16 sectors report losses greater than
0.5%. Figure 15 shows that the impact varies even within the top-5 most
impacted sectors under the Compound Cyber-Physical Threat. The “State and
local government enterprises” face 7% losses under the Compound Cyber-Physical
Threat with $\beta=2$ and 34% losses when $\gamma=3$. The top-4 sectors are
the same for both scenarios and include also the “Federal government
enterprises (GFE)” and “Amusements, gambling, and recreation industries”, with
average losses of 1% when $\beta=2$ and 5% when $\gamma=3$. The “Educational
services” sector completes the top-5 with 3% for $\gamma=3$ and the ranks 5th
with 1% for $\beta=2$.

Figure 14: Economic losses measured as decrease of daily GDP compared to
Baseline GDP. Doubling the cyberattacker resources aggravates economy-wide
losses by 10 times, while enhanced transmission vulnerability can lead to a
60-fold increase of economy-wide losses. Figure 15: Top-5 most impacted
sectors under the Compound Cyberphysical Threat scenario with $\gamma=3$ (up)
and $\beta=2$ (down). The most affected sector is the “State and local
government enterprises (GSLE)” in both cases.

##  5 Discussion

Our results reveal which electric power infrastructure components are most
vulnerable to cyberattacks under extreme weather stressors. Depending on the
budget of the attacker, unserved load in the NYISO increases from 0 to 3% as
the budget doubles, however the impact is highly localized. Less
interconnected NYISO zones are heavily targeted, leading to disproportionally
more unserved load compared to more interconnected zones with greater local
electricity supply, which remains relatively unaffected.

###  5.1 Electric Power Network Resilience and Economy-Wide Response

#####  The cyberattacker targets nodes of the NYISO electric power network
with low network connectivity.

The Long Island and New York City zones have the second and third greatest
installed capacity. Moreover, the aggregate Long Island and New York City
demand absorbs all of Long Island and New York City electricity production
respectively. The two zones are similar in that they connect with other NYISO
zones via relatively low-capacity interconnections. Therefore, compromising
generation and transmission capacity in these two nodes translates directly to
unserved load. Note that excess transmission capacity in the Long Island zone
implies that the cyberattacker needs to compromise a larger portion of NYISO
transmission capacity, which is costlier compared to compromising generation
capacity, to cause unserved load. Thus, between the two zones with low
connectivity, the cyberattacker compromises the Long Island zone first.

#####  For each targeted node, the cyberattacker compromises generation
capacity before compromising power flow capacity.

In Section 4, we observed that the first targeted zone was Long Island. For
example, under default cyberattacker costs and budget, i.e., $\beta=1$ and
$\gamma=1$, at 17:00, the cyberattacker compromises 20% of Long Island
generation capacity. As the attacker budget doubles, the cyberattacker
compromises more generation capacity. In the case of Long Island and New York
City, the cyberattacker prefers compromising generation capacity before
transmission capacity for two reasons. First, compromising generation capacity
requires less resources. In Section 2 we assume that power generation follow
less rigorous cybersecurity protocols, compared to the transmission system,
hence are more vulnerable to cyberattacks. For that, a cyberattacker needs
less resources to compromise a MW of generation capacity compared to
transmission capacity. Second, the low connectivity of both zones with the
other NYISO zones implies that the two zones have less access to excess
capacity within the NYISO grid. Given their low connectivity and that the Long
Island and New York City zones cover almost all of their load from local
production, compromising generation capacity has a direct impact on unserved
load and the cyberattacker targets the Long Island zone first.

##### Increasing the vulnerability of transmission components increases
unserved load exponentially, while increasing the cyberattacker resources
availability increases unserved load linearly.

We model increasing vulnerability of transmission components by decreasing the
cyberattacker transmission cost, which allows the cyberattacker to compromise
more transmission capacity for the same budget. For $\gamma=3$, the cost of
compromising transmission capacity is so low, that the cyberattacker can
compromise most of the NYISO transmission capacity. On average, the
cyberattacker compromises 82% of transmission capacity in the 24-hour period
for $\gamma=3$. When $\gamma=3$, the cyberattacker can compromise more than
70% of transmission capacity, which leaves most zones largely dependent on
their domestic generation. When transmission components vulnerability is
lower, the attacker can compromise less transmission capacity, due to their
resource constraint, and excess transmission capacity allows the NYISO to
absorb the impact. However, under high vulnerability, significant transmission
capacity is compromised, and additional losses in transmission capacity lead
to unserved load in multiple zones. For example, compromising transmission
between Hudson Valley (G) and Dunwoodie (I) disconnects most major upstate
generators from supplying New York City and Long Island demand. On the other
hand, increasing the budget under default cyberattacker costs leads to
targeting generation components. However, excess generation capacity in the
NYISO grid and the ability to substitute zonal generation with imports from
other NYISO zones leads to an increase of 19.3 GWh in unserved load, affecting
more than 594,000 customers, for double the cyberattacker resources.

#####  A Compound Cyber-Physical Threat compounds also economy-wide losses.

The lack of regional electricity supply constrains the activity of other
sectors. Across all economic sectors, the Heatwave and Cyberattack do not
impact other sectors, while a Compound Cyber-Physical Threat can decrease
economic activity by $2 million. We observe the largest impact for the largest
operational disruption when $\gamma=3$, which leads to losses of $120 million
or 2.7% of daily New York State GDP. Beyond the “Utilities” sector, we
identify the “Federal government enterprises”, “State and local enterprises”
and “Transit and ground passenger transportation” sectors as particularly
vulnerable to disruptions in electricity supply in New York State. These
enterprises provide essential services for local well-being, e.g.,
transportation. Therefore, disruptions in the activity of these three sectors
can impact a broad range of population groups. The assumptions underpinning a
CGE, discussed in the next Section, constrain the quantification of economy-
wide losses. Our analysis does not aim to quantify economy-wide losses or the
VOLL, which requires extensive ex-post surveys. Our goal is to provide an ex-
ante computation of the lower bound, meaning that the impact of the Compound
Cyber-Physical Threat to New York State GDP and individual sectors will be at
least as high as in our analysis. For example, across scenarios, we derive
losses between $1,054-$1,078 per MWh in the 24-hour period. For comparison,
the estimated losses for other than compound cyber-physical stressors, e.g.,
26-hour long 1977 New York City blackout and 48-hour long 2003 Northeastern
blackout, were $3,405 and $9,300 per MWh, respectively [1]. The response of
the economy is the result of the microeconomic principles underpinning CGEs,
hence the results serve as an indicator of sectors which are more likely to
face disruptions.

###  5.2 Limitations

The limitations of our study are inherent to the nature of the Compound Cyber-
Physical Threat against critical infrastructures and the computational
capabilities of existing methodologies.

First, actual data on features of the transmission and distribution system are
confidential [9]. Hence, our results are sensitive to electric power network
features, which in our case include transmission capacities, voltage angle
difference limits, and generation capacities.

Second, our results are sensitive to the assumptions on the cyberattacker’s
rationale. In this paper, we assume scenarios that are rapid and have
immediate, localized impact on consumers. Alternatively, a cyberattacker may
opt for longer-term strategies to exploit less severe events. For example, a
cyberattacker may want to maximize the operating cost of the NYISO. In this
case, the cyberattacker may attempt to inflict less severe damage such that
they go unnoticed. Moreover, under a scenario of NYISO cost maximization, the
cyberattacker may opt to compromise low-cost resources first in order to force
the NYISO to dispatch more expensive fuel-fired generators into the fuel mix
more often. When factoring in interdependent infrastructures, the
cyberattacker may decide to compromise power supply to zones with critical
telecommunications components or wastewater treatment facilities. Our
framework allows the flexibility to adjust the cyberattacker objectives
$\tilde{C}(\cdot)$. Moreover, the MCP reformulation of the lower-level problem
can be enhanced to account for interdependencies with other critical sectors,
e.g., the water/ wastewater sector.

Third, our results are sensitive to the elasticities of substitution between
production inputs of CGE sectors. In Section 3.4 we mentioned that CGEs have
yearly time steps and thus input substitutions happen at a longer time-scale
compared to the 24-hour time-frame studied here. However, it may be more
difficult to substitute between inputs in the short-term, e.g., between
different energy sources. In addition, post-disruption elasticities of
substitution are not available, even for individual sectors. Hence, our
analysis is not suitable for the quantification of economic losses. Instead,
we aim to identify vulnerable sectors and their response, and provide an ex-
post calculation of a lower bound on economic losses which is consistent in
the orders of magnitude.

##  6 Conclusions

###  6.1 Summary and Contributions

This paper introduces the novel Compound Cyber-Physical Threat, arising when a
cyberattack targets electric power network components when they are
critically-stressed due to extreme weather events. Our work contributes to the
identification and management of electric power network vulnerabilities, and
economy-wide losses. We find that when generators are more easily compromised,
the cyberattcker targets NYISO zones with low network connectivity. While a
Cyberattack under normal operating conditions or a Heatwave barely disrupt the
NYISO electric power network, a Compound Cyber-Physical Threat exacerbates
regional unserved electricity demand under the individual threats by more than
three times, leading to 13% unserved load in a 24-hour summer period. When
transmission lines are easier to compromise, the cyberattacker targets NYISO
zones with low network connectivity and a lack of generation capacity. The
Compound Cyber-Physical Threat exacerbates also economy-wide losses across
sectors. The individual Heatwave and Cyberattack scenarios result in an
insignificant GDP reduction. However, depending on the cyberattacker cost of
each strategy, GDP reduction ranges between 0-2.7% of daily GDP.

In the context of cyberattacks timed with extreme weather events, we introduce
a framework that accounts for the operational capabilities of a cyberattacker
and test the sensitivity of our framework against the cyberattacker resource
availability. In addition, our framework couples the cyberattacker strategy
with the response of all sectors in the economy. We also contribute to the
derivation of a lower bound of economy-wide losses at the state level. We show
that our estimation is consistent in orders of magnitude with estimates of
losses in past blackouts. Although the quantification of economy-wide losses
is beyond the capabilities of our framework, further refinement of the
proposed methodology can contribute to estimating a “reference point” for the
VOLL under Compound Cyber-Physical Threats. Therefore, further development of
our framework and realistic data can assist system operators and federal,
state and local agencies in infrastructure planning against extreme weather
events and cyberattacks.

###  6.2 Future Research

This paper focuses on heatwaves, but our framework is applicable to emerging,
billion-dollar extreme weather events, including hurricanes, storms, and
wildfires. Hurricanes, storms, and wildfires follow a trajectory, thus
disproportionately affecting infrastructures across the trajectory compared to
the rest of the system. Therefore, their impact is highly localized.
Understanding the regional impact of extreme weather events, also across
communities with varying sociodemographic characteristics, requires a
methodology with greater spatial resolution. Hence, future research should
model both the distribution and the transmission system in a single model.
Moreover, the analysis of distribution-level impacts requires a county-level
CGE dedicated to extreme events. Such CGE requires recalibrating the sectoral
inputs substitution elasticities to account for the short-term effects during
and following the extreme event.

Second, this paper focuses on electric power networks, due to their critical
contribution to the operation of most sectors. However, a disruption in the
electricity sector can propagate to multiple critical infrastructures through
their interdependencies and further compound the impact of an extreme event on
regional infrastructures and economic activity. Therefore, going beyond
operational disruptions and accounting for the public health impacts of
Compound Cyber-Physical Threats requires understanding the propagation of
failures in the water/wastewater, telecommunications, and transportation
infrastructures.

##  7 Acknowledgements

The authors are grateful to Dr. Maxwell Brown and Dr. Jon Becker for their
generous advice on the implementation of WiNDC.

## References

  * Abidi and Martins, [2015] Abidi, M. M. and Martins, N. (2015).  Impact of Power System Blackouts.  Panel Session at the IEEE PES GM 2015 “Cascading Failures: Advanced Methodologies, Restoration and Industry Perspectives”.
  * Acharya et al., [2020] Acharya, S., Dvorkin, Y., and Karri, R. (2020).  Public plug-in electric vehicles + grid data: Is a new cyberattack vector viable?  IEEE Transactions on Smart Grid, 11(6):5099–5113.
  * Ashok et al., [2015] Ashok, A., Wang, P., Brown, M., and Govindarasu, M. (2015).  Experimental evaluation of cyber attacks on automatic generation control using a cps security testbed.  In 2015 IEEE Power & Energy Society General Meeting, pages 1–5.
  * Balistreri et al., [2020] Balistreri, E. J., Zhang, W., and Beghin, J. (2020).  The State-level Burden of the Trade War: Interactions between the Market Facilitation Program and Tariffs.  Center for agricultural and rural development (card) publications, Center for Agricultural and Rural Development (CARD) at Iowa State University.
  * BEA, [2022] BEA (2022).  Sagdp1 gross domestic product (gdp) summary, annual by states.  Accessed August 31, 2022.
  * Becker and Carbone, [2020] Becker, J. and Carbone, J. (2020).  Modeling Sectoral Labor Transitions with WiNDC.
  * Bessec and Fouquau, [2008] Bessec, M. and Fouquau, J. (2008).  The non-linear link between electricity consumption and temperature in europe: A threshold panel approach.  Energy Economics, 30(5):2705–2721.
  * Biehl et al., [2017] Biehl, E., Buzogany, S., Huang, A., Chodur, G., and Neff, R. (2017).  Baltimore food system resilience advisory report.  Technical report, Johns Hopkins Center for A Livable Future.
  * Birchfield et al., [2017] Birchfield, A. B., Gegner, K. M., Xu, T., Shetye, K. S., and Overbye, T. J. (2017).  Statistical considerations in the creation of realistic synthetic power grids for geomagnetic disturbance studies.  IEEE Transactions on Power Systems, 32(2):1502–1510.
  * Boyd and Vandenberghe, [2004] Boyd, S. and Vandenberghe, L. (2004).  Convex optimization.  Cambridge university press.
  * Capitanescu and Wehenkel, [2008] Capitanescu, F. and Wehenkel, L. (2008).  A new iterative approach to the corrective security-constrained optimal power flow problem.  IEEE Transactions on Power Systems, 23(4):1533–1541.
  * Castillo et al., [2019] Castillo, A., Arguello, B., Cruz, G., and Swiler, L. (2019).  Cyber-physical emulation and optimization of worst-case cyber attacks on the power grid.  In 2019 Resilience Week (RWS), volume 1, pages 14–18.
  * Ceferino et al., [2021] Ceferino, L., Lin, N., and Xi, D. (2021).  Bayesian Updating of Solar Panel Fragility Curves and Implications of Higher Panel Strength for Solar Generation Resilience.  engrXiv.
  * Ceferino et al., [2020] Ceferino, L., Liu, C., Alisjahbana, I., Patel, S., Sun, T., Kiremidjian, A., and Rajagopal, R. (2020).  Earthquake resilience of distributed energy resources.  In 17th World Conference on Earthquake Engineering, Tokyo, Japan.
  * CI, [2022] CI (2022).  2021 h2 healthcare data breach report.  Technical report, Critical Insight.
  * CJBS, [2015] CJBS (2015).  World cities risk 2015-2025.  Technical report, Cambridge Judge Business School, University of Cambridge.  Cambridge Judge Business School, University of Cambridge with the generous support of LLOYD’S.
  * Cohen et al., [2018] Cohen, J., Moeltner, K., Reichl, J., and Schmidthaler, M. (2018).  Effect of global warming on willingness to pay for uninterrupted electricity supply in European nations.  Nature Energy, 3(1):37–45.
  * Colombo et al., [2017] Colombo, A. W., Karnouskos, S., Kaynak, O., Shi, Y., and Yin, S. (2017).  Industrial cyberphysical systems: A backbone of the fourth industrial revolution.  IEEE Industrial Electronics Magazine, 11(1):6–16.
  * Cong et al., [2020] Cong, S., Nock, D., Qiu, Y. L., and Xing, B. (2020).  Unveiling hidden energy poverty using the energy equity gap.  Nature Communications, 13:2456.
  * Dreyer et al., [2018] Dreyer, P., Jones, T. M., Klima, K., Oberholtzer, J., Strong, A., Welburn, J. W., and Winkelman, Z. (2018).  Estimating the Global Cost of Cyber Risk: Methodology and Examples.  RAND Corporation, Santa Monica, CA.
  * Eckstrom et al., [2022] Eckstrom, S., Murphy, G., Ye, E., Acharya, S., Mieth, R., , and Dvorkin, Y. (2022).  Outing Power Outages: Real-time and Predictive Socio-demographic Analytics for New York City.  Accepted for publication in the 2022 IEEE Power and Energy Society General Meeting, Denver CO.
  * Eldridge et al., [2018] Eldridge, B., O’Neill, R., and Castillo, A. (2018).  An improved method for the dcopf with losses.  IEEE Transactions on Power Systems, 33(4):3779–3788.
  * Eling et al., [2022] Eling, M., Elvedi, M., and Falco, G. (2022).  The economic impact of extreme cyber risk scenarios.  North American Actuarial Journal, 0(0):1–15.
  * EPRI, [2015] EPRI (2015).  Analysis of selected electric sector high risk failure scenarios.  Technical report, Electric Power Research Institute.  National Electric Sector Cybersecurity Organization Resource (NESCOR).
  * ERCOT, [2013] ERCOT (2013).  Estimating the value of lost load.  Technical report, Electric Reliability Council of Texas, Inc.  Briefing paper prepared for the Electric Reliability Council of Texas, Inc. by London Economics International LLC.
  * FBI, [2022] FBI (2022).  Ransomware Attacks on Agricultural Cooperatives Potentially Timed to Critical Seasons.  Federal Bureau of Investigation Private Industry Notification 20220420-001.
  * Feng et al., [2018] Feng, M., Mitchell, J., Pang, J., Shen, X., and Waechter, A. (2018).  Complementarity formulations of $l0$-norm optimization.  Pacific Journal of Optimization, 14(2):273–305.
  * FERC, [2020] FERC (2020).  Audit of new york independent system operator. docket no. pa19-1-000.  Published July 7, 2020. Accessed August 31, 2022.
  * Gabriel et al., [2012] Gabriel, S. A., Conejo, A. J., Fuller, J. D., Hobbs, B. F., and Ruiz, C. (2012).  Complementarity Modeling in Energy Markets.  Springer.
  * Granato and Polacek, [2019] Granato, A. and Polacek, A. (2019).  The growth and challenges of cyber insurance.  Chicago Fed Letter, No. 426, 2019.
  * HMO, [2018] HMO (2018).  Houston prepares for critical infrastructure threats involving cyberattacks.
  * Huppmann and Siddiqui, [2018] Huppmann, D. and Siddiqui, S. (2018).  An exact solution method for binary equilibrium problems with compensation and the power market uplift problem.  European Journal of Operational Research, 266(2):622–638.
  * IECC, [2021] IECC (2021).  The state of cyber report 2017-2021.  Technical report, Indiana’s Executive Council on Cybersecurity.
  * II et al., [2022] II, J. L. B., Hagen, A., and Berg, R. (2022).  Hurricane ida (al092021).  Technical report, National Hurricane Center.  National Hurricane Center Tropical Cyclone Report, 26 August–1 September 2021.
  * Jones and O’Neill, [2016] Jones, B. and O’Neill, B. C. (2016).  Spatially explicit global population scenarios consistent with the shared socioeconomic pathways.  Environmental Research Letters, 11(8):084003.
  * Karangelos and Wehenkel, [2022] Karangelos, E. and Wehenkel, L. (2022).  Cyber–physical risk modeling with imperfect cyber-attackers.  Electric Power Systems Research, 211:108437.
  * Khan et al., [2022] Khan, H. A. U., Kim, J., and Dvorkin, Y. (2022).  Risk-informed participation in t&d markets.  Electric Power Systems Research, 202:107624.
  * Khaw et al., [2021] Khaw, Y. M., Abiri Jahromi, A., Arani, M. F. M., Sanner, S., Kundur, D., and Kassouf, M. (2021).  A deep learning-based cyberattack detection system for transmission protective relays.  IEEE Transactions on Smart Grid, 12(3):2554–2565.
  * Klinenberg, [1999] Klinenberg, E. (1999).  Denaturalizing disaster: A social autopsy of the 1995 chicago heat wave.  Theory and Society, 28(2):239–295.
  * Krishnamurthy et al., [2019] Krishnamurthy, P., Salehghaffari, H., Duraisamy, S., Karri, R., and Khorrami, F. (2019).  Stealthy rootkits in smart grid controllers.  In 2019 IEEE 37th International Conference on Computer Design (ICCD), pages 20–28.
  * Latto et al., [2021] Latto, A., Hagen, A., and Berg, R. (2021).  Hurricane isaias (al092020).  Technical report, National Hurricane Center.  National Hurricane Center Tropical Cyclone Report, 30 July–4 August 2020\.
  * Leyffer et al., [2006] Leyffer, S., López-Calva, G., and Nocedal, J. (2006).  Interior methods for mathematical programs with complementarity constraints.  SIAM Journal on Optimization, 17(1):52–77.
  * Leyffer and Munson, [2010] Leyffer, S. and Munson, T. (2010).  Solving multi-leader–common-follower games.  Optimization Methods and Software, 25(4):601–623.
  * Liang et al., [2016] Liang, J., Sankar, L., and Kosut, O. (2016).  Vulnerability analysis and consequences of false data injection attack on power system state estimation.  IEEE Transactions on Power Systems, 31(5):3864–3872.
  * Luo et al., [1996] Luo, Z.-Q., Pang, J.-S., and Ralph, D. (1996).  Mathematical Programs with Equilibrium Constraints.  Cambridge University Press.
  * Markhvida et al., [2020] Markhvida, M., Walsh, B., Hallegatte, S., and Baker, J. (2020).  Quantification of disaster impacts through household well-being losses.  Nature Sustainability, 3:538–547.
  * McLaughlin et al., [2016] McLaughlin, S., Konstantinou, C., Wang, X., Davi, L., Sadeghi, A.-R., Maniatakos, M., and Karri, R. (2016).  The cybersecurity landscape in industrial control systems.  Proceedings of the IEEE, 104(5):1039–1057.
  * Moynihan et al., [2022] Moynihan, E., Avraam, C., Siddiqui, S., and Neff, R. (2022).  Optimization based modeling for the food supply chain’s resilience to outbreaks.  Frontiers in Sustainable Food Systems, 6.
  * NYC, [2013] NYC (2013).  A stronger, more resilient new york.  Technical report, NYC Special Initiative for Rebuilding & Resiliency.
  * [50] NYISO (2019a).  Load & capacity data.  Technical report, New York Independent System Operator, Inc.  Gold Book. A report by The New York Independent System Operator, Inc.
  * [51] NYISO (2019b).  Load & capacity data.  Technical report, New York Independent System Operator, Inc.  Gold Book. A report by The New York Independent System Operator, Inc.
  * OIES, [2022] OIES (2022).  Natural gas demand in europe: The impacts of covid-19 and other influences in 2020.  Technical report, Oxford Institute for Energy Studies.
  * Patel et al., [2021] Patel, S., Ceferino, L., Liu, C., Kiremidjian, A., and Rajagopal, R. (2021).  The disaster resilience value of shared rooftop solar systems in residential communities.  Earthquake Spectra, 37(4):2638–2661.
  * Puri et al., [2013] Puri, K., Dietachmayer, G., Steinle, P., Dix, M., Rikus, L., Logan, L., Naughton, M., Tingwell, C., Xiao, Y., Barras, V., Bermous, I., Bowen, R., Deschamps, L., Franklin, C., Fraser, J., Glowacki, T., Harris, B., Lee, J., Le, T., and Engel, C. (2013).  Implementation of the initial access numerical weather prediction system.  Australian Meteorological and Oceanographic Journal, 63:265–284.
  * Ralph, [2008] Ralph, D. (2008).  Mathematical programs with complementarity constraints in traffic and telecommunications networks.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 366(1872):1973–1987.
  * Ramamurthy et al., [2017] Ramamurthy, P., Li, D., and Bou-Zeid, E. (2017).  High-resolution simulation of heatwave events in New York City.  Theoretical and Applied Climatology, 128(1-2):89–102.
  * Romitti and Wing, [2022] Romitti, Y. and Wing, I. (2022).  Heterogeneous climate change impacts on electricity demand in world cities circa mid-century.  Scientific Reports, 12:4280.
  * Rose and Guha, [2004] Rose, A. and Guha, G.-S. (2004).  Computable General Equilibrium Modeling of Electric Utility Lifeline Losses from Earthquakes, pages 119–141.  Springer Berlin Heidelberg, Berlin, Heidelberg.
  * Rose et al., [2007] Rose, A., Oladosu, G., and Liao, S.-Y. (2007).  Business interruption impacts of a terrorist attack on the electric power system of los angeles: Customer resilience to a total blackout.  Risk Analysis, 27(3):513–531.
  * Rutherford and Schreiber, [2019] Rutherford, T. F. and Schreiber, A. (2019).  Tools for Open Source, Subnational CGE Modeling with an Illustrative Analysis of Carbon Leakage.  Journal of Global Economic Analysis, 4(2):1–66.
  * Sankaranarayanan et al., [2018] Sankaranarayanan, S., Feijoo, F., and Siddiqui, S. (2018).  Sensitivity and covariance in stochastic complementarity problems with an application to north american natural gas markets.  European Journal of Operational Research, 268(1):25–36.
  * Scheel and Scholtes, [2000] Scheel, H. and Scholtes, S. (2000).  Mathematical programs with complementarity constraints: Stationarity, optimality, and sensitivity.  Mathematics of Operations Research, 25(1):1–22.
  * Scott et al., [2018] Scott, A. A., Waugh, D. W., and Zaitchik, B. F. (2018).  Reduced urban heat island intensity under warmer conditions.  Environmental Research Letters, 13(6):064003.
  * Siddiqui and Gabriel, [2013] Siddiqui, S. and Gabriel, S. (2013).  An SOS1-Based Approach for Solving MPECs with a Natural Gas Market Application.  Networks and Spatial Economics, 13(2):205–227.
  * Smith, [2021] Smith, A. (2021).  2020 U.S. billion-dollar weather and climate disasters in historical context.  Technical report, National Oceanic and Atmospheric Administration.
  * Sullivan et al., [2009] Sullivan, M., Mercurio, M., and Schellenberg, J. (2009).  Estimated value of service reliability for electric utility customers in the united states.  Technical report, Lawrence - Berkeley National Laboratory.
  * USCB, [2021] USCB (2021).  U.S. Census 2020.
  * USCPSOTF, [2004] USCPSOTF (2004).  Final report on the august 14, 2003 blackout in the united states and canada: Causes and recommendations.  Technical report, U.S.-Canada Power System Outage Task Force.
  * USDC, [2013] USDC (2013).  Economic impact of hurricane sandy potential economic activity lost and gained in new jersey and new york.  Technical report, United States Department of Commerce.  Economics and Statistics Administration Office of the Chief Economist.
  * USDHS, [2015] USDHS (2015).  NCCIC/ICS-CERT Year in Review.  Technical report, U.S. Department of Homeland Security.
  * USTR, [2020] USTR (2020).  Economic and Trade Agreement between the Government of the United States of America and the Government of the People’S Republic of China.
  * Wang and Karri, [2016] Wang, X. and Karri, R. (2016).  Reusing hardware performance counters to detect and identify kernel control-flow modifying rootkits.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 35(3):485–498.
  * Wang et al., [2015] Wang, X., Konstantinou, C., Maniatakos, M., and Karri, R. (2015).  Confirm: Detecting firmware modifications in embedded systems using hardware performance counters.  In 2015 IEEE/ACM International Conference on Computer-Aided Design (ICCAD), pages 544–551.
  * Wang et al., [2016] Wang, X., Konstantinou, C., Maniatakos, M., Karri, R., Lee, S., Robison, P., Stergiou, P., and Kim, S. (2016).  Malicious firmware detection with hardware performance counters.  IEEE Transactions on Multi-Scale Computing Systems, 2(3):160–173.
  * WEF, [2022] WEF (2022).  Why the energy sector’s latest cyberattack in Europe matters.  World Economic Forum.
  * Wogrin and Gayme, [2015] Wogrin, S. and Gayme, D. F. (2015).  Optimizing storage siting, sizing, and technology portfolios in transmission-constrained networks.  IEEE Transactions on Power Systems, 30(6):3304–3313.
  * Xu and Hobbs, [2019] Xu, Q. and Hobbs, B. F. (2019).  Value of model enhancements: quantifying the benefit of improved transmission planning models.  IET Generation, Transmission & Distribution, 13(13):2836–2845.
  * Yuan et al., [2011] Yuan, Y., Li, Z., and Ren, K. (2011).  Modeling load redistribution attacks in power systems.  IEEE Transactions on Smart Grid, 2(2):382–390.

## Supplementary Material: Abbreviations and Nomenclature

Abbreviations

DC-OPF | Direct Current Optimal Power Flow
---|---
MCP | Mixed Complementarity Problem
MILP | Mixed Integer Linear Program
MPEC | Mathematical Program with Equilibrium Constraints
MPCC | Mathematical Program with Complementarity Constraints
ICC | Information communication and control
Table 2:  Abbreviations.

Nomenclature

Definition | Description
---|---
$\mathcal{N}$ |  Set of electric power network nodes of Manhattan
$\mathcal{E}$ |  Set of electric power network edges of Manhattan
$\mathcal{S}$ |  Set of weather seasons
$\mathcal{H}$ |  Set of 1-hour time steps
$\mathcal{G}$ |  Set of power generation technologies by node
$\mathcal{Y}$ |  Set of feasible space of all DC-OPF variables
$\mathcal{Y}_{1}$ |  Set of feasible space of all non-negative DC-OPF variables
$\mathcal{Y}_{2}$ |  Set of feasible space of all free DC-OPF variables
$\mathcal{Y}_{1sh}^{(i)}$ |  Set of feasible space of coordinate $i$ of set $\mathcal{Y}_{1}$, $i=1,2,\dots,5$
$\mathcal{Y}_{2sh}^{(i)}$ |  Set of feasible space of coordinate $i$ of set $\mathcal{Y}_{2}$, $i=1,2,\dots,8$
$N=|\mathcal{N}|$ |  Parameter of cardinality of set $\mathcal{N}$
$E=|\mathcal{E}|$ |  Parameter of cardinality of set $\mathcal{E}$
$H=|\mathcal{H}|$ |  Parameter of cardinality of set $\mathcal{H}$
$G=|\mathcal{G}|$ |  Parameter of cardinality of set $\mathcal{T}$
$Y=|\mathcal{Y}|$ |  Parameter of cardinality of set $\mathcal{Y}$
$Y_{1}=|\mathcal{Y}_{1}|$ |  Parameter of cardinality of set $\mathcal{Y}_{1}$
$Y_{2}=|\mathcal{Y}_{2}|$ |  Parameter of cardinality of set $\mathcal{Y}_{2}$
$A\in M_{E\times N}\left(\mathbb{R}\right)$ |  Incidence matrix of the electric power network
$b\in\mathbb{R}^{E}$ |  Parameter of susceptance of transmission lines
$B=\text{diag}\left\\{b\right\\}\in M_{E\times E}\left(\mathbb{R}\right)$ |  Susceptance Matrix of transmission lines
$e_{i}\in\mathbb{R}^{N}$ |  Vector of all zeros except for entry $i$ value of “1”
$z^{g}_{sh}\in\mathbb{R}^{G}$ |  Upper Level Variable of compromised power generation capacity
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (MWh)
$z^{\theta}_{sh}\in\mathbb{R}^{N}$ |  Upper Level Variable of compromised angle difference limit
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (rad)
$z^{f}_{sh}\in\mathbb{R}^{E}$ |  Upper Level Variable of compromised power flow capacity in all edges
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (MWh)
$\tilde{b}\in\mathbb{R}_{\geq 0}$ |  Parameter of budget for all resources of the cyberattacker
$\tilde{c}^{g}_{sh}\in\mathbb{R}^{G}$ |  Parameters of cost of compromising power generation capacity
| (proportional to budget $\tilde{b}$)
$\tilde{c}^{\theta}_{sh}\in\mathbb{R}^{E}$ |  Parameters of cost of compromising voltage angle differences
| (proportional to budget $\tilde{b}$)
$\tilde{c}^{f}_{sh}\in\mathbb{R}^{E}$ |  Parameters of cost of compromising power flow capacity
| (proportional to budget $\tilde{b}$)

Nomenclature (continued)

Definition | Description
---|---
$g_{sh}\in\mathbb{R}^{G}$ |  Lower-Level Primal Variable of power generation
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (MWh)
$\theta_{sh}\in\mathbb{R}^{N}$ |  Lower-Level Primal Variable of angle difference between nodes
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (rad)
$f_{sh}\in\mathbb{R}^{E}$ |  Lower-Level Primal Variable of power flow in all edges
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (MWh)
$u_{sh}\in\mathbb{R}^{N}_{\geq 0}$ |  Lower-Level Primal Variable of unserved electricty demand
| in all nodes and in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$
(MWh)
$y^{g}_{sh}\in\mathbb{R}^{G}$ |  Lower-Level Primal Variable of investment
| in power generation capacity in season $s\in\mathcal{S}$,
$y^{f}_{sh}\in\mathbb{R}^{E}$ |  Lower-Level Primal Variable of investment
| in power flow capacity in season $s\in\mathcal{S}$, during hour
$h\in\mathcal{H}$ (MW)
$\underline{g}_{sh},\bar{g}_{sh}\in\mathbb{R}^{G}$ |  Parameters of lower and upper power generation limits
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (MWh)
$\underline{f}_{sh},\bar{f}_{sh}\in\mathbb{R}^{E}$ |  Parameters of lower and upper power flow limits
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (MWh)
$\underline{\theta}_{sh},\bar{\theta}_{sh}\in\mathbb{R}^{E}$ |  Parameters of lower and upper angle difference limits across all edges
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (rad)
$c^{g,l}_{sh}\in\mathbb{R}^{G}_{\geq 0}$ |  Parameters of linear term of generation cost function
$c^{g,q}_{sh}\in\mathbb{R}^{G}_{\geq 0}$ |  Parameters of quadratic term of generation cost function
$c^{f,l}_{sh}\in\mathbb{R}^{G}_{\geq 0}$ |  Parameters of linear term of flow cost function
$c^{f,q}_{sh}\in\mathbb{R}^{G}_{\geq 0}$ |  Parameters of quadratic term of flow cost function
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (1000$/MWh)
$c^{u}_{sh}\in\mathbb{R}^{N}_{\geq 0}$ |  Parameter of value of lost load
| in season $s\in\mathcal{S}$, during hour $t\in\mathcal{T}$ (1000$/MWh)
$d_{sh}\in\mathbb{R}^{N}_{\geq 0}$ |  Parameter of total electricity demand
$d^{w}_{sh}\in\mathbb{R}^{N}_{\geq 0}$ |  Parameter of electricity demand from the water/wastewater sector
| in season $s\in\mathcal{S}$, during hour $h\in\mathcal{H}$ (MWh)
$d^{o}_{sh}\in\mathbb{R}^{N}_{\geq 0}$ |  Parameter of electricity demand from all sectors excluding
| the water/wastewater sector in season $s\in\mathcal{S}$,
| during hour $h\in\mathcal{H}$ (MWh)
$\pi^{d}_{sh}\in\mathbb{R}^{N}$ |  Lower-Level Dual Variable of market-clearing constraint (1f)
| (1000$/MWh)
$\pi^{f}_{sh}\in\mathbb{R}^{E}$ |  Lower-Level Dual Variable of power flow voltage
| angle constraint (1g) (1000$/MWh)
$\underline{\rho}^{g}_{sh},\bar{\rho}^{g}_{sh}\in\mathbb{R}^{G}_{\geq 0}$ |  Lower-Level Dual Variables of power generation constraints (1000$/MW)
$\underline{\rho}^{f}_{sh},\bar{\rho}^{f}_{sh}\in\mathbb{R}^{E}_{\geq 0}$ |  Lower-Level Dual Variables of power flow constraints (1000$/MWh)
$\underline{\rho}^{\theta}_{sh},\bar{\rho}^{\theta}_{sh}\in\mathbb{R}^{E}_{\geq 0}$ |  Lower-Level Dual Variables of voltage angle constraints (1000$/rad)
$\delta_{sh}\in\mathbb{R}$ |  Lower-Level Dual Variable of reference angle constraint (1000$/rad)
$C:\mathcal{Y}\rightarrow\mathbb{R}$ |  Function of power system operating cost (1000$)
$\tilde{C}:\mathcal{Y}\times\prod\limits_{sh}\mathbb{R}^{G}\times\prod\limits_{sh}\mathbb{R}^{E}\times\prod\limits_{sh}\mathbb{R}^{E}\rightarrow\mathbb{R}$ |  Function of cyberattacker targeted disruption
Table 3:  Nomenclature.