Robust Design of Interdependent Networks Considering Intra-Network Support Flow

Cyber-physical systems comprise multiple interdependent networks supporting each other. This paper considers a combination of power grid and computer network as a typical instance of the interdependent networks. A node in an interdependent network fails if it cannot receive the required volume of support from the nodes in the other interdependent network. Meanwhile, each node supporting the other network is required to receive a sufficient volume of intra-network flow named “support flow”, which corresponds to the flow of electric power in the power grid and the flow of control messages in the computer network. Each sink node of the support flow must terminate a sufficient volume of support flow generated at source nodes in the same network in order to achieve the support of the other network. This means that the support flow is one of the essential elements to make the interdependent networks operable. The conditions related to the support flow increases the possibility of a small initial node failure invoking a cascading failure with a wide range. This paper proposes a robust design method for the interdependent networks to maximize the robustness against the cascading failures while assuring satisfaction of the constraints on the support flow. The method can derive the interdependency to minimize the impact of the cascading failures induced from a given set of possible initial node failures when the configuration of the respective interdependent networks is specified beforehand. The robust design problem is formulated using a mixed integer linear programming model and two kinds of metaheuristic methods are proposed to solve the problem efficiently even if a great number of initial node failures are specified. Several characteristics of the risk of the cascading failures are revealed from the results of the robust design. Simulation experiments demonstrate the effectiveness of the proposed metaheuristic methods.


Robust Design of Interdependent Networks
Considering Intra-Network Support Flow Nagao Ogino , Member, IEEE Abstract-Cyber-physical systems comprise multiple interdependent networks supporting each other.This paper considers a combination of power grid and computer network as a typical instance of the interdependent networks.A node in an interdependent network fails if it cannot receive the required volume of support from the nodes in the other interdependent network.Meanwhile, each node supporting the other network is required to receive a sufficient volume of intra-network flow named "support flow", which corresponds to the flow of electric power in the power grid and the flow of control messages in the computer network.Each sink node of the support flow must terminate a sufficient volume of support flow generated at source nodes in the same network in order to achieve the support of the other network.This means that the support flow is one of the essential elements to make the interdependent networks operable.The conditions related to the support flow increases the possibility of a small initial node failure invoking a cascading failure with a wide range.This paper proposes a robust design method for the interdependent networks to maximize the robustness against the cascading failures while assuring satisfaction of the constraints on the support flow.The method can derive the interdependency to minimize the impact of the cascading failures induced from a given set of possible initial node failures when the configuration of the respective interdependent networks is specified beforehand.The robust design problem is formulated using a mixed integer linear programming model and two kinds of metaheuristic methods are proposed to solve the problem efficiently even if a great number of initial node failures are specified.Several characteristics of the risk of the cascading failures are revealed from the results of the robust design.Simulation experiments demonstrate the effectiveness of the proposed metaheuristic methods.
Index Terms-Interdependent networks, robust design, optimization model, metaheuristic method, support flow.

I. INTRODUCTION
S MART cities can be achieved on a basis of cyber-physical systems where a computer network is integrated with various physical networks such as a power grid and a transportation network [1], [2], [3].Multiple integrated networks composing the cyber-physical system rely on and support each other [4].For example, each node in the computer network is supplied with the electric power from the power grid.Meanwhile, each node in the power grid is operated efficiently using the computer network.However, such interdependency between the multiple networks causes a risk of cascading failure with a wide range since an initial node failure within one network recursively affects the other networks [5], [6], [7].When a node fails initially in an interdependent network, nodes supported by the failed node also incur damage in the other interdependent network.This may cause additional node failures in the interdependent network where the initial node failure has occurred.Hence, the cascading failures appearing in the interdependent networks have been extensively studied from a variety of perspectives [8], [9], [10].The dynamics of the cascading failures and the vulnerability in the interdependent networks have been revealed theoretically and experimentally [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25].Several methods for designing the interdependent networks with the robustness against the cascading failures have been proposed [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37].Recovery schemes from the cascading failures have also been proposed for achieving the resilient interdependent networks [38], [39], [40], [41].
A combination of power grid and computer network has been intensively studied as a typical instance of the interdependent networks.Each distribution node in the power grid must be supplied with sufficient electric power from the generation nodes in order to provide the computer network with the required electric power [17], [22].Likewise, each edge node operating the power grid needs to receive a sufficient amount of control messages from the server nodes via the computer network.This paper names the intra-network flow of the electric power and the control messages "support flow".The support flow is generated at source nodes of the support flow and terminated at sink nodes of the support flow in the same network.Each sink node of the support flow must be able to terminate a sufficient volume of support flow while achieving the required support for the other network.Nevertheless, no existing design method for the interdependent networks has recognized the above condition on the support flow as one of the requirements for the operable interdependent networks [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37].Hence, this paper proposes a robust design method for the interdependent networks that guarantees satisfaction of the constraints imposed on the support flow.The proposed method can derive the interdependency between the networks to minimize the significant risk of the cascading failures when the configuration of the respective networks and a set of possible initial node failures are specified beforehand.
The main contributions of this paper are as follows.
1) The robust design problem for the interdependent networks accompanied with the constraints on the support flow is formulated using a MILP (Mixed Integer Linear Programming) model, where β-CVaR (β-Conditional Value at Risk) of cascading failure cost is utilized as the risk measure.
2) Two kinds of metaheuristic methods are proposed to achieve the robust design of interdependent networks efficiently even when an enormous number of initial node failures must be taken into account.
3) The characteristics of the achievable β-CVaR are clarified in relation to the interdependent networks with the spatial random topology through exhaustive simulation experiments.
4) The effectiveness of the proposed metaheuristic methods is verified using the simulation experiments.
The remainder of this paper is organized as follows.Section II explains related work as the background to this paper.Section III models the considered interdependent networks.Section IV formulates the robust design problem of the interdependent networks.Section V proposes two kinds of metaheuristic methods for designing the robust interdependent networks efficiently under a great number of initial node failures.Simulation experiments on the robust design are executed in Section VI.Finally, Section VII concludes this paper.

II. RELATED WORK
The cascading failures in the interdependent networks have been tackled from various aspects [8], [9], [10].Vulnerable nodes whose failures lead to the largest cascading failures were identified under the assumption of different forms of interdependency [11], [12], [13].However, two works [11] and [12] only focused on the interdependency between the multiple networks and disregarded the topology of each of the networks.A graph theoretic metric for measuring the cascade effect on the interdependent networks and significance metrics to assess the impact of the individual nodes on the spreading of failures have been proposed [14], [15], [16].A realistic model of the interdependent networks was proposed to study the evolution of the cascading failures [17].The robustness of interdependent networks has been analyzed on a basis of the percolation theory [18], [19], [20], [21].In the percolation theory, each node is regarded as being in a failure state when the node is not included in the largest connected component of each interdependent network.In particular, work [19] has analyzed the interdependent networks with the spatial random topology, which will be evaluated in this paper.A cascading failure mechanism considering the load on each node was evaluated using the simulation experiments [22].The dynamics of the cascading failures has been analyzed while taking the propagation speed into account [23].A cascading failure model considering the load on the nodes and the mode of the interdependency was evaluated for the interdependent networks with dependency groups [24].Vulnerability of various attacks has also been assessed experimentally [25].
Several design methods for the interdependent networks have been proposed to derive the optimum interdependency.The interdependency strategy against the random attacks has been derived theoretically and experimentally [26].Two metaheuristic methods for deriving the interdependency between the components in the cyber physical systems were proposed [27].An efficient method for restructuring the interdependency has been proposed to improve the survivability of the interdependent networks [28].However, the above three works only considered the cases where the topology of the individual networks is unknown.The interdependency to maximize the size of initial node failures that violate the connectivity in each of the networks has also been revealed [29].
The optimum interdependency was characterized by a networked Markov chain model for representing the propagation mechanism of the cascading failures [30].A design method for the robust interdependent networks has been derived according to the budget for structuring the interdependency [31].A design method has been proposed on a basis of the percolation theory under the assumption of different forms of interdependency and attacks [32].A design method based on metrics of the robustness reflecting the different propagation mechanisms of the cascading failures has been proposed [33].This method considered that a node also fails if its current load including the load offloaded from the failed nodes exceeds its capacity.A memetic algorithm for optimizing the interdependency was also proposed to enhance the robustness against the malicious attacks [34].Recently, an analytical method to deploy the optimum backup power supply for the communication network has been proposed [35].The above design methods regard all the outside nodes of the largest connected component of each network as being in a failure state on the basis of the percolation theory [30], [31], [32], [33], [34], [35].A heuristic method has minimized the required number of interconnections for preventing the cascading failures invoked from single node failures [36].The method guaranteed the connectivity of each distribution node with one of the generation nodes in the power grid.A dynamic game approach for designing the robust interdependent networks has also been proposed [37].
Several methods for recovering from the cascading failures have been proposed to improve the resiliency of the interdependent networks.The optimum scheduling for recovering the interdependent cyber-physical systems from a largescale failure has been derived [38].The dynamics of failure propagation and self-healing ability acquired by the cyber network was analyzed [39].A robustification method was proposed to rewire a small number of links in real-time during the progress of the cascading failures [40].A cascade mitigation method and a recovery method were proposed for the uncertain cascading failures [41].Algorithms for computing the most reliable paths and the most reliable pairs of paths in each interdependent network have also been proposed [42].
From the above literature review, no existing design method has recognized the intra-network support flow as one of the components for the interdependent networks to be operable [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37].This paper proposes a robust design method for the interdependent networks, where each sink node supporting the other network is assured to terminate a sufficient volume of support flow generated at the source nodes in the same network.In contrast to the existing recovery methods, the proposed design method enhances the robustness of the interdependent networks prior to the occurrence of the cascading failures.

III. MODEL OF INTERDEPENDENT NETWORKS
This paper considers a pair of interdependent networks supporting each other.In particular, a pair of power grid and computer network is supposed as a typical instance of the interdependent networks.The power grid transmits the electric power from the generation nodes to the distribution nodes, which provide each node in the computer network with the required electrical power.Meanwhile, each node composing the power grid is controlled by the edge nodes to which the required amount of control messages is transferred from the server nodes through the computer network.Hereinafter, the above intra-network flow of the electric power and the control messages is together named "support flow".In the power grid, source nodes and sink nodes of the support flow correspond to the generation nodes and the distribution nodes respectively.Meanwhile, they indicate the server nodes and the edge nodes in the computer network.Each sink node in a network is required to terminate the support flow whose volume is equivalent to the total support volume that the sink node provides for the other network.
From the above consideration, the following conditions must be met when the cascading failure has converged.
Condition 1) Each node in a network must receive the required volume of support from the working sink nodes in the other network.Otherwise, the node fails.
Condition 2) Each sink node must be connected with at least one of the working source nodes in the same network.Otherwise, the sink node fails.
Condition 3) Each working sink node in a network must be able to terminate a sufficient volume of support flow in order to support the working nodes in the other network.
Condition 1 is essential for the operable interdependent networks.If a sink node is disconnected from all the working source nodes in the same network, the sink node can be regarded as in a failure state since it loses its capability of supporting the other network.Condition 2 represents this constraint imposed on the sink nodes.The largest connected component of each network may include no source node of the support flow.Hence, it is not considered as a criterion of the node failure whether nodes are included in the largest connected component or not.Constraints on the capacity of the nodes and links for handling the support flow can be included in the model.The capacity of the nodes and links gives an influence on Condition 3 and subsequently affects Condition 1. Naturally, the failed node cannot forward the support flow generated at the working source nodes.Each node may be a source node of the support flow and a sink node of the support flow simultaneously.Furthermore, each node may be neither a source node nor a sink node.The above three conditions increase the possibility of an initial node failure within a small region invoking a cascading failure with a wide range.Hence, it is crucially important to minimize the significant impact of the cascading failures.Figure 1 shows an example of two interdependent networks supporting each other.Each node in network 0 and network 1 is a source node of the support flow or a sink node of the support flow.Two nodes in an identical network are connected using an intra-network link.The interdependency between network 0 and network 1 is represented by directional internetwork links from the sink nodes in one network toward the nodes supported by the sink nodes in the other network.In case of Fig. 1, an initial failure on node m6 invokes a cascading failure on all the nodes in both the networks.When sink node m6 fails, sink node n5 fails since it can receive no support from the working sink nodes in network 0 (Condition 1).When sink node n5 fails, sink node n6 fails and can terminate no support flow since it is disconnected from source nodes n1 and n2 (Conditions 2 and 3).Next, sink node m4 fails since it can receive no support from failed sink node n6 (Condition 1).Furthermore, sink node m5 also fails and can terminate no support flow since it is disconnected from source nodes m1 and m2 due to the failure in sink node m4 (Conditions 2 and 3).Consequently, sink nodes n3 and n4 fail since they can receive no support from failed sink node m5 (Condition 1).Then, all the remaining nodes m1, m2, and m3 in network 0 fail since they can receive no support from the working sink nodes in network 1 (Condition 1).Finally, all the remaining nodes n1 and n2 in network 1 also fail for the same reason (Condition 1).

A. Concept of Robust Design
The robust design of interdependent networks achieves the desirable interdependency between two networks when the topology of each of the two networks including the placement of source nodes and sink nodes of the support flow is specified beforehand.Aiming to maximize the robustness of the interdependent networks, a given number of sink nodes to support each node in the other network are selected from the candidate sink nodes specified beforehand.The robust design of the interdependent networks minimizes the risk of the cascading failures when a set of possible initial node failures Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and probability of their occurrence are given.A possible initial node failure may be a multiple-node failure.As the risk measure, β-CVaR (β-Conditional Value at Risk) of cascading failure cost is considered to minimize the impact of the cascading failures in the worst cases [43], [44], [45].Here, the cascading failure cost indicates the magnitude of the damage induced from each of the cascading failures and is calculated as the sum of the costs of the node failures involved in the considered cascading failure.The β-CVaR of the cascading failure cost is defined as the expected cascading failure cost when the cascading failure cost exceeds β-VaR (β-Value at Risk), which indicates the probabilistic upper bounds for the cascading failure cost.The cascading failure cost is less than the β-VaR with probability β whose value is specified as probability level.This means that the β-VaR is equivalent to 100β-percentile of the cascading failure cost.The β-CVaR always has an attractive mathematical property, i.e., coherency including the subadditivity, in comparison with the β-VaR [43].Hence, the robust design of the interdependent networks adopts the β-CVaR of the cascading failure cost as the risk measure to be minimized.
Figure 2 illustrates the β-CVaR of the cascading failure cost.Figure 2(a) shows the considered interdependent networks.Nodes m1 and n1 are the source nodes and the sink nodes simultaneously.All the other nodes are merely the sink nodes.In Fig. 2, all the single-node failures are considered as the initial node failures and each initial node failure is notated by the node that fails initially.Probability P t that each initial node failure t(∈ M ) occurs is given as shown in Fig. 2(a).
Here, notation M indicates the set of nodes comprising both the interdependent networks.Probability P t is normalized such that the total sum of probabilities ( t∈M P t ) is 1.0.Since the failure cost is assumed to be 1.0 in all the nodes, the cascading failure cost indicates the final number of filed nodes after the cascading failure has converged.For example, node n3 fails according to Condition 1 when an initial node failure m3 occurs.Consequently, nodes m4 and n4 also fail and can terminate no support flow according to Conditions 2 and 3.This means that the cascading failure cost induced from initial node failure m3 is 4.0.In Fig. 2(b), all the cascading failures are aligned from the right side to the left side in the ascending order of their costs.The vertical axis indicates cost e(t) of each cascading failure invoked by an initial node failure t(∈ M ).The horizontal axis indicates the cumulated probability P t of probability P t .In Fig. 2(b), the β-CVaR of the cascading failure cost indicates the average cascading failure cost in the shaded area.When the value of β is set at 0.75 in Fig. 2(b), the values of β-VaR and β-CVaR are 6.0 and 6.8 respectively.
When cascading failure cost e(t) and normalized probability of initial node failure occurrence P t are given, the value of β-CVaR (φ β ) can be derived from the following mathematical programming model including auxiliary variable d(t) [44].
Here, notation T indicates the set of initial node failures.In the model, the value of variable α β is identical to the β-VaR of the cascading failure cost when the objective function is minimized.

B. Formulation of Robust Design
The robust design problem of the interdependent networks can be formulated using a MILP (Mixed Integer Linear Programming) model on an augmented network.In the augmented network, an augmented node is added to each of the interdependent networks.Furthermore, the augmented node is connected with each of the source nodes of the support flow in the same network using an augmented directional link.Each intra-network link is represented by a pair of links with opposite directions.Several candidates sink nodes are specified for each node in the other network beforehand.The candidate sink nodes are indicated by candidate directional inter-network links from the candidate sink nodes to the node that may be supported by the candidate sink nodes.The optimum internetwork links can be selected from the candidate inter-network links by solving the MILP model.Figure 3 depicts an example of the augmented network for the interdependent networks shown in Fig. 1.
The MILP model for the robust design of interdependent networks can be formulated as follows.Table I summarizes notations for constants and sets used in the MILP model.Notation M i indicates the set of nodes comprising network i and excluding the augmented node.In contrast, notation E i indicates the set of links including the augmented link for network i. Set of nodes A i (m) includes the augmented node  for network i when node m is a source node.When the natural disaster in a certain area causes an initial node failure, nodes within the area fail initially.Multiple nodes or no node may be included in the area.This means that set of nodes F t specifying initial node failure t may include multiple nodes in both the networks or may be empty.Probability of initial node failure t occurring is indicated by notation P t and normalized so that the total sum of probabilities ( t∈T P t ) can be adjusted to 1.0.Failure cost Q i m represents the magnitude of the damage induced from the failure at node m in network i.  failure cost and can be expressed using the auxiliary variables, as shown in the previous subsection. (1) 2) Constraint on the Number of Selectable Sink Nodes: If each node can be redundantly supported by multiple sink nodes, the robustness of the node is improved against the cascading failures.Constraint on the number of sink nodes selectable for node m is given as follows.
3) Constraints Representing Condition 1 for Each Node (Part 1): Constraints indicating that each working node must be supported by at least one working sink node are expressed as follows, when initial node failure t occurs.3a) Constraint indicating that sink node n selected to support node m fails.
The values of variables g(t) i m,n are zero whenever possible to reduce the value of the objective function.Thus, they are one only when the values of both variables f (t) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The above two constraints can be replaced by the following constraint when k i m = 1.3c) Constraint indicating that node m fails if a sink node n selected to support node m fails.
4) Constraints Representing Condition 1 for Each Node (Part 2): Constraints indicating that each working node must be provided with the required support volume from the working sink nodes are expressed as follows, when initial node failure t occurs.4a) Constraint indicating that the support volume received from failed sink node n is zero.
4c) Constraint indicating that the support volume received from unselected sink node n is zero.
) Constraint for ensuring that working node m in network i can receive the required support volume from the selected working sink nodes in network 1 − i.
5) Constraints Representing Condition 2 for Each Sink Node: Constraints to guarantee the connectivity of each working sink node with one of the working source nodes are shown under the occurrence of initial node failure t.Each working sink node must be able to terminate the positive connectivity flow generated at the augmented node in the same network.5a) 5b) Constraints indicating that the connectivity flow passing through failed node m is zero.
5c) Constraint indicating that the connectivity flow is certainly terminated at working sink node m c .
5d) Constraint indicating the preservation of connectivity flow at working non-sink node m.
5e) Constraint indicating that a volume of connectivity flow starting from the augmented node is equivalent to the total volume of connectivity flow terminated at the working sink nodes.
6) Constraints Representing Condition 3 for Support Flow: Constraints on the support flow are given as follows, when initial node failure t occurs.In the power grid, the support flow is also required to satisfy the Kirchhoff's law [22].However, one of the feasible solutions of the MILP model achieving the same volume of support flow generated at each source node and terminated at each sink node and achieving the same values of binary variables f (t) i m satisfies the Kirchhoff's law certainly when the capacity of each node and each link to forward the support flow is unrestricted.The values of binary variables f (t) i m indicate the network topology after the convergence of the cascading failure.Since the capacity of each node and each link to forward the support flow is not considered in the formulation, the value of objective function in the formulated MILP model is not affected by the Kirchhoff's law.Hence, the MILP model include no constraint related to the Kirchhoff's law.6a) 6b) Constraints indicating that a volume of support flow passing through failed node m is zero.
6c) Constraint indicating that a volume of support flow terminated at working sink node m c in network 1 − i is equivalent to the total support volume provided by working sink node m c for network i.
6d) Constraint indicating the preservation of support flow volume at working non-sink node m.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
6e) Constraint indicating that a volume of support flow starting from the augmented node for network i is equivalent to the total support volume provided by all the working sink nodes in network i.
6f) Constraint on the capacity of source node m s for generating the support flow is given as follows.
) Constraint on Each Initial Node Failure: Constraint defining initial node failure t is given as follows.Each initial node failure t may include failed nodes in both the networks.
8) Constraint on the Cascading Failure Cost: Constraint for indicating the cascading failure cost caused by initial node failure t is expressed as follows.
9) Constraint on the Auxiliary Variables: As shown in Section IV-A, constraints on the auxiliary variables to compute the value of β-CVaR are given as follows.
The robust design of interdependent networks can be achieved by solving the above MILP model.The optimum interdependency to minimize the largest risk of the cascading failures is expressed by the values of binary variables x i m,n in the derived solution.Although the above MILP model includes the constraint on the capacity of each source node, constraint on the capacity of each sink node to terminate the support flow can be added to the model as a part of the constraints imposed on the support flow.Furthermore, constraint on the capacity of each node and each link to forward the support flow can be added to the MILP model.However, the model is required to consider the Kirchhoff's rule as additional constraints imposed on the support flow in case of the power grid [22].The budget constraint for structuring the interdependency can be added to the MILP model by specifying the cost of inter-network links beforehand.The MILP model can be solved exactly using the branch and bound method [46].

V. METAHEURISTIC METHODS FOR ROBUST DESIGN
The required computational effort for exactly solving the formulated MILP model increases due to an increase in the scale of the considered interdependent networks.Furthermore, the required computational effort increases rapidly as the number of possible initial node failures increases.Hence, two kinds of metaheuristic methods are proposed to obtain a suboptimal solution of the robust design problem with a small computational effort even when a great number of initial node failures must be taken into account.The two methods are chosen on a basis of the analysis of the structure of the formulated MILP model.

A. Lagrangian Relaxation Method Accompanied With Lagrangian Decomposition
All the constraints other than constraint (2) in the MILP model are only imposed on the individual initial node failures.Hence, the Lagrangian relaxation problem derived from the original MILP model due to the relaxation of a small number of constraints can be decomposed into multiple small subproblems related to the respective initial node failures [47].Although the number of subproblems increases due to an increase in the number of initial node failures, each subproblem can be solved with a small computational effort.This means that the Lagrangian relaxation problem can be solved easily even if a great number of initial node failures are specified beforehand.The Lagrangian relaxation method solves the Lagrangian relaxation problem repeatedly while updating the Lagrange multipliers.The Lagrange multipliers for maximizing the solution of the Lagrangian relaxation problem can be obtained by solving the Lagrangian dual problem.The Lagrangian dual problem can be solved instantly since it only includes the real variables.According to the weak duality theorem, the maximized solution of the Lagrangian relaxation problem can well approximate the optimal solution of the original problem [47].Figure 4 is an illustrative flow chart for the Lagrangian relaxation method.The procedure can be finished according to the current gap between the two solutions of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the Lagrangian relaxation problem and the Lagrangian dual problem.However, the accuracy of the solution achieved in a given computational time is evaluated in comparison to the other methods.Thus, the procedure finishes when the elapsed computational time (time) reaches the given allowed maximum time (Time).The followings are explanation on the detailed procedure for each step in Fig. 4.

1) Solving the Lagrangian Relaxation Problem (Step 1):
The Lagrangian relaxation problem is separated into multiple subproblems each of which computes the interdependency to minimize the cascading failure cost induced from an initial node failure.Since binary variable x i m,n in the original MILP model is independent from each initial node failure, the following binary variable is newly defined for the Lagrangian relaxation problem instead of binary variable x i m,n .x (t) i m,n : Binary variable indicating whether sink node n in network 1 -i (i = 0, 1) is selected to support node m in network i (= 1) or not (= 0) when initial node failure t occurs.
Furthermore, the following constraint is added to indicate that the value of binary variable x (t) i m,n must be identical between all the initial node failures.
The above constraint is relaxed in the Lagrangian relaxation problem since the problem considers the optimum interdependency for each initial node failure separately.Constraint (9a) is also relaxed in the Lagrangian relaxation problem since the objective function of the problem should include cascading failure cost e(t) to be minimized.
Hence, the Lagrange multipliers in the Lagrangian relaxation problem are defined as follows.

Minimize Obj
The above objective function can be transformed into the following objective function and three constraints on the Lagrange multipliers, as shown in the Appendix.

Minimize
The above constraints (13a) through (13c) are taken into account when the Lagrangian dual problem is solved.
The objective function ( 12) is formulated as the sum of the terms associated with the individual initial node failures.Hence, the minimization of the objective function ( 12) is equivalent to the minimization of each term related to an initial node failure.The Lagrangian relaxation problem is required to consider constraints (2) through ( 8) and (9b) in the original MILP model.All the above constraints can be classified into different groups, each of which is related to an initial node failure when binary variable x i m,n is replaced with binary variable x (t) i m,n .Consequently, the Lagrangian relaxation problem can be decomposed into multiple subproblems where the individual terms of the objective function (12) are minimized under the constraints included in each of the groups related to the respective initial node failures.Each small-scale subproblem can be solved exactly using the branch and bound method [46].All the Lagrange multipliers are constant in the Lagrangian relaxation problem.For example, all the Lagrange multipliers are set at zero at the beginning.
2) Solving the Lagrangian Dual Problem (Step 2): The Lagrangian dual problem computes the Lagrange multipliers to maximize the solution of the Lagrangian relaxation problem.When all the Lagrange multipliers γ(t) and μ(t) i m,n are simply represented by notations γ and μ respectively and the solution of the Lagrangian relaxation problem is indicated by notation y(γ, μ), function y(γ, μ) is convex upward.Hence, the value of function y(γ, μ) must be less than or equal to a tangent plane of function y(γ, μ) obtained from a solution of the Lagrangian relaxation problem.This means that each solution of the Lagrangian relaxation problem conducts an additional constraint related to the cutting plane for the Lagrangian dual problem.While the Lagrangian relaxation problem is solved repeatedly using the updated Lagrange multipliers, the Lagrangian dual problem is solved more accurately owing to an increase in the number of constraints.
The following notations are used for representing the conducted constraints related to the cutting plane.
LN: Number of the Lagrangian relaxation problems that have been solved repeatedly.
x ln (t) Then, the constraints related to the cutting plane are expressed as follows.
The Lagrangian dual problem is approximately solved by the bundle method aiming for the stable convergence [48].The objective function is represented by a quadratic function in the bundle method.Hence, the Lagrangian dual problem can be formulated using a quadratic programming model.The model can be solved efficiently with the use of the AMDD (Alternating Direction Method of Multipliers) [49].The Lagrange multipliers derived from the model are utilized in the subsequent Lagrangian relaxation problem.

3) Computation of Feasible Solutions (Step 3):
The solution of the Lagrangian relaxation problem may be unfeasible since a part of the constraints in the original MILP model are relaxed.Hence, a feasible solution of the original MILP model must be computed from the latest solution of the Lagrangian relaxation problem.The unfeasible solution of the Lagrangian relaxation problem may specify different interdependency for the individual initial node failures.The specified interdependency is more significant when the corresponding initial node failure induces a larger cascading failure cost.Thus, a restricted number (f sol ) of feasible solutions are derived from the interdependency specified for the initial node failures that cause larger cascading failure costs.Finally, one feasible solution achieving the smallest β-CVaR of cascading failure cost is selected as the final solution in the Lagrangian relaxation method.The value of parameter f sol is determined on a basis of the tradeoff relationship between the accuracy of the final solution and the additional computational effort.The β-CVaR of cascading failure cost in each feasible solution can be obtained from the original MILP model taking the following additional constraint into account.
Here, notation t 0 indicates a critical initial node failure selected for deriving a feasible solution.

B. Neighbor Local Search Method
When the number of possible initial node failures increases, the number of binary variables relying on the respective initial node failures increases.This increase in the number of binary variables mainly causes a large volume of computation.Thus, consideration to only a part of the binary variables relating with the critical initial node failures can reduce the required computational volume significantly without a great loss of the accuracy in the solution.The neighbor local search method  Since the accuracy of the solution attained in a given computational time is assessed in the following section, the procedure finishes instantly after the elapsed computational time (time) reaches the allowed maximum time (Time).
First, a trivial solution is set as an initial feasible solution of the MILP model.For example, the solution where all the binary variables f (t) i m are one is set as the initial feasible solution.A given number of initial node failures are selected and only the values of binary variables f (t) i m related to the selected initial node failures are changeable in the subsequent MILP models.All the critical initial node failures that have contributed to the β-CVaR of cascading failure cost in the previous solution are selected preferentially.When more initial node failures can be selected, additional initial node failures are randomly selected from the remaining ones aiming Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
at the extension of the range for seeking better solutions of the original MILP model.The modified MILP models where the unselected binary variables f (t) i m are regarded as the constants can be solved exactly using the branch and bound method [46].The modified MILP models have objective function (1) and include the following constraint in addition to constraints (2) through (9b) in the original MILP model.
In the above constraint, notation f 0 (t) i m indicates the value of binary variable f (t) i m in the previous solution and notation U indicates the set of initial node failures where the values of binary variables f (t) i m are changeable in the subsequent modified MILP model to be solved.The unselected binary variables f (t) i m (t ∈ T − U ) are fixed at values f 0 (t) i m due to constraint (16) in the subsequent modified MILP model.The number of changeable binary variables (|U|) is determined on a basis of the tradeoff relationship between the degree of the improvement in the updated feasible solutions and the required computational effort for solving the modified MILP models.
The solution of the subsequent MILP model is better than or equivalent to those of the preceding MILP models.However, the solutions of the modified MILP models have possibility of converging to a local optimal solution in the neighbor local search method.If the β-CVaR of cascading failure cost is not improved in a given number of continuous solutions (Iter), the procedure finishes instantly prior to the elapse of the allowed maximum time specified beforehand.

A. Simulation Settings
Simulation experiments related to the robust design of small-scale interdependent networks are executed to clarify the characteristics of the achievable β-CVaR of the cascading failure cost.The robust design of small-scale networks is achieved by solving the formulated MILP model exactly.The robust design of large-scale interdependent networks is also attempted by solving the MILP model approximately using the two metaheuristic methods.Simulation results on the robust design demonstrate the effectiveness of the proposed metaheuristic methods.Tables III and IV summarize the settings for the simulation experiments in small-scale networks and large-scale networks, respectively.All the considered networks have the connected spatial random topology in order to reflect the practical topology of the interdependent networks [8], [9], [10], [19].Each network is configured on a square region with a coordinate system.When the total number of nodes and the total number of links are specified, each connected spatial random network is configured on a basis of the Euclidean distance defined on the square region as follows [50].First, all the nodes are placed randomly on the square region.Next, the minimum spanning tree is constructed and pairs of nodes not still connected directly using a link and nearest to each other are sequentially connected using the remaining links.
A metropolitan area is considered as the square region on which each interdependent network is configured.Hence, the numbers of nodes and links in each network are assumed as  shown in Tables III and IV.As shown in Fig. 6, the square region is partitioned into multiple square sub-regions by dividing the region into an equal number of sub-fields vertically and horizontally.The number of square sub-regions is decided so that one or two nodes could be included in each sub-region on average.Each initial node failure represents node failures occurring within each of the square sub-regions.Hence, each initial node failure includes one or two failed nodes on average.However, a part of initial node failures may include no failed node or more than two failed nodes.For instance, the local storm and the flood damage within each sub-region cause an initial node failure.The numbers of candidate sink nodes and selectable sink nodes are uniform in all the nodes.The required support volume is set at 1.0 in all the nodes.The probability of each initial node failure occurring is set at an identical value.Since the failure cost is set at 1.0 in all the nodes, the cascading failure cost indicates the final number of failed nodes after the cascading failure has converged.The probability level (β) is set at 0.90 so that several critical initial node failures can always be assessed using the β-CVaR of the cascading failure cost.
All the nodes in the considered networks can be the sink nodes of the support flow.This means that all the nodes in a network have the capability of supporting the nodes in the other network.A given number of sink nodes closer to each node in the other network are considered as the candidate sink nodes for the node in the other network.Three kinds of methods are considered for placing the source nodes of the support flow in each network.A given number of nodes nearest to any one of the four corners of the square region are selected as the source nodes in "Corner" method.The generation nodes or substation nodes in the power grid may be placed using the "Corner" method when the electric power is supplied from the outside of the metropolitan area.In "Center" method, a given number of nodes nearest to the center of the square region are selected as the source nodes.The server nodes in the computer network may be deployed using the "Center" method when the data center is placed in the central office of the telecommunication company.In the end, a given number of nodes are randomly selected in "Random" method.All the nodes selected as the source nodes are the sink nodes simultaneously.
A CPU with a clock rate of 3.20 GHz and memory with a capacity of 32.0 GB are used as the experimental environment.The branch and bound method included in a Japanese commercial solver is adopted as the strict method [51].The AMDD for solving the Lagrangian dual problem is executed using the OSQP solver [52].The simulation programs for the two metaheuristic methods are written in Python language.

B. Simulation Results of Strict Design in Small-Scale Networks
Two sets of 20 pairs of interdependent networks, each of which comprises 16 nodes, are evaluated.The interdependent networks included in one set comprise 24 links, and those in the other set comprise 32 links.The number of initial node failures (|T|) taken into account is 25 in this subsection.The robust design of each pair of interdependent networks is achieved using the branch and bound method as one of the strict methods [46].The average values obtained from each set of the 20 pairs of interdependent networks are shown as the simulation results.The optimal solutions can be computed within 3600sec in all the evaluation scenarios.When the number of source nodes increases, the cascading failure cost generally decreases due to an improvement in the connectivity of each sink node with one of the working source nodes.The β-CVaR of cascading failure cost in the "Center" method is smaller than that in the "Corner" method.In the "Center" method, each sink node can connect with one of the working source nodes more easily since each source node has a larger node degree.In particular, the "Center" method is the most advantageous when the number of source nodes is only one.However, the difference between the "Center" method and the "Corner" method decreases when the number of source nodes increases.The "Center" method has a risk of one initial node failure invoking multiple simultaneous failures on most of the concentrated source nodes even if the number of source nodes increases.The β-CVaR of cascading failure cost shows intermediate characteristics between the "Corner" and the "Center" methods when the "Corner" method is adopted in one network and the "Center" method is adopted in the other network.When the number of source nodes is more than one, the "Random method" can minimize the β-CVaR of cascading failure cost since the source nodes are evenly distributed in the entire network.
2) Node Degree in Selected Sink Nodes: Figure 8 shows the average node degree in the selected sink nodes.The number of links (|L|) in each interdependent network is 24.The capacity of each source node is unrestricted in Fig. 8(a).Meanwhile, the required capacity is equally distributed into each of the source nodes in Fig. 8(b).In each figure, the average node degree is shown using the solid line when the number of selectable sink nodes (k) is only one, and using the broken line when the number of selectable sink nodes (k) is two.In general, the average node degree in the selected sink nodes is greater than that in the entire network.The sink nodes with a larger node degree are selected preferentially from the perspective of the connectivity with the working source nodes.This also causes that the average node degree in the selected sink nodes approaches that in the entire network due to an increase in the number of selected sink nodes.On the other hand, a given number of sink nodes closer to the source nodes tend to be selected as those for supporting each node actually.Thus, the average node degree in the "Corner" method is smaller than that in the "Center" method.In particular, the Fig. 8. Average node degree in selected sink nodes.average node degree in the "Corner" method is smaller than that in the entire network when the number of source nodes is one or two.The difference in the average node degree is reduced between the placement methods of the source nodes when the number of source nodes increases.This is because the placement of the source nodes approaches the uniform arrangement in all the placement methods due to an increase in the number of source nodes.

C. Simulation Results in Large-Scale Networks Accompanied With Numerous Initial Node Failures
The robust design of interdependent networks is executed in terms of a set of 20 pairs of interdependent networks comprising 32 nodes and 48 links.The total number of initial node failures (|T|) taken into account increases to 49 in this subsection.The robust design is executed using the Lagrangian relaxation method and the neighbor local search method.The robust design is also attempted using the branch and bound method [44].In the Lagrangian relaxation method, 16 feasible solutions are derived from the final solution of the Lagrangian relaxation problem aiming to obtain the best feasible solution certainly with the least additional computational effort.The values of binary variables related to only 16 initial node failures are changeable in the modified MILP models that the neighbor local search method solves repeatedly.This number of initial node failures represents the maximum number in the MILP model for which the branch and bound method can compute the optimal solution.This means that the neighbor local search method prioritizes the improvement of the feasible solutions achieved by the subsequent modified MILP models.In this subsection, the "Random" method is adopted for the placement of the source nodes since the method achieves better robustness when the number of source nodes is more than one.The allowed maximum computational time (Time) is set at 3600sec since the optimal solutions for all the small-scale interdependent networks can be computed within 3600sec in the previous subsection.
1) Limitation of Strict Method: Figure 9 shows the ratio of the MILP models for which the branch and bound method can obtain a feasible solution and the optimal solution within the allowed maximum computational time (3600sec).The ratio of a feasible solution being acquired is shown using the solid line and the ratio that the optimal solution can be computed is shown using the broken line.The capacity of source nodes is unrestricted in Fig. 9(a) and restricted in the same manner as the previous subsection in Fig. 9(b).As shown in Fig. 9, a feasible solution can only be obtained for a part of the MILP models and the optimal solution cannot be computed in most of the MILP models when one sink node is selectable to support each node (k = 1).A feasible solution can be obtained for most of the MILP models when the number of selectable sink nodes (k) is more than one.However, the optimal solution can be computed merely in a part of the MILP models when the capacity of source nodes is restricted.Unlike the branch and bound method, the Lagrangian relaxation method and the neighbor local search method can always obtain a feasible solution.
Figure 10 shows the average β-CVaR of cascading failure cost obtained from the part of the MILP models for which the branch and bound method can acquire a feasible solution.For comparison, Fig. 10 also shows the average β-CVaR derived from the identical MILP models using the Lagrangian relaxation method and the neighbor local search method.The average β-CVaR is shown using the solid line, broken line, and chain line when the number of selectable sink nodes (k) is one, two, and three, respectively.The capacity of source nodes is unrestricted in Fig. 10(a) and restricted in Fig. 10(b).As shown in Fig. 10, the average β-CVaR of cascading failure cost in the branch and bound method is larger than those in the two metaheuristic methods except when the capacity of source nodes is unrestricted and the number of selectable sink nodes is more than one.This reveals that the feasible solution derived using the branch and bound method is often inferior to those in the two metaheuristic methods even if the branch and bound method can acquire a feasible solution.When the capacity of source nodes is unrestricted and the number of selectable sink nodes is more than one, the branch and bound method can acquire the optimal solution in most of the MILP models.However, the two metaheuristic methods can also obtain the optimal solution in most of the MILP models and the average β-CVaR of cascading failure cost is almost identical between the branch and bound method and the two metaheuristic methods in this scenario.
2) Comparison of Two Metaheuristic Methods: Figure 11 shows the average β-CVaR of cascading failure cost obtained using the Lagrangian relaxation method and the neighbor local search method.Figure 11 shows the average β-CVaR derived from all the MILP models formulated for the 20 pairs of interdependent networks.The average β-CVaR is shown using the solid line, broken line, and chain line when the number of selectable sink nodes (k) is one, two, and three, respectively.As shown in Fig. 11(a), the average β-CVaR of cascading failure cost in the neighbor local search method is smaller than that in the Lagrangian relaxation method when the capacity of source nodes is unrestricted and the number of selectable sink nodes (k) is one.When the number of selectable nodes (k) is more than one, both the methods can acquire the optimal solution and the average β-CVaR are identical between the two methods.As shown in Fig. 11(b), the Lagrangian relaxation method is more advantageous due to an increase in the number of source nodes when the capacity of source nodes is restricted and the number of selectable sink nodes (k) is only one.When the number of selectable nodes (k) is more than one, the Lagrangian relaxation method can always acquire more accurate feasible solution.Although the approach to the global optimal solution is guaranteed in the Lagrangian relaxation method, the neighbor local search method is in danger of being trapped into a local optimal solution.From the simulation results, the neighbor local search method is easily trapped when a volume of support flow traversing each route is restricted and the diversity of the routes for the support flow increases.In those cases, the Lagrangian relaxation method can maximize the robustness of the interdependent networks accompanied with an enormous number of initial node failures.

VII. CONCLUSION
This paper proposed a robust design method for the interdependent networks to minimize the largest impact of the cascading failures invoked by a set of possible initial node failures.The robust design problem taking the constraints on the support flow into account was formulated using a MILP model, and two kinds of metaheuristic methods were proposed for solving the problem efficiently even when a great number of initial node failures are specified.The robust design was executed in relation to the interdependent networks with the connected spatial random topology.From the results of the robust design, it was clarified that the size of the cascading failures can be reduced by placing a source node of the support flow at the center of the network and arranging the multiple source nodes evenly in the entire network.In general, sink nodes with a larger degree and closer to the source nodes should be selected to support the other network.The proposed metaheuristic methods can acquire the suboptimal solutions even when the strict method cannot obtain a feasible solution due to an increase in the number of initial node failures.Even if the strict method can acquire a feasible solution, the solution may be inferior to ones that the two metaheuristic methods can obtain.Further evaluation is required for accessing the impacts of the capacity restriction on each node and each link.The influence from the budget limitation to achieve the interdependency also needs to be evaluated.Since the proposed metaheuristic methods need to solve the small-scale MILP models iteratively, they may be unable to solve the robust design problem with a huge size.Finding an approximation algorithm to achieve the robust design in a polynomial time remains for further study.If the following constraints on the Lagrange multipliers are not satisfied, the objective function can be reduced unlimitedly since no upper bound is specified for variables α β and d(t).Thus, the following constraints should be satisfied so that the objective function to be minimized can be meaningful.
When constraints (A1) and (A2) are satisfied, the above objective function is minimized by setting the values of variables α β and d(t) at zero.Hence, the objective function is expressed as follows.In the above objective function, the value of expression within parentheses is invariant even if an arbitrary constant is added to all the Lagrange multipliers included in a set {μ(t) i m, n |t ∈ T }.Hence, the following constraint can be assumed.
Consequently, the objective function of the Lagrangian relaxation problem is converted into the following objective function (A4) and constraints (A1) through (A3).former senior researcher of KDDI Research Inc. for their encouragement throughout the study.
Figure 5 is an illustrative flow chart for the neighbor local search method.

1 )
Minimized β-CVaR of Cascading Failure Cost: Figure 7 shows the minimized β-CVaR of cascading failure cost obtained from the accurate solution of the MILP model.The number of links (|L|) in each interdependent network is 24 in Figs.7(a) and 7(b), and 32 in Figs.7(c) and 7(d).In Figs.7(a) and 7(c), the capacity of each source node is unrestricted.This means that constraint (6f) is excluded in the MILP model to be solved.Meanwhile, the required capacity for generating the support flow is equally distributed into each of the source nodes in Figs.7(b) and 7(d).In each figure, the β-CVaR of cascading failure cost is shown using the solid line when the number of selectable sink nodes (k) is only one.Meanwhile, the β-CVaR of cascading failure cost is shown using the broken line when the number of selectable sink nodes (k) is two.Naturally, the cascading failure cost can be reduced due to an increase in the number of sink nodes selected for each node.The β-CVaR of cascading failure costs in Figs.7(c) and 7(d) somewhat decrease in comparison with those in Figs.7(a) and 7(b) respectively.This is because the connectivity of each sink node with the working source nodes is improved due to an increase in the average node degree in the interdependent networks.The β-CVaR of cascading failure costs in Figs.7(b) and 7(d) increase compared with those in Figs.7(a) and 7(c) respectively.This indicates that each sink node is required to terminate a sufficient volume of support flow in addition to its connectivity with one of the working source nodes.The connectivity of each sink node with the working source nodes with a sufficient residual capacity is easily broken in Figs.7(b) and 7(d).

Fig. 10 .
Fig. 10.Average β-CVaR of cascading failure cost obtained from a part of MILP models.

Fig. 11 .
Fig. 11.Average β-CVaR of cascading failure cost obtained from all the MILP models.

= α β 1
APPENDIX TRANSFORMATION OF OBJECTIVE FUNCTION IN LAGRANGIAN RELAXATION PROBLEMThe objective function of the Lagrangian relaxation problem can be transformed as follows.Obj lr = α β + t)(e(t) − α β − d(t)) t∈T )e(t) − i∈{0,1} m∈M i n∈N i S (m) μ(t) i m,n x (t) i m,n ⎞ ⎟ ⎠ (A4)ACKNOWLEDGMENTThe author would like to gratitude Dr. Nakamura, president & CEO, Dr. Otani, former executive director, and Dr. Kitahara, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I NOTATIONS
FOR CONSTANTS AND SETS

TABLE II NOTATIONS
FOR VARIABLES

TABLE III SIMULATION
SETTINGS FOR SMALL-SCALE NETWORKS

TABLE IV SIMULATION
SETTINGS FOR LARGE-SCALE NETWORKS