Distance-based critical node detection and its application to network vulnerability assessment

Abstract

A popular strategy for network vulnerability assessment is to identify critical nodes, whose deletion maximally degrades the connectivity of the original network. Decision makers usually have information about the network structure they intend to obtain, but do not know the number of nodes to target. In this context, we study a distance-based critical node problem with unknown budget, known as the \beta-distance-based vertex disruptor problem. To solve it, we first derive three integer linear programming formulations, i.e., recursive, triangular connectivity-based, and reduced path-based ones. They are then solved using the CPLEX solver. To approximately solve large instances that an exact solver fails to solve, we propose a simple and effective integrated strategic oscillation search that combines hill climbing and strategic oscillation. It examines a large neighborhood by alternating between a tabu-enhanced destruction procedure and a random order-based improving construction procedure. Extensive experiments on both real-world and synthetic benchmark instances demonstrate the advantages of the proposed formulations and heuristic. In particular, the reduced path-based formulation outperforms both recursive and triangular connectivity-based ones. The proposed heuristic is significantly faster than exact and state-of-the-art algorithms. Finally, we perform a case study on air transportation networks to gain insights into the impacts of COVID-19.

Graphical abstract

Keywords

network vulnerability assessment / critical node detection / integer programming / hill climbing / strategic oscillation

1 Introduction

Network (or graph) data structures are widely used to describe numerous types of interactive systems, such as transportation, communication, electrical, supply chain, biological, and social networks. For such networks, decision makers are interested in evaluating their vulnerability and robustness properties to avoid or thwart adversarial attacks, natural disasters, disease epidemics, and random failures (Zhou et al., 2021; Liu et al., 2023b; Wang et al., 2025). An effective strategy for assessing network vulnerability is based on the concept of “critical nodes,” which are defined as a subset of nodes whose removal maximally degrades the connectivity of the original network, or maximally fragments the network according to some predefined measures. The resulting optimization problem is known as a critical node detection problem (CNDP) (Arulselvan et al., 2009; Zhou et al., 2024). In some cases, the information about the number of nodes to be deleted (i.e., k) is known, and the goal is to optimize the connectivity metric under a budgetary constraint that limits the number of nodes that can be removed from the original graph. This type of CNDP is known as k-vertex CNDP (Arulselvan et al., 2009; Zhou et al., 2024).

However, in many cases, we do not know the number of nodes that should be deleted, but we do have information about the network structure we need to obtain. The relevant problems are studied as \beta-connectivity CNDPs (Dinh et al., 2010; Dinh and Thai, 2015; Furini et al., 2020; Liu et al., 2023a). Given a graph G=(V,E) and a connectivity measure \sigma, the \beta-connectivity CNDPs are to minimize the number of deleted nodes, such that the connectivity measure \sigma of the residual graph is bounded by a given threshold \beta (). The \beta-connectivity CNDP represents a category of CNDPs, including a \beta-vertex disruptor problem (\beta-VDP) (Dinh et al., 2010; Dinh and Thai, 2015), \alpha-separator problem (Zhou et al., 2023b), cardinality constrained critical node problem (CC-CNP) (Arulselvan et al., 2011; Zhou et al., 2019; Liu et al., 2023a), and component cardinality-constrained critical node problem (3C-CNP) (Lalou et al., 2016).

The connectivity of a network is usually evaluated by the number of node pairs that are connected via some paths. For example, the well-known critical node problem (CNP) (Arulselvan et al., 2009) aims to minimize a pairwise connectivity measure of a residual network via the removal of a limited subset of nodes from the original network. According to Salemi and Buchanan (2022), the requirement that a path exists between two nodes is preferably accompanied by the requirement that the path length should be short. This is especially true for social networks and airline networks. In this vein, Veremyev et al. (2015) present a distance-based critical node problem (DCNP), which generalizes CNP by using a distance-based connectivity measure to evaluate the network that remains after a set of nodes is removed from the original one. Both CNP and DCNP are computationally challenging NP-hard problems (Veremyev et al., 2015).

Due to the theoretical and practical significance of DCNP, considerable effort has been devoted to solving it. Existing algorithms for it can be divided into exact and heuristic algorithms. To solve DCNP optimally, Veremyev et al., (2015) present two exact solution approaches. However, they are time-consuming on graphs with a large diameter and non-unit edge lengths. To overcome this difficulty, Hooshmand et al. (2020) present a new integer programming (IP) model with a simple structure and solve it by an efficient Bender decomposition algorithm. Recently, Salemi and Buchanan (2022) introduce two new IP models (i.e., thin and path-like), and compare them with the existing one (Veremyev et al., 2015). Alozie et al. (2021) present a new path-based model for DCNP by using separation heuristics and valid inequalities that exploit the structure of the problems to enhance the models. Aringhieri et al. (2019) present dynamic programming algorithms for DCNP over special graph classes such as paths, trees, and series-parallel graphs. To solve DCNP approximately, Alozie et al. (2022) present a centrality-based heuristic algorithm for DCNP. Recently, Zhou et al. (2023a) propose a fast tri-individual memetic search approach for DCNP, which maintains a small population of only three individuals during its search process.

Existing studies on CNDPs have primarily focused on k-vertex CNDPs and their distance-based variants (Arulselvan et al., 2009; Baggio et al., 2021; Zhou et al., 2024). However, there has been relatively little research devoted to solving \beta-connectivity CNDPs. To date, we are aware of only four studies that have specifically focused on them. Lalou et al. (2016) study a 3C-CNP that seeks to find a minimal set of nodes, whose removal is constrained by requiring the size of each connected component in the residual graph to satisfy a given bound. Furini et al. (2020) study a k-vertex cut problem that aims to determine the minimum weight subset of nodes whose removal produces a graph consisting of at least k connected components. Zhou et al. (2023b) solve an \alpha-separator problem that attempts to find a minimum node separator such that any connected components in the residual graph contains at most \lceil \alpha \cdot |V|\rceil nodes. A comprehensive survey on CNDPs can be found in Lalou et al. (2018).

Notably, these studies do not take into account the distance-based connectivity measure that depends on the lengths of pairwise distances between nodes within the remaining connected components. To the best of our knowledge, no study has focused on the distance-based \beta-connectivity CNDPs. However, in many real-world applications such as communication, transportation, and social networks, the separation of two nodes is considered effective only when their distance is sufficiently great.

To close this gap, we study a distance-based critical node problem with an unknown budget constraint named \beta-distance-based vertex disruptor problem (\beta-DVDP), which extends a \beta-vertex disruptor problem (\beta-VDP) by considering the distance-based connectivity measure. Both \beta-DVDP and \beta-VDP belong to the category of \beta-connectivity CNDPs. Table 1 presents a detailed comparison of different CNDP variants. We aim to make the following main contributions.

1) Methodologically, we introduce, for the first time, a distance-based critical node problem with an unknown budget constraint, denoted as \beta-DVDP, which is a distance-based variant of the \beta-VDP. This problem is motivated by practical applications in vulnerability assessments of various networks. We first derive three integer linear programming formations, i.e., recursive, triangular connectivity-based and reduced path-based ones. For solving larger instances, we propose an integrated strategic oscillation search (ISOS) approach inspired by hill climbing and strategic oscillation. ISOS iteratively improves an initial starting solution to obtain a high-quality one by examining a large neighborhood based on a tabu-enhanced destruction procedure and a random order-based improving construction procedure.

2) Empirically, we conduct extensive computational experiments with real-world and synthetic benchmark instances that reveal interesting insights and demonstrate the advantages and limitations of existing exact and heuristic algorithms. Among other outcomes, our experiments establish that the reduced path-based formulation outperforms its two peers. They also show that exact IP solvers can handle only very small sized instances and demonstrate that ISOS is significantly faster than them. Finally, we perform a vulnerability assessment using ISOS to understand the impact of the CVOID-19 pandemic on air transportation networks.

The rest of this paper is organized as follows. Section 2 describes \beta-VDP and its distance-based variant \beta-DVDP. Section 3 introduces three proposed IP formulations for \beta-DVDP. Section 4 presents ISOS. Computational results on both exact and heuristic algorithms are reported in Section 5. To show the practical value of our approach, a case study on the vulnerability assessment of an air transportation network is performed in Section 6. Finally, conclusions are made in Section 7.

2 Problem description

The well-known CNP (Arulselvan et al., 2009; Zhou et al., 2024; Tan et al., 2025) and its distance-based version DCNP (Veremyev et al., 2015; Salemi and Buchanan, 2022; Zhou et al., 2023a) are k-vertex CNDPs. They usefully model a category of practical applications where we have information about the number of nodes to target but do not know how to optimally achieve our goal. However, in some cases, we have information about a network structure we need to obtain, but we may lack knowledge about the specific nodes to target. In this situation, as previously noted, the relevant problems become \beta-connectivity CNDPs. The \beta-vertex disruptor problem and its distance-based variant (i.e., the \beta-DVDP), are two \beta-connectivity CNDPs and constitute important graph-theoretical optimization problems that underlie the network vulnerability assessment.

2.1 \beta-vertex disruptor problem

The challenge of network vulnerability assessment embodied in \beta-VDP (Dinh et al., 2010) is how to find a minimum set of nodes whose removal degrades pairwise connectivity to a desired degree. Given a strongly connected graph G=(V,E) with |V|=n nodes and |E|=m edges, and a fixed number , \beta-VDP aims to find a subset S\subseteq V of minimum size that can be removed from G such that the total pairwise connectivity in the remaining graph G[V\mathrm{\setminus }S] is no more than \beta \cdot (\genfrac{}{}{0ex}{}{n}{2}). Such a set S is called a \beta-vertex disruptor.

(1){min}_{S\subseteq V}\sum _{i=1}^n{x}_i\text{，}

(2)\mathrm{s}.\mathrm{t}.\sum _{i=1}^T(\genfrac{}{}{0ex}{}{|C_i|}{2})\le \beta \cdot (\genfrac{}{}{0ex}{}{n}{2})\text{，}

where x_i=1 if node i\in S, and otherwise 0, T is the total number of connected components C_i in the residual graph G[V\mathrm{\setminus }S].

2.2 \beta-distance-based VDP

\beta-DVDP further extends \beta-VDP by considering the distance-based connectivity measure. In many real-world applications—such as communication, transportation, and social networks—two nodes are considered effectively separated only when their distance exceeds a sufficiently large value. The objective is to minimize the number of deleted nodes so that the distance-based connectivity metric \sigma is bounded by a given threshold \beta, i.e.,

(3){min}_{S\subseteq V}\sum _{i=1}^n{w_i\cdot x}_i,

(4)

where x_i=1 if node i\in S, and 0 otherwise, w_i\in R_{+} is the weight of node i\in V, function \psi :Z_{+}\cup \{+\mathrm{\infty }\}\to R represents a distance-based connectivity metric, and c_{ij}\in {\mathbb{R}}_{+} is the cost associated with the connection between nodes i and j, which can represent the physical lengths of each pair of nodes (i,j). The distance d(i,j) is the length of the shortest path between iand j in G[V\mathrm{\setminus }S] if they are connected, and otherwise +\mathrm{\infty }. Finally, N presents the total distance-based connectivity of the original graph G, i.e.,

(5)

The distance-based connectivity measure \psi (\cdot ) is a function of the pairwise distances between nodes in the remaining graph (Veremyev et al. 2015), rather than just considering whether nodes are connected or not in CNDP. The definition of “distance-based connectivity” in a network varies depending on practical applications, resulting in different variants of CNDP. There are many distance-based connectivity measures available in the literature, such as minimizing the total number of pairs of nodes connected by a hop distance of at most \theta (a positive integer), minimizing the Harary index (which is equivalent to minimizing the graph’s efficiency), minimizing the sum of power functions of distances in a graph, maximizing the generalized Wiener index, and maximizing the shortest path between two given nodes in a graph.

In this work, we consider the general case of \beta-DVDP, i.e., minimizing the number of node pairs connected by a path of length at most θ. Correspondingly, the distance-based connectivity measure is defined as follows:

(6)\psi (d(i,j))=\{\begin{array}{l}1,d(i,j)\le \theta ,\\ 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\end{array}

where \theta represents a maximum distance limit. This definition implies that a pair of nodes can only be considered as truly connected when the hop distance between them does not exceed \theta. Additionally, we focus on the unweighted case, i.e., c_{ij}=1, \mathrm{\forall }i,j\in V and . Note that when \theta \ge n-1, \psi (d(i,j)) only considers whether there is a path between i and j, in which case \beta-DVDP reduces to the standard \beta-VDP (Dinh et al. 2010). Both \beta-VDP and \beta-DVDP belong to the category of \beta-connectivity CNDPs.

3 Integer programming formulations

This section presents three formulations for \beta-DVDP—recursive, triangular connectivity-based, and reduced path-based ones, which are inspired by the existing DCNP studies (Veremyev et al. 2015, Hooshmand et al. 2020, Salemi and Buchanan 2022). Note that \beta-DVDP belongs to the class of \beta-connectivity CNDPs; whereas DCNP falls into the category of k-vertex CNDPs.

3.1 Recursive (REC) model

For any node i\in V, we define a decision variable x_i as

(7)\begin{array}{c}{x}_i=\{\begin{array}{l}1,\text{if; node}i\text{is deleted}\\ 0,\text{otherwise}\end{array}\end{array}

Similarly, we define the connectivity variables y_{ij}^{\theta } as follows:

(8)\begin{array}{c}{y}_{ij}^{\theta }=\{\begin{array}{l}1,\text{if}i\text{and}j\text{are linked by a path with}d(i,j)\le \theta \\ 0,\text{otherwise}\end{array}\end{array}

Following Veremyev et al. (2015), the REC model is provided as follows:

\beta -{\mathbf{D}\mathbf{V}\mathbf{D}\mathbf{P}}_{\mathbf{R}\mathbf{E}\mathbf{C}}:

(9)min\sum _{i\in V}{x}_i,

(10)\text{s.t.}{y}_{ij}^1+x_i+x_j\ge 1,\mathrm{\forall }\left(i,j\right)\in E,

(11)y_{ij}^{\theta }=y_{ij}^1,\mathrm{\forall }\left(i,j\right)\in E,\theta \in \{2,\dots ,k\},

(12)y_{ij}^{\theta }+x_i\le 1,\mathrm{\forall }i,j\in V,\theta \in \{1,\dots ,k\},

(13)y_{ij}^{\theta }+x_j\le 1,\mathrm{\forall }i,j\in V,\theta \in \{1,\dots ,k\},

(14)y_{ij}^{\theta }\le \sum _{t:(i,t)\in E}{y}_{tj}^{\theta -1},\mathrm{\forall }\left(i,j\right)\notin E,\theta \in \{2,\dots ,k\},

(15)y_{ij}^{\theta }\ge y_{tj}^{\theta -1}-x_i,\mathrm{\forall }\left(i,t\right)\in E,\theta \in \{2,\dots ,k\},

(16)\sum _{i,j\in V}\left(y_{ij}^1+\sum _{\theta =2}^k\left(y_{ij}^{\theta }-y_{ij}^{\theta -1}\right)\right)\le \beta \cdot N,

(17)x_i\in \{0,1\},\mathrm{\forall }i\in V,

(18)y_{ij}^{\theta }\in \{0,1\},\mathrm{\forall }i,j\in V,\theta \in \{1,\dots ,k\}.

Constraints (10)–(11) ensure that y_{ij}^{\theta }=1 for adjacent node pairs (i,j), if neither i nor j is deleted and there exists a path of length at most \theta between them in the remaining graph, where \theta \in \{2,...,k\} and k denotes the length of the longest possible path in the residual graph. Constraints (12)–(13) enforce y_{ij}^{\theta } to be zero if either i or j is deleted. Constraints (14)-(15) ensure that there is a path of length at most \theta between non-adjacent nodes i and j if i is not deleted and there is a path of length at most \theta -1 between some non-deleted node t and j for in the neighborhood of node i. (16) limits the distance-based connectivity of the residual graph to a specified proportion \beta of the initial distance-based connectivity of the original graph. Constraints (17)–(18) are binary restrictions on the decision variables. The following improvements can be made to enhance \beta-{\text{DVDP}}_{\text{REC}}.

(1) Since y_{ij}^l=y_{ji}^l, it is sufficient to define only one of these variables, such as those with , and consequently reduce the number of variables and constraints. In this case, a slight modification of the constraints in the recursive model is necessary, which involves the replacement of each occurrence of y_{ij}^{\theta } with by y_{ij}^{\theta }. Additionally, it is necessary to replace (14) with Constraints (19), and (15) with (20)–(21), i.e.,

(19)

(20)

(21)

(2) Variable y_{ij}^l with d(i,j)>\theta in G before the deletion of the nodes can be omitted from the model (or fixed to zero), thereby reducing the number of variables.

(3) Let N_1 denote all nodes in V with degree 1 in G. These nodes are also known as the leaf nodes of G. Then for any leaf node i\in N_1, its associated variable x_i=0 if its stem j (i.e., (i,j)\in E) is not in N_1. Intuitively, it is preferable to delete the stem of i instead of i\in N_1.

(4) The integrality restrictions for variables y_{ij}^l\in \{0,1\} can be relaxed to be continuous in [0,1].

3.2 Triangular connectivity-based (TCB) model

We define a set of nonnegative continuous variables , which represent the length of the shortest path between i and j in G[V\mathrm{\setminus }S] if i and j are connected by a path with d(i,j)\le \theta in G[V\mathrm{\setminus }S]. Otherwise, u_{ij} is equal to M, which can be set to \theta +1 when distance is hop-based. We additionally define a set of binary variables {\mathit{\text{z}}}_{ij}(i,j\in V: ) that equals zero if u_{ij}\le \theta, and 1 otherwise.

Following Hooshmand et al. (2020), a TCB model is as follows:

\beta -{\mathbf{D}\mathbf{V}\mathbf{D}\mathbf{P}}_{\mathbf{T}\mathbf{C}\mathbf{B}}:

(22)min\sum _{i\in V}{x}_i,

(23)\text{s.t.}{x}_i=0,\mathrm{\forall }i\in N_1,

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)x_i\in \{0,1\},\mathrm{\forall }i\in V,

(32){\mathit{\text{z}}}_{ij}\in \{0,1\},\mathrm{\forall }i,j\in V,

(33)u_{ij}\ge 0,\mathrm{\forall }i,j\in V.

Constraint (23) ensure that the leaf nodes are not deleted. Constraint (24) ensure that for each edge (i,j)\in E, if neither i nor j is deleted, the length of the shortest path between i and j is at most equal to 1. Constraint (25) establish the upper bound limit M for the variable u_{ij}. Constraints (26)–(28) establish the triangle inequality for each triple (i,j,k) with . More precisely, Constraint (26) state that for three nodes i,j, and k in the graph, the shortest path between i and k is not longer than the sum of the shortest path between i and j and the shortest path between j and k. Similarly, Constraint (27) and (28) enforce the triangle inequality on the shortest path between j and k, and between i and j, respectively. Constraint (29) together with Constraint (30) ensure that {\mathit{\text{z}}}_{ij} takes a value of 1 whenever u_{ij}>\theta. Constraint (30) bounds the distance-based connectivity of the residual graph. Constraints (31)–(33) show that the decision variables are binary.

In the above model, the values of variables u_{ij} are determined by using triangle inequalities. As suggested by Hooshmand et al. (2020), for any pair of nodes i and j, if the length

of the shortest path between them is less than their distance, their shortest one must first pass through a node in the neighborhood of i before connecting i and j. Therefore, Constraints (26)–(28) can be replaced by Constraints (34)–(36).

(34)

(35)

(36)

3.3 Reduced path-based (RPB) model

Considering the characteristics of a distance-based connectivity metric, we find that shorter paths contribute more to the objective function than longer paths. Since only paths of length at most \theta are important, we restrict consideration to paths whose lengths do not exceed this threshold. The reduced path-based model is based on the minimal length bounded vertex connectors.

Definition 1 (Minimal Length-\theta i,j-Connector). A subset of vertices P\subseteq V that contains i and j is called a length-\theta i,j-connector in G=(V,E) if d_{G[P]}(i,j)\le \theta, where d_{G[P]}(i,j) is the length of the shortest path between i and j in G[P]=(P,(P\times P)\cap E) that is the subgraph induced by P in G=(V,E). If no proper subset of P is a length- \theta i,j-connector, then P is a minimal length- \theta i,j-connector.

Given a pair of nodes (i,j) with d_{G[P]}(i,j)\le \theta, let P_{\theta }(i,j) denote the set of all candidate minimal length-\theta i,j-connectors. Based on the concept of a minimal length-\theta i,j-connector, we propose a new IP model using a new set of connectivity variables y_{ij}:

(37)\begin{array}{c}{y}_{ij}=\{\begin{array}{l}1,\text{if}i\text{and}j\text{are linked by a path with}d(i,j)\le \theta ,\\ 0,\text{otherwise}.\end{array}\end{array}

When \theta is large, the above RPB model can result in too many constraints in (39) as their count grows exponentially with \theta. Adding them to the formulation may make the problem intractable. However, as only a small number of constraints are active in the formulation, they can be added as “lazy constraints” and introduced into the formulation only when they are violated by the current solution, thereby reducing the overall number of constraints. To decide if (39) is violated, we utilize the separation routines (Salemi and Buchanan, 2022), which solve a shortest path problem and run in polynomial time.

Using the former set of node deletion variables (7) along with the new connectivity variables (37), we h ave the RPB model for \beta-DVDP:

\beta -{\mathbf{D}\mathbf{V}\mathbf{D}\mathbf{P}}_{\mathbf{R}\mathbf{P}\mathbf{B}}:

(38)min\sum _{i\in V}{x}_i,

(39)

(40)

(41)x_i\in \{0,1\},\mathrm{\forall }i\in V,

(42)

Constraint (39) ensures that node pairs i and j are \theta-distance disconnected if at least one vertex for each minimal length-\theta i,j-connector connecting i and j is deleted. Constraint (40) bounds the distance-based connectivity in the residual graph to a specified proportion \beta of the initial distance-based connectivity of the original graph, while Constraints (41)–(42) ensure binary decision variables.

When \theta is small, Constraint (39) can be satisfied by enumerating all possible minimal length- \theta i,j-connectors that connect nodes i and j for each pair of nodes (i,j). In this case, the following family of valid inequalities can be incorporated to strengthen the model.

(43)

where \cap P_{\theta }(i,j) is the intersection of all minimal length-\theta i,j-connectors. Constraint (43) states that if all the minimum length-\theta connectors between nodes i and j pass through node r, then the removal of r from the graph disconnects nodes i and j.

When θ is large, the above RPB model can result in too many constraints in Constraint (39) as their count grows exponentially with θ. Adding them to the formulation may make the problem intractable. However, as only a small number of constraints are active in the formulation, they can be added as “lazy constraints” and introduced into the formulation only when they are violated by the current solution, thereby reducing the overall number of constraints. In order to decide if Constraint (39) is violated, we utilize the separation routines (Salemi and Buchanan, 2022), which solve a shortest path problem and run in polynomial time.

Finally, the RPB model can be enhanced by omitting the definition of variables y_{ij} for node pairs (i,j) with d(i,j)>\theta in the original graph from the model, which leads to a significant decrease in the number of variables. Furthermore, the integrality restrictions for variables y_{ij}\in \{0,1\} can also be relaxed to allow y_{ij} to take continuous values in [0,1].

4 Integrated strategic oscillation search

Exact solvers can achieve the optimal performance on small problem instances, but do not work well for medium and large-scale instances. In this section, we present a new heuristic approach called integrated strategic oscillation search (ISOS) for solving \beta-DVDP. We focus on its unweighted case, i.e., w_{ij}=1, \mathrm{\forall }(i,j)\in E. Without loss of generality, f(\cdot ) is assumed to be a nonnegative function, i.e., f:Z_{>0}\cup \{+\mathrm{\infty }\}\to R_{+}, and f(+\mathrm{\infty })=0. We start with the solution representation and evaluation, followed by the general framework of ISOS and detailed introduction to each of its key algorithmic component.

4.1 Solution representation and evaluation

Given a graph G=(V,E) and threshold \beta (), a solution can be represented by a subset of the removed nodes, and any subset S\subseteq V can be considered as a feasible solution if its removal from G leads to G[V\mathrm{\setminus }S] with distance-based connectivity metric bounded by \beta. To measure the total connectivity of G[V\mathrm{\setminus }S], we employ in this work a distance-based connectivity metric , where d_{G[V\mathrm{\setminus }S]}(i,j) denotes the distance between i and j in G[V\mathrm{\setminus }S]. Therefore, S\subset V is a feasible solution if and only if \phi (S)\le \beta \cdot N, where . Then the quality of S is evaluated by its size f(S)=|S|. Given two feasible solutions S and S^{\prime }, S is better than S^{\prime } if . Consequently, the task of solving \beta-DVDP is to find a feasible solution of minimum size.

4.2 General framework

ISOS combines the fundamental concepts of both hill climbing (HC) (Kirkpatrick et al., 1983) and strategic oscillation (SO) (Glover, 1989). HC is a fast and greedy local search that can converge prematurely and become trapped in local optima, whereas SO may be viewed as a diversification strategy that alternates between a destruction procedure and construction one.

● HC (Kirkpatrick et al., 1983) is one of the simplest local search algorithms. It only accepts candidate solutions with the better cost than the current one, which leads the search to quickly being trapped in a local optimum. A better solution (in a longer computation time) can be achieved via alternative techniques, which allow the limited acceptance of non-improving moves. Motivated by this, several effective HC variants have been proposed in the literature, such as late acceptance HC (Burke and Bykov, 2017) and step counting HC search (Bykov and Petrovic, 2016).

● SO (Glover, 1989) is closely linked to tabu search (Zuo et al., 2019; Zhang et al., 2023) as one of the earliest proposed strategies for diversifying the search. It integrates a construction procedure to “add” elements and a destruction once to “drop” elements with respect to a critical event. It has been successfully applied to solve many NP-hard problems, such as multidimensional knapsack problem (Glover and Kochenberger, 1996) and quadratic minimum spanning tree problem (Lozano et al., 2014).

Definition 2 (Non-improving and Non-worsening (NINW) Solution). Given a feasible starting solution S, a feasible solution S^{\prime } with an objective function value equaling that of S, i.e., f(S^{\prime })=f(S), is considered as an NINW solution.

ISOS combines the benefits of both HC and SO by incorporating destruction and construction procedures to augment HC. During each of its iteration, SO is employed to generate highly promising neighboring solutions. This is achieved by using a destruction procedure to add a node to the current solution and a construction procedure to remove one or more nodes from the solution. Its use for neighborhood exploration offers three key benefits compared to an explicitly defined neighborhood. First, it allows the algorithm to potentially generate high-quality solutions by using complex moves that operate on multiple nodes simultaneously. Second, it enables the algorithm to avoid exploring the entire neighborhood induced by these complex moves, which would be very large and time-consuming to evaluate. This can greatly speed up the search process. Third, it creates a more flexible search behavior as the number of removed nodes depends on the current search context and is not predefined. This flexibility can potentially lead to the discovery of more promising solutions.

Then, ISOS employs a simulated HC strategy to accept the solution produced by SO. Unlike HC that only accepts improving solutions, ISOS can accept both improving and NINW solutions. This acceptance criterion enables ISOS to more effectively escape from local optima by accepting some non-improving neighboring solutions (i.e., NINW solutions) as part of its strategy. Note that SO in general allow for worsening solutions as well as NINW solutions.

The general scheme of our ISOS consists of three main algorithmic components: solution initialization, tabu-enhanced destruction, and random order-based improving construction as realized in Algorithm 1. The method starts from an initial solution generated by a solution initialization procedure and then iteratively improves it to a high-quality solution by examining a neighborhood based on a tabu-enhanced destruction procedure and a random order-based improving construction procedure. If the resulting solution produced by SO is an improving or NINW solution, it is accepted by the acceptance criterion. Otherwise, the search process stops and restart with a new initial solution. ISOS repeats the above process until a given stopping condition (e.g., the time limit) is satisfied. Figure 1 shows the flowchart of ISOS to solve-DVDP.

4.3 Solution initialization

ISOS starts its search from an initial solution. As betweenness centrality (Brandes, 2001) is widely used to evaluate the importance of a node, we build an initial solution based on k-betweenness centrality (Jiang et al., 2009).

Definition 3 (k-Betweenness Centrality). Let d(s,t) denote the length of the shortest path between nodes s and t. We define {\sigma }_{s,t}^k as the number of paths between nodes s and t of length at most d(s,t)+k, and {\sigma }_{s,t}^k(v) as the number of these paths of length at most d(s,t)+k that pass-through node v. The k-betweenness centrality is given by

(44){\phi }_k(v)={\sum }_{s\ne v\ne t\in V}\frac{{\sigma }_{s,t}^k\left(v\right)}{{\sigma }_{s,t}^k}.

Initially, all nodes in G are sorted according to their k-betweenness centrality values in descending order and stored in a node sequence L, which can be finished in O(|E||V|+|V|\mathrm{l}\mathrm{o}\mathrm{g}|V|). Starting from an empty set, at each step, a node of L with the maximum betweenness centrality value is added into S with a given removal probability \tau (). Once a node is successfully added into S, it is removed from L. The process is repeated until \phi (S)\le \beta \cdot |N|. Correspondingly, a feasible solution consisting of all removed nodes is obtained. Note that |N| represents the total number of valid node pairs in the original G, which can be effectively computed by building k-depth breadth first search (BFS) trees for each node. It requires O(b^k) time to build all trees, where b denotes the branching factor of a tree.

4.4 Tabu-enhanced destruction procedure

The destruction procedure aims to dismantle the structure of the current solution. The degree of destruction is a key consideration when implementing this procedure. If only a small portion of a solution is dismantled, ISOS may have difficulty to explore search space. If a very large portion of a solution is dismantled, then ISOS degrades into repeated re-optimization. It can be time consuming or yield a poor solution depending on how a partial solution is repaired.

This procedure typically incorporates an element of stochasticity to ensure that different parts of a solution are dismantled in each invocation. To this end, we propose a tabu-enhanced destruction procedure that incorporates the basic idea of the tabu search mechanism (Glover, 1989) into the general destruction procedure. Specifically, to dismantle a solution, a best allowed node is removed from G[V\mathrm{\setminus }S] and added to S (i.e., S\leftarrow S\cup \{v\}). Note that we aim to minimize the number of vertices in S, adding a vertex to S leads to a deteriorating solution. The removed node v should not only be a non-tabu node, but also have the maximum k-betweenness centrality value. Once v is removed, it cannot be reassigned to the remaining graph again in the next \delta +Random(\mu ) iterations, where \delta and \mu >1 are two constants (i.e., a tabu constant and a random constant), and Random(\mu ) is a function to return a random number in \{0,1,\dots ,\mu -1\}. Note that a tabu node whose removal leads to a solution better than the best solution found so far (termed an aspiration criterion) is always accepted.

4.5 Random order-based improving construction procedure

Once a dismantled solution S^{\prime } is obtained, a construction procedure is applied to rebuild it. There is a lot of freedom in choosing a construction procedure, which may be optimal or heuristic.

We use a random order-based improving construction procedure to repair a dismantled solution S^{\prime }. The construction at each step maintains or improves the objective function value of the starting solution while generating a more complete solution. At the beginning, all nodes of S^{\prime } are randomly sorted, then each node u is checked to determine whether it can be put back into the remaining graph G[V\mathrm{\setminus }{S}^{\prime }]. The node u can be removed from S^{\prime } and added to S if \phi (S^{\prime }\mathrm{\setminus }\{u\})\le \beta \cdot N is satisfied. Using this procedure, ISOS can obtain a better solution S than the dismantled solution S^{\prime } (i.e., f(S^{\prime })-f(S)=\mathrm{\Delta }f\ge 1), as both the tabu-enhanced destruction procedure and the random order-based improving construction procedure operate on feasible solutions. In particular, when \mathrm{\Delta }f>1, ISOS can obtain a feasible solution that improves upon the starting solution; otherwise, ISOS finds a feasible solution with the same objective value as the starting one.

4.6 Connections with related studies

Table 2 presents a comparison of ISOS with related algorithms. Although the proposed ISOS method incorporates fundamental concepts from both hill climbing (HC) and strategic oscillation (SO), it distinguishes itself from these algorithms in the following three aspects.

1) In contrast to HC which explicitly defines the neighborhood, the neighborhood of ISOS is defined implicitly by a destruction procedure and a construction one. The former aims to dismantle a part of the current solution while the latter rebuilds the dismantled one.

2) Employing the idea of a construction/destruction type of SO, where a construction procedure is used to “add” elements and a destruction procedure is applied to “drop” elements with respect to the critical event, our SO strategy comes in two forms in problems where a feasibility boundary exists: a 2-sided strategy that crosses back and forth between feasible and infeasible search space, and a 1-sided strategy that remains on one side of a boundary between these regions. Our ISOS adopts a 1-sided strategy that employs destruction and construction procedures to perform neighborhood search in feasible solution space rather than crossing between feasible one and infeasible one.

3) ISOS adopts a large neighborhood strategy, which shares an idea in common with large neighborhood search (LNS) (Shaw, 1998) and is also emphasized as a feature of SO. The main idea behind LNS heuristic is that the large neighborhood allows the heuristic to navigate in solution space easily, even if an instance is tightly constrained, in contrast to a small neighborhood search heuristic like HC. ISOS distinguishes itself from LNS with the following features. a) our tabu-enhanced destruction procedure incorporates the basic idea of a tabu search mechanism into the destruction procedure; b) our random order-based improving construction procedure only adds one node at each time instead of multiple nodes; c) unlike LNS, ISOS creates a more flexible search behavior by operating on a variable number of nodes; d) ISOS always accepts a new neighboring solution because it is either an improving or NINW solution.

4.7 Complexity analysis

To analyze the computational complexity of ISOS, we consider each algorithmic component of Algorithm 1. ISOS starts its search from an initial solution S, which is generated by a solution initialization procedure in O(|E||V|+|V|log|V|+|S|b^k). At each iteration, ISOS sequentially executes a tabu-enhanced destruction procedure and a random order-based improving construction procedure to obtain a neighboring solution. The former requires O(|E||V|+|V|log|V|+b^k) time, while the latter takes O(|S|b^k) time. Therefore, the total time complexity of ISOS is O(|E||V|+|V|log|V|+|S|b^k) at each iteration.

5 Computational experiments

5.1 Benchmark instances

Our computational experiments are conducted on a set of benchmark instances, which have been widely used to evaluate the algorithms for DCNPs (Veremyev et al., 2015; Alozie et al., 2021; Zhou et al., 2023a). They can be divided into two categories:

● Real-world benchmark consists of two groups of 22 instances: R1 and R2. R1 is composed of 11 real-world networks selected from the Pajek and UCINET data sets. R2 includes 11 large real-world instances, which are selected from the original CNP instances.

● Synthetic benchmark is also composed of two groups of 33 instances: S1 and S2. The former includes 21 instances generated by using NetworkX random graph generators, while the latter consists of two original CNP instances and 10 new instances. Note that these new instances have the same sizes as the original benchmark instances in Ventresca (2012).

Following Alozie et al. (2022), we consider each instance of both real-world and synthetic benchmarks under two different budgetary constraints, i.e., \beta =0.25 and \beta =0.50.

5.2 Experimental settings

Our algorithms are implemented in C++ for CPLEX version 12.8.0. and compiled with gcc 8.1.0 and the flag “-O0.” All experiments are conducted on a standard PC equipped with an AMD Ryzen 7 5900U CPU, clocked at 1.9 GHz, and 16 GB of RAM operating under the Windows 10 system. In our experiments, we set the hop distance as three, i.e., \theta =3, which is an appropriate hop distance for most of benchmark instances, as suggested in Veremyev et al. (2015) and Salemi and Buchanan (2022). Small values of \theta ensure that paths are counted only when they are short, which is what distinguishes \beta-DVDP from \beta-VDP.The time limit for each run is set to one hour.

Table 3 provides the settings for several parameters required to obtain the results of ISOS. Since the performance of ISOS is sensitive to its parameter values, we employ the well-known automatic parameter configuration tool irace (Lopez-Ibanezet al., 2016) to determine suitable parameter values for ISOS. Since these benchmark instances have different scales and \beta values, our automatic configuration experiments are conducted on eight instances randomly selected from both real-world and synthetic benchmarks. The tuning process is given a budget of 100 runs for each instance with a time budget of one hour per run. The range and corresponding final value of each parameter are provided in the third and fourth column of Table 3, respectively. To make a fair comparison, we also adjust the parameter values of simulated annealing by the irace tool based on the same benchmark instances.

5.3 Comparison of three proposed formulations

In this section, we conduct a systematic comparison of the three proposed formulations, \beta-DVDP_REC, \beta-DVDP_TCB, and \beta-DVDP_RPB, using both real-world and synthetic benchmark instances. To solve them on each instance, we run CPLEX 12.8.0. until it either reaches one hour or runs out of memory. Detailed comparative results for the three proposed formulations under \beta =0.25 and \beta = 0.50 are presented in Tables 4 and 5, respectively. From these results, it is clear that \beta-DVDP_RPB outperforms \beta-DVDP_REC and \beta-DVDP_TCB. Table 6 provides a summary of the performance comparison results based on these two tables, grouped according to the four classes of instances. The first column shows the class with the number of instances (in brackets). The first block with six columns (2-7) shows the results given \beta =0.50. In particular, columns 2-3, 4-5 and 6-7 show the number of solved instances and the computation time (in seconds) for the three compared formulations, respectively. The second block with four columns (8-13) presents the results under \beta =0.25.

As shown in Table 6, \beta-DVDP_RPB formulation demonstrates overall superior performance compared with \beta-DVDP_REC and \beta-DVDP_TCB ones across the four sets of benchmark instances. On the R1 set, with \beta-DVDP_RPB, we are able to solve all instances to optimality within one hour. In contrast, all three formulations perform poorly on the R2 set, as we cannot solve any instances to optimality. For the S1 benchmark, \beta-DVDP_RPB shows slightly better performance than its two peers (3 versus 1 and 0 when \beta =0.25, and 3 versus 0 and 0 when \beta =0.50). A similar pattern is observed for the S2 benchmark. These findings confirm the overall superiority of \beta-DVDP_RPB over the other two. This superiority is reasonable, as it can be strengthened by omitting the definition of variables y_{ij} for node pairs (i,j) with d(i,j)>\theta in the original graph from the model, leading to a substantial reduction in the number of variables and thus improving computational efficiency.

5.4 Comparison between ISOS and baseline algorithms

In this section, we experimentally compare the performance of ISOS with two baseline algorithms (i.e., the RPB formulation solved by CPLEX and the simulated annealing (SA) algorithm) on both real-world and synthetic benchmarks. To make a fair comparison, we run them on the same computing platform with the same experimental conditions.

5.4.1 Results on real-world benchmarks

Tables 7 and 8 present the performance comparisons between ISOS and baseline algorithms on real-world benchmarks with \beta =0.25 and \beta =0.50, respectively. In them, column 1 presents the name of the instance set (Instance). Column 2 reports the best lower bound (LB) and upper bound (UB) found by the CPLEX solver at termination and the corresponding gap between LB and UP is shown in in Column 3. Column 4 presents the computation time (t) in seconds. Columns 5-7 and 8-10 list the results of SA and ISOS, respectively, including the best solution value (\hat{f}) found during 10 runs, the average solution value (\overline{f}) and the average computation time (\overline{t}) across the 10 runs. At their bottom, we give the number of instances on which ISOS reaches a better (#Wins), equal (#Ties), and worse (#Losses) results in terms of each indicator than the corresponding methods. At the last row, we also provide the p-values of the Wilcoxon signed ranks test.

From Table 7, we observe that ISOS can quickly find the optimal solutions for all 11 real-world instances of R1. In particular, it is significantly better than SA in terms of \overline{f} at a significance level of 0.05. The 11 real-world instances of R2 cannot be optimally solved by CPLEX in one hour. However, ISOS demonstrates excellent performance by improving the lower and upper bounds of all 11 instances. At a significance level of 0.05, ISOS significantly outperforms SA in terms of both \hat{f} and \overline{f}.

Next, Table 8 shows the comparative results for ISOS and SA under β = 0.50, where both ISOS and SA achieve the best \hat{f} values at least once on all 11 real-world instances of R1. However, ISOS shows more consistent performance than SA, by achieving the best \hat{f} values with a 100% success rate. ISOS also performs excellent performance on the 11 real-world instances of R2, which cannot be solved by CPLEX, and it improves the lower and upper bounds of many instances. At a significance level of 0.05, ISOS is significantly better than SA in terms of both \hat{f} and \overline{f}.

5.4.2 Results on synthetic benchmarks

Tables 9 and 10 present the performance comparisons between ISOS and the baseline algorithms on synthetic benchmarks under \beta =0.25 and \beta =0.50, respectively. From Table 9, we see that ISOS improves the lower and upper bounds of many instances. At a significance level of 0.05, ISOS significantly better than SA in terms of both \hat{f} and \overline{f} on both benchmarks S1 and S2.

Table 10 shows similar performance comparisons for \beta =0.50. In particular, ISOS improves the lower and upper bounds of many instances and significantly outperforms SA in terms of both \hat{f} and \overline{f} on benchmark S1 at a significance level of 0.05. For benchmark S2, however, there is no significant performance difference between the methods.

Finally, Table 11 summarizes the performance comparisons between ISOS and SA on different benchmarks and under different \beta values. From it, we observe that ISOS achieves better or equal performance than SA on all instances except for the ER100 instance with \beta =0.50, where ISOS and SA obtain the same \hat{f} value, and ISOS is slightly worse than SA in terms of \overline{f}. At a significance level of 0.05, ISOS is significantly better than SA in terms of both \hat{f} and \overline{f} on benchmarks R2 (with both \beta =0.25 and \beta =0.50), S1 (with both \beta =0.25 and \beta =0.50), and S2 (with \beta =0.25).

5.5 Experimental analysis

To gain a deeper understanding of ISOS, we perform two groups of additional experiments: 1) to compare the running performance of ISOS and SA; and 2) to evaluate the benefit of the tabu mechanism used in the destructive procedure. Our following experiments are conducted on eight representative instances, which are selected from both real-world and synthetic benchmarks.

5.5.1 Running profiles

To obtain a deeper insight into the performance of ISOS and SA, we experimentally compare their running performance curves, by defining a function t\mapsto f(t), where f(t)presents the best solution value found at the computational time t. Figure 2 presents the detailed performance curves between ISOS and SA on the eight instances, where the x-axis denotes the computation time t in seconds and the y-axis indicates the best solution value f(t) found at t.

From Fig. 2, we observe that ISOS consistently outperforms SA by finding smaller solution values on all eight instances. For the powergrid instance, we can see that ISOS quickly finds a better solution than SA, and reaches a steady value of 420 at 200 s. For the WS500 instance, we observe that ISOS shows a worse performance than SA at the beginning, but eventually outperforms SA at the end of the search.

5.5.2 Benefit of Tabu-enhanced destructive procedure

To show the benefit of the tabu mechanism integrated into the destructive procedure, we experimentally compare ISOS with a variant called ISOS* which is obtained by disabling the tabu mechanism in the tabu-enhanced destructive procedure. Detailed comparative results between ISOS and ISOS* are summarized in Table 12.

In Table 12, columns 1–2 present the names of the instances (Instance) and its \beta value, columns 3–5 give the results of ISOS*, including the best solution value (\hat{f}) found during the 10 runs, the average solution value (\overline{f}) and the average computation time (\overline{t}) across the 10 runs. Correspondingly, columns 6-8 list the results of ISOS. Note that the better values of each performance indicator are in bold. At the bottom of the table, we count the number of instances on which ISOS researches a better (#Wins), equal (#Ties), and worse (#Losses) results in terms of both \hat{f} and \overline{f} than ISOS*. From Table 12, we observe that ISOS demonstrates better performance than ISOS*. In particular, it outperforms ISOS* in terms of \hat{f} on four instances and achieves the same results on the remaining four instances. The same observation can also be made for performance indicator \overline{f}. These results confirm the effectiveness of the tabu-enhanced destructive procedure.

6 Case study

In this section, we present a practical application of network vulnerability assessment by investigating the influence of COVID-19 on the air transportation network. As an ongoing global pandemic of coronavirus disease, COVID-19 has had a significant impact on many industries, of which air travel has been one of the hardest hit, with many flight cancellations and airport closures as a consequence. By analyzing the characteristics of air transportation networks, we aim to gain insights into the effects of COVID-19 on the air transportation industry. Using official flight data, we build four similar US air transportation networks—USAir2018, USAir2019, USAir2020 and USAir2021—where nodes represent airports and edges indicate existing flight connections. All node weights are set to 1, and we take \theta =3. The resulting network instances are available on the above-mentioned GitHub website.

In our experiment, we consider \beta =0.25. For each network, we execute ISOS under the time limit \hat{t}=30 minutes, and then record the best solution. All elements of the best solution can be considered as critical airports, and the remaining nodes in the air transportation network are non-critical airports. To investigate the impact of COVID-19 on the air transportation network, we compare all identified critical airports in US air transportation networks from 2018 to 2021. Table 13 summarize the statistical information of the critical airports identified. The first column is the airport code, columns 2-3 present the results of USAir2018, including the number of passengers served and a flag (“\surd”) to indicate whether an airport is critical. Correspondingly, columns 4-5, 6-7, and 8-9 show the results of USAir2019, USAir2020 and USAir2021, respectively. The last column reports the frequency with which each airport is identified as a critical node. In particular, a frequency of 4 indicates that the airport is consistently identified as being critical from 2018 to 2021, and can therefore be regarded as a major hub. At the end of Table 12, we show the total number of critical airports together with their construction numbers (four times | three times | two times | one time) for each network. The acronyms and abbreviations of the critical airports are provided in Appendix A, as shown in Table A1. From Table 13, three interesting observations can be made.

First, notable changes in passenger counts are observed between 2018 and 2021. Specifically, passenger volumes experienced only a minor shift from 2018 to 2019, followed by a sharp decline from 2019 to 2020, and then a rapid recovery from 2020 to 2021. This pattern is consistent across all airports. This changes likely reflect the onset of the COVID-19 pandemic at the end of 2019, which resulted in widespread flight cancellations and airport closures.

Second, notable changes are also observed in the number of critical airports. More specifically, 89, 106, 86, and 108 critical airports are identified in 2018, 2019, 2020, and 2021, respectively. The increase of 17 critical airports from 2018 to 2019 (from 89 to 106) suggests a relatively stable and expanding air transportation network prior to the COVID-19 pandemic. In contrast, the number of critical airports decreases by 20 in 2020 compared with 2019 (from 106 to 86), reflecting the severe disruption caused by widespread flight cancellations and temporary airport closures. The subsequent rebound in 2021, with an increase of 22 critical airports (from 86 to 108), indicates a substantial recovery of network connectivity as air travel gradually resumed, leading to the re-emergence of several airports as critical nodes.

Finally, 69 airports can be regarded as major hubs in the US air transportation network, as they are consistently identified as critical airports from 2018 to 2021. Notably, while the set of critical airports exhibits some variation over time, the primary impact of COVID-19 is reflected in changes in passenger volumes rather than in the structural roles of these major hubs.

The US air transportation networks are illustrated in Fig. 3, where critical and non-critical airports are represented by red and blue dots, respectively. As expected, the \beta-DVDP solutions partition the network into multiple connected components. Comparing the results across years, only minor changes are observed between USAir2018 and USAir2019, indicating a relatively stable network structure during this period. In contrast, a substantial restricting occurs between USAir2019 and USAir2020, coinciding with the onset of the COVID-19 pandemic. During this time, many flights were cancelled and numerous airports experienced temporary closures, leading to a significant reduction in connectivity. For example, the number of critical airports dropped from 106 in 2019 to 86 in 2020, and several previously connected components became fragmented. These patterns demonstrate that while the overall set of major hubs remains largely stable, the pandemic caused pronounced disruptions in network connectivity and passenger flows, highlighting the vulnerability of the air transportation system to large-scale shocks. decision-makers should pay close attention to the identified critical airports and prioritize their protection to maintain the basic connectivity of the network.

Based on the case study results, the proposed method offers practical guidance for decision-makers. By identifying critical airports and vulnerable network structures, it helps emergency management agencies prioritize response efforts and develop targeted contingency plans for high-impact disruptions. The quantitative assessment of node importance and performance degradation also supports more efficient resource allocation, enabling limited repair, security, and backup resources to be deployed where they can generate the greatest resilience improvement. Furthermore, the findings provide valuable insights for long-term infrastructure protection and investment planning, assisting decision-makers in strengthening key facilities to preserve overall network connectivity and robustness.

7 Conclusions and future work

In this work, we study a distance-based critical node detection problem with an unknown budget, called the \beta-distance-based vertex disruptor problem (\beta-DVDP), which is an important graph-theoretical optimization problem in network vulnerability assessment. We first derive three linear integer programming models, i.e., recursive, triangular connectivity-based and reduced path-based ones. We then solve small \beta-DVDP instances with the exact solver. To handle larger instances, we propose a simple and effective integrated strategic oscillation search (ISOS) approach that extends the hill-climbing strategy with strategic oscillation. ISOS iteratively improves an initial solution to yield a high-quality solution by examining a large neighborhood based on a tabu-enhanced destruction and a random order-based improving construction procedure.

Extensive experiments on both real-world and synthetic instances show the reduced path-based formulation outperforms the recursive and triangular connectivity-based models. ISOS is significantly faster than CPLEX. To verify and disclose the nature of the significant impact of the COVID-19 pandemic on the air traffic industry, we employ ISOS to perform the network vulnerability assessment on the US air transportation networks from 2018 to 2021. The proposed models and algorithms provide decision makers with valuable tools for evaluating network vulnerability and robustness, helping them prevent or mitigate adversarial attacks, natural disasters, disease outbreaks, and random failures. Our future research intends to address the following two issues. First, we plan to integrate ISOS with genetic algorithm and ant colony optimization. Second, we aim to explore its application to the vulnerability assessments of other real-world networks, such as supply chain and financial networks.

8 Appendix A

References

[1] Alozie G U, Arulselvan A, Akartunalı K, Pasiliao E L Jr, (2021). Efficient methods for the distance-based critical node detection problem in complex networks. Computers & Operations Research, 131: 105254

[2] Alozie G U, Arulselvan A, Akartunalı K, Pasiliao E L Jr, (2022). A heuristic approach for the distance-based critical node detection problem in complex networks. Journal of the Operational Research Society, 73( 6): 1347–1361

[3] Aringhieri R, Grosso A, Hosteins P, Scatamacchia R, (2019). Polynomial and pseudo-polynomial time algorithms for different classes of the distance critical node problem. Discrete Applied Mathematics, 253: 103–121

[4] Arulselvan A, Commander C W, Elefteriadou L, Pardalos P M, (2009). Detecting critical nodes in sparse graphs. Computers & Operations Research, 36( 7): 2193–2200

[5] Arulselvan A, Commander C W, Shylo O, Pardalos P M (2011). Cardinality-constrained critical node detection problem, 79–91 (New York, NY: Springer New York)

[6] Baggio A, Carvalho M, Lodi A, Tramontani A, (2021). Multilevel approaches for the critical node problem. Operations Research, 69( 2): 486–508

[7] Brandes U, (2001). A faster algorithm for betweenness centrality. Journal of Mathematical Sociology, 25( 2): 163–177

[8] Burke E K, Bykov Y, (2017). The late acceptance hill-climbing heuristic. European Journal of Operational Research, 258( 1): 70–78

[9] Bykov Y, Petrovic S, (2016). A step counting hill climbing algorithm applied to university examination timetabling. Journal of Scheduling, 19( 4): 479–492

[10] Dinh T N, Thai M T, (2015). Network under joint node and link attacks: Vulnerability assessment methods and analysis. IEEE/ACM Transactions on Networking, 23( 3): 1001–1011

[11] Dinh T N, Xuan Y, Thai M T, Park E K, Znati T (2010). On approximation of new optimization methods for assessing network vulnerability. Proceedings of the 29th Conference on Information Communications 2678–2686

[12] Furini F, Ljubic I, Malaguti E, Paronuzzi P, (2020). On integer and bilevel formulations for the k-vertex cut problem. Mathematical Programming Computation, 12( 2): 133–164

[13] Glover F, (1989). Tabu Search-Part I. ORSA Journal on Computing, 1( 3): 190–206

[14] Glover F, Kochenberger G A (1996). Critical event tabu search for multidimensional knapsack problems, 407–427 (Boston, MA: Springer US)

[15] Hooshmand F, Mirarabrazi F, MirHassani S A, (2020). Efficient benders decomposition for distance-based critical node detection problem. Omega, 93: 102037

[16] Jiang K, Ediger D, Bader D A (2009). Generalizing k-betweenness centrality using short paths and a parallel multithreaded implementation. 2009 International Conference on Parallel Processing, 542–549 (IEEE)

[17] Kirkpatrick S, Gelatt C D Jr, Vecchi M P, (1983). Optimization by simulated annealing. Science, 220( 4598): 671–680

[18] Lalou M, Tahraoui M, Kheddouci H, (2016). Component-cardinality-constrained critical node problem in graphs. Discrete Applied Mathematics, 210: 150–163

[19] Lalou M, Tahraoui M, Kheddouci H, (2018). The critical node detection problem in networks: A survey. Computer Science Review, 28: 92–117

[20] Liu C, Ge S, Zhang Y, (2023a). Identifying the cardinality-constrained critical nodes with a hybrid evolutionary algorithm. Information Sciences, 642: 119140

[21] Liu J, Xi Y, Wang J, (2023b). Resilience strategies for sustainable supply chains under budget constraints in the post COVID-19 era. Frontiers of Engineering Management, 10( 1): 143–157

[22] López-Ibáñez M, Dubois-Lacoste J, Pérez Cáceres L, Birattari M, Stützle T, (2016). The irace package: Iterated racing for automatic algorithm configuration. Operations Research Perspectives, 3: 43–58

[23] Lozano M, Glover F, García-Marítnez C, Rodríguez F J, Marít R, (2014). Tabu search with strategic oscillation for the quadratic minimum spanning tree. IIE Transactions, 46( 4): 414–428

[24] Salemi H, Buchanan A, (2022). Solving the distance-based critical node problem. INFORMS Journal on Computing, 34( 3): 1309–1326

[25] Shaw P Using constraint programming and local search methods to solve vehicle routing problems. Principles and Practice of Constraint Programming: 4th International Conference, CP98 Pisa, Italy, October 26–30, 1998 Proceedings 1998, 417–431 (Springer)

[26] Tan X, Zhou Y, Zhou M, Fu Z H, (2025). Learning to detect critical nodes in sparse graphs via feature importance awareness. IEEE Transactions on Automation Science and Engineering, 22: 3772–3782

[27] Ventresca M, (2012). Global search algorithms using a combinatorial unranking-based problem representation for the critical node detection problem. Computers & Operations Research, 39( 11): 2763–2775

[28] Veremyev A, Prokopyev O A, Pasiliao E L, (2015). Critical nodes for distance-based connectivity and related problems in graphs. Networks, 66( 3): 170–195

[29] Wang F, Li K, Chen X, Zhang W, (2025). Assessing supply chain risks for chip industry with LDA and multi-layer bayesian network method. Frontiers of Engineering Management, 12( 4): 1037–1057

[30] Zhang J, Liu H, Zu P, Zhao M, Wang C, Albeshri A, Abusorrah A, Zhou M, (2023). Using tabu search to avoid concave obstacles for source location. IEEE Transactions on Intelligent Transportation Systems, 24( 11): 11720–11732

[31] Zhou Y, Hao J K, Glover F, (2019). Memetic search for identifying critical nodes in sparse graphs. IEEE Transactions on Cybernetics, 49( 10): 3699–3712

[32] Zhou Y, Kundu T, Qin W, Goh M, Sheu J B, (2021). Vulnerability of the worldwide air transportation network to global catastrophes such as COVID-19. Transportation Research Part E, Logistics and Transportation Review, 154: 102469

[33] Zhou Y, Li J, Hao J K, Glover F, (2024). Detecting critical nodes in sparse graphs via reduce-solve-combine memetic search. INFORMS Journal on Computing, 36( 1): 39–60

[34] Zhou Y, Wang G, Hao J K, Geng N, Jiang Z, (2023a). A fast tri-individual memetic search approach for the distance-based critical node problem. European Journal of Operational Research, 308( 2): 540–554

[35] Zhou Y, Zhang X, Geng N, Jiang Z, Wang S, Zhou M, (2023b). Frequent itemset-driven search for finding minimal node separators and its application to air transportation network analysis. IEEE Transactions on Intelligent Transportation Systems, 24( 8): 8348–8360

[36] Zuo X, Li B, Huang X, Zhou M, Cheng C, Zhao X, Liu Z, (2019). Optimizing hospital emergency department layout via multi-objective tabu search. IEEE Transactions on Automation Science and Engineering, 16( 3): 1137–1147

RIGHTS & PERMISSIONS

The Author(s)