# Vulnerability Assessment of Charging Stations in the Electrified Road Network

Ke Zhang, Yi-Ping Fang, Hongping Wang, Yufeng Zhuang, Enrico Zio

► To cite this version:

Ke Zhang, Yi-Ping Fang, Hongping Wang, Yufeng Zhuang, Enrico Zio. Vulnerability Assessment of Charging Stations in the Electrified Road Network. IEEE Transactions on Intelligent Transportation Systems, 2025, pp.1-12. 10.1109/TITS.2025.3548586. hal-05050521

HAL Id: hal-05050521

https://hal.science/hal-05050521v1

Submitted on 29 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Vulnerability Assessment of Charging Stations in the Electrified Road Network

Ke Zhang<sup>1,0</sup>, Yi-Ping Fang<sup>1,0</sup>, Member, IEEE, Hongping Wang<sup>1,0</sup>, Member, IEEE, Yufeng Zhuang, and Enrico Zio<sup>1,0</sup>, Fellow, IEEE

**Abstract**—As the adoption of electric vehicles (EVs) within electrified road networks (ERNs) continues to grow, the criticality of fast-charging stations (FCSs) to support this growth and alleviate range anxiety for EV users is increasingly evident. Due to the complex operating environment of FCSs, their reliable, safe operation faces significant challenges, from exposure to natural disasters, deliberate attacks, and technical failures. This leads to the need to analyze the vulnerability of FCSs within ERNs. This paper provides a mathematical framework to assess the vulnerability of the ERN to intentional attacks and identify critical FCSs within ERNs. We propose a System Optimal Dynamic Mixed Traffic Flow Assignment Model (SODTA) game-theoretical for an ERN where electric and fuel vehicles coexist. Building upon this model, an attacker-defender (AD) framework is established for assessing the vulnerability of FCSs within the ERN. We transform this AD problem into a mixed integer linear programming (MILP) one to solve it efficiently. Two indices are proposed to evaluate the vulnerability and the reallocation of charging loads within the ERN. Finally, we apply the proposed framework to conduct a vulnerability assessment of the ERN in North Carolina, USA. The key findings are as follows: 1) As the EV penetration increases, the studied ERN becomes more vulnerable to attacks. 2) With limited attack resources, mitigating performance losses over time is possible through charging load redistribution. 3) Variation in overall vulnerability levels with increasing attack resources reveals distinct patterns, ranging from mitigable to unmitigable vulnerability.

**Index Terms**—Electric vehicles (EVs), electrified road network (ERNs), vulnerability assessment, fast-charging stations (FCSs), defender-attacker game model, system optimal dynamic mixed traffic flow assignment.

## I. INTRODUCTION

THE introduction of electric vehicles (EVs) in electrified road networks (ERNs) is gaining momentum [1]. Despite

Received 5 December 2023; revised 1 July 2024 and 30 November 2024; accepted 1 March 2025. This work was supported in part by the Research Innovation Fund for College Students of Beijing University of Posts and Telecommunications under Grant 202412011, in part by the 2025 Beijing National Science Foundation-Xiaomi Innovation Joint Fund under Grant L243024, and in part by the Fundamental Research Funds for the Central Universities under Grant 2023RC37. The work of Yi-Ping Fang was supported by the Agence Nationale de la Recherche (ANR) under Grant ANR-22-CE39-0003. The Associate Editor for this article was Y. Kim. (Corresponding author: Hongping Wang.)

Ke Zhang, Hongping Wang, and Yufeng Zhuang are with the School of Intelligent Engineering and Automation, Beijing University of Posts and Telecommunications, Beijing 100876, China (e-mail: hongping.wang@bupt.edu.cn).

Yi-Ping Fang is with the Laboratoire Génie Industriel, CentraleSupélec, Université Paris-Saclay, 91190 Gif-sur-Yvette, France.

Enrico Zio is with Mines Paris-PSL, Centre de Recherche sur les Risques et les Crises (CRC), 06560 Sophia Antipolis, France, and also with the Energy Department, Politecnico di Milano, 20133 Milan, Italy.

This article has supplementary downloadable material available at https://doi.org/10.1109/TITS.2025.3548586, provided by the authors.
Digital Object Identifier 10.1109/TITS.2025.3548586

the increasing number of charging stations (CSs), it is anticipated that there will be a rise in the number of EVs to be served per CS [1]. Then, reliable and safe operation of CSs is needed to meet the expectations of EV users of uninterrupted functionality of the ERNs. Nonetheless, intentional cyber-attacks [2], extreme weather events [3], [4], and unusual operation conditions pose potential threats to the operation of CSs and the ERNs which accommodate them. Yet, the vulnerability that CSs introduce to the ERN is a relatively understated aspect.

In general, the definition of vulnerability may differ depending on the particular context. In the context of network systems (like the ERNs), vulnerability can be categorized into various types, such as connectivity vulnerability, accessibility vulnerability [5], and maximum flow-based vulnerability [6]. These concepts of vulnerability may also be contingent upon the system's performance metrics, topological, such as the shortest path between source and sink nodes, and physical, such as the traffic or power flow. To analyze the vulnerability of the ERN from a system-level perspective, this paper adopts a broad definition of vulnerability as articulated in [7], i.e., the transportation system's susceptibility to significant disruptions caused by extreme events.

Then, different network vulnerability assessment methods exist, which can be broadly categorized into three types: 1) static topology-based methods [5], [8], [9], 2) simulation-based methods [10], [11], [12], and 3) optimization-based methods [13], [14], [15]. In static topology-based methods, the system is typically abstracted into a static network. Then, network theory-based methods are used to assess the network's vulnerability and nodes' importance by various connectivity, closeness, and betweenness-based indices. For example, in [5], multiple topological indices, including normalized edge betweenness centrality, were utilized to measure the impact of extreme floods on the road network. In [8], the vulnerability of the Eurasian maritime transportation system was analyzed using three centrality measures, which were integrated via joint entropy to assess the importance of ports. However, due to the static nature of the topology-based methods, they have limitations in representing the spatio-temporal characteristics of dynamic systems such as ERNs.

Simulation-based approaches involve simulating the traffic or power flows within the system. For the vulnerability assessment, components within the system, such as roads or electric lines, are stochastically or consecutively removed to analyze various future scenarios [16], [17], [18]. System performance is, then, evaluated based on the system's response. In [4], several indicators based on CS utilization were used to evaluate
the impact of flooding on CSs, in terms of average, extreme, and variance utilization, and in terms of the closest distance between EVs and available CSs. However, it is important to note that in order to reduce computational time, the simulating of the flow in the network after disturbances may not precisely represent vulnerability reality, thus potentially leading to inaccurate assessment.

Optimization-based vulnerability assessment methods utilize a two-layer optimization known as the Attacker-Defender model (AD), also referred as the bi-level interdiction problem [19], [20]. This model is able to pinpoint the most vulnerable links or nodes within the system and quantify the loss in system performance under worst-case scenarios. It offers the flexibility to comprehensively account for various factors, including network topology, optimal responses after disturbances, and the spatio-temporal characteristics of the system. Thanks to the advancement of efficient algorithms [21], [22], [23], [24], this method has gained popularity for assessing the vulnerability of complex systems. For instance, [25] used this model to assess the vulnerability of a multi-mode equilibrium traffic network. Reference [26] added a defender layer to the AD model, resulting in the defender-attacker-defender model (DAD). The inner-most two layers were used to identify the critical lines and nodes in the power system, whereas the outer-most layer was used to devise optimal defense strategies that reduce the impact of worst-case attacks and enhance the system's resilience.

Recently, the vulnerability of systems comprising EVs has emerged as a significant area of research. Previous studies have predominantly analyzed the vulnerability of these systems in the context of power grids. For instance, [27] investigated the vulnerability of power grids with EVs using simulation methods, approximating the traffic flow of EVs as virtual power flow. Similarly, [14] developed a vulnerability assessment method based on the AD model for coupled transportation and multi-energy networks. In their study, traffic flow is calculated using a static user equilibrium model and is also simplified into virtual flow, employing a methodology akin to that of [27], [28], and [29] discussed existing flaws within the EV charging system and examined the vulnerabilities of the grid introduced by substantial EV loads manipulated by malicious attackers. However, in these studies, the intricacies of the road network level, including EV behavior, dynamic charging processes, and traffic flow, are often overlooked or simplified. While extensive research has focused on the lines within grid and road segments, FCSs remain relatively understudied as objects of vulnerability assessment. To address this research gap, this study assesses the vulnerability of road networks with EVs, i.e., ERNs, especially from a transportation systems perspective. The primary objectives are: 1) to identify vulnerable FCSs within the ERN whose malfunctions result in significant performance losses and 2) to evaluate the impact of these malfunctions on ERNs using a quantifiable vulnerability index. To this aim, this paper proposes a game theory framework to assess the vulnerability of an ERN, adopting a worst-case and disruption-agnostic perspective.

The main contributions of this paper are listed as follows:

* This work fills the research gap in analyzing the vulnerability of FCS in the ERN by establishing an AD model,

* The mixed SODTA model proposed in [30] is improved by simplifying the charging constraints of EVs.

* The proposed AD model is first transformed into an equivalent mixed-integer linear programming (MILP) model, and an accelerated Benders decomposition algorithm is further employed to enhance the computational efficiency.

The rest of this paper is arranged as follows. Section II presents the vulnerability assessment framework, including the ERN model in Section II-A, the AD model in Section II-B, along with its solution method in Section. II-C, and the metrics used to perform the vulnerability assessment proposed in Section II-D. Section III thus takes the use of the presented framework for the vulnerability assessment of an experimental network. Section IV gives a summary of the work.

## II. VULNERABILITY ASSESSMENT FRAMEWORK

### A. Electrified Road Network Modeling

Assessing the vulnerability of the system using optimization-based methods requires an operational model for ERNs. The difference between modeling an ERN as opposed to a conventional road network lies in the representation of the EVs' charging process and energy update. It is well established that there are two methods to assign traffic flow of conventional vehicles: system optimal (SO) and user equilibrium (UE) [31]. The fundamental principle can be extended to accommodate mixed traffic flow involving EVs and gasoline vehicles (GVs). Apart from the ways for assigning traffic flow, modeling the ERN can also be categorized as static and dynamic. For instance, [32] employed a static user equilibrium approach to model the ERN, representing CSs using queuing theory, with the assumption that EVs charged a predetermined amount of energy at each CS. Given the dynamic characteristics of FCSs and EVs, it is more appropriate to model the ERN dynamically. To name but a few, [33] introduces an advanced decentralized framework for electricity trading among connected EVs in parking lots. It details an intelligent contract mechanism for maximizing EVs' profitability through a stochastic bidding process. In [30], a dynamic system optimal model was employed to describe the ERN. This model is grounded in the cell transmission model originally proposed and substantiated in line with the hydrodynamic traffic flow theory in [34], and later formulated into linear programming by [35]. The computational efficiency of the cell transmission model was improved through the link transmission model [36], which was also formulated into linear programming, as detailed in [37]. The link transmission model was further extended to mixed traffic flow scenarios in [38]. It is worth noting that [38] highlighted certain redundant variables within the model used to describe the charging status of EVs, and revising the constraints might further improve the model's solution efficiency. The present study addresses this issue by eliminating these redundant variables and revising the corresponding constraints. Then, this refined model is employed for the vulnerability analysis of an ERN.

The ERN is represented by $G(N, A)$. The set of links is denoted by $A$, and the set of nodes by $N$. The nodes are classified into two types: general nodes $N_G$ and source/sink
nodes $\mathcal{N}_{SR}$. General nodes are intersections on the highway connecting links, whereas source/sink nodes represent cities in the ERN and serve as both the origin and destination of front paths on the ERN. The four kinds of links in the model are the source, general, charging, and sink links. A source link $a \in \mathcal{A}_R$ is where the mixed traffic flow enters the ERN, whereas a sink link $a \in \mathcal{A}_S$ is where it exits upon reaching its destination. Each pair of source and sink links is connected to a source/sink node. Charging links $a \in \mathcal{A}_C$ represent FCSs at intersections and provide electricity for EVs. Source, sink, and charging links are dummy links. The general link $a \in \mathcal{A}_G$ represents homogeneous highways allowing GVs and EVs to transit. Charging links can be represented as a circle with a self-loop structure, which allows EVs to recharge on the road.

We use a discrete set $\mathcal{E} = \{1, 2, \dots, E\}$ to describe the state of charge (SOC) of the EVs. $E$ is the maximum energy level (EL) of EVs. The time horizon is also discretized into equal length $\delta$, and $\mathcal{T} = \{t|t = 0, 1, \dots, T_{max}\}$ is the set for all discretized time steps. The length of the time step corresponds to real time intervals. Shorter time steps facilitate a more detailed analysis of the system's dynamics; however, they also demand greater computational resources. $UG_a^s(t)$ and $VG_a^s(t)$ represent the cumulative number of GVs that have entered and left link $a$ to destination $s$ by the end of time step $t$, respectively. $UE_a^{s,e}(t)$ and $VE_a^{s,e}(t)$ represent the cumulative number of EVs with energy level $e$ that enter and leave link $a$ to destination $s$ by the end of the time step $t$.

The attributes of each link are inflow, outflow, and carrying capacity. The inflow, outflow, and the number of cars on link $a$ in time step $t$ should not go beyond the inflow, outflow, and carry capacity of the link, respectively:

$$ \sum_{s \in \mathcal{N}_{SR}} [VG_a^s(t) - VG_a^s(t-1)] + \sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [VE_a^{s,e}(t) - VE_a^{s,e}(t-1)] \le f_a^O(t), \forall a \in \mathcal{A} \setminus \{\mathcal{A}_C\}, \forall t, s \quad (1) $$

$$ \sum_{s \in \mathcal{N}_{SR}} [UG_a^s(t) - UG_a^s(t-1)] + \sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [UE_a^{s,e}(t) - UE_a^{s,e}(t-1)] \le f_a^I(t), \forall a \in \mathcal{A} \setminus \{\mathcal{A}_C\}, \forall t, s \quad (2) $$

$$ \sum_{s \in \mathcal{N}_{SR}} [VG_a^s(t) - VG_a^s(t-\beta_a)] + \sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [UE_a^{s,e}(t) - VE_a^{s,e}(t-\beta_a)] \le N_a(t), \forall a \in \mathcal{A} \setminus \{\mathcal{A}_C\}, \forall t, s \quad (3) $$

where $f_a^O(t)$, $f_a^I(t)$, $N_a(t)$ denote the outflow capacity of link $a$ inflow capacity $a$, and the maximum number of cars that can be carried by link $a$ at the time step $t$, respectively. $\beta_a$ is the travel time required by the backward shock wave from the exit to the entry of link $a$. GVs and EVs need to spend at least $v_a$ time steps to traverse link $a$, where $v_a$ is the free-flow travel time on link $a$. Therefore, the cumulative number of GVs/EVs that enter link $a$ before the time step $t - v_a$ is less than or equal to the cumulative number of GVs/EVs that leave link $a$ before time step $t$. This condition can be expressed as:

$$ \sum_{s \in \mathcal{N}_{SR}} VG_a^s(t) \le \sum_{s \in \mathcal{N}_{SR}} UG_a^s(t - v_a), \forall a \in \mathcal{A} \setminus \{\mathcal{A}_C\}, \forall t \quad (4) $$

For EVs, (5a) guarantees that the outflow of EVs should be less than or equal to the inflow before $v_a$, and it also ensures that the ELs of the outflow of EVs are updated; (5b) ensures that the ELs of all EVs should be less than $E$:

$$ VE_a^{s,e-\rho_a}(t) \le UE_a^{s,e}(t - v_a), \quad \forall a \in \mathcal{A} \setminus \{\mathcal{A}_C\}, e \in \mathcal{E} \cap \{e > \rho_a\}, \forall s, t \quad (5a) $$

$$ UE_a^{s,e}(t) = 0, \forall a \in \mathcal{A} \setminus \{\mathcal{A}_C\}, e \in \mathcal{E} \cap \{e \le \rho_a\}, \forall s, t \quad (5b) $$

where $\rho_a$ represents the energy consumption of the travel link $a$, which can be calculated by the following equation: $\rho_a = L_a/(\delta \cdot v_f)$, and $L_a$ denotes the physical length of link $a$.

The traffic demand is met by letting the cumulative inflow of GVs/EVs at the source link be equal to the cumulative demand, as follows:

$$ UG_a^s(t) = DG_a^s(t), \forall a \in \mathcal{A}_R, \forall s, t \quad (6a) $$

$$ U E_a^{s,e}(t) = D E_a^{s,e}(t), \forall a \in \mathcal{A}_R, e \in \mathcal{E}, \forall s, t \quad (6b) $$

where $DG_a^s(t)/DE_a^{s,e}(t)$ represents the cumulative traffic demand for GVs and EVs, respectively, at the entrance of the source link $a$ by the end of the time step $t$.

To ensure traffic flow conservation, we make the outflow at ordinary nodes equal to their inflow, as shown below:

$$ \sum_{a \in B(i)} VG_a^s(t) = \sum_{b \in A(i)} UG_a^s(t), \forall i \in \mathcal{N}_G, \forall s, t \quad (7a) $$

$$ \sum_{a \in B(i)} V E_a^{s,e}(t) = \sum_{b \in A(i)} U E_a^{s,e}(t), \forall i \in \mathcal{N}_G, \forall e \in \mathcal{E}, \forall s, t \quad (7b) $$

where $A(i)/B(i)$ represents the set of links whose tail/head is at node $i$.

For EVs, the existing occupancy on the charging link $a$ should be less than or equal to the number of charging piles on that link:

$$ \sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [UE_a^{s,e}(t) - VE_a^{s,e}(t)] \le NC_a(t), \forall a \in \mathcal{A}_C, \forall t \quad (8) $$

where $NC_a(t)$ is the actual number of charging piles for charging link $a$ at time step $t$.

The update process of the EV's EL at the FCS is described as follows:

$$ V E_a^{s,e}(t) = 0, \forall a \in \mathcal{A}_C, \forall e \in \{e < \alpha_a^t\}, \forall s, t \quad (9a) $$

$$ V E_a^{s,e}(t) \le U E_a^{s,e-\alpha_a^t}(t-1), \quad \forall a \in \mathcal{A}_C, \forall e \in \{\alpha_a^t \le e < E\}, \forall s, t \quad (9b) $$

$$ V E_a^{s,E}(t) \le \sum_{l=1}^{\alpha_a^t} U E_a^{s,E-l}(t-1), \forall a \in \mathcal{A}_C, \forall s, t \quad (9c) $$

where $\alpha_a^t$ represents the number of ELs that the EVs can recharge in charging link $a$ within time $\delta$. When an EV is charged at the FCS for one time step, its EL cannot be less than $\alpha_a^t$, corresponding to (9a). After the EV is charged, its EL should be less than the maximum EL, as shown in (9c). (9b) describes the charging process of EVs whose ELs are not in
the boundary condition. The cumulative number of EVs and
GVs are nonnegative and nondecreasing:

$$
VG_a^s(t) - VG_a^s(t-1) \ge 0, \forall a \in \mathcal{A}, \forall s, t \quad (10a)
$$

$$
VE_{a}^{s,e}(t) - VE_{a}^{s,e}(t-1) \ge 0, \forall a \in \mathcal{A}, \forall e \in \mathcal{E}, \forall s, t \quad (10b)
$$

$$
UG_a^s(t) - UG_a^s(t-1) \ge 0, \forall a \in \mathcal{A}, \forall s, t \quad (10c)
$$

$$
UE_{a}^{s,e}(t) - UE_{a}^{s,e}(t-1) \ge 0, \forall a \in \mathcal{A}, \forall e \in \mathcal{E}, \forall s, t \quad (10d)
$$

The number of cumulative EVs and GVs entering each link and leaving each link at time step 0 is initialized to 0 by the constraints below:

$$
UG_a^s(0) = VG_a^s(0) = 0, \forall a \in \mathcal{A}, \forall s \quad (11a)
$$

$$
UE_{a}^{s,e}(0) = VE_{a}^{s,e}(0) = 0, \forall a \in \mathcal{A}, \forall e \in \mathcal{E}, \forall s \quad (11b)
$$

The constraints above provide a set of basic feasible solu-
tions for traffic flow assignment:

$$
\Omega = \{\mathbf{x}|\text{Subject to (1) - (11)}\} \quad (12)
$$

where **x** represents the vector of a mixed flow pattern,
**x** = [[**U**G<sub>a</sub><sup>s</sup>(t), **V**G<sub>a</sub><sup>s</sup>(t), a ∈ **A**/{**A**<sub>C</sub>}, s ∈ **N**<sub>SR</sub>, t ∈ **T**], [**U**E<sub>a</sub><sup>s,e</sup>(t), **U**E<sub>a</sub><sup>s,e</sup>(t), a ∈ **A**, s ∈ **N**<sub>SR</sub>, e ∈ **E**, t ∈ **T**]].

The objective of the SODTA is to minimize the total travel time of both GVs and EVs in the ERN. In the link-based SODTA model, the total travel time is determined via the time steps when each vehicle stops on each link. The SODTA model of the ERN is formulated as follows:

$$
\begin{equation}
\begin{split}
\min_{x \in \Omega} \bar{\omega} = & \sum_{s \in \mathcal{N}_{SR}} \sum_{t \in \mathcal{T}} \sum_{a \in \mathcal{A} / \{\mathcal{A}_C, \mathcal{A}_S\}} [UG_a^s(t) - VG_a^s(t)] \\
& + \sum_{s \in \mathcal{N}_{SR}} \sum_{t \in \mathcal{T}} \sum_{a \in \mathcal{A} / \mathcal{A}_S} \sum_{e \in \mathcal{E}} [UE_a^{s,e}(t) - VE_a^{s,e}(t)]
\end{split}
\tag{13}
\end{equation}
$$

where $\bar{\omega}$ is the objective function of the SODTA problem,
representing the total system travel time.

B. Attacker-Defender Model

We make the following basic assumptions: 1) In the ERN,
we only consider the damage of FCSs since their damage tends
to have severe consequences. 2) If an FCS is attacked, none of
its charging piles can provide any service; therefore, we can
use a set of binary variables to indicate whether an FCS has
been attacked. Under these assumptions, we can rewrite the
constraint of the number of charging piles in the charging
links (8) as follows:

$$
\sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [UE_a^{s,e}(t) - VE_a^{s,e}(t)] \le \mu_a \cdot NC_a(t), \\
\forall a \in \mathcal{A}_C, \forall t \tag{14}
$$

where $\mu_a, \forall a \in \mathcal{A}_C$ represents the state of FCS $a$ and it is a binary variable $\mu_a$ equal to 0 if FCS $a$ is attacked, and 1 otherwise; $\boldsymbol{\mu} \in \mathcal{M}$ is the attack strategy for the FCSs in the entire ERN; $\mathcal{M} = \{\boldsymbol{\mu} | \sum_{a \in \mathcal{A}_C} (1 - \mu_a) \le n, \boldsymbol{\mu} \in \{0, 1\}^{|\mathcal{A}_C|}\}$ is the set of all attack strategies for the ERN, where $n$ represents the attack resource level, i.e., the number of FCSs that can be attacked at most. The feasible domain of the modified SODTA,

being the inner level problem of the AD model, can be attained
via substituting (8) with (14):

$$
\Omega(\boldsymbol{\mu}) = \{\boldsymbol{x} | \text{Subject to (1) -- (7), (9) -- (11), (14)}\} \quad (15)
$$

Then, the inner level modified SODTA problem aims to
minimize the total travel time of the ERN under an unknown
attack strategy **μ**, and the outer problem is to select an attack
strategy **μ** from the set of attack strategies **M** to maximize
the objective function. This AD model can be formulated as
follows:

$$
\begin{equation}
\begin{split}
& \max_{\mu \in \mathcal{M}} \min_{x \in \Omega(\mu)} \bar{\omega}(\mu) \\
& = \sum_{s \in \mathcal{N}_{SR}} \sum_{t \in \mathcal{T}} \sum_{a \in \mathcal{A} / \{\mathcal{A}_C, \mathcal{A}_S\}} [UG_a^s(t) - VG_a^s(t)] \\
& \quad + \sum_{s \in \mathcal{N}_{SR}} \sum_{t \in \mathcal{T}} \sum_{a \in \mathcal{A} / \mathcal{A}_S} \sum_{e \in \mathcal{E}} [UE_a^{s,e}(t) - VE_a^{s,e}(t)]
\end{split}
\tag{16}
\end{equation}
$$

By solving the aforementioned AD model, we can identify
the worst-case attack strategy **μ** and the vulnerable FCSs,
thereby achieving our first objective of identifying vulnerable
FCSs within the ERN.

C. Solution Method for the Attacker-Defender Model

The bi-level mixed integer programming problem (16) is
solved using the following method. First, we reformulate the
original problem into a single-level mixed integer problem
with a nonlinear term in the objective function (17). Then,
a linearization technique is used to linearize it in the objec-
tive function of the reformulated single-level problem, which
renders this problem an MILP. Lastly, an accelerated Benders
decomposition is applied to solve this MILP [39], [40]. The
main procedures of the proposed method are illustrated in
Fig. 1. The inner layer problem is an LP, and with the strong
duality theorem, we can obtain the dual inner level of the bi-
level problem:

$$
\begin{equation}
\begin{split}
& \max_{\mu \in \mathcal{M}} \max_{y \in \Omega^{-1}} \vartheta(\mu) \\
& = \sum_{t \in \mathcal{T}} \left\{ \sum_{a \in \mathcal{A}_R} \sum_{s \in \mathcal{N}_{SR}} [g \lambda_a^s(t) \cdot DG_a^s(t)] \right. \\
& \qquad + \sum_{a \in \mathcal{A}_R} \sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [e \lambda_a^{s,e}(t) \cdot DE_a^{s,e}(t)] \\
& \qquad + \sum_{a \in \mathcal{A} / \{\mathcal{A}_C\}} [\pi_a^1(t) \cdot f_a^O(t) + \pi_a^2(t) \cdot f_a^I(t)] \\
& \qquad + \pi_a^3(t) \cdot N_a(t)] + \sum_{a \in \mathcal{A}_C} [\pi_a^4(t) \cdot \mu_a \cdot NC_a(t)]
\end{split}
\tag{17}
\end{equation}
$$

where $g\lambda_a^s(t)$ is the dual variable of the constraint of satisfying
GV traffic demand, corresponding to (6a) and is a free con-
tinuous variable; $e\lambda_a^{s,e}(t)$ is the dual variable of the constraint
of satisfying EV traffic demand, corresponding to (6b), and
is a free continuous variable; $\pi_a^1(t)$, $\pi_a^2(t)$, $\pi_a^3(t)$, $\pi_a^4(t)$ are
the dual variables for the maximum outflow constraint of
road traffic flow in (1), the maximum inflow constraint of
road traffic flow in (2), the capacity constraint of road traffic
flow in (3), and the number of charging piles in charging
stations constraint in (8), respectively; they are all non-positive
<figure><img src="image_2.png" /><figcaption>Fig. 1. Flowchart of the solution technique for the attacker-defender problem.</figcaption></figure>

continuous variables. $y$ is the vector of the dual variables of the original problem. It can be observed that the objective function contains a nonlinear term consisting of binary and continuous variables, which can be linearized by introducing a new set of variables $\sigma_a(t)$ and a set of linearization constraints:

$$ \sigma_a(t) \ge \pi_a^4(t), \forall a \in \mathcal{A}_C, \forall t \quad (18a) $$

$$ \sigma_a(t) \ge M \cdot \mu_a, \forall a \in \mathcal{A}_C, \forall t \quad (18b) $$

$$ \sigma_a(t) \le \pi_a^4(t) + M \cdot (\mu_a - 1), \forall a \in \mathcal{A}_C, \forall t \quad (18c) $$

$$ \sigma_a(t) \ge M, \forall a \in \mathcal{A}_C, \forall t \quad (18d) $$

where $M$ is a predefined negative constant that represents the lower bound of the variable $\pi_a^4(t)$. As the dual variable of the constraint in (8), the economic interpretation of $\pi_a^4(t)$ is the impact on the objective function when the right-hand parameter $NC_a(t)$ of (8) changes by one unit. This impact always decreases the total travel time (objective function) of the system, so the variable takes a non-positive value. Based on the actual meaning of this variable, we can consider $\pi_a^4(t) \in [M, 0]$. When $\mu_a = 0, \forall a \in \mathcal{A}_C$, (18b) ensures that $\sigma_a(t) \ge 0, \forall a \in \mathcal{A}_C, \forall t$, and since $\sigma_a(t)$ is a non-positive continuous variable, we have $\sigma_a(t) \le 0, \forall a \in \mathcal{A}_C, \forall t$. Therefore, $\sigma_a(t) = 0, \forall a \in \mathcal{A}_C, \forall t$. When $\mu_a = 1, \forall a \in \mathcal{A}_C$, (18a) ensures that $\sigma_a(t) \ge \pi_a^4(t), \forall a \in \mathcal{A}_C, \forall t$, and (18c) ensures that $\sigma_a(t) \le \pi_a^4(t), \forall a \in \mathcal{A}_C, \forall t$. Therefore, $\sigma_a(t) = \pi_a^4(t), \forall a \in \mathcal{A}_C, \forall t$. After adding the set of linearization constraints (18), the nonlinear term $\sum_{a \in \mathcal{A}_C} [\pi_a^4(t) \cdot \mu_a \cdot NC_a(t)]$ can be equivalently replaced by the linear term $\sum_{a \in \mathcal{A}_C} [\sigma_a(t) \cdot NC_a(t)]$. Additionally, the objective function and the objective of the two levels are the same. This bi-level problem can be, therefore, transformed into a single-level MILP:

$$ \begin{aligned} \max_{(y, \mu) \in \Omega_L^{-1}} & \sum_{t \in T} \left\{ \sum_{a \in \mathcal{A}_R} \sum_{s \in \mathcal{N}_{SR}} [g \lambda_a^s(t) \cdot D G_a^s(t)] \right. \\ & + \sum_{a \in \mathcal{A}_R} \sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [e \lambda_a^{s,e}(t) \cdot D E_a^{s,e}(t)] \\ & + \sum_{a \in \mathcal{A}/\{\mathcal{A}_C\}} [\pi_a^1(t) \cdot f_a^O(t) + \pi_a^2(t) \cdot f_a^I(t) \\ & + \pi_a^3(t) \cdot N_a(t)] + \sum_{a \in \mathcal{A}_C} [\sigma_a(t) \cdot NC_a(t)] \end{aligned} \quad (19) $$

where $\Omega_L^{-1}$ is feasible solution set of the MILP. The MILP can also be written in compact form:

$$ \begin{array}{ll} \max & c^T y \\ \text{s.t.} & Ay \le b - B\mu \\ & \sum_{a \in \mathcal{A}_C} (1 - \mu_a) \le n \\ & \mu \in \{0, 1\}^{|\mathcal{A}_C|} \end{array} \quad (20) $$

where $A, B$ are the coefficient matrices of the variable $x, \mu$, respectively; $b$ is the vector of right-side parameters in the constraints; $c$ is the cost coefficient vector.

The resulting MILP (19) is often of large-scale in the context of a real-world ERN with numerous nodes and links. For its efficient solution as mentioned earlier, an accelerated Benders Decomposition (BD) is employed. BD is a well-established approach applicable to many mathematical programming problems with discrete variables [41], [42], [43]. When employing BD in the traditional way, the main problem (MP) is typically supplied with either an optimality or feasibility cut during each iteration to reduce the gap. However, it is widely admitted that the classical BD exhibits slow convergence and may require a substantial number of iterations to attain the optimal solution [44], [45]. To overcome these limitations, the BD is accelerated by generating multiple cuts using multiple solutions of MP. This has been shown effective in reducing BD's solution time [44], [46].

In the accelerated BD, the MILP (19) is first decomposed into SP and MP. The attack strategy $\mu$ solved from the MP is passed to the SP as a known parameter. Since the SP is a pure linear programming, we can obtain the dual of the SP as below:

$$ \begin{array}{ll} \min_{\bar{x} \in \Omega_{DSP}} & \bar{\omega}(\mu) = (b - B\mu)^T \bar{x} \\ \text{s.t.} & A^T \bar{x} \le c \end{array} \quad (21) $$

where $\bar{x}$ is the vector of the variables in the dual SP. In fact, the form of dual SP is similar to the SODTA problem (13) and the variable $\bar{x}$ consists of $x$ in (13) and the dual variables of the linearization constraints (18). The set of constraints for dual SP $\Omega_{DSP}$ together with $\Omega_L^{-1}$ and $\Omega_L^{-1}$ is omitted in the main text for brevity and is provided in [47]. The objective function $(b - B\mu)^T \bar{x}$ denotes $\bar{\omega}(\mu)$.

The MP together with its optimality cut can be formulated as follows:

$$ \max_{\mu \in \mathcal{M}} q \quad (22) $$

$$ \text{s.t.} \quad (b - B\mu)^T \bar{x}_o^{(i)} \ge q \quad (23a) $$

$$ (b - B\mu)^T \bar{x}_s^{(i)} \ge q \quad (23b) $$

where $q$ is an auxiliary variable representing the value of the objective function in the dual SP. For the MP, the decision variable is $\mu$, and $\bar{x}$ is a fixed parameter which is provided by the solution of the dual SP. The index $i$ indicates the number of cuts that have been added to the MP. By solving the MP, we can obtain attack strategy $\mu$ as an input for the dual SP.
**Algorithm 1** Accelerated Benders Decomposition With Multiple Cuts for Problem (19)

1: Initialize $UB \leftarrow +\infty, LB \leftarrow 0$

2: $\forall \mu_o^{(0)}, \mu_s^{(0)} \in \mathcal{M} \leftarrow \text{Any initial attack strategies}$

3: $\epsilon \leftarrow$ a very small convergence tolerance (e.g., 0.001), $i \leftarrow 0$ the current iteration number, $I_M \leftarrow$ the maximum iteration number

4: **while** $UB - LB/LB \ge \epsilon$ and $i \le I_M$ **do**

5: Solve the dual SP (21) with the optimal attack strategy $\mu_o^{(i)}$ obtained in (22)

6: Solve the dual SP (21) with the sub-optimal attack strategy $\mu_s^{(i)}$ obtained in (22)

7: Update the lower bound $LB \leftarrow \max\{LB, \bar{\omega}(\mu_o^{(i)}), \bar{\omega}(\mu_s^{(i)})\}$

8: **if** $\bar{\omega}(\mu_o^{(i)}) < LB$ or $\bar{\omega}(\mu_s^{(i)}) < LB$ **then**

9: Add constraint (24) of $\bar{\omega}(\mu_o^{(i)})/\bar{\omega}(\mu_s^{(i)})$ to the MP

10: **else**

11: Skip

12: **end if**

13: Add cuts $(b - B\mu)^T \bar{x}_o^{(i)} \ge q$ and $(b - B\mu)^T \bar{x}_s^{(i)} \ge q$ to the MP

14: Solve the MP (22)

15: Update the upper $UB \leftarrow \min\{UB, q\}$

16: Obtain the optimal attack strategy $\mu_o^{(i+1)}$ and sub-optimal attack strategy $\mu_s^{(i+1)}$ in current MP

17: Update the iteration number $i \leftarrow i + 1$

18: **end while**

Then, a cut can be generated each time we solve the dual SP with given attack strategy $\mu$. The optimal solution and sub-optimal solution of the MP are used here to generate multiple cuts for MP. Only optimality cuts are generated for the MP, because the dual SP is always feasible. The interaction process between the MP and the dual SP, as well as the solution algorithm, are provided in the algorithm (1). At the beginning of the algorithm, we first initialize the upper and lower bounds of the objective function as $UB = +\infty$ and $LB = 0$, respectively. Then, we randomly select two initial attack strategies from $\mathcal{M}$ as the initial optimal and sub-optimal solutions. We set the iteration information of the algorithm, initialize the current iteration number $i$ to 0, set the upper and lower bound convergence tolerance $\epsilon$ (usually set to 0.001), and set the maximum number of iterations $I_M$.

In this algorithm, the dual SP is solved twice every iteration with the optimal solution $\mu_o^{(i)}$ and the sub-optimal solution $\mu_s^{(i)}$ of the MP. Not only will the solution of both dual SPs be taken advantage of to generate two optimality cuts for MP, but the objective function values of both dual SPs will be compared to each other and the current lower bound $LB$. If the objective function value corresponding to one attack strategy $\mu^l$ is lower than the objective function value of another attack strategy or the lower bound, a cut is added to avoid this inferior attack strategy being obtained again in the MP, which occasionally happens in practice:

$$ \sum_{\forall a \in \mathcal{A}_C, \mu_a^l = 0} (1 - \mu_a) \le n - 1 \quad (24) $$

This constraint restricts the attack FCSs in the inferior attack strategy so that they will not be attacked simultaneously again, i.e., this inferior attack strategy is excluded in the set of attack strategies $\mathcal{M}$.

## D. Performance Indices

To assess the impact of attacks on the ERN and fulfill the second objective of this work, we introduce several indices to evaluate system performance under attacks. These indices encompass overall performance, vulnerability index, load and utilization rate of FCSs, average travel time for GV and EV drivers, and queuing time at FCSs. The overall system performance is measured by the cumulative traffic demand satisfaction rate, which is formulated as follows:

$$ P(t) = \frac{\sum_{s \in \mathcal{N}_{SR}} \sum_{a \in \mathcal{A}_S} U G_a^s(t) + \sum_{s \in \mathcal{N}_{SR}} \sum_{a \in \mathcal{A}_S} \sum_{e \in \mathcal{E}} U E_a^{s,e}(t)}{\sum_{s \in \mathcal{N}_{SR}} \sum_{a \in \mathcal{A}_R} D G_a^s(t) + \sum_{s \in \mathcal{N}_{SR}} \sum_{a \in \mathcal{A}_R} \sum_{e \in \mathcal{E}} D E_a^{s,e}(t)} \quad (25) $$

A higher $P(t)$ represents a higher proportion of traffic flow reaching its destination (i.e., the numerator in (25)) in the total traffic demand up to time step $t$ (i.e., the denominator in (25)). The performance index $P(t)$ is bounded between 0 and 1, making it easy to interpret its meaning: a value of 1 indicates complete satisfaction of traffic demand, whereas 0 signifies none met and complete breakdown of ERN. This performance index is widely used in literature, e.g. [26], [30], and [38]. The system vulnerability index (VI) can be further formulated as follows according to the system performance index and the definition of vulnerability (i.e., the sensitivity to disruptions [7]) adopted in this paper:

$$ VI(t) = \frac{P(t) - P_{\mu}(t)}{P(t)} \quad (26) $$

where $P_{\mu}(t)$ indicates the system performance index under attack strategy $\mu$.

Note that the system's VI is also a cumulative value, representing the performance loss up to time step $t$. When $VI(t) = 0$ ($P(t) = P_{\mu}(t)$), it means that the attack has had no detrimental impact on the system's performance through time step $t$; conversely, when $VI(t) = 1$ ($P_{\mu}(t) = 0$), it indicates a complete collapse of the ERN under the attack. Thus, a larger value of $VI(t)$ corresponds to greater vulnerability of the ERN to the attack strategy $\mu$ at time step $t$.

Another index is proposed to investigate the redistribution of FCS's utilization after the attack:

$$ U_a(t) = \frac{\sum_{s \in \mathcal{N}_{SR}} \sum_{e \in \mathcal{E}} [U E_a^{s,e}(t) - V E_a^{s,e}(t)]}{N C_a(t)}, \quad \forall a \in \mathcal{A}_C, \forall t \quad (27a) $$

$$ \Delta U_a(t) = U_a^{\mu}(t) - U_a(t), \quad \forall a \in \mathcal{A}_C, \forall \mu \in \mathcal{M}, \forall t \quad (27b) $$

where $U_a(t)$, $U_a^{\mu}(t)$ represent the utilization rate of FCS during time step $t$ under the normal circumstance and attack strategy $\mu$, respectively. $\Delta U_a(t)$ denotes the FCS utilization variation at time step $t$ after attack strategy $\mu$.

The charging load of each FCS and the whole system can be calculated as:

$$ \mathcal{L}_a = \sum_{t \in \mathcal{T}} \sum_{e \in \mathcal{E}} \alpha_a^t \cdot [U E_a^{s,e}(t) - V E_a^{s,e}(t)], \quad \forall a \in \mathcal{A}_C \quad (28a) $$

$$ \mathcal{L} = \sum_{a \in \mathcal{A}_C} \mathcal{L}_a \quad (28b) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Flowchart of the vulnerability assessment framework.</figcaption></figure>

The average travel time (ATT) of EVs and HVs with different destinations is used to assess the impact of attack strategy $\mu$ on the different kinds of vehicles with different destinations. The calculation of ATT and the changed ATT ratio involves the following equation:

$$ GATT_s = \frac{\sum_{a \in A \setminus \{A_S, A_C\}} \sum_{t \in T} [UG_a^s(t) - VG_a^s(t)] \cdot \delta}{\sum_{a \in A_R} DG_a^s(T_{max})}, \quad (29a) $$

$$ EATT_s = \frac{\sum_{a \in A \setminus A_S} \sum_{t \in T} \sum_{e \in \mathcal{E}} [UE_a^{s,e}(t) - VE_a^{s,e}(t)] \cdot \delta}{\sum_{a \in A_R} \sum_{e \in \mathcal{E}} DE_a^{s,e}(T_{max})}, \quad (29b) $$

where the numerator in (29) represents the total travel time of all HVs/EVs with destination $s$. The denominator in (29) represents the accumulated demand of HVs/EVs with destination $s$ up to the final time step $T_{max}$. The ATT of EVs and HVs with destination $s$ are denoted as $EATT_s$ and $GATT_s$ separately.

We employ the changed ratio of ATT of different groups of vehicles, i.e., the vehicles of the different kinds with different destinations, to assess how these vehicles of different groups are impacted under attack strategy $\mu$:

$$ \Delta GATT_s^\mu = \frac{GATT_s^\mu - GATT_s}{GATT_s}, \quad \forall s \in \mathcal{N}_{SR}, \forall \mu \in \mathcal{M} \quad (30a) $$

$$ \Delta EATT_s^\mu = \frac{EATT_s^\mu - EATT_s}{EATT_s}, \quad \forall s \in \mathcal{N}_{SR}, \forall \mu \in \mathcal{M} \quad (30b) $$

where $GATT_s^\mu/EATT_s^\mu$, $\Delta GATT_s^\mu/\Delta EATT_s^\mu$ denote the ATT of HVs/EVs under attack strategy $\mu$ and the changed ratio of the ATT of HVs/EVs under attack strategy $\mu$, respectively. These two indexes reflect the vulnerability of traffic conditions of EVs and HVs from an individual and static perspective. While the $VI(t)$ can reflect the vulnerability of traffic conditions from a systematic and dynamic perspective.

Our SODTA model does not include a queuing link at the FCS. Instead, EVs wait for charging on the general links $a \in A_G$. Therefore, the average queuing time (AQT) for EVs at FCSs can be determined by calculating the weighted difference

<figure><img src="image_4.png" /><figcaption>Fig. 3. The topology structure of the experimental ERN [38].</figcaption></figure>

in traverse time between EVs and HVs on the general links:

$$ AQT = \sum_{a \in A_G} \times \frac{\sum_{e \in \mathcal{E}} \sum_{s \in \mathcal{N}_s} \{VE_a^{s,e}(T_{max}) \cdot [e\tau_a(T_{max}) - g\tau_a(T_{max})]\}}{\sum_{e \in \mathcal{E}} \sum_{s \in \mathcal{N}_s} VE_a^{s,e}(T_{max})} \quad (31) $$

where $e\tau_a(t)$, $g\tau_a(t)$ indicate the actual traverse time of EVs and HVs at general link $a$ at time step $t$, respectively. $e\tau_a(t)$ and $g\tau_a(t)$ are derived from the cumulative vehicle numbers, i.e., the solution of the optimization model. The methodology for calculating the link traverse time is similar to [36], [48], and [49]. For brevity, the method is detailed in our Appendix [47].

### III. NUMERICAL EXAMPLE

#### A. Data Description

As a numerical case study of the proposed vulnerability assessment approach, the ERN in [38] was used, as shown in Fig. 3. The data is available in the appendix [47]. The length of time step $\delta$ is set to 6 minutes. All experiments were conducted on a computer equipped with an AMD Ryzen 7 6800H 3.20GHz-CPU and 16GB of RAM. IBM ILOG CPLEX 12.10.0 was utilized as the solver for all LP and MILP problems.

#### B. Experimental Results

1) One Fast Charging Station Attacked Scenario: Our investigation started focusing on the scenario where a single
<figure><img src="image_5.png" /><figcaption>Fig. 4. Evolution of vulnerability index under different EV penetration.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Charging stations utilization rate variety $\Delta U_a(t)$ under 50% EV penetration.</figcaption></figure>

FCS is under attack, i.e., $n = 1$. The main procedure of the assessment is given in Fig. 2. To determine the mixed traffic flow assignment and charging load distribution under normal situations, we solved the SODTA problem (12). Subsequently, we resolved the AD model’s equivalent formulation (19), identifying the optimal attack strategy and system responses in the worst-case scenario. The FCS that is attacked in the worst-case scenario is FCS #302 under an EV penetration of 25%, 50%, 75%, and 100%. The method outlined in Section II-D was employed to calculate the system VI, as illustrated in Fig. 4.

It is noticeable that as the EV penetration level increases, system VI also increases, indicating an escalating vulnerability of the system. Additionally, the evolution of system VI under four different levels of EV penetration follows a consistent pattern: the worst-case attack takes effect at time step 8; VI reaches its peak at either time step 8 or 10; subsequent to the peak, it can be observed that the system gradually absorbs the impact and recovers from it; by the last time step, the system performance nearly returns to a state of normality, particularly when EV penetration rate is 25%.

Fig. 5 provides an intuitive depiction of the redistribution of the charging load within the ERN under EV penetration rates 50%. It explains why the system can recover from the worst-case attack. It is shown in this Figure that: the charging load initially allocated to the attacked FCS #302 is redirected to the FCSs that are topologically close, primarily

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of VI under different attack resources and 50% EV penetration.</figcaption></figure>

<table><caption>TABLE I
THE CHANGED RATIO OF ATT OF GVS/EVS WITH DIFFERENT DESTINATION AFTER WORST-CASE ATTACK</caption><thead><tr><th>EV penetration</th><th>Vehicle Type</th><th>201</th><th>202</th><th>203</th><th>204</th><th>205</th><th>206</th><th>207</th><th>208</th></tr></thead><tbody><tr><td rowspan="2">0%</td><td>GV</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td></tr><tr><td>GV</td><td>0.71%</td><td>0.10%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.27%</td><td>0.00%</td><td>0.00%</td></tr><tr><td rowspan="2">25%</td><td>GV</td><td>0.00%</td><td>-0.43%</td><td>35.17%</td><td>-1.77%</td><td>5.50%</td><td>-0.80%</td><td>0.00%</td><td>14.82%</td></tr><tr><td>GV</td><td>0.00%</td><td>-0.30%</td><td>0.15%</td><td>0.41%</td><td>0.00%</td><td>-0.27%</td><td>0.00%</td><td>0.00%</td></tr><tr><td rowspan="2">50%</td><td>EV</td><td>0.00%</td><td>0.25%</td><td>40.81%</td><td>0.00%</td><td>13.91%</td><td>-1.78%</td><td>-0.94%</td><td>24.54%</td></tr><tr><td>EV</td><td>0.00%</td><td>0.10%</td><td>0.00%</td><td>0.00%</td><td>-0.25%</td><td>-2.16%</td><td>0.31%</td><td>0.00%</td></tr><tr><td rowspan="2">75%</td><td>EV</td><td>0.00%</td><td>0.25%</td><td>45.76%</td><td>-1.21%</td><td>23.09%</td><td>-1.27%</td><td>0.00%</td><td>28.94%</td></tr><tr><td>GV</td><td>~</td><td>~</td><td>~</td><td>~</td><td>~</td><td>~</td><td>~</td><td>~</td></tr><tr><td rowspan="2">100%</td><td>EV</td><td>0.00%</td><td>-0.32%</td><td>77.11%</td><td>2.02%</td><td>31.30%</td><td>-0.76%</td><td>-0.14%</td><td>33.59%</td></tr></tbody></table>

FCS #304. At time step 8, for FCS #304 we can observe a significant increase in utilization rate, which coincides with the rise in the system’s VI curve as shown in Fig. 4. This observation suggests that the EVs, originally intended to reach their destinations at time step 9, had to divert to FCS #304 for recharging due to the damage inflicted upon FCS #302. Once the EVs finish recharging in FCS #304, they resume their trip to their respective destination. This phenomenon explains why, in the subsequent time steps in Fig. 4, the overall ERN managed to withstand the attack on FCS #302. This highlights the importance of FCSs closely located to the vulnerable FCS: much nearby FCSs share the charging load in case the vulnerable FCS is damaged, ultimately reducing the overall vulnerability of the system.

2) *Multiple Fast Charging Stations Attacked Scenario:* In this part, we study the impact of worst-case attacks on the traveling of EVs and GVs. The change in ATT of GVs/EVs with different destinations after the attacks is listed in Tab. I. Depending on the extent to which the vehicles are affected, they can be divided into two categories. The first category consists of GVs and some EVs with specific destinations whose ATT is barely affected by the attack. For example, all GVs and EVs destined for nodes 201, 202, 204, 206, and 208 fall into this category. The second category, which comprises affected EVs, can be further divided into two sub-categories: 1) Affected EVs that can still reach their destinations. These EVs are able to take detours to other functional FCSs. This represents a mitigable vulnerability at the system level. 2) Affected EVs that cannot reach their destinations. Because the closest FCS is damaged, and their remaining electrical charge is insufficient to reach any other FCSs within the ERN. This represents an unmitigable vulnerability at the system level.
<table><caption>TABLE II
OPTIMAL ATTACK STRATEGIES AND SYSTEM VULNERABILITY INDEX</caption><thead><tr><th>n</th><th>μ<sub>a</sub> = 0</th><th>VI<sup>1</sup></th><th>Total Charging Load (ELs)</th></tr></thead><tbody><tr><td>0</td><td></td><td>0</td><td>3930</td></tr><tr><td>1</td><td>302</td><td>0.011</td><td>4200</td></tr><tr><td>2</td><td>302, 304</td><td>0.123</td><td>2676</td></tr><tr><td>3</td><td>302, 304, 307</td><td>0.182</td><td>2592</td></tr><tr><td>4</td><td>302, 304, 306, 307</td><td>0.201</td><td>1278</td></tr><tr><td>5</td><td>301, 302, 303, 304, 308</td><td>0.217</td><td>2061</td></tr><tr><td>6</td><td>301, 302, 303, 304, 305, 308</td><td>0.248</td><td>1473</td></tr><tr><td>7</td><td>301, 302, 303, 304, 306, 307, 308</td><td>0.291</td><td>582</td></tr></tbody></table>

As the attack resources increase, the proportions of these three categories of vehicles will change, leading to changes in the system's VI. In the following cases, we study scenarios where multiple FCSs are attacked. Fig. 6 displays the VI curve of the system with 50% EV penetration under different attack resources. We can observe that with an increase in the attack resources, the system's VI curve exhibits two different trends: 1) Transition from mitigable to unmitigable vulnerability: initially, the VI curve experiences a decrease at a specific time step; however, with increases in attack resources, it starts to rise. For example, when the attack resources increase from 1 to 2, the VI curve, which drops at time step 12, instead registers an increase. This occurs because some EVs previously able to reach their destinations by taking detours for charging are transformed into EVs that cannot reach any functional FCS. In essence, a portion of mitigable vulnerability morphs into unmitigable vulnerability. 2) Overall upward shift in VI curve: following the increase in attack resources, the VI curve experiences a collective upward shift. This trend is exemplified by the shift in the VI curve when moving from 3 to 4 attacked FCSs. As the attack resources increase from 3 to 4, the VI values for all time steps increase by approximately 0.015. This phenomenon can be explained by a portion of initially unaffected EVs now being transformed into EVs unable to reach any available FCSs for recharging and, consequently, unable to reach their destinations. In essence, new instances of unmitigable vulnerabilities emerge within the system.

In Tab. II, we present the optimal attack strategies under different attack resources, along with the corresponding VI scores for each scenario. These strategies are associated with the identification of FCSs by their average utilization rate and total charging load over the considered time steps. Three fundamental observations emerge: 1) Attack targets do not align with utilization or load: it is notable that FCSs subjected to attacks in the worst-case scenarios are not the ones with the highest average utilization rate or the heaviest charging load. For instance, when the attack resources $n = 1$, the optimal attack strategy is to impair FCS #302, which is not the FCS with the highest average utilization rate (which is FCS #308) or with the highest charging load (which is FCS #301). 2) Subset dynamics in attack strategies: generally,

<table><caption>TABLE III
TABLE OF CHARGING STATION IDs RANKED BY THEIR AVERAGE UTILIZATION RATE AND CHARGING LOAD</caption><thead><tr><th>Rank by Average Utilization Rate</th><th>Charging Station ID</th><th>Average Utilization Rate</th><th>Rank by Charging Load</th><th>Charging Station ID</th><th>Charging Load (Unit: ELs)</th></tr></thead><tbody><tr><td>1</td><td>308</td><td>15.58%</td><td>1</td><td>301</td><td>738</td></tr><tr><td>2</td><td>301</td><td>10.25%</td><td>2</td><td>307</td><td>720</td></tr><tr><td>3</td><td>307</td><td>10.00%</td><td>3</td><td>304</td><td>612</td></tr><tr><td>4</td><td>304</td><td>8.50%</td><td>4</td><td>303</td><td>576</td></tr><tr><td>5</td><td>303</td><td>8.00%</td><td>5</td><td>302</td><td>576</td></tr><tr><td>6</td><td>302</td><td>5.33%</td><td>6</td><td>308</td><td>561</td></tr><tr><td>7</td><td>306</td><td>3.83%</td><td>7</td><td>306</td><td>138</td></tr><tr><td>8</td><td>305</td><td>0</td><td>8</td><td>305</td><td>0</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 7. EVs' average queuing time under different EV penetration and attack resources.</figcaption></figure>

the optimal attack strategy for targeting $n - 1$ FCSs is the subset for the optimal attack strategy for attacking $n$ FCSs, with a notable exception when $n$ increases from 4 to 5. Obviously, $\{302, 304, 306, 307\} \not\subseteq \{301, 302, 303, 304, 308\}$. 3) Complexity of combinatorial attack strategies: the optimal strategies of attacking multiple FCSs do not equal a simple superposition of attacks on the most vulnerable FCS. For a single FCS attack, the best strategy is to attack FCS #302, as shown in Tab. II, whereas a sub-optimal strategy is to attack FCS #307. However, when targeting 2 FCSs, the optimal strategy shifts to attacking FCS #302 and #304. These findings highlight the significance of evaluating the ERN's vulnerability and identifying critical FCSs from a systemic perspective.

The vulnerability of EVs' queuing time to malicious attacks is shown in Fig. 7. From the perspective of EV penetration, the average queuing time (AQT) tends to rise as EV penetration increases, which underscores the importance of expanding the scale and size of CS to ensure the service quality of EVs at FCSs as EV penetration continues to increase. From the perspective of attack resources, we can perceive that: 1) For each level of EV penetration, the queuing time first increases as the number of attack resources increases, reaching a peak when 4 attack resources are deployed. This initial rise and peak can be attributed to the increasing disruption of the charging infrastructure, which reduces the available charging capacity and forces more EVs to wait longer for charging. 2) After this point, the queuing time declines as more attack resources are added. This counter-intuitive decrease can be explained by the limited number of FCSs survived in the ERN and fewer opportunities for EVs to charge and queue. 3) Ultimately,

<footnote n="1">¹ refers to the system VI up to the final time step, i.e., VI(20)</footnote>
<table><caption>TABLE IV
SOLUTION TIME OF THE PROPOSED AD WITH DIFFERENT METHODS UNDER 50% EV PENETRATION</caption><thead><tr><th rowspan="2">Instance</th><th colspan="2">Enumerate</th><th colspan="2">Solve the MILP<br>without Decomposition</th><th colspan="2">Traditional Benders<br>Decomposition</th><th colspan="2">Accelerated Benders<br>Decomposition</th><th rowspan="2">|M|</th></tr><tr><th>n</th><th>Time (sec)</th><th>Obj.</th><th>Time (sec)</th><th>Obj.</th><th>Time (sec)</th><th>Obj.</th><th>Time (sec)</th><th>Obj.</th></tr></thead><tbody><tr><td>1</td><td>239.55 ± 12.25<sup>1</sup></td><td>18820</td><td>N/A<sup>2</sup></td><td>18218<sup>3</sup></td><td>195.75 ± 7.71*</td><td>18820</td><td>221.23 ± 9.44</td><td>18820</td><td>8</td></tr><tr><td>2</td><td>859.49 ± 48.71</td><td>20850</td><td>N/A</td><td>19122</td><td>639.92 ± 23.12</td><td>20850</td><td>609.17 ± 22.87*</td><td>20850</td><td>28</td></tr><tr><td>3</td><td>1659.01 ± 77.62</td><td>21238</td><td>N/A</td><td>19520</td><td>1168.36 ± 41.11</td><td>21238</td><td>1019.75 ± 31.54*</td><td>21238</td><td>56</td></tr><tr><td>4</td><td>2069.37 ± 107.94</td><td>21668</td><td>N/A</td><td>21604</td><td>1455.06 ± 52.53</td><td>21668</td><td>1137.22 ± 40.77*</td><td>21668</td><td>70</td></tr><tr><td>5</td><td>1671.58 ± 83.24</td><td>23238</td><td>N/A</td><td>21686</td><td>1207.29 ± 42.57</td><td>23238</td><td>1009.96 ± 29.97*</td><td>23238</td><td>56</td></tr><tr><td>6</td><td>853.65 ± 41.78</td><td>23606</td><td>N/A</td><td>21738</td><td>649.0 ± 24.38</td><td>23606</td><td>611.33 ± 21.85*</td><td>23606</td><td>28</td></tr><tr><td>7</td><td>242.52 ± 11.79</td><td>24044</td><td>N/A</td><td>22490</td><td>187.43 ± 7.51*</td><td>24044</td><td>222.46 ± 9.96</td><td>24044</td><td>8</td></tr></tbody></table>

when the attack resources reach 8, the queuing time drops to 0, as all FCSs in the ERN become malfunctioned.

### C. Solution Efficiency

To verify the effectiveness and accuracy of our proposed accelerated Benders Decomposition, we utilized several methods to solve the model. These methods include: 1) Directly solve the two-layer AD model (16) without reformulation (enumerate). Due to the intractable two-layer structure, it is inviable to directly solve the AD model with a solver. Hence, we enumerate over the set of attack strategies $M$ to yield the solution of the AD model. This process may be time-consuming, but it guarantees that we get the optimal solution as we search throughout the entire set $M$. 2) Directly solve the equivalent reformulated MILP (19) without decomposition. 3) Decompose the MILP (19) and solve it via traditional Benders Decomposition. 4) Decompose the MILP (19) and solve it via proposed accelerated Benders Decomposition. Using the aforementioned methods solving the proposed model repetitively for 10 times, we can yield the result listed in Tab.IV. The optimization problem (16) of ERN displayed in Fig.3 comprises a significantly large number of variables (397, 153) and constraints (684, 057). The problem reduced by CPLEX has a number of non-zeros in the lower triangle of $A \times A^T = 2,064,690$. Solving such a problem with fixed integer variables takes about 30 seconds using CPLEX. In comparison to the relevant literature involving the AD model [25], [26], our approach demonstrates high complexity and large-scale trait. This is because our model is sophisticated at detailing the specific features of ERNs. The enumerate method consistently took the longest time across all instances, with the time generally proportional to the size of the set of attack strategies $|M|$, as the enumerate method searches all over the set $M$. Surprisingly, it is rather time-consuming to directly solve the reformulated MILP (19), whose solution time exceeds the set time limit. This is likely because the large scale, high dimensionality, and a substantial number of binary variables of the reformulated MILP (19) is tricky for the solver. Nonetheless, by decomposing the reformulated MILP and employing Benders Decomposition, we observe improved computational time, which emphasizes the superiority of reformulating and simplifying the proposed AD model. Accelerated Benders Decomposition outperforms the traditional one except for the

instances where $n = 1$ and $n = 7$, especially when the size of the set of attack strategies is large (like $n = 4$). The correctness of accelerated Benders Decomposition can be validated by the observation that the objective function obtained by accelerated Benders Decomposition is consistent with the one yielded by the enumerate method.

## IV. CONCLUSION

This paper proposes a novel framework to identify critical fast charging stations (FCSs) within an electrified road network (ERN) and assess the ERN's vulnerability under FCS failures, a consideration often overlooked in existing literature. Based on a system optimal dynamic mixed traffic flow assignment (SODTA), we formulated an attack-defender (AD) model and reformulated it into an equivalent mixed integer linear programming (MILP). This MILP is further solved via an accelerated Benders decomposition, incorporating a multiple-cut strategy to enhance convergence. Three metrics were proposed for assessing the vulnerability, focusing on the charging load redistribution and the impacts on various vehicle groups. The framework's effectiveness is demonstrated through an example regarding an ERN in North Carolina. Key findings include that: 1) When the attack budget is limited, the vulnerability of ERN can be partially mitigated by itself via optimally assigning the mixed traffic flow and redistributing the charging load in the ERN. 2) When the attack budget increases, the vulnerability index (VI) curve changes in two distinct manners. 3) The importance of FCSs is not necessarily correlated with their historical usage, suggesting that relying solely on usage data may not reliably identify critical FCSs. 4) The queuing time of EVs at FCSs becomes longer as the EV penetration increases. The attack strategy $\mu$ obtained from solving the AD model guides us in hardening the FCS within the ERN to reduce the system vulnerability and to enhance its resilience from a preparedness perspective. The redistributing mixed traffic flow $x$ holds promise for its use in managing responses to given FCS failure events.

While the SODTA model does have the capability to represent the dynamic features of the FCS and EVs with a certain level of accuracy, the computational efficiency of this model could be further improved. On the one hand, the computational efficiency of the proposed method can be improved by increasing the length of the time step or

<footnote n="1">¹ The error is quantified by the standard deviation of the time taken to solve the problem repeatedly for 10 times.</footnote>

<footnote n="2">² Not available, because this method exceeds the time limit 5000 seconds we set.</footnote>

<footnote n="3">³ This listed the objective function obtained when the solver terminate at the time limit (maybe not optimal).</footnote>

<footnote>* Method with the highest speed under given instance.</footnote>
decreasing the size of the set of energy levels to reduce the number of variables. On the other hand, the gap tolerance could be set larger to reduce computation time. If we have a limiting requirement in the computational time and are not too sensitive to the solution accuracy, the two strategies mentioned above can be applied. To enhance the computational efficiency without loss of accuracy, the unique block-angular structure inherent in the SODTA model could be utilized. The Dantzig-Wolfe decomposition [50] can be attempted to decompose the SODTA into a MP and two sub-problems. Parallel computation could be further used to solve these two sub-problems simultaneously.

In future research, the following aspects will be considered: 1) Improving the solution time of the SODTA model through reformulation and decomposition. 2) Extending the framework to three layers to incorporate considerations of optimal design to protect ERNs. Adaptive robust optimization [51] and stochastic optimization [52] can be explored for this. 3) Investigating the interplay between power grid topology and traffic network dynamics, which is not yet considered in this study, to understand their collective impact on overall system vulnerability.

## REFERENCES

[1] E. M. Bibra et al., "Global EV outlook 2022: Securing supplies for an electric future," Int. Energy Agency, Paris, France, Rep., 2022. [Online]. Available: https://iea.blob.core.windows.net/assets/ad8fb04c-4f75-42fc-973a-6e54c8a4449a/GlobalElectricVehicleOutlook2022.pdf

[2] A. Sanghvi and T. Markel, "Cybersecurity for electric vehicle fast-charging infrastructure," in Proc. IEEE Transp. Electrific. Conf. Expo (ITEC), Jun. 2021, pp. 573-576.

[3] J. Zhao, H. Zhang, Q. Lu, C. Xu, and G. Yang, "Lightning protection of electric vehicle charging infrastructure," in Proc. 33rd Int. Conf. Lightning Protection (ICLP), Sep. 2016, pp. 1-4.

[4] G. Raman, G. Raman, and J. C.-H. Peng, "Resilience of urban public electric vehicle charging infrastructure to flooding," Nature Commun., vol. 13, no. 1, p. 3213, Jun. 2022.

[5] T. Papilloud and M. Keiler, "Vulnerability patterns of road network to extreme floods based on accessibility measures," Transp. Res., D, Transp. Environ., vol. 100, Nov. 2021, Art. no. 103045.

[6] J. Lu, D. Tang, L. Hong, J. Lu, and C. Wei, "Maximum flow-based vulnerability assessment of coupled power-traffic systems," in Proc. 4th Int. Conf. Electr. Eng. Control Technol. (CEECT), Dec. 2022, pp. 642-648.

[7] K. Berdica, "An introduction to road vulnerability: What has been done, is done and should be done," Transp. Policy, vol. 9, no. 2, pp. 117-127, Apr. 2002.

[8] T. Wen, Q. Gao, Y.-W. Chen, and K. H. Cheong, "Exploring the vulnerability of transportation networks by entropy: A case study of Asia-Europe maritime transportation network," Rel. Eng. Syst. Saf., vol. 226, Oct. 2022, Art. no. 108578.

[9] X. Yang, A. Chen, B. Ning, and T. Tang, "Measuring route diversity for urban rail transit networks: A case study of the Beijing metro network," IEEE Trans. Intell. Transp. Syst., vol. 18, no. 2, pp. 259-268, Feb. 2017.

[10] M. Pirani, S. Baldi, and K. H. Johansson, "Impact of network topology on the resilience of vehicle platoons," IEEE Trans. Intell. Transp. Syst., vol. 23, no. 9, pp. 15166-15177, Sep. 2022.

[11] H. Wang, Y.-P. Fang, and E. Zio, "Risk assessment of an electrical power system considering the influence of traffic congestion on a hypothetical scenario of electrified transportation system in new York state," IEEE Trans. Intell. Transp. Syst., vol. 22, no. 1, pp. 142-155, Jan. 2021.

[12] D. Said, M. Elloumi, and L. Khoukhi, "Cyber-attack on P2P energy transaction between connected electric vehicles: A false data injection detection based machine learning model," IEEE Access, vol. 10, pp. 63640-63647, 2022.

[13] K. H. Thompson and H. T. Tran, "Operational perspectives into the resilience of the U.S. air transportation network against intelligent attacks," IEEE Trans. Intell. Transp. Syst., vol. 21, no. 4, pp. 1503-1513, Apr. 2020.

[14] Y. Tao, J. Qiu, S. Lai, X. Sun, and J. Zhao, "Vulnerability assessment of coupled transportation and multi-energy networks considering electric and hydrogen vehicles," IEEE Trans. Intell. Transp. Syst., vol. 24, no. 11, pp. 12614-12626, Nov. 2023.

[15] X. Chen and Y. Deng, "Evidential software risk assessment model on ordered frame of discernment," Expert Syst. Appl., vol. 250, Sep. 2024, Art. no. 123786.

[16] Q. Zhou, W. Pedrycz, and Y. Deng, "Order-2 probabilistic information fusion on random permutation set," IEEE Trans. Knowl. Data Eng., vol. 37, no. 2, pp. 837-850, Feb. 2025, doi: 10.1109/TKDE.2024.3484009.

[17] J. Deng, Y. Deng, and J.-B. Yang, "Random permutation set reasoning," IEEE Trans. Pattern Anal. Mach. Intell., vol. 46, no. 12, pp. 10246-10258, Dec. 2024, doi: 10.1109/TPAMI.2024.348349.

[18] Y. Wang, Z. Li, and Y. Deng, "A new orthogonal sum in random permutation set," Fuzzy Sets Syst., vol. 490, Aug. 2024, Art. no. 109034.

[19] J. C. Smith and Y. Song, "A survey of network interdiction models and algorithms," Eur. J. Oper. Res., vol. 283, no. 3, pp. 797-811, Jun. 2020.

[20] H. Chen et al., "A novel game-based belief rule base," Expert Syst. Appl., vol. 254, Nov. 2024, Art. no. 124374.

[21] B. Zeng and L. Zhao, "Solving two-stage robust optimization problems using a column-and-constraint generation method," Oper. Res. Lett., vol. 41, no. 5, pp. 457-461, Sep. 2013.

[22] R. Fakhry, E. Hassini, M. Ezzeldin, and W. El-Dakhakhni, "Tri-level mixed-binary linear programming: Solution approaches and application in defending critical infrastructure," Eur. J. Oper. Res., vol. 298, no. 3, pp. 1114-1131, May 2022.

[23] B. Zeng and Y. An, "Solving bilevel mixed integer program by reformulations and decomposition," Optim. Online, vol. 165, pp. 1-34, Jun. 2014.

[24] T. Zhao, Z. Li, and Y. Deng, "Linearity in Deng entropy," Chaos, Solitons Fractals, vol. 178, Jan. 2024, Art. no. 114388.

[25] H. Wen-Jun and Z. Xi-Zhao, "A multi-mode two-level transportation model based on network vulnerability," in Proc. 5th Int. Conf. Transp. Inf. Saf. (ICTIS), Jul. 2019, pp. 104-110.

[26] M. Ouyang and Y. Fang, "A mathematical framework to optimize critical infrastructure resilience against intentional attacks," Comput.-Aided Civil Infrastruct. Eng., vol. 32, no. 11, pp. 909-929, Nov. 2017.

[27] N. Liu, X. Hu, L. Ma, and X. Yu, "Vulnerability assessment for coupled network consisting of power grid and EV traffic network," IEEE Trans. Smart Grid, vol. 13, no. 1, pp. 589-598, Jan. 2022.

[28] M. A. Sayed, R. Atallah, C. Assi, and M. Debbabi, "Electric vehicle attack impact on power grid operation," Int. J. Electr. Power Energy Syst., vol. 137, May 2022, Art. no. 107784.

[29] M. A. Sayed, M. Ghafouri, R. Atallah, M. Debbabi, and C. Assi, "Grid chaos: An uncertainty-conscious robust dynamic EV load-altering attack strategy on power grid stability," Appl. Energy, vol. 363, Jun. 2024, Art. no. 122972.

[30] H. Wang, A. F. Abdin, Y. Fang, and E. Zio, "Resilience assessment of electrified road networks subject to charging station failures," Comput.-Aided Civil Infrastruct. Eng., vol. 37, no. 3, pp. 300-316, Mar. 2022.

[31] J. G. Wardrop, "Road paper. Some theoretical aspects of road traffic research," Proc. Inst. Civil Eng., vol. 1, no. 3, pp. 325-362, May 1952.

[32] W. Wei, L. Wu, J. Wang, and S. Mei, "Network equilibrium of coupled transportation and power distribution systems," IEEE Trans. Smart Grid, vol. 9, no. 6, pp. 6764-6779, Nov. 2018.

[33] D. Said, "A decentralized electricity trading framework (DETF) for connected EVs: A blockchain and machine learning for profit margin optimization," IEEE Trans. Ind. Informat., vol. 17, no. 10, pp. 6594-6602, Oct. 2021.

[34] C. F. Daganzo, "The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory," Transp. Res. B, Methodol., vol. 28, no. 4, pp. 269-287, Aug. 1994.

[35] A. K. Ziliaskopoulos, "A linear programming model for the single destination system optimum dynamic traffic assignment problem," Transp. Sci., vol. 34, no. 1, pp. 37-49, Feb. 2000.

[36] I. Yperman, "The link transmission model for dynamic network loading," Ph.D. dissertation, Dept. Civil Eng., Katholieke Universiteit Leuven, Leuven, Belgium, 2007.

[37] J. Long and W. Y. Szeto, "Link-based system optimum dynamic traffic assignment problems in general networks," Oper. Res., vol. 67, no. 1, pp. 167-182, Jan. 2019.

[38] H. Wang, Y.-P. Fang, and E. Zio, "Resilience-oriented optimal post-disruption reconfiguration for coupled traffic-power systems," Rel. Eng. Syst. Saf., vol. 222, Jun. 2022, Art. no. 108408.
[39] T. L. Magnanti and R. T. Wong, "Accelerating benders decomposition: Algorithmic enhancement and model selection criteria," *Oper. Res.*, vol. 29, no. 3, pp. 464–484, Jun. 1981.

[40] S. Lumbreras and A. Ramos, "How to solve the transmission expansion planning problem faster: Acceleration techniques applied to benders' decomposition," *IET Gener., Transmiss. Distrib.*, vol. 10, no. 10, pp. 2351–2359, Jul. 2016.

[41] K. N. Reddy, A. Kumar, A. Choudhary, and T. C. E. Cheng, "Multi-period green reverse logistics network design: An improved benders-decomposition-based heuristic approach," *Eur. J. Oper. Res.*, vol. 303, no. 2, pp. 735–752, Dec. 2022.

[42] J. Yin, F. Pu, L. Yang, A. D'Ariano, and Z. Wang, "Integrated optimization of rolling stock allocation and train timetables for urban rail transit networks: A benders decomposition approach," *Transp. Res., B, Methodol.*, vol. 176, Oct. 2023, Art. no. 102815.

[43] F. Alkaabneh, A. Diabat, and H. O. Gao, "Benders decomposition for the inventory vehicle routing problem with perishable products and environmental costs," *Comput. Oper. Res.*, vol. 113, Jan. 2020, Art. no. 104751.

[44] N. B. Asl and S. A. MirHassani, "Accelerating benders decomposition: Multiple cuts via multiple solutions," *J. Combinat. Optim.*, vol. 37, no. 3, pp. 806–826, Apr. 2019.

[45] C. I. Fábián and Z. Szoke, "Solving two-stage stochastic programming problems with level decomposition," *Comput. Manag. Sci.*, vol. 4, pp. 313–353, Jul. 2007. [Online]. Available: https://api.semanticscholar.org/CorpusID:1689822

[46] Z. Zhao, L. Fan, and Z. Han, "Optimal data center energy management with hybrid quantum-classical multi-cuts benders' decomposition method," *IEEE Trans. Sustain. Energy*, vol. 15, no. 2, pp. 847–858, Apr. 2024.

[47] K. Zhang, Y.-P. Fang, H. Wang, Y. Zhuang, and E. Zio. (2023). *Appendix of Vulnerability Assessment of Charging Stations in the Electrified Road Network*. [Online]. Available: https://github.com/moriyasukowa/Vulnerability-Assessment-of-Charging-Stations-in-the-Electrified-Road-Network-Appendix

[48] L. Han, S. Ukkusuri, and K. Doan, "Complementarity formulations for the cell transmission model based dynamic user equilibrium with departure time choice, elastic demand and user heterogeneity," *Transp. Res. B, Methodol.*, vol. 45, no. 10, pp. 1749–1767, Dec. 2011.

[49] S. V. Ukkusuri, L. Han, and K. Doan, "Dynamic user equilibrium with a path based cell transmission model for general traffic networks," *Transp. Res. B, Methodol.*, vol. 46, no. 10, pp. 1657–1684, Dec. 2012.

[50] G. B. Dantzig and P. Wolfe, "The decomposition algorithm for linear programs," *Econometrica*, vol. 29, no. 4, pp. 767–778, Oct. 1961.

[51] Y.-P. Fang and E. Zio, "An adaptive robust framework for the optimization of the resilience of interdependent infrastructures under natural hazards," *Eur. J. Oper. Res.*, vol. 276, no. 3, pp. 1119–1136, 2019.

[52] K. Wang, Z. Xu, Y. Liu, and Y. Fang, "Resilience enhancement for multistate interdependent infrastructure networks: From a preparedness perspective," *IEEE Trans. Rel.*, vol. 72, no. 1, pp. 190–203, Mar. 2023.

<figure><img src="image_9.png" /></figure>

**Hongping Wang** (Member, IEEE) received the Ph.D. degree in computer science from Centrale-Supélec, Université Paris-Saclay, France, in 2021. She is currently a Lecturer with the School of Intelligent Engineering and Automation, Beijing University of Posts and Telecommunications, Beijing, China. Her research interests concentrate on the modeling and analyzing of complex networked systems, particularly focusing on systems' vulnerability and risk analysis.

<figure><img src="image_10.png" /></figure>

**Yufeng Zhuang** received the Ph.D. degree in signal and information processing engineering from Beijing University of Posts and Telecommunications in 2010. He is currently a Full Professor with the School of Intelligent Engineering and Automation, Beijing University of Posts and Telecommunications. He has led and participated in over 40 research projects, including several national science and technology support programs, such as National Key Basic Research and Development Plans (973 projects), National Natural Science Fund,

and Beijing Natural Science Fund. He has published over 30 papers in prestigious international journals and applied for or holds more than ten patents. His main research directions include intelligent robotics technology, logistics, informatization and IoT technology and key technologies and applications for smart logistics and flexible supply chain.

**Ke Zhang** is currently pursuing the degree with Beijing University of Posts and Telecommunications. His research interests include modeling and vulnerability assessment of transportation systems and optimization under uncertainty.

<figure><img src="image_11.png" /></figure>

**Yi-Ping Fang** (Member, IEEE) received the Ph.D. degree in industrial engineering from the École Centrale Paris, Paris, France, in 2015. He is currently a Full Professor with the Chair of Risk and Resilience of Complex Systems, Laboratoire Génie Industriel, CentraleSupélec, Université Paris-Saclay, Paris. He was the Post-Doctoral Research Fellow with ETH, Zürich, Switzerland, from 2015 to 2017. His research interests include the study and development of advanced computational methods for risk, reliability, and resilience analytics of critical

cyber-physical systems (including smart grids, intelligent transportation, and 5G-and-beyond-systems), stochastic and robust optimization, risk and decision analysis, and machine learning.

<figure><img src="image_12.png" /></figure>

**Enrico Zio** (Fellow, IEEE) received the M.Sc. degree in nuclear engineering from the Politecnico di Milano, Italy, in 1991, the M.Sc. degree in mechanical engineering from the University of California at Los Angeles, CA, USA, in 1995, and the Ph.D. degree in nuclear engineering from Politecnico di Milano, in 1996, and the Ph.D. degree in probabilistic risk assessment from Massachusetts Institute of Technology, Cambridge, MA, USA, in 1998. He is currently a Full Professor with Mines Paris-PSL, Centre de Recherche sur les Risques et les Crises

(CRC), Sophia Antipolis, France; a Full Professor and the President of the Alumni Association, Politecnico di Milano; an Eminent Scholar with Kyung Hee University, Republic of Korea; Distinguished Guest Professor with Tsinghua University, Beijing, China; an Adjunct Professor with the City University of Hong Kong, Beihang University, and Wuhan University, China; and the Co-Director of the Center for Reliability and Safety of Critical Infrastructures (CRESCI) and the Sino-French Laboratory of Risk Science and Engineering (RISE), Beihang University. His research interests include the modeling of the failure-repair-maintenance behavior of components and complex systems, for the modeling and analysis of their reliability, maintainability, safety, vulnerability, and resilience characteristics, by monte carlo simulation methods, artificial intelligence techniques, and optimization heuristics.