#  Cyber Insurance Against Cyberattacks
on Electric Vehicle Charging Stations

Samrat Acharya, _Student Member, IEEE_ , Robert Mieth, _Member, IEEE_ ,
Charalambos Konstantinou, _Senior Member, IEEE_ , Ramesh Karri, _Fellow, IEEE_
, and Yury Dvorkin _Member, IEEE_

###### Abstract

Even with state-of-the-art defense mechanisms, cyberattacks in the electric
power distribution sector are commonplace. Particularly alarming are load-
altering (demand-side) cyberattacks launched through high-wattage assets,
which are not continuously monitored by electric power utilities. Electric
Vehicle Charging Stations (EVCSs) are among such high-wattage assets and,
therefore, cyber insurance can be an effective mechanism to protect EVCSs from
economic losses caused by cyberattacks. This paper presents a data-driven
cyber insurance design model for public EVCSs. Under some mildly restrictive
assumptions, we derive an optimal cyber insurance premium. Then, we robustify
this optimal premium against uncertainty in data and investigate the risk of
insuring the EVCSs using Conditional Value-at-Risk. A case study with data
from EVCSs in Manhattan, New York illustrates our results.

## Nomenclature

Sets and Indexes

$\mathcal{S}$

Set of typical days, indexed by $s$.

$\mathcal{T}$

Set of hours in a day, indexed by $t$.

$\mathcal{B}$

Set of buses, indexed by $b$.

$\mathcal{I}$

Set of generators, indexed by $i$.

$\mathcal{I}_{b}$

Set of generators connected to bus $b$, $\mathcal{I}_{b}\in\mathcal{I}$.

$\mathcal{L}$

Set of power lines, indexed by $l$.

l(o), l(r)

Indices of sending, receiving buses of line $l$.

Parameters

$C_{i}^{g}$

Energy price offered by generator $i$, $\$/MWh$.

$\overline{G}_{i}$

Maximum power output by generator $i$, $MW$.

$\overline{F}_{l}$

Flow limit on power line $l$, $MW$.

$z_{l}$

Reactance of line $l$, $\mathrm{\SIUnitSymbolOhm}$.

$d_{t}^{s}$

EVCS charging demand at time $t$ for typical day $s$, $MW$.

$d^{d,s}_{b,t}$

Nodal power demand at bus $b$ at time $t$ for typical day $s$, $MW$.

$r$

Insurance profit loading factor.

$\phi^{s}$

Likelihood of typical day $s$.

$\gamma$

EVCS cyberattack risk sharing factor.

$\kappa$

Coefficient of increment in insurance premium per cyberattack on EVCS.

$\mathbb{P}(A)$

Probability of a cyberattack on EVCS during the insurance period.

$A_{h}$

Number of cyberattacks on EVCS in the past.

$\rho$

Penalty to EVCS due to cyberattack, $\$/kW$.

Variables

x

Insurance premium for insurance period.

$\lambda_{t}^{c}$

EV charging price offered by an EVCS at time $t$, $\$/kWh$ .

$\lambda_{t,b}^{u}$

Electricity price offered by a power utility at time $t$ at bus $b$, $\$/kWh$
.

$g_{i,t}^{s}$

Power output of generator $i$ at time $t$ on typical day $s$, $MW$.

$\theta_{l,t}^{s}$

Voltage phase angle at line end $l$ at time $t$ on typical day $s$, $rad$.

$f_{l,t}^{s}$

Power flow in line $l$ at time $t$ on typical day $s$, $MW$.

##  I Introduction

Smart grids have become increasingly prone to cyber intrusions. For example,
in 2018, the U.S. Department of Homeland Security issued 1211% more security
threats for entities in the U.S. electricity sector relative to 2010 [1].
Similarly, IBM reported a 13% increase in financial losses from cyber breaches
in the electricity sector from 2019 to 2020 [2]. While there are several
initiatives to mitigate such threats, their complete elimination is unlikely
due to the complex cybersecurity terrain shaped by both known and undiscovered
vulnerabilities. In particular, concerns arise in case of so-called load-
altering (demand-side) cyberattacks, which leverage countless access points
enabled by demand-side Information and Communication Technologies (ICT) to
compromise operations of high-wattage appliances and processes (e.g. Electric
Vehicles (EVs), demand-side management and microgrid solutions, grid-
interactive building systems, etc.) [3, 4, 5, 6, 7, 8]. These appliances are
not directly monitored by the utilities, which obstructs continuous validation
of their trustworthiness. As a result, there is no means to ensure robustness
of these appliances against cyberattacks and to promote cybersecurity
compliance and hygiene among electricity consumers. To bridge this gap, this
paper designs a cyber insurance mechanism that promotes cybersecurity in smart
grids and hedges financial risks imposed by load-altering attacks. The
proposed mechanism is compatible with available real-world data and we
demonstrate its effectiveness for Electric Vehicle Charging Stations (EVCSs)
using data from Manhattan, New York (NY).

Cyber insurances are designed to cover losses and liabilities caused by cyber
incidents such as data breaches, business interruption, health repercussions
and equipment damage [9]. The global cyber insurance market is expected to
grow from $4.8 billion in 2018 to $28.6 billion in 2026, with a compound
average annual growth rate of 24.9% [10]. At the same time, the contribution
of the energy sector to this market remains thin relative to the banking and
financial service, Information Technology (IT), and healthcare sectors.
However, with the proliferation of high-impact cyberattacks on the energy
sector, such as the Colonial Pipeline attack in May 2021, cyber insurances
have gained prominent attention [11]. Nevertheless, present risk modeling and
insurance evaluation techniques treat the electricity sector similarly to the
IT sector and do not account for physical aspects of electricity delivery
[10]. Such considerations may lead to erroneous insurance policies and reduce
viability of these mechanisms. In fact, a lack of domain-specific risk
modeling techniques has already affected the bottom line of the cyber
insurance industry. For example, the loss ratio of the US cyber insurers
increased from 47% in 2019 to 73% in 2020 [12].

Following the growing demand for cyber insurances, several recent studies
analyzed methods for designing cyber insurance premiums [13, 14]. Romanosky et
al. [13] studied 235 cyber insurance policies underwritten between 2007 and
2017 in New York, Pennsylvania, and California and outlined core insurance
design principles and the range of their design parameters (e.g. how to
characterize cyber incidents, differentiate risks between customers, and price
these risks). We detail these foundations in Section II. Woods et al. [14]
studied real-world data from 26 cyber insurers in California and summarized
premium differences across various insurance characteristics (e.g. coverage
ratio, risk exposure, revenue of the insured entity). The results from [14]
demonstrate that risk exposure from cyber incidents are best modeled using
Gamma, Log-Normal, Polynomial, Pareto, and Weibull distributions, which
departs from typical exponential distribution assumptions used to deal with
adversarial uncertainties in smart grids (e.g. unplanned outages [15]). Taking
a more theoretical perspective, Zhang et al. [16] proposed a non-cooperative
game between the cyber insurer and the insured entity to design the optimal
premium, where the upper-level models an insurer and the lower-level models an
insured entity, which play a zero-sum game with an attacker. The important
outcome of [16] is a demonstration that cyber insurance products encourage
insured entities to adopt security mechanisms that will potentially reduce the
success rate of cyberattacks and subsequently insurance premiums. The studies
in [13, 14, 16] consider generic cyberattack risks, which do not leverage a
domain-specific context. By contrast, this paper focuses specifically on
designing a data-driven optimization for computing cyber insurance premiums
informed by domain-specific insights from real-world operations of electric
power distribution systems with EVCSs, which are shown as a potent load-
altering attack vector [17, 18, 19, 3, 20].

Recent studies in [21, 22, 23] have considered cyber insurance mechanisms for
power grids, and [24] specifically focuses on EVs. In [21], Maynard and
Beecroft qualitatively studied the need for cyber insurances in the U.S. power
grids using the Erebos trojan attack, which can incur financial loss up to
$223.8 billion. Lau et al. in [22] demonstrated that insurance premiums for
transmission companies depend on the allocation of defense resources across
the network and can be reduced up to $\approx 50\%$, if optimized with domain-
specific insights from transmission network operations. In [22], power grid
cyberattacks are modeled using a Semi-Markov Process (SMP), which makes it
possible to realistically capture stochastic transitions between “nominal” and
“attacked” system states. The case study uses the synthetic IEEE Reliability
Test System (RTS)-96 test bed. Similarly, using another artificial IEEE
benchmark, Yang et al. demonstrated in [23] the effects of generator ramping
rates, bus loading, and intrusion probabilities on cyber insurance premiums.
The results in [23] suggest that premiums increase as the likelihood of
intrusion, which nonlinearly depends on bus loading, increases. Honag et al.
in [24] demonstrated that a cyber insured EV minimizes its gross charging
cost, even under Denial-of-Service (DoS) attacks during charging. However, the
studies in [21, 22, 23, 24] are limited in two notable aspects. First, they
assume that the insurer knows insurance policy parameters, e.g. the
probability of an attack on the insured entity, with perfect accuracy. Second,
attack scenarios (e.g. Erebos trojan in [21]) or attack implementations deal
with artificial data sets, which are customized to show profound attack
effects.

Departing from [16, 13, 14, 21, 22, 23, 24], this paper develops a cyber
insurance mechanism that focuses on a specific and realistic load-altering
attack vector stemming from electric vehicle charging and uses real-world
charging and power grid data from Manhattan, NY. The paper makes the following
contributions:

  1. 1.

We develop an optimization framework to design an insurance premium with (i) a
predetermined (regulated) electricity tariff which mirrors the current
practice, and (ii) a dynamic electricity tariff settled using the Distribution
Locational Marginal Price (DLMP) framework [25], which aims to represent
electric power distribution of the future. The first case uses a Bi-Level (BL)
game between the insurer and the EVCS, while the second one uses a Tri-Level
(TL) game including the insurer, the EVCS, and the distribution system.

  2. 2.

In the case of the predetermined tariff, we derive an analytical solution
determining the optimal insurance premium and analyze its dependence on
external factors, i.e. EVCS attack probability, EVCS risk-attitude, and
insurance profit loading factor. Under the DLMP tariff, we obtain a numerical
solution using a computationally tractable column-and-cut (C&CG) algorithm.

  3. 3.

In both cases, we employ a SMP informed by real-world data to estimate the
probability of the attack on EVCSs. This allows for an accurate
parameterization of unknown external factors that inform insurance decisions.
Further, the premium optimization is robustified using Conditional-Value-At-
Risk (CVaR), which internalizes a risk attitude of the EVCS.

  4. 4.

Finally, we use real-world data from Manhattan, NY to evaluate the role of
cyber insurances on the financial impact of cyberattacks.

##  II Cyber Insurance Design Principles

This section describes foundational principles of cyber insurance design as
applicable for EVCSs and electric power distribution systems. Consider an EVCS
acquiring a cyber insurance for a given time period. The insurer collects a
monthly premium (payment) from the EVCS. In return, the insurer fully or
partially covers the EVCS financial loss caused by a cyberattack on the EVCS.
This loss can be first-party and third-party. First-party losses refer to the
loss directly incurred by the insured entity, e.g. from damage caused by
cyberattacks on physical and cyber assets of the EVCS and/or a service
interruption. Third-party losses, or liability, refer to payouts and legal
costs that the insured entity needs to pay if it is found legally liable. For
example, an EVCS may be liable for damages or injuries of its employees,
customers, EVs and the power grid caused by a cyberattack.

One of the key challenges of the optimal cyber insurance design is to assess
and evaluate uncertain first- and third-party losses with high accuracy. Both
the insurer and insured entity benefit from increasing the accuracy of loss
assessments and, therefore, often enter proprietary, non-public agreements
defining various aspects of loss coverage. As a result of these proprietary
agreements, they also often obfuscate their specific methods for assessing
attack risks and optimal insurance premiums. However, some regulatory
authorities require that insurers submit their policies for review by
authorities [13]. For example, U.S. insurers file their issued policies with
the National Association of Insurance Commissioners (NAIC), a non-for-profit
organization supervising insurance policies. These authorities check if the
insurance policies are: (i) affordable to the insured entity, (ii) adequate
for the insurers to run their business, and (iii) uniformly priced across
entities with similar risk profiles (e.g. financial institutions face similar
cyberattack risks). Following [13], we summarize the following types of
insurance policies:

  1. 1.

_Flat rate_ : The insurer estimates the cyberattack probability of the insured
entity and the expected financial loss for a given time period. Then, the flat
rate premium is computed to trade-off desired profit and risk exposure levels
of the insurer. This premium policy, for instance, can be applied uniformly
for similar small businesses.

  2. 2.

_Flat rate with hazard groups_ : In addition to the flat rate policy above,
this policy labels cyberattack risks of the insured entities into low, medium,
and high risks. Such risks are qualitatively analyzed based on the sensitivity
of the cyber space of the insured entities, i.e. the volume of their networks
and connected devices. Then, each hazard group is assigned its own flat-rate
policy.

  3. 3.

_Base rate_ : In addition to the cybersecurity factors described above, this
policy also accounts for the asset value or revenue of the insured entity and
industry-dependent policy factors (e.g. cyberattack frequency and historical
insurance claims). Furthermore, the base rate calculation often incorporates
the use of the security mechanisms of the insured entity to reduce cyberattack
risks (e.g. multi-factor authentications and data encryption).

In this paper, we use the base rate premium design principle for assessing the
cyber insurance of the EVCSs. In contrast to the two flat-rate policies, the
base-rate policy allows for the insurer to incorporate the risk-attitude of
the insured entity (e.g. the portion of loss coverage, attack history, and
defense measures) and design a data-rich insurance policy.

##  III Insurance Design Models

###  III-A Predetermined Tariff

(a) (b)

Figure 1: Structure of optimization models for EVCS cyber insurance: (a) Bi-
Level (BL) for the predetermined electricity tariff and (b) Tri-Level (TL) for
the DLMP based electricity tariff.

Currently, many EVCSs operate in power grids with a predetermined (regulated)
electricity tariff. For example, Con Edison of New York established a
customized volumetric tariff for EV charging [26]. Under this tariff,
interactions between the insurer and EVCS can be modeled using a non-
cooperative game structure in Fig. 1(a), where the Upper-Level (UL) models the
insurer and the Lower-Level (LL) is the insured EVCS. Mathematically this game
can be expressed as the following BL optimization:

$\displaystyle\min_{\Xi^{I}}$ $\displaystyle~{}x$ (1a) $\displaystyle s.t.$
$\displaystyle~{}CL(d,\lambda_{t}^{c,*})=\underbrace{\frac{\gamma(1+\kappa
A_{h})\mathbb{P}(A)}{1-r}}_{\text{Policy factors
$\Pi$}}\underbrace{\sum_{s=1}^{S}\phi^{s}\sum_{t=1}^{T}d_{t}^{s}\lambda_{t}^{c,*}}_{\text{EVCS
revenue}},$ (1b) $\displaystyle x\geq CL,$ (1c)
$\displaystyle\boldsymbol{\lambda^{c,*}}=\operatorname*{arg\,min}_{\Xi^{E}}~{}\\{||{\boldsymbol{\lambda^{c}}}||_{2}^{2}:\quad\boldsymbol{\lambda^{c}}=[\lambda_{t}^{c},\ldots,\lambda_{T}^{c}]^{\top}\\},$
(2a) $\displaystyle
s.t.~{}C(\lambda^{c},\lambda_{t}^{u*},\hat{x}^{*})\\!=\\!\overbrace{(1-\mathbb{P}(A))\\!\sum_{s=1}^{S}\phi^{s}\\!\sum_{t=1}^{T}d_{t}^{s}(\lambda_{t}^{u,*}-\lambda_{t}^{c})}^{\text{Net
cost during unattacked time}}$
$\displaystyle+\underbrace{\mathbb{P}(A)\sum_{s=1}^{S}\phi^{s}\sum_{t=1}^{T}(d_{t}^{s}(\rho-\gamma\lambda_{t}^{c}))}_{\text{Net
cost during
attack}}+\underbrace{\sum_{s=1}^{S}\phi^{s}\sum_{t=1}^{T}\hat{x}^{*}d_{t}^{s}}_{\text{Cost
of insurance}}\leq 0.$ (2b)

where the insurer and EVCS are modeled in Eqs. (1a)-(1c) and Eqs. (2a)-(2b)
with decisions sets ${\Xi^{I}}=\\{x\geq 0\\}$ and
${\Xi^{E}}=\\{\lambda_{t}^{c}\geq 0\\}$. The insurer in Eq. (1) seeks to
determine the lowest insurance premium $(x)$ for a given period that recovers
the expected cyber loss $(CL)$ caused by a cyberattack on the EVCS. The
expected loss in Eq. (1b) is computed as the product of the expected EVCS
revenue and industry-dependent policy factors
$\Pi=\\{\mathbb{P}(A),r,\gamma,\kappa,A_{h}\\}$. The expected revenue is
computed over a given set of typical operational days (indexed by
$s\in\mathcal{S}$ and with the likelihood of $\phi^{s}$), which are provided
by the EVCS to the insurer. Policy factors $\Pi$ are computed by the insurer
and include an attack probability during the insured period $(\mathbb{P}(A))$,
a profit loading factor $(r)$, a risk sharing factor $(\gamma\in[0,1])$, a
premium increment factor for historical cyber incidents $(\kappa)$, and a
count of historical cyber incidents on the EVCS $(A_{h})$. Parameter
$\mathbb{P}(A)$ captures the cyberattack risk of an EVCS. Accordingly, a
greater value of $r$ implies a greater confidence level in estimating the
cyber losses and, thus, potentially, lower profit margins of the insurer, and
vice versa. The risk sharing factor makes it possible to arbitrage risk
sharing between the insurer and EVCS, where $\gamma=1$ indicates a complete
transfer of the cyber risk to the insurer and $\gamma=0$ indicates no risk
sharing, i.e., no insurance. Parameter $\kappa\geq 0$ imposes a penalty for
$A_{h}$, while $\kappa=0$ indicates no marginal increment for $A_{h}$.

In turn, the EVCS in Eq. (2) minimizes the hourly EV charging rate
$(\lambda_{t}^{c})$ for $t\in\mathcal{T}$ hours that recovers the expected
EVCS operating cost $(C)$ during the insured period. The expected operating
cost in Eq. (2b) includes three components: (i) net cost during the unattacked
time (probability is $1-\mathbb{P}(A)$), (ii) net cost during attacks
(probability is $\mathbb{P}(A)$), and (iii) insurance premium cost. During the
unattacked period (i.e., the first term in Eq. (2b)), the EVCS operates
normally. Thus, the net EVCS cost is calculated as the EVCS revenue collected
from charging EVs minus the the electricity cost paid to the power utility
based on the electricity tariff $(\lambda_{t}^{u,*})$. During the attack on
EVCS (i.e., the second term in Eq. (2b)), the EVCS does not operate. Thus, the
EVCS not only loses revenue from EV charging during the attack period, but
also loses the confidence of EV drivers to use the attacked EVCS in the
future. These losses are captured by
$\mathbb{P}(A)\sum_{s}\phi^{s}\sum_{t}d_{t}^{s}\rho$, where $\rho(\$/kW)$ is
the penalty incurred by the attack. However, because the EVCS is insured, it
is entitled to receive the insurance payment equal to $\gamma$ times the lost
revenue. Thus, the net cost under attack is the loss in revenue and goodwill
minus the insurance claim. Irrespective of the state of the attack on EVCS,
the EVCS pays the insurance premium, as shown in the third term of Eq. (2b).
Notably, we denote $\sum_{s}\phi^{s}\sum_{t}\hat{x}^{*}d_{t}^{s}=x^{*}$ in Eq.
(2b). All three terms of Eq. (2b) are stochastic over $s\in\mathcal{S}$
typical days with likelihood $\phi^{s}$. Next, using Eqs. (1a)-(2b), we state
and prove:

######  Theorem 1.

Consider the model in Eqs. (1)-(2) and let $\lambda_{t}^{u,*},~{}\forall
t\in\mathcal{T},$ be a predetermined electricity tariff for the insurance
period (e.g. [26]). Then, the optimal cyber insurance premium follows as:

$\displaystyle x^{*}\\!\\!=$
$\displaystyle\frac{C}{[(\mathbb{P}(A)(\\!\gamma\\!-\\!1\\!)\\!+\\!1\\!)\\!-\\!C]}\\!\\!\bigg{[}\\!\mathbb{P}(A)\rho\\!\\!\sum_{t}\\!\\!D_{t}\\!\\!+\\!\\!(1\\!-\\!\mathbb{P}(A))\\!\\!\sum\\!\\!D_{t}\lambda_{t}^{u,\\!*}\\!\bigg{]},$
(3)

where parameters $C=\frac{\mathbb{P}(A)\gamma(1+\kappa A_{h})}{1-p}$ and
$D_{t}=\sum_{s}\phi^{s}d_{t}^{s}$.

Proof: As per Eq. (2b), the EVCS breaks even when $\lambda_{t}^{c}$ is such
that Eq. (2b) is binding. Thus, the Lagrangian function of the LL problem is:

$\displaystyle\mathcal{L}=$
$\displaystyle\sum_{t}(\lambda_{t}^{c})^{2}+\omega\bigg{[}(1-\mathbb{P}(A))\sum_{t=1}^{T}\bigg{(}D_{t}(\lambda_{t}^{u,*}-\lambda_{t}^{c})\bigg{)}$
$\displaystyle+\sum_{t}^{T}(\hat{x}^{*}D_{t})+\mathbb{P}(A)\sum_{t=1}^{T}\bigg{(}D_{t}(\rho-\gamma\lambda_{t}^{c}\bigg{)}\bigg{]},$
(4)

where $\omega$ is a Lagrangian dual variable associated with the binding
constraint in Eq. (2b). Next, the first-order optimally conditions with
respect to $\lambda_{t}^{c},\forall t\in\mathcal{T},$ is:

$\displaystyle\frac{\partial{L}}{\partial{\lambda^{c}_{t}}}$
$\displaystyle=2\lambda_{t}^{c}-\omega\bigg{(}\mathbb{P}(A)\gamma
D_{t}+(1-\mathbb{P}(A))D_{t}\bigg{)}=0.$ (5)

Solving Eq. (5) for $\lambda^{c}_{t}$ leads to:

$\displaystyle\lambda^{c,*}_{t}=\frac{1}{2}\omega\bigg{[}(\mathbb{P}(A)(\gamma-1)+1)D_{t}\bigg{]},\quad\forall
t\in 1:T.$ (6)

In turn, $\omega$ can be expressed from setting
$\frac{\partial{L}}{\partial{w}}=0$ and comparing the result with Eq. (6),
which returns:

$\displaystyle\omega=$
$\displaystyle\frac{2}{\sum_{t}D_{t}^{2}[\mathbb{P}(A)(\gamma-1)+1]^{2}}\bigg{[}\bigg{(}\mathbb{P}(A)\rho+\lambda_{t}^{u,*}(1-\mathbb{P}(A))$
$\displaystyle+\hat{x}^{*}\bigg{)}\sum_{t}D_{t}\bigg{]}.$ (7)

Next, the expression for $\lambda_{t}^{c,*}$ in Eq. (6) is used to solve the
UL optimization in Eq. (1). When solving the optimization in Eq. (1), Eq. (1c)
becomes binding to attain the least-cost objective value and, thus, the
Lagrangian function of the LL problem in Eq. (1) follows as:

$\displaystyle\mathcal{L}$
$\displaystyle=x+\delta\bigg{[}C\sum_{t}D_{t}\lambda_{t}^{c,*}-x\bigg{]},$ (8)

where $\delta$ is a Lagrangian dual variable of Eq. (1c). Using
$\lambda_{t}^{c*}$ in Eq. (6) and then Eq. (6) in Eq. (8), the first-order
optimality condition with respect to $\delta$ is:

$\displaystyle\frac{\partial{L}}{\partial{\delta}}=$
$\displaystyle\frac{1}{2}\omega
C\sum_{t}\bigg{[}D_{t}\bigg{(}(\mathbb{P}(A)(\gamma-1)+1)D_{t}\bigg{)}\bigg{]}-x=0.$
(9)

Solving Eq. (9) for $x^{*}$ gives Eq. (3). ∎

Remark 1 (The role of EVCS demand forecast error): The cyber insurance premium
in Eq. (3), which is obtained from Eqs. (1)-(2), is directly proportional to
the expected EVCS demand forecast, i.e., $\sum_{s}\phi^{s}d_{t}^{s}=D_{t}$.
Thus, the EVCS operator can reduce its premium by strategically reporting a
lower value of the EVCS demand forecast. However, we assume this premium
reduction is very unlikely because the commercial EVCSs report their demand
forecast to multiple entities (e.g. power grid operator, insurer, public
users), and hence, the insurer may cross-validate the reported forecasts.

###  III-B DLMP-based Tariff

As opposed to the predetermined electricity tariff underlying the analytical
solution in Theorem 1, electric power distribution systems are expected to
gradually transition towards dynamic and more spatially and temporally
granular tariffs (e.g. by means of the DLMPs [25]). To account for this
eventuality, the optimization in Eqs. (1a)-(2b) is recast as the TL
optimization model shown in Fig. 1(b). In this TL optimization, the UL and
Middle-Level (ML) problems are the insurer and EVCS with the roles explained
above, while the LL is a model of the distribution system operator used to
compute DLMPs. Mathematically, the UL and ML optimizations are given by Eq.
(10a) and Eq. (10b), respectively, while the LL optimization is given in Eqs.
(10c)-(10g). The resulting TL optimization follows then as:

$\displaystyle\text{Eqs.~{}
\eqref{eq:objective_upper_bi}-\eqref{eq:premium_bound}}\hskip
36.98857pt:\text{UL problem}$ (10a) $\displaystyle\text{Eqs.~{}
\eqref{eq:objective_lower_bi}-\eqref{eq:EVCS_cost}}\hskip 36.98857pt:\text{ML
problem}$ (10b)
$\displaystyle\boldsymbol{\lambda^{u,*}}=\operatorname*{arg\,min}_{\Xi^{G}}~{}\\{C^{s}_{LL}:=\sum_{t}\sum_{i}C_{i}^{g}g_{i,t}^{s},~{}\forall
s\in\mathcal{S}\\},$ (10c) $\displaystyle\quad s.t.~{}0\leq
g_{i,t}^{s}\leq\overline{G}_{i}:(\underline{\alpha}_{i,t}^{s},\overline{\alpha}_{i,t}^{s}),~{}\forall
s\in\mathcal{S},~{}\forall t\in\mathcal{T},~{}\forall i\in\mathcal{I},$ (10d)
$\displaystyle\quad
f_{l,t}^{s}z_{l}=\theta_{o(l),t}^{s}-\theta_{r(l),t}^{s}:(\xi_{l,t}^{s}),~{}\forall
s\in\mathcal{S},~{}\forall t\in\mathcal{T},~{}\forall l\in\mathcal{L},$ (10e)
$\displaystyle\quad-\overline{F}_{l}\leq
f_{l,t}^{s}\leq\overline{F}_{l}:(\underline{\delta}_{l,t}^{s},\overline{\delta}_{l,t}^{s}),\forall
s\in\mathcal{S},~{}\forall t\in\mathcal{T},~{}\forall l\in\mathcal{L},$ (10f)
$\displaystyle\quad d_{b,t}^{d,s}=\sum_{i\in
I_{b}}g_{i,t}^{s}+\sum_{l|r(l)=b}f_{l,t}^{s}-\sum_{l|o(l)=b}f_{l,t}^{s}-d_{b,t}^{s}:(\lambda_{b,t}^{s}),$
$\displaystyle\hskip 39.83385pt~{}\forall s\in\mathcal{S},~{}\forall
b\in\mathcal{B},~{}\forall t\in\mathcal{T},$ (10g)

where ${\Xi^{G}}=\\{g_{i,t}^{s},\theta_{n,t}^{s},f_{l,t}^{s}\geq 0\\}$ and
$\boldsymbol{\lambda^{u}}=[\lambda_{b,t}^{s},\ldots,\lambda_{b,T}^{s}]^{\top}$.
The objective of the LL optimization in Eq. (10c) is to minimize the expected
generation cost $(C_{LL})$ in the electric power distribution system. Eq.
(10d) enforces power generation limits on the generators, Eq. (10e) models DC
power flows, Eq. (10f) enforces flow capacity limits, and Eq. (10g) ensures
the nodal power balance. The dual variable of Eq. (10g), $\lambda_{b,t}^{s}$,
is adopted as a DLMP at bus $b$ for time $t$ and day $s$. Note that we use DC
power flow as a proxy to compute DLMPs, but other power flow formulations,
e.g. LinDistFlow [27], can also be adopted to compute additional DLMP
components.

The TL optimization in Eq. (10a)-(10g) does not have an analytical solution
and can only be solved numerically, e.g. by means of the C&CG algorithm, [28,
29], described in Section V-B. However, the use of this algorithm requires a
reformulation of the TL optimization as an equivalent BL problem. Since the LL
problem is convex and linear, the strong duality condition holds and the
single-level equivalent of the ML and LL can be obtained using the Dual LL
(DLL) formulation given as ($\forall s\in\mathcal{S}$):

$\displaystyle\min_{\Xi^{DLL}}C_{DLL}^{s}:=\sum_{t\in\mathcal{T}}\bigg{(}\sum_{i\in\mathcal{I}}\overline{\alpha}_{i,t}^{s}\overline{G}_{i}+\sum_{i\in\mathcal{L}}(\overline{\delta}_{l,t}^{s}+\underline{\delta}_{l,t}^{s})\overline{F}_{l}$
$\displaystyle\hskip
56.9055pt-\sum_{b\in\mathcal{B}}(d_{b,t}^{d,s}+d_{b,t}^{s})\lambda_{b,t}^{s}\bigg{)}$
(11a) $\displaystyle
s.t.\quad\overline{\alpha}_{i,t}^{s}-\underline{\alpha}_{i,t}^{s}-\lambda_{b(i),t}^{s}=-C_{i}^{g},\quad\forall
i\in\mathcal{I},~{}\forall t\in\mathcal{T},$ (11b)
$\displaystyle\xi_{l,t}^{s}+\overline{\delta}_{l,t}^{s}-\underline{\delta}_{l,t}^{s}+\lambda_{o(l),t}^{s}-\lambda_{r(l),t}^{s}=0,~{}\forall
l\in\mathcal{L},~{}\forall t\in\mathcal{T},$ (11c) $\displaystyle\sum_{l\mid
r(l)=b}\frac{\xi_{l,t}^{s}}{z_{l}}-\sum_{l\mid
o(l)=b}\frac{\xi_{l,t}^{s}}{z_{l}}=0,\quad\forall l\in\mathcal{L},~{}\forall
t\in\mathcal{T},$ (11d)

where
${\Xi}^{DLL}=\\{\overline{\alpha}_{i,t}^{s},\underline{\alpha}_{i,t}^{s},\overline{\delta}_{l,t}^{s},\underline{\delta}_{l,t}^{s}\geq
0,\xi_{l,t}^{s},\lambda_{b,t}^{s}:\text{free}\\}$. Given Eq. (11)-(11d), the
TL model is equivalent to:

$\displaystyle\text{Eq.~{}\eqref{eq:objective_upper_bi}}\hskip
66.86397pt:\text{UL problem}$ (12a) $\displaystyle\text{s.t~{}
Eqs.~{}\eqref{eq:loss_revenue}-\eqref{eq:premium_bound}},$
$\displaystyle\lambda_{t}^{c}\in\operatorname*{arg\,min}\\{$ (12b)
$\displaystyle\text{Eqs.~{}\eqref{eq:objective_lower_bi}-\eqref{eq:EVCS_cost}},\hskip
41.25648pt:\text{ML problem}$ (12c)
$\displaystyle\text{Eqs.~{}\eqref{eq:gen_const}-\eqref{eq:load_const},}\hskip
29.87538pt:\text{LL primal constraints}$ (12d)
$\displaystyle\text{Eqs.~{}\eqref{eq:dll_constraint1}-\eqref{eq:dll_constraint3}},\hskip
31.29802pt:\text{LL dual constraints}$ (12e) $\displaystyle
O_{LL}^{s}=O_{DLL}^{s},~{}\forall s\in\mathcal{S}\hskip 5.69054pt:\text{LL
strong duality}\\}.$ (12f)

The BL equivalent in Eqs. (12a)-(12f) is then solved using the C&CG algorithm
as explained in Section V-B.

##  IV Estimating $\mathbb{P}(A)$ using a Semi-Markov Process

GIDFC$H_{GI}$$H_{ID}$$H_{IF}$$H_{DC}$$H_{CG}$$H_{FG}$ Figure 2: Semi-Markov
Process (SMP) for cyberattacks on EVCSs. $H_{(\cdot)}$ defines the Cumulative
Distribution Function (CDF) of the state transition time.

The optimization problems in Eqs. (1a)-(2b) and Eqs. (12a)-(12f) both require
an estimation of insurance policy factors including $\mathbb{P}(A)$ as
explained in Section 2. These factors can be parameterized using the current
insurance industry practice and using the SMP, which accommodates non-
exponential distributions of Markovian state transitions. We consider five
states that represent cyberattacks on the EVCS and its recovery [30, 22]: Good
(G), Intrusion (I), Detection (D), Containment (C), and Failure (F). During
normal operations, the EVCS is in state $G$. State $I$ indicates that the EVCS
and/or its server are infected with malware (see our review of attack routes
in [3]). After the infection, the attack is either detected in state $D$,
which returns the EVCS to state $G$ via state $C$, or remains undetected,
which leads to the failure of the EVCS operation given by state $F$. Once the
EVCS detects the attack, the EVCS contains the cause of the attack to avoid
its spread at state $C$ before returning the EVCS to state $G$. Following the
success of the attack at state $F$, the EVCS recovers to state $G$.

Let $\mathcal{S}=\\{G,I,D,C,F\\}$ be a state space of the Markov chain
representing the effect of cyberattacks on the EVCS and let
$X_{n}\in\mathcal{S},~{}n\in\mathcal{N}=\\{0,1,2,\ldots\\}$ denote the states
of the cyberattack after $n$ transitions in the Markov chain. The state
transitions occur at jump times $T_{n}$. Thus, the SMP is a sequence of two-
dimensional random variables $(X_{n},T_{n})$, where the attack states
$G,I,D,C,\text{and}~{}F$ are stochastic over time. The SMP evolves over
continuous time $t\in[T_{n},T_{n+1})$ and is formalized as:

$\displaystyle\mathbb{P}(X_{n+1}=j,\tau_{n+1}\leq
t\mid(X_{0},T_{0}),\ldots,(X_{n}=i,T_{n}))$
$\displaystyle=\mathbb{P}(X_{n+1}=j,\tau_{n+1}\leq t\mid X_{n}=i)$ (13)
$\displaystyle=\mathbb{P}_{ij}H_{ij}(t),\quad\forall
n\in\mathcal{N},~{}\forall i,j\in\mathcal{S},$

where $\tau_{n+1}=T_{n+1}-T_{n}$ is the sojourn time, i.e., the time spent in
a given state between transitions $n$ and $n+1$, $\mathbb{P}_{ij}$ is the
probability of transitioning from state $i$ to $j$, and $H_{ij}(t)$ is the
Cumulative Distribution Function (CDF) of the transition time between states
$i$ and $j$. Below we provide a two-stage method to solve the SMP based on
[31, 30], where the first stage returns the transition probability of the
embedded Markov chain of the SMP and the second stage returns the steady-state
probabilities of the SMP states. Notably, the embedded Markov chain ignores
the time spent at a state and only accounts for the transition probability
between the states.

####  IV-1 Stage 1

The kernel matrix denoted as $K(t)$ is given as:

$K(t)=\begin{bmatrix}0&k_{GI}&0&0&0\\\ 0&0&k_{ID}&0&k_{IF}\\\
0&0&0&k_{DC}&0\\\ k_{CG}&0&0&0&0\\\ k_{FG}&0&0&0&0\\\ \end{bmatrix}.$ (14)

Assuming a Weibull distribution of state transitions given as
$H_{ij}(t,\alpha_{ij},\beta_{ij})=1-\exp{-(\frac{t}{\alpha_{ij}})^{\beta_{ij}}}$,
the elements of $K(t)$ in Eq. (14) are determined as:

$\begin{split}&k_{GI}(t)=H_{GI}(t),\quad
k_{ID}(t)=\int_{0}^{t}\bar{H}_{IF}(t)dH_{ID}(t)\\\
&k_{IF}(t)=\int_{0}^{t}\bar{H}_{ID}(t)dH_{IF}(t),\quad k_{DC}(t)=H_{DC}(t),\\\
&k_{CG}(t)=H_{CG}(t),\quad k_{FG}(t)=H_{FG}(t),\end{split}$ (15)

where $\bar{H}_{ij}=1-H_{ij}$. The Weibull parameters of state transitions
(i.e., $\alpha,\beta)$ are obtained by fitting empirical data from historical
cyberattacks (e.g. [32]). Due to the limited data which are publicly available
in regards to cyberattacks on EVCS, we model the distributions using data from
cyberattacks in the ICT industry, which is consistence with the real-world
practice [33, 14].

The steady-state probability matrix of the embedded Markov chain is determined
for $t\to\infty$. Thus, we state the kernel matrix in Eq. (14) as
$K(t\to\infty)=1$, which means that the probability of transitioning between
states approaches 1 in infinite time. Since the embedded Markov chain is
discrete-time, applying the steady-state condition to the kernel matrix, leads
to the vector of steady-state probabilities of the embedded Markov chain given
as $\boldsymbol{p}=[p_{G}~{}p_{I}~{}p_{D}~{}p_{C}~{}p_{F}]$:

$\boldsymbol{p}=\boldsymbol{p}K^{\top}(t\to\infty);\quad\boldsymbol{p}\boldsymbol{e}^{\top}=1,$
(16)

where $\boldsymbol{e}=[1~{}1~{}1~{}1~{}1].$

####  IV-2 Stage 2

This stage calculates sojourn time $T_{s}$ spent in state $s$ before
transitioning to another state:

$\begin{split}&T_{G}=\int_{0}^{\infty}\bar{H}_{GI}(t)dt,\quad
T_{I}=\int_{0}^{\infty}\bar{H}_{ID}(t)\bar{H}_{IF}dt,\\\
&T_{D}=\int_{0}^{\infty}\bar{H}_{DC}(t)dt,\quad
T_{C}=\int_{0}^{\infty}\bar{H}_{CG}(t)dt,\\\
&T_{F}=\int_{0}^{\infty}\bar{H}_{FG}(t)dt.\end{split}$ (17)

Let $\boldsymbol{T}=[T_{G}~{}T_{I}~{}T_{D}~{}T_{C}~{}T_{F}]^{\top}$. Using the
values of $\boldsymbol{p}$ in Eq. (16) and $\boldsymbol{T}$ in Eq. (17), the
steady-state probability of transitioning to failure state $F$, e.g. the
probability of attack success $\mathbb{P}(A)$, is:

$\mathbb{P}(A)=\mathbb{P}_{F}=\frac{p_{F}T_{F}}{\boldsymbol{p}\boldsymbol{T}^{\top}}.$
(18)

Since the estimation of $\mathbb{P}(A)$ invokes dealing with uncertainty,
Section V robustifies the optimizations in Section III against the uncertainty
in $\mathbb{P}(A)$.

##  V Distributionally Robust and Risk- Averse Insurance Premium Design

This section extends the optimizations in Section III to account for the
uncertainty of insurance policy factors $\Pi$ and risk-averse EVCS attitudes.

###  V-A Robust and Risk-Averse Bi-Level Optimization Model

Inspired by the analytical expression in Eq. (3) of Theorem 1, we seek to
robustify this result against uncertain policy factors $\mathbb{P}(A)$, $r$,
and $\kappa$. To this end, we first assume that any parameter estimation error
is contained in a box-constrained confidence interval as practiced by the
industry [13, 12], where the box limits can be elicited using expert
knowledge, industry practices, and/or statistical analysis. Second, we revisit
Eqs. (2a)-(2b) to include a CVaR metric as a means of considering EVCS risk
attitudes.

####  V-A1 Robust and Risk-Averse EVCS Optimization Problem

Parameter $\mathbb{P}(A)$ in the EVCS optimization in Eq. (2) is subject to
uncertainty due to data insufficiency and noise (e.g. see the SMP in Section
IV). Therefore, this uncertainty can be modeled as
$\mathbb{P}(A):=\overline{\mathbb{P}}(A)\pm\epsilon_{p}$, where
$\overline{\mathbb{P}}(A)$ is the estimated value of the steady-state attack
probability in Section IV and $\epsilon_{p}$ is an user-defined estimation
error. Note that the subscript p in $(\cdot)_{p}$ refers to the parameter
$\mathbb{P}(A)$. Accordingly, the box-constrained confidence interval is given
as $\mathbb{P}(A)\in[\underline{\Gamma}_{p},~{}\overline{\Gamma}_{p}]$, where
$\overline{\Gamma}_{p}=\overline{\mathbb{P}}(A)+\epsilon_{p}$ and
$\underline{\Gamma}_{p}=\overline{\mathbb{P}}(A)-\epsilon_{p}$. For example,
assume $\epsilon_{p}$ follows a normal distribution with zero mean
$\mathbb{E}({\epsilon_{p})}=0$ and variance
$\text{Var}(\epsilon_{p})=\sigma_{p}^{2}$. Then, the box constraint can be
defined as :

$\underline{\Gamma}_{p}=\overline{\mathbb{P}}(A)-t_{(1-\xi/2)}\frac{\sigma_{p}}{\sqrt{N}},~{}\overline{\Gamma}_{p}=\overline{\mathbb{P}}(A)+t_{(1-\xi/2)}\frac{\sigma_{p}}{\sqrt{N}},$
(19)

where $\overline{\mathbb{P}}(A)$ is an empirical mean modeled by the
$t$-distribution, $t_{(1-\xi/2)}$ is the ${(1-\xi/2)}$-quantile of the
$t$-distribution computed from $N$ observations of ${\mathbb{P}}(A)$.

Using Eq. (19), the risk-averse formulation of the EVCS optimization problem
in Eq. (2) can be recast in a robust manner as follows:

$\displaystyle\min_{\lambda^{c}}$
$\displaystyle~{}||{\boldsymbol{\lambda^{c}}}||_{2}^{2};\quad\quad\boldsymbol{\lambda^{c}}=[\lambda_{t}^{c},\ldots,\lambda_{T}^{c}]^{\top},$
(20a) $\displaystyle
s.t.~{}\text{CVaR}_{\alpha}(C(d^{s},\phi^{s},\lambda^{c}))\leq 0,$ (20b)

where $\text{CVaR}_{\alpha}$ of the EVCS cost $C(d^{s},\phi^{s},\lambda^{c})$,
with probability distribution $\phi^{s}$, is defined as the mean of
$\alpha-$tail distribution of $C(d^{s},\phi^{s},\lambda^{c})$ [34, 35], that
is:

$\displaystyle\text{CVaR}_{\alpha}(C(d^{s},\phi^{s},\lambda^{c}))=\sup_{\tilde{\phi}\in
Q_{\alpha}}\sum_{s\in S}\tilde{\phi}^{s}\tilde{C}(d^{s},\lambda^{c}),$ (21a)
where: $\displaystyle
Q_{\alpha}=\\{\tilde{\phi}:\sum_{s}\tilde{\phi}^{s}=1,~{}\tilde{\phi}^{s}\leq\frac{1}{\alpha}\phi^{s}\\},$
(21b) $\displaystyle\tilde{C}(d^{s},\lambda^{c})=$
$\displaystyle(1-\mathbb{P}(A))\sum_{t=1}^{T}d_{t}^{s}(\lambda_{t}^{u,*}-\lambda_{t}^{c})+\sum_{t=1}^{T}\hat{x}^{*}d_{t}^{s}$
$\displaystyle+\mathbb{P}(A)\sum_{t=1}^{T}(d_{t}^{s}(\rho-\gamma\lambda_{t}^{c})).$
(21c) Since $\mathbb{P}(A)$ holds a linear relationship with the EVCS cost in
Eq. (21c), the upper bound of $\mathbb{P}(A)$ in Eq. (19) gives the
$\alpha$-worst-case ($\alpha$-WC) estimation of the EVCS cost. Thus, Eq. (21c)
is recast as: $\displaystyle\tilde{C}^{WC}(d^{s},\lambda^{c})=$
$\displaystyle(1-\overline{\Gamma}_{p})\sum_{t=1}^{T}d_{t}^{s}(\lambda_{t}^{u*}-\lambda_{t}^{c})+\sum_{t=1}^{T}\hat{x}^{*}d_{t}^{s}$
$\displaystyle+\overline{\Gamma}_{p}\sum_{t=1}^{T}(d_{t}^{s}(\rho-\gamma\lambda_{t}^{c})).$
(21d)

Thus, using Eq. (21d) in the EVCS optimization in Eq. (20), the EVCS operator
designs a minimum $\lambda_{t}^{c,WC},~{}\forall t\in\mathcal{T}$, that at
least recovers the $\alpha$-WC EVCS cost. In this setting, the risk-averse
EVCS optimization problem in Eq. (20) follows:

$\displaystyle\min_{\lambda^{c,WC},\zeta,v}$
$\displaystyle~{}||{\boldsymbol{\lambda^{c,WC}}}||_{2}^{2};\quad\quad\boldsymbol{\lambda^{c,WC}}=[\lambda_{t}^{c,WC},\ldots,\lambda_{T}^{c,WC}]^{\top}$
(22a) $\displaystyle s.t.~{}v+\sum_{s}\phi^{s}\zeta_{s}\leq 0,:(\eta)$ (22b)
$\displaystyle\alpha\zeta^{s}\geq\tilde{C}^{WC}(d^{s},\lambda^{c})-v,:(\varphi^{s}),~{}\forall
s\in 1:S,$ (22c) $\displaystyle\zeta^{s}\geq 0,:(\mu^{s})~{}\forall s\in 1:S,$
(22d) $\displaystyle\lambda_{t}^{c,WC}\geq 0,:(\beta_{t})~{}\forall t\in 1:T,$
(22e)

where $v$ is an auxiliary variable used to recast Eq. (21a) as an minimization
problem and $\eta$, $\varphi^{s}$, $\beta_{t}$, and $\mu^{s}$ are the
Lagrangian dual variables of Eq. (22b)-(22e). Next, using the KKT conditions
of the optimization in Eq. (22):

$\displaystyle v+\sum_{s}\phi^{s}\zeta^{s}\leq 0$ (23a)
$\displaystyle-\alpha\zeta^{s}-v+\tilde{C}^{WC}(d^{s},\lambda^{c})\leq 0,\quad
s\in 1:S,$ (23b) $\displaystyle-\zeta^{s}\leq 0,~{}\varphi^{s}\geq 0,\quad
s\in 1:S,$ (23c) $\displaystyle\eta\geq 0,~{}\mu\geq 0,\beta_{t}\geq 0,\quad
t\in 1:T,$ (23d) $\displaystyle\eta(v+\sum_{s}\phi^{s}\zeta^{s})=0,$ (23e)
$\displaystyle\varphi^{s}(-\alpha\zeta^{s}-v+\tilde{C}^{WC}(d^{s},\lambda^{c}))=0,\quad
s\in 1:S,$ (23f) $\displaystyle\mu^{s}(-\zeta^{s})=0,\quad s\in 1:S,$ (23g)
$\displaystyle\eta\phi^{s}-\alpha\varphi^{s}-\mu^{s}=0,\quad s\in 1:S,$ (23h)
$\displaystyle\beta_{t}(-\lambda_{t}^{c,WC})=0,\quad t\in 1:T,$ (23i)
$\displaystyle\eta-\sum_{s}\varphi^{s}=0,$ (23j) $\displaystyle
2\lambda_{t}^{c,WC}-[\overline{\Gamma}_{p}(\gamma-1)+1]\sum_{s}(\varphi^{s}d_{t}^{s})-\beta_{t}=0,\quad
t\in 1:T.$ (23k)

we finally obtain:

$\displaystyle(1-\alpha)\sum_{s}\varphi^{s}=\sum_{s}\mu^{s}.$ (24)

This relationship captures the trade-off solved by the CVaR between adding
scenarios to the set of worst-cases and adjusting their relative probability.
For example, if $\alpha=1$, then Eq. (24) yields $\sum_{s}\mu^{s}=0$, which
implies $\mu^{s}=0,\forall s\in\mathcal{S},$ and $\zeta^{s}\geq 0,\forall
s\in\mathcal{S}$ as per Eq. (23g). In other words, for $\alpha=1$ all
scenarios are considered with their original probability $\phi^{s}$ such that
$\sum_{s}\phi^{s}\zeta_{s}$ recovers EVCS cost in expectation as in Section
III-A. On the other hand, if $\alpha=0$, Eq. (22) only considers the single WC
scenario as $\zeta_{s}=0$ and $v\geq\tilde{C}^{WC}(d^{s},\lambda^{c}),\forall
s\in\mathcal{S}$, as per Eq. (22d). Since $\sum_{s}\varphi_{s}=\eta$ in Eq.
(23j), $\varphi^{s}$ is a proxy for the probability of typical days. Thus, we
define a relative probability of typical day $s$ as
$\hat{\varphi}_{s}=\frac{\varphi_{s}}{\eta},~{}\forall s\in S$. Given
$\hat{\varphi}_{s}$, Eq. (22) determines ${\lambda_{t}}^{c,WC},\forall
t\in\mathcal{T},$ based on the $\alpha$-WC EVCS cost.

####  V-A2 Robust Insurer Optimization Problem

Given the value of $\lambda_{t}^{c,WC}$ and EVCS demand profiles
$d_{t}^{s}~{}\forall s\in\mathcal{S},~{}\forall t\in\mathcal{T}$, we robustify
the insurer optimization by introducing box constraints on uncertain
parameters $\mathbb{P}(A)$, $\kappa$, and $r$. Note that the insurer uses the
box constraint for $\mathbb{P}(A)$ as described in Eq. (19). Consider the box
constraints for $\kappa$ and $r$ as
$[\underline{\kappa},~{}\overline{\kappa}]$ and
$[\underline{r},~{}\overline{r}]$. Since $\kappa$ and $r$ are industry-
specific parameters, the range for the box constraint is also industry-
specific [13]. Accordingly, the insurer optimization problem in Eq. (1) is
recast as:

$\displaystyle\min$ $\displaystyle~{}x^{WC}$ (25a) $\displaystyle s.t.$
$\displaystyle~{}CL=\underbrace{\frac{\gamma(1+\overline{\kappa}A_{h})\overline{\Gamma}_{p}}{1-\overline{r}}}_{\text{policy
factors}}\underbrace{\bigg{(}\sum_{s=1}^{S}\phi^{s}\sum_{t=1}^{T}d_{t}^{s}\lambda_{t}^{c,WC}\bigg{)}}_{\text{EVCS
revenue}}$ (25b) $\displaystyle x\geq CL.$ (25c)

Note that the risk-averse EVCS has $\hat{\varphi}^{s}$ as the relative
probability of typical days. However, the insurer in Eq. (25) may not have the
same information on $\hat{\varphi}^{s}$ as the EVCS (see Remark 1 in Section
III). For example, the insurer optimizes the premium using the expected EV
charging demand profiles during the insured period, but the EVCS optimization
uses WC typical days instead. Therefore, an exact analytical expression for
the insurance premium cannot be obtained in this case, but it can be
determined numerically. Since the risk-averse LL problem in Eq. (22) is
convex, we convert the robust and risk-averse BL optimizations, insurer in Eq.
(25) and EVCS in Eq. (22), to the following single-level equivalent [36]:

$\displaystyle\text{Eq.~{}}\eqref{eq:ul_robust},\quad\text{: UL problem}$
(26a) $\displaystyle\text{Eq.~{}}\eqref{eq:EVCS_CVAR_primal},\quad\text{: LL
primal problem}$ (26b) $\displaystyle\text{Eq.~{}}\eqref{eq:KKT},\quad\text{:
KKT conditions for ML Problem}$ (26c)

The insurance premium for $\alpha$-WC EVCS cost is thus obtained by
numerically solving the single-level equivalent in Eq. (26).

###  V-B Robust and Risk-Averse Tri-Level Optimization Model

This section recasts the original TL optimization in case of the robust
insurer and EVCS presented in Sections V-A1 and V-A2. Thus, the TL
optimization in Section III-B results in:

$\displaystyle\text{Eq.~{}}\eqref{eq:objective_upper_bi_robust}\hskip
62.59605pt:\text{UL problem}$ (27a) $\displaystyle
s.t~{}\text{Eqs.~{}}\eqref{eq:loss_revenue_robust}-\eqref{eq:premium_bound_robust}$
$\displaystyle\lambda_{t}^{c,WC}\in\operatorname*{arg\,min}\\{$ (27b)
$\displaystyle\text{Eqs.~{}}\eqref{eq:EVCS_CVAR_primal_obj}-\eqref{eq:primal_cvar_constraint_last}\hskip
25.60747pt:\text{ML problem}$ (27c)
$\displaystyle\text{Eqs.~{}}\eqref{eq:gen_const}-\eqref{eq:load_const}\hskip
24.18483pt:\text{LL primal constraints}$ (27d)
$\displaystyle\text{Eqs.~{}}\eqref{eq:dll_constraint1}-\eqref{eq:dll_constraint3}\hskip
24.18483pt:\text{LL dual constraints}$ (27e) $\displaystyle
O_{LL}^{s}=O_{DLL}^{s},~{}\forall s\in\mathcal{S}\hskip 5.69054pt:\text{LL
strong duality}\\}.$ (27f)

We use the C&CG algorithm to solve the TL optimization in Eq. (27) as applied
in [28, 29]. The C&CG algorithm is based on the decomposition of the TL
optimization in Eq. (27) into the principal problem and subproblems:

####  V-B1 Principal problem

$\displaystyle\min_{\Xi^{MP}}~{}$ $\displaystyle
x^{WC}+||{\boldsymbol{\lambda^{c,WC}}}||_{2}^{2};$ (28a) $\displaystyle
s.t.~{}\text{Eqs.~{}}\eqref{eq:loss_revenue_robust}-\eqref{eq:premium_bound_robust}\hskip
7.11317pt\text{: UL constraints}$ (28b)
$\displaystyle\text{Eqs.~{}}\eqref{eq:primal_cvar_constraint1}-\eqref{eq:primal_cvar_constraint_last}\hskip
24.18483pt\text{: ML constraints}$ (28c)
$\displaystyle\text{Eqs.~{}}\eqref{eq:gen_const}-\eqref{eq:load_const}\hskip
24.18483pt\text{: LL primal constraints}$ (28d)
$\displaystyle\text{Eq.~{}}\eqref{eq:dll_constraint1}-\eqref{eq:dll_constraint3}\hskip
25.60747pt\text{: LL dual constraints}$ (28e) $\displaystyle
O_{LL}^{s}=O_{DLL}^{s},~{}\forall s\in\mathcal{S}\hskip 2.84526pt\text{: LL
strong duality}\\},$ (28f) with the C&CG constraints updated at every
iteration $k$ as:

$\displaystyle||{\boldsymbol{\lambda^{c,WC}}}||_{2}^{2}\geq||{\boldsymbol{\lambda^{c,WC,(j-1),*}}}||_{2}^{2},~{}j=2,\dots,k,$
(28g)
$\displaystyle\lambda_{t}^{c,WC}\geq\lambda_{t}^{c,WC,(j-1),*},~{}\forall
t\in\mathcal{T},~{}j=2,\dots,k,$ (28h)

####  V-B2 Subproblem

$\displaystyle\min_{\Xi^{SP}}$
$\displaystyle~{}||{\boldsymbol{\lambda^{c,WC}}}||_{2}^{2}$ (29a)
$\displaystyle s.t.~{}v+\sum_{s}\phi^{s}\zeta_{s}\leq 0,$ (29b)
$\displaystyle\alpha\zeta^{s}\geq-\tilde{C}^{WC}(d^{s},\lambda^{c},x^{(k),*})-v,~{}\forall
s\in 1:S$ (29c) $\displaystyle\zeta^{s}\geq 0,~{}\forall s\in 1:S$ (29d)
$\displaystyle\lambda_{t}^{c,WC}\geq 0,~{}\forall t\in 1:T$ (29e)
$\displaystyle\text{Eqs.~{}}\eqref{eq:gen_const}-\eqref{eq:load_const}\hskip
24.18483pt:\text{LL primal constraints}$ (29f)
$\displaystyle\text{Eqs.~{}}\eqref{eq:dll_constraint1}-\eqref{eq:dll_constraint3}\hskip
24.18483pt:\text{LL dual constraints}$ (29g) $\displaystyle
O_{LL}^{s}=O_{DLL}^{s},~{}\forall s\in\mathcal{S}\hskip 5.69054pt:\text{LL
strong duality}\\}.$ (29h)

where $\Xi^{MP}=\\{\Xi^{I}\cup\Xi^{E}\cup\Xi^{G}\\}$ and
$\Xi^{SP}=\Xi^{E}\cup\Xi^{G}$. The upper and lower bounds, for the convergence
of the C&CG algorithm, are calculated as the values of
$x^{WC}+||{\boldsymbol{\lambda^{c,WC}}}||_{2}^{2}$ in the subproblems and the
principal problem, respectively. For any iteration $k$, the value of
$x^{WC,(k),*}$ obtained from the principal problem is passed to the
subproblems, which returns the value of $\lambda_{t}^{WC,c*,j},\forall j\leq
k$. Note that at the first iteration, the principal problem is solved without
considering C&CG constraints in Eqs. (28g) and (28h). The principal problem
and subproblems are iterated until the lower and upper bounds converge within
an user-defined tolerance. The additional cut in Eq. (28h) is introduced in
the ML objective to boost the convergence of the C&CG algorithm.

##  VI Case Study

The case study is performed using power grid and commercial EVCS data from
Manhattan, NY. Fig. 3 displays weekly EVCS demand profiles collected as
explained in [20]. Also based on [20], the power grid model is a 7-bus network
with 4 load buses (#3, 4, 5, and 6) and generation resources located at buses
#1, 2, 5, and 7. The non-EVCS demand profiles at each bus are built based on
load distribution curves for New York City published by the New York
Independent System Operator (NYISO). This distribution network is connected to
NYISO at buses #1 and 7, where we set LMPs using data from [37]. Additionally,
generation resources at buses #2 and 5 have incremental production costs of
$15/MWh and $10/MWh.

We now define the data for insurance policy factors $\Pi$. The box-constrained
confidence interval for the uncertain insurance policy factors are defined as
$r=[0.25~{}0.35]$, $\kappa=[0.2~{}0.3]$, and
$\mathbb{P}(A)=[0.03582~{}0.04378]$. The box limits for $r$ and $\kappa$ are
constructed using industry practices across insurers [13, 12], while the
limits for $\mathbb{P}(A)$ are constructed around the value of $\mathbb{P}(A)$
obtained using the SMP in Section IV. Notably, the box limits of
$\mathbb{P}(A)$ consider $\epsilon=\pm 10\%$ in $\mathbb{P}(A)$ estimated by
the SMP. To compute $\mathbb{P}(A)$, Table I shows the data-driven Weibull
distribution parameters for the SMP state transitions and the data source for
modeling the transition distributions. Since there is limited data on EVCSs-
related cyberattacks recorded publicly, we use cyberattack data from IT
industry available from SANS Institute [32], the Canadian Institute for
Cybersecurity [38], and empirical studies [39, 33, 30]. Similarly, Table II
shows the sojourn time and the steady-state probabilities of each SMP state.

The case study is carried out on a 64-bit personal computer with a 2.2 GHz
Intel Core i7 processor and 8 GB RAM. The implementation uses the Ipopt and
Gurobi solvers. All the simulation instances presented below are solved within
1 hour.

TABLE I: Weibull distribution parameters for SMP state transitions.

State Transition | Parameter | Data | Source
---|---|---|---
Shape $\beta$ | Scale $\alpha$
G-I | 2.0675 | 18.8178 | Real | [38]
I-D | 1.9293 | 16.0712 | Real | [32]
D-C | 1.5698 | 18.4858 | Real | [32]
C-G | 1.3816 | 15.7033 | Real | [32]
I-F | 0.7000 | 400.000 | Approx. | [30]
F-G | 0.6783 | 13.4487 | Approx. | [39]

TABLE II: Sojourn time and steady state probabilities of SMP states.

State $s$ | Entity
---|---
Sojourn Time $T_{s}$ (hour) | Probability $\mathbb{P}_{s}$
G | 9.1016 | 0.2010
I | 13.3431 | 0.2943
D | 11.7754 | 0.2366
C | 11.3659 | 0.2283
F | 19.8271 | 0.03980

###  VI-A EVCS With A Predetermined Tariff

This section studies the performance of the premium obtained using the
analytical result in Theorem 1 and its robust and risk-averse reformulation
described in Section V. Alongside, we also study the sensitivity of the
insurance premium to policy factors $\Pi$. To carry out these analyses,
predetermined electricity tariffs in the power grid are obtained by solving
the optimization problem in Eqs. (10c)-(10g). The value of the predetermined
electricity tariffs are shown in Fig. 4(a). Note that since the predetermined
electricity tariff at bus #5 is same as for other buses except bus #4, we
study only buses #4 and #5. It is also observed that the predetermined
electricity tariff $(\lambda^{u})$ in Fig. 4(a) is independent of the EVCS
demand, which follows from the dissimilarity of the hourly profiles for
$\lambda^{u}$ and the EV charging price $(\lambda^{c})$ in Fig. 4(a). Fig.
4(b) presents the values of $\hat{\lambda}^{u}$, $\hat{\lambda}^{c}$, and
$\hat{x}$ at buses #4 and 5, where $\hat{\lambda}^{u}$ and $\hat{\lambda}^{c}$
are the average values of $\lambda^{u}$ and $\lambda^{c}$ in Fig. 4(a) and
$\hat{x}$ is the insurance premium per KWh EVCS demand calculated using Eq.
(3). Note that $\sum_{s}\phi^{s}\sum_{t}\hat{x}d_{t}^{s}=x$. Since
$\hat{\lambda}^{u}$ at bus #4 is greater than at bus #5, as Fig. 4(b) shows,
$\hat{\lambda}^{c}$ and $\hat{x}$ at bus #4 are also greater than those at bus
#5.

Figure 3: Hourly power consumption of EVCSs in Manhattan, NY. Figure 4:
Performance of EVCS at buses $\\#$4 and 5: (a) Hourly profiles of the
predetermined electricity tariff ($\lambda^{u}$) and the EV charging price
offered by the EVCS ($\lambda^{c}$) and (b) The insurance premium per kWh of
the EVCS demand ($\hat{x}$), the 24-hour average electricity tariff
($\hat{\lambda}^{u}$), and the 24-hour average EV charging price
($\hat{\lambda}^{c}$).  Figure 5: Sensitivity of the cyber insurance premium
per kWh EVCS demand $(\hat{x})$ based on the expected EVCS cost at bus #4 with
respect to changes in the policy parameters: (a) Changes in $r$,
$\mathbb{P}(A)$, and $\gamma$ and (b) Changes in $\kappa$ and $A_{h}$. The
expected premium is calculated with $r=0.3$, $\mathbb{P}(A)=0.03980$,
$\gamma=1$, $\kappa=0.25$, $A_{h}=1$, and $\rho=3$.

Fig. 5 illustrates the performance of the premium as a function of policy
factors $\Pi$ at bus #4. To analyze the sensitivity of these factors, we study
the analytically calculated expected insurance premium in Eq. (3), which uses
the mean value of $\mathbb{P}(A)$, $r$, and $\kappa$ obtained from their box-
constrained confidence interval defined above, i.e., $\mathbb{P}(A)=0.03989$,
$r=0.3$, and $\kappa=0.25$. Furthermore, we assume that the EVCS transfers
$100\%$ of its cyber risk to the insurer, i.e., $\gamma=1$ and the EVCS has
encountered one cyberattack in the past, i.e., $A_{h}=1$. In this setting, as
Fig 5(a) shows, the insurance design is more sensitive to $\mathbb{P}(A)$ than
to other insurance policy factors. Similarly, the premium increases at a
greater rate when the insurer is less confident about the values of the policy
factors and the revenue reported by the EVCS, i.e. premium $x$ increases
faster after the insurer increases profit loading factor $r$ above $50\%$ as
seen in Fig. 5(a). On the other hand, changes in risk transfer factor $\gamma$
affect the premium linearly.

Figure 6: Changes in (a) the average EV charging price $(\hat{\lambda}^{c})$
and (b) the cyber insurance premium per kWh EVCS power $(\hat{x})$ for
different levels of the EVCS risk-aversion measured by $\alpha$. The EV
charging price and premium are bounded by the uncertainty in policy parameters
$\mathbb{P}(A)$, $r$, and $\kappa$.

Similarly, Fig 5(b) demonstrates how the premium changes with the cyberattack
history of the EVCS. A close-to-linear relationship is established between
number of EVCS cyberattacks $A_{h}$ and premium $x$. Notably, an increase in
the value of $\kappa$ leads to a faster premium increase.

In addition to the premium obtained from expected cost in Fig. 5, Fig. 6
presents the bounds of premium and EV charging price for a risk-averse EVCS.
Note that the bounds are determined numerically using Eq. (26). The lower and
the upper bounds of $\hat{x}$ and $\lambda^{c}$ in Fig. 6 use the lower and
upper limits of the box-constrained confidence interval for policy parameters.
That is, $\mathbb{P}(A)=0.03582$, $r=0.25$, and $\kappa=0.2$ for lower bounds,
while $\mathbb{P}(A)=0.04378$, $r=0.35$, and $\kappa=0.3$ for upper bounds. In
this setting, as shown in Fig. 6(a), as the EVCS becomes more risk-averse
(i.e., a decrease in $\alpha$), it chooses the “worse-case” typical days. As a
result, the EVCS optimization in Eq. (20) chooses a greater value of
$\lambda^{c}$ such that the EVCS net cost in Eq. (20b) is at most zero.
Following the increased value of $\lambda^{c}$, as shown in Fig. 6(b), the
optimization in Eq. (27) also returns a higher value of $\hat{x}$. Moreover,
as shown in Fig. 6, the bounds on $\hat{x}$ and $\lambda^{c}$ widen as the
EVCS becomes more risk-averse (i.e., the value of user-defined parameter
$\alpha$ reduces).

###  VI-B EVCS With A DLMP-based Tariff

This section investigates the relationship between insurance premium $\hat{x}$
and dynamic electricity tariff $\lambda^{u}$ set by the utility using the
DLMPs. Since the current penetration of EVCS demand as in Fig. 3 is small
relative to the system load, we analyze the effect of $\lambda^{u}$ on the
premium design by increasing EVCS demands as indicated in Fig. 7 and Tables
III and IV. Fig. 7 present the sensitivity of the insurance premium, the EV
charging price offered by the EVCS, and the electricity tariff at buses #4 and
5 to an increase in the EVCS demand. As Fig. 7 shows $\hat{x}$,
$\hat{\lambda}^{c}$, and $\hat{\lambda}^{u}$ start to change faster with a
smaller increase in EVCS demand (e.g. 100x) at bus #4 than at bus #5. Thus, we
focus on bus #4. Table III presents the changes in the bounds of
$\hat{\lambda}^{c}$, designed by a risk-averse and robust EVCS, with an
increase in the EVCS demand. As Table III shows the bounds on
$\hat{\lambda}^{c}$ widen when the EVCS becomes more risk-averse (i.e., the
value of user-defined parameter $\alpha$ reduces). Also, the results of Table
III demonstrate that the EV charging price offered by the EVCS increases
because of the increased EVCS demand in the distribution system, which is
followed by the increase in the electricity tariff. Similarly, Table IV shows
the changes in the insurance premium with an increase in the EVCS demand.
Similar to the $\hat{\lambda}^{c}$ in Table III, the premium in Table IV
increases as the EVCS power capacity increases and the EVCS opts for a more
risk-averse insurance design.

Figure 7: Changes in $\hat{x}$, $\hat{\lambda}^{u}$, and $\hat{\lambda}^{c}$
with an increase in EVCS demand by 100x and 1000x at (a) bus #4 and (b) bus
#5. TABLE III: $\hat{\lambda}^{c}$ with an increase in EVCS demand at bus #4
(¢/kWh)

d | $\alpha=1$ | $\alpha=0.5$ | $\alpha=0$
---|---|---|---
LB | EX | UB | LB | EX | UB | LB | EX | UB
1x | 26.74 | 28.16 | 29.71 | 28.57 | 30.16 | 21.92 | 29.90 | 31.70 | 33.72
100x | 26.95 | 28.37 | 29.93 | 28.80 | 30.40 | 32.16 | 30.18 | 31.98 | 34.01
400x | 27.26 | 28.69 | 30.24 | 29.15 | 30.75 | 32.51 | 30.52 | 32.32 | 34.36
800x | 27.68 | 29.11 | 30.67 | 29.69 | 31.29 | 33.07 | 31.18 | 32.99 | 35.04
1000x | 27.98 | 29.41 | 30.98 | 30.08 | 31.69 | 33.47 | 31.66 | 33.49 | 35.55

  * •

LB: Lower Bound, EX: Expected, UB: Upper Bound, d: EVCS Demand

TABLE IV: $\hat{x}$ with an increases in EVCS demand at bus #4 (¢/kWh)

d | $\alpha=1$ | $\alpha=0.5$ | $\alpha=0$
---|---|---|---
LB | EX | UB | LB | EX | UB | LB | EX | UB
1x | 1.801 | 2.351 | 3.057 | 2.062 | 2.699 | 3.519 | 2.519 | 3.310 | 4.338
100x | 1.815 | 2.369 | 3.089 | 2.079 | 2.720 | 3.546 | 2.542 | 3.340 | 4.375
400x | 1.836 | 2.395 | 3.111 | 2.105 | 2.754 | 3.588 | 2.570 | 3.375 | 4.420
800x | 1.863 | 2.372 | 3.155 | 2.144 | 2.803 | 3.649 | 2.625 | 3.445 | 4.508
1000x | 1.884 | 2.409 | 3.187 | 2.172 | 2.838 | 3.694 | 2.666 | 3.497 | 4.573

  * •

LB: Lower Bound, EX: Expected, UB: Upper Bound, d: EVCS Demand

##  VII Conclusion

This paper develops models for optimizing cyber insurances against
cyberattacks on commercial EVCSs, which internalize operational dynamics of
the EVCSs and the power distribution company. The insurance premium for EVCSs
operations under the predetermined electricity tariff is modeled using a bi-
level optimization problem, while the premium for EVCSs operating under the
DLMPs is modeled using a tri-level optimization problem. Furthermore, the
models are robustified against parameter uncertainty in insurance
characteristics and risk attitudes of the EVCS. The analytical results
demonstrate that the insurance design is very sensitive to the probability of
attacks on EVCSs in comparison to other design parameters. The paper presents
how the insurance premium bounds for a range of design parameters. Further,
numerical results show that the premium and the EV charging price increase for
more risk-aware EVCSs.

## References

  * [1] F. Nick _et al._ , “Cybersecurity threats escalate in the energy sector,” Energy Monitor. [Online]. Available: https://tinyurl.com/an8mpn7z
  * [2] IBM Security, “Cost of a data breach report,” IBM Security, Armonk, NY, USA, July 2020. [Online]. Available: https://tinyurl.com/3hmjnv5t
  * [3] S. Acharya _et al._ , “Cybersecurity of smart electric vehicle charging: A power grid perspective,” _IEEE Access_ , vol. 8, 2020.
  * [4] S. Soltan _et al._ , “Blackiot: Iot botnet of high wattage devices can disrupt the power grid,” in _USENIX Security 18_ , 2018, pp. 15–32.
  * [5] S. Amini _et al._ , “Dynamic load altering attacks against power system stability,” _IEEE Trans. Smart Grid_ , vol. 9, no. 4, 2016.
  * [6] G. Raman _et al._ , “Manipulating residents’ behavior to attack the urban power distribution system,” _IEEE Trans. Ind. Inf._ , vol. 15, no. 10, 2019.
  * [7] S. Acharya _et al._ , “Causative cyberattacks on online learning-based automated demand response systems,” _IEEE Trans. Smart Grid_ , 2021.
  * [8] J. Ospina _et al._ , “On the feasibility of load-changing attacks in power systems during the covid-19 pandemic,” _IEEE Access_ , vol. 9, 2021.
  * [9] Cybersecurity & Infrastructure Security Agency. ”Cybersecurity Insurance,”. [Online]. Available: https://tinyurl.com/p56h2dm2
  * [10] G. Aarti _et al._ , “Cyber insurance market outlook - 2026,” Applied Market Research, March 2020. [Online]. Available: https://tinyurl.com/jnp7ak3r
  * [11] L. Sanicola, “After colonial attack, energy companies rush to secure cyber insurance.” [Online]. Available: https://tinyurl.com/e9zf6dde
  * [12] Fitch, “Sharply rising cyber insurance claims signal further risk challenges,” Fitch. [Online]. Available: https://tinyurl.com/46t27evu
  * [13] S. Romanosky _et al._ , “Content analysis of cyber insurance policies: how do carriers price cyber risk?” _J. Cyber._ , vol. 5, no. 1, p. tyz002, 2019\.
  * [14] D. Woods _et al._ , “The county fair cyber loss distribution: Drawing inferences from insurance prices,” in _W. Eco. Inf. Sec._ , 2019.
  * [15] R. N. Allan _et al._ , _Reliability evaluation of power systems_. Springer Science & Business Media, 2013.
  * [16] R. Zhang _et al._ , “A bl game approach to attack-aware cyber insurance of computer networks,” _IEEE J. Sel. Areas Comm._ , vol. 35, no. 3, 2017.
  * [17] K. Harnett _et al._ , “DoE/DHS/DoT volpe technical meeting on electric vehicle and charging station cybersecurity report,” Tech. Rep., 2018.
  * [18] O. G. M. Khan _et al._ , “Impact of electric vehicles botnets on the power grid,” in _2019 IEEE EPEC_. IEEE, 2019, pp. 1–5.
  * [19] M. Basnet _et al._ , “Exploring cybersecurity issues in 5g enabled electric vehicle charging station with deep learning,” _arXiv:2104.08553_ , 2021.
  * [20] S. Acharya _et al._ , “Public plug-in electric vehicles+ grid data: Is a new cyberattack vector viable?” _IEEE Trans. Smart Grid_ , 2020.
  * [21] T. Maynard _et al._ , “Business blackout: The insurance implications of a cyber-attack on the us power grid,” _Lloyd’s of London._ , vol. 15, 2015\.
  * [22] P. Lau _et al._ , “A cybersecurity insurance model for power system reliability considering optimal defense resource allocation,” _IEEE Trans. Smart Grid_ , vol. 11, no. 5, pp. 4403–4414, 2020.
  * [23] Z. Yang _et al._ , “Premium calculation for insurance businesses based on cyber risks in IP-Based power substations,” _IEEE Access_ , vol. 8, 2020\.
  * [24] D. T. Hoang _et al._ , “Charging and discharging of PEVs in V2G systems: A cyber insurance-based model,” _IEEE Access_ , vol. 5, 2017.
  * [25] M. Caramanis _et al._ , “Co-optimization of power and reserves in dynamic T&D power markets with nondispatchable renewable generation and distributed energy resources,” _Proc. IEEE_ , vol. 104, no. 4, 2016.
  * [26] Consolidated Edison. ”Electric Vehicle”, accessed march 20, 2021. [Online]. Available: https://tinyurl.com/yvjx2e8a
  * [27] M. E. Baran _et al._ , “Optimal capacitor placement on radial distribution systems,” _IEEE Trans. Pow. Del._ , vol. 4, no. 1, 1989.
  * [28] B. Zeng _et al._ , “Solving two-stage robust optimization problems using a column-and-constraint generation method,” _Oper. Res. Lett._ , 2013.
  * [29] C. Ruiz and A. J. Conejo, “Robust transmission expansion planning,” _Eu. J. Oper. Res._ , vol. 242, no. 2, pp. 390–401, 2015.
  * [30] Q. Liu _et al._ , “Probabilistic modeling and analysis of sequential cyber-attacks,” _Engineering Reports_ , vol. 1, no. 4, 2019.
  * [31] G. Kumar _et al._ , “Availability analysis of repairable mechanical systems using analytical semi-markov approach,” _Quality Engineering_ , vol. 25, no. 2, 2013.
  * [32] B. Matt, “SANS 2019 incident response survey: It’s time for a change,” SANS Institute, USA. [Online]. Available: https://tinyurl.com/ysrfnxz6
  * [33] B. Schroeder _et al._ , “A large-scale study of failures in high-performance computing systems,” _IEEE Trans. Depen. & Secure Comp._, vol. 7, 2009.
  * [34] R. T. Rockafellar, “Coherent approaches to risk in optimization under uncertainty,” in _OR Tools and Applications_. Informs, 2007.
  * [35] C. A. Hans _et al._ , “Risk-averse model predictive operation control of islanded microgrids,” _IEEE Trans. Ctrl. Sys. Tech._ , vol. 28, no. 6, 2019\.
  * [36] Z. H. Gümüş and C. A. Floudas, “Global optimization of nonlinear bilevel programming problems,” _J. Glo. Opt._ , vol. 20, no. 1, 2001.
  * [37] NYISO. Accessed May 1, 2021. [Online]. Available: https://www.nyiso.com/real-time-dashboard
  * [38] Canadian Institute for Cybersecurity. ”CSE-CIC-IDS2018 on AWS,”. [Online]. Available: https://tinyurl.com/zbcv7rfd
  * [39] U. Franke, H. Holm, and J. König, “The distribution of time to recovery of enterprise it services,” _IEEE Trans. Rel._ , vol. 63, no. 4, 2014.