The multi-terminal maximum-flow network-interdiction problem

European Journal of Operational Research 211 (2011) 241–251 Contents lists available at ScienceDirect European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor Discrete Optimization The multi-terminal maximum-flow network-interdiction problem _ Ibrahim Akgün a,⇑, Barbaros Ç. Tansel a, R. Kevin Wood b a b Department of Industrial Engineering, Bilkent University, Bilkent 06800, Ankara, Turkey Operations Research Department, Naval Postgraduate School, Monterey, CA 93943, USA a r t i c l e i n f o Article history: Received 23 May 2007 Accepted 19 December 2010 Available online 24 December 2010 Keywords: OR in military Integer programming Network flows Network interdiction a b s t r a c t This paper defines and studies the multi-terminal maximum-flow network-interdiction problem (MTNIP) in which a network user attempts to maximize flow in a network among K P 3 pre-specified node groups while an interdictor uses limited resources to interdict network arcs to minimize this maximum flow. The paper proposes an exact (MTNIP-E) and an approximating model (MPNIM) to solve this NP-hard problem and presents computational results to compare the models. MTNIP-E is obtained by first formulating MTNIP as bi-level min–max program and then converting it into a mixed integer program where the flow is explicitly minimized. MPNIM is binary-integer program that does not minimize the flow directly. It partitions the node set into disjoint subsets such that each node group is in a different subset and minimizes the sum of the arc capacities crossing between different subsets. Computational results show that MPNIM can solve all instances in a few seconds while MTNIP-E cannot solve about one third of the problems in 24 hour. The optimal objective function values of both models are equal to each other for some problems while they differ from each other as much as 46.2% in the worst case. However, when the postinterdiction flow capacity incurred by the solution of MPNIM is computed and compared to the objective value of MTNIP-E, the largest difference is only 7.90% implying that MPNIM may be a very good approximation to MTNIP-E. Ó 2010 Elsevier B.V. All rights reserved. 1. Introduction This paper investigates what we call the multi-terminal maximum-flow network-interdiction problem (MTNIP) for which no previous work exists. In MTNIP, there are two opponents, a network user/defender and an interdictor/attacker. The network user wishes to maximize flow among K P 3 node groups in an undirected network while the interdictor tries to minimize the network user’s maximum flow by using limited interdiction resources (e.g., aerial sorties, missiles) to destroy the arcs of the network. Consider two hostile forces AT and DF where AT is the attacker/interdictor on a communication (transportation) network and DF is the defender/network user. AT’s interest is to minimize DF’s inter-force communication (transportation) capabilities by attacking a subset of DF’s communication (transportation) lines. The set of lines that can be attacked is limited by the availability of AT’s interdiction resources. Locations of DF’s forces may or may not be precisely known to AT. If they are precisely known, they are taken to be source and sink locations that mutually exchange information (materials). In the remaining case, we assume that AT has sufficient information to confine these locations to K node groups that are taken to be source and sink groups among which information (material) exchange takes place. AT’s problem is to identify a set of arcs whose deletion from the network limits DF’s ability to transfer flows (signals/materials) between exact or possible sources and sinks while DF aims to maximize flow through the intact part of the network. This problem can be modeled as a bilevel min–max problem where the inner maximization is a flow maximization problem given that a subset of arcs is interdicted while the outer minimization involves the minimization of the maximum objective value of the inner maximization over the set of binary vectors each satisfying the upper bound on the interdiction resource. Each binary vector specifies which arcs to be interdicted and which ones to be left intact. The resulting bi-level min–max problem is what we refer to as MTNIP. The mathematical formulation will be given in Section 2. MTNIP is a useful model to analyze possible courses of action to protect critical infrastructures against possible terrorist attacks. Such critical infrastructures may include telecommunication lines, power lines, subways, highways, energy delivery lines (e.g., natural gas, petroleum), and the like. ⇑ Corresponding author. Tel.: +90 312 290 1262; fax: +90 312 266 4054. _ Akgün). E-mail address: iakgun@bilkent.edu.tr (I. 0377-2217/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2010.12.011 242 _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. MTNIP is a generalization of the maximum-flow network-interdiction problem (MFNIP) (e.g., Wood, 1993). The problem context in MFNIP is the same as the one in MTNIP except that the interdictor tries to minimize the maximum flow from a source node s to a sink node t instead of among three or more groups of nodes. That is, MFNIP is a special case of MTNIP with K = 2. MFNIP is proven to be NP-hard even if a single unit of resource is required to interdict each arc (Wood, 1993) and hence MTNIP is also NP-hard. The goal of this paper is to extend the single-commodity structure of MFNIP to a multi-commodity structure so that more realistic and general problem settings can be handled. Even though MTNIP is new to the literature, MFNIP is well-studied. The notable contributions are Wollmer (1964, 1970a), McMasters and Mustin (1970), Helmbold (1971), Ghare et al. (1971), Lubore et al. (1971), Wood (1993), Cormican et al. (1998), and Whiteman (1999). Almost all studies prior to Wood (1993) are specific to the application and are not extendible to more general contexts. Wood (1993) is the first to adopt mathematical programming methods. He develops a min–max formulation of MFNIP and then converts it to an integer-programming model. Cormican et al. (1998) study a stochastic variation of MFNIP. Whiteman (1999) adapts Wood (1993)’s model to select target sets. Another category of network-interdiction problem is that of maximizing the shortest path (MXSP) in which a set of arcs is disabled to maximize the length of a shortest path between s and t through the usable portion of the network. Notable contributions are Fulkerson and Harding (1977), Golden (1978), Israeli (1999), and Israeli and Wood (2002). Lim and Smith (2007) study a multi-commodity network-interdiction problem where the network user makes profit by delivering multiple commodities to certain destinations while the interdictor tries to minimize the network user’s profit by destroying arcs. The authors develop two models and a partitioning algorithm along with a heuristic procedure for the partial and complete interdiction of arcs, respectively. This study is closer to ours than others due to the multi-commodity flow structure; however, the two problems are structurally different. Other studies similar in spirit to MTNIP but applied to different fields are as follows. Wollmer (1970b) and Washburn and Wood (1994) develop game-theoretic network-interdiction models. Assimakopoulos (1987) suggests an interdiction model for preventing hospital infections. Anandalingam and Apprey (1991) investigate conflict resolution problems. Church et al. (2004) study the interdiction of supply and emergency response facilities. Salmeron et al. (2004) study the disruptions to electric power grids. Brown et al. (2005) describe a model for planning the pre-positioning of defensive missile interceptors. Brown et al. (2006) apply optimization models to make critical infrastructure more resilient against attacks. Scaparra and Church (2008) study the problem of allocating protective resources among the facilities of a system. Smith et al. (2007) examine the problem of fortifying a network to defend against attacks in the context of survivable network design (e.g., Alevras et al., 1998; Myung et al., 1999; Ouveysi and Wirth, 1999). Desai and Sen (2010) consider the problem of designing reliable networks that satisfy several constraints while simultaneously allocating multiple resources to mitigate the arc failure probabilities such that the total cost of network design and resource allocation is minimized. In the remainder of this paper, we (1) develop an exact formulation for MTNIP, (2) develop an approximating formulation for MTNIP, (3) present computational results, and (4) conclude with further research directions. 2. Exact formulation of MTNIP MTNIP is defined on a capacitated, undirected network G = (N, A) with node set N and arc set A consisting of unordered pairs of distinct nodes. Flow on (i, j) 2 A can move from i to j or from j to i. The total flow on (i, j) 2 A, defined by the sum of flows from i to j and from j to i, is restricted by a positive integral capacity uij. The network user aims to maximize total flow among K P 3 disjoint, pre-specified node groups N 01 # N; . . . ; N 0K # N where each node S group acts both as a source and a sink. We define N 0 ¼ k N 0k to be special nodes and N  N0 to be regular nodes. It is also natural to as0 sume that jN j  jNj. The total flow among K node groups is taken to be the sum of K single-commodity flows distinguished by their source groups and restricted by joint capacity constraints. The kth single-commodity flow originates in nodes in N 0k and is delivered through the arcs of the network to nodes in N 0  N 0k . That is, node group N 0k is a source for commodity k and a sink for any other commodity k0 – k. Maximization of the sum of K single-commodity flows is equivalent to maximizing the total flow routed between K (K-1)
pairs of node groups N 0k ; N 0k with k – k0 . The network user’s resulting problem is a multi-commodity maximum-flow problem (Costa et al., 2005). In MTNIP, an interdictor aims to minimize the maximum flow achievable by the network user by destroying arcs. We assume that the interdictor uses a single type of interdiction resource with a total of R units. Interdicting an arc (i, j) 2 A requires rij > 0 units of the resource. Partial interdiction of an arc is not allowed, i.e. an arc is either interdicted or not interdicted. In MTNIP, the network user and the interdictor engage in a two-step, sequential decision-making process: the attacker first allocates limited interdiction resources to destroy arcs so that the maximum flow achievable by the network user is minimized and then the network user maximizes flow through the network given the interdiction decisions of the attacker. In this sense, the interdictor is the leader and the network user is the follower. This leader–follower relationship is similar to the one in a static Stackelberg game (Siman and Cruz, 1973) except that a more general Stackelberg game continues in alternating plays between the leader and the follower. Such a game can be expressed mathematically as a bi-level programming problem (Dempe, 2002). In accordance with this, the interdictor’s problem MTNIP is modeled as a bi-level min–max program. Later, we convert it into a mixed-integer linear program (MIP). In the following subsections, the network user’s and interdictor’s problems are modeled, respectively. In the rest of the paper, zP and Z LP P will represent the optimal objective function values for P and for the linear programming relaxation of P, respectively. Similarly, xP and xLP P will represent the optimal solutions of P and the LP relaxation of P, respectively. 2.1. The formulation of the network user’s problem and its dual The network user’s problem is modeled as a multi-commodity maximum flow problem (MXF). Let yijk and yjik be the amounts of flow, respectively, from node i to node j and from node j to i on arc (i, j) for which the source node is any node in N 0k . _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. 243 2.1.1. Model MXF: Network user’s multi-commodity maximum-flow model z ¼ max y X k¼1;...;K 0 @ X s:t: X k¼1;...;K ði;jÞ2A:j2N 0k yijk  j:ði;jÞ2A[ðj;iÞ2A X yijk þ ði;jÞ2A:i2N 0k X X 1 yjik A; yjik ¼ 0; ð1Þ k ¼ 1; . . . ; K; i 2 N  N0 : aik ; ð2Þ j:ði;jÞ2A[ðj;iÞ2A
yijk þ yjik 6 uij ; ði; jÞ 2 A : bij ; yijk P 0; yjik P 0; ð3Þ k ¼ 1; . . . ; K; ði; jÞ 2 A: ð4Þ N 0k MXF is the network user’s multi-commodity maximum-flow model to maximize flow among node groups when there is no interdiction. The objective function (1) maximizes the sum of flows originating from each node group N 0k . Constraints (2) are flow-balance constraints for regular nodes. The arc-capacity constraints (3) restrict the amount of flow on each arc to the arc’s nominal capacity. Constraints (4) are non-negativity constraints. We note that MXF can also be modeled by creating a super source connected to N 0k , a super sink connected to N 0  N 0k , and by maximizing the sum of flows on return arcs from super sinks to super sources. Even though this is the more common approach in most network flow formulations (e.g., Ahuja et al., 1993), we prefer to maximize the sum of flows leaving each node group N 0k in our formulation. This eliminates the flow-balance constraints for nodes in N0 . We assume in our formulation that no flow of commodity k occurs within the node group N 0k and in arcs leading from nodes outside of 0 N k into nodes in N 0k . That is, yijk  yjik  0 for i; j 2 N 0k ; k ¼ 1; . . . ; K, and yijk  0 for j 2 N 0k ; i 2 N  N 0k ; k ¼ 1; . . . ; K. We further assume that nodes in N 0k are neither sources nor transshipment nodes for any other commodity k0 – k. Accordingly, yijk  0 for i 2 N 0k0 ; j 2 N  N 0k for each k0 – k. These assumptions can be incorporated into the model by preprocessing the data regarding the network structure. We may define, for example, a three-dimensional matrix A whose rows and columns are associated with the nodes of the network and whose layers are associated with the commodities so that the entry aijk takes on the value of 1 if flow is allowed from node i to node j for commodity k and 0 otherwise. Then, set (1) aijk = 0 for i; j 2 N 0k ; k ¼ 1; . . . ; K, (2) ajik = 0 for i 2 N 0k , j is an element of N  N 0k ; k ¼ 1; . . . ; K, and (3) aijk0 ¼ 0 for i 2 N 0k , 0 j is an element of N  N 0k ; k –k; k ¼ 1; . . . ; K. Next, we give the dual problem D-MXF associated with MXF and derive some results about the properties of the dual variables that will be used later. 2.1.2. Model D-MXF: The dual of the multi-commodity maximum flow model MXF z ¼ min a;b s:t: X uij bij ; ð5Þ ði;jÞ2A  aik þ ajk þ bij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i; j 2 N  N0 ; ð6Þ  ajk þ aik þ bij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i; j 2 N  N0 ; ð7Þ aik þ bij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 N  N0 ; j 2 N0  N0k ; ajk þ bij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; j 2 N  N0 ; i 2 N0  N0k ;  ajk þ bij P 1; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 N 0k ; j 2 N  N0 ;  aik þ bij P 1; k ¼ 1; . . . ; K; ði; jÞ 2 A; j 2 N0k ; i 2 N  N0 ; ð11Þ k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 N 0k ; j 2 N  N0k ð12Þ bij P 1; 0 and j 2 N0k ; i 2 N0  N0k ; ð8Þ ð9Þ ð10Þ aik ; free k ¼ 1; . . . ; K; i 2 N  N ; ð13Þ bij P 0; ð14Þ ði; jÞ 2 A: In D-MXF, aik and bij are dual variables for constraints (2) and (3), respectively. bij can be viewed as a distance label on the arc (i, j) and aik as the potential corresponding to commodity k on node i. Thus, the dual problem is an assignment of potentials to non-terminal/non-special nodes (a zero potential is assigned to terminal/special nodes) and non-negative distance labels to arcs. We observe that there is an optimal solution to D-MXF such that 1 6 aik 6 0, "i 2 N, k = 1, . . . , K and 0 6 bij 6 1, "(i, j) 2 A. This is justified by observing that the coefficient of bij is positive in the objective function so that making each bij as small as possible as permitted by the constraints does not cause a loss of optimality. Let A1, A2, A3 and A4 be the sets of arcs defined for the constraint pairs (6)–(11), and for  (12), respectively. That is, A1 ¼ fði; jÞ 2 A : i; j 2 N  N 0 g; A2 ¼ ði; jÞ 2 A : i 2 N  N 0 ; j 2 N 0  N 0k ; A3 ¼ fði; jÞ 2 A : i 2 N k0 ; j 2 N  N 0 g, and A4 ¼ fði; jÞ 2 A : i 2 N k0 ; j 2 N 0  N k0 and j 2 N k0 ; i 2 N 0  N 0k0 g. Clearly, A1, A2, A3 and A4 partition the arc set A into four disjoint subsets. Observe that no two variables aik and aik0 with the same node index i but different commodity indices k and k0 appear in the same constraint. Accordingly, the restriction of a variable aik to the interval [1, 0] does not affect any other aik for k – k0 . Constraints (6) and (7) imply that, for each arc (i, j) 2 A1, the variable bij is bounded below by the maximum of aik  ajk and aik + ajk for k = 1, . . . , K. Hence, bij is bounded below by the maximum over k of these bounds. The restriction of the variables aik and ajk to the interval [1, 0] for these arcs implies that the lower bound on bij enforced by constraints (6) and (7) is at most 1. Accordingly, restricting bij to the interval [0, 1] for such arcs maintains feasibility without causing a loss of optimality. Similarly, constraints (8) and (9) imply that bij is bounded below by the maximum of  aik and ajk for (i, j) 2 A2 and k = 1, . . . , K. With the restriction of the variables aik and ajk to the interval [1, 0] for these arcs, the implied lower bound is at most 1. Hence, we may again restrict bij to the interval [0, 1] for the arc group A2. For arcs (i, j) 2 A3, the constraints (10) and (11) imply bij is bounded below by the maximum of 1 + aik and 1 + ajk for k = 1, . . . , K. Restriction of the variables aikand ajk to the interval [1, 0] for these arcs implies that the maximum of these lower bounds is again at most 1. Hence, we may restrict bij to the interval [0, 1] for these arcs as well. Constraints (12) imply bij is bounded below by 1 for arcs in A4. Restriction of bij to the interval [0, 1] for these arcs does not _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. 244 cause loss of optimality and yields bij = 1 in an optimal solution for these arcs. This concludes the justification of our initial claim that there is an optimal solution to D-MXF such that all aikvalues are restricted to the interval [1, 0] and all bij values are restricted to the interval [0, 1]. This result will be useful in converting the bi-level programming formulation of MTNIP into a mixed integer linear program (MIP). Due to the duality relationship between maximum flow and minimum cut problems, D-MXF is closely related to the NP-hard minimum multi-way (multi-terminal) cut problem (MMCP). See, for example, Dahlhaus et al. (1994) and Costa et al. (2005). In MMCP, the purpose is to find a set of arcs with minimum total capacity whose removal from G puts each terminal (specific node groups) in a different connected component Gk = (Nk, Ak),k = 1, . . . , K. When K = 2, the well-known maximum-flow minimum-cut theorem (Ford and Fulkerson, 1956) holds and the set of saturated arcs in a maximum flow identifies also a minimum cut for the dual problem. Hence, the minimum two-way cut in MMCP is directly available as an optimal solution to D-MXF for K = 2. On the other hand, the maximum flow among K P 3 node groups need not be integral and does not in general give a multi-way minimum cut solution. In this case, the strong duality holds for MXF and D-MXF (e.g., Garg et al., 1996; Costa et al., 2005) and this implies that an integral solution xDMXF gives the minimum multi-way cut for MMCP. 2.2. The formulation of the interdictor’s problem In this subsection, the interdictor’s problem is modeled as a bi-level, min–max program and then converted into a MIP. The interdictor’s decision variable xij takes on the value of 1 if arc (i, j) is interdicted and 0 otherwise. 2.2.1. Model MTNIP-BI: MTNIP Formulation as a Bi-level Program X  z ¼ min max x2X y k¼1;...;K 0 @ X yijk þ ði;jÞ2A:i2N0k X ði;jÞ2A:i2N 0k 1 yjik A; ð15Þ s:t: Constraints ð2Þ; ð4Þ; and X

yijk þ yjik 6 uij 1  xij ði; jÞ 2 A : hij ; ð16Þ k¼1;...;K where X¼ ( x 2 f0; 1g jAj : X ði;jÞ2A ) rij xij 6 R : ð17Þ MTNIP-BI is the interdictor’s model to minimize the maximum flow achievable in MXF. For fixed x, the inner maximization is the network user’s maximum-flow model. The objective function (15) minimizes the maximum flow among the subsets N 0k . Constraints (16) set P k¼1;...;k ðyijk þ yjik Þ to zero when xij = 1 and to uij when xij = 0. Constraints (17) limit the expenditure of interdiction resource and require interdiction variables to be binary. MTNIP-BI is impossible to solve with standard optimization software. It may be possible to solve it by developing specialized decomposition techniques as offered by Israeli and Wood (2002). However, we prefer a simpler method that allows us to convert MTNIP-BI into a MIP and then to solve it directly by using standard software. Our method consists of (1) taking the dual of the inner maximization by fixing x temporarily and then releasing x to obtain a mixedinteger nonlinear ‘‘min–min’’ model, which is simply a minimization model, and (2) linearizing the nonlinear model to get a MIP. 2.2.2. Model MTNIP-MINP: MTNIP formulation as a mixed integer nonlinear program z ¼ min min x s:t: a;h X uij ð1  xij Þhij ; ð18Þ ði;jÞ2A  aik þ ajk þ hij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i; j 2 N  N0 ; ð19Þ  ajk þ aik þ hij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i; j 2 N  N0 ; ð20Þ 0 0 N  N0k ; N0  N0k ; 0 aik þ hij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 N  N ; j 2 ajk þ hij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; j 2 N  N0 ; i 2  ajk þ hij P 1; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 N0k ; j 2 N  N ;  aik þ hij P 1; k ¼ 1; . . . ; K; ði; jÞ 2 A; j 2 N0k ; i 2 N  N0 ; hij P 1; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 X r ij xij 6 R; N0k ; 0 j2N  N0k and j 2 ð21Þ ð22Þ ð23Þ ð24Þ N0k ; 0 i2N  N0k ; ð25Þ ð26Þ ði;jÞ2A xij 2 f0; 1g ði; jÞ 2 A; ð27Þ 0 aik ; free k ¼ 1; . . . ; K; i 2 N  N ; ð28Þ hij P 0; ð29Þ ði; jÞ 2 A; Let MXF (x) be the version of MXF with upper bounds uij in MXF replaced by the upper bounds uij(1  xij) and let D-MXF (x) be the dual of MXF (x). That is, MXF (x) is the inner maximization in MTNIP-BI defined by (15), (2), (4), and (16) and D-MXF (x) is the inner minimization in MTNIP-MINP defined by (18)–(25), (28), and (29). The dual variables aik and hij in MTNIP-MINP correspond to constraints (2) and _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. 245 (16), respectively. Observe that hij plays in D-MXF (x) the role of bij in D-MXF so that the restriction of hij to the interval [0, 1] does not cause a loss of optimality in D-MXF (x). This restriction allows us to use the linearization that replaces (1  xij)hij with gij P 0 and adding the set of constraints gij P hij  xij. This yields the following MIP. 2.2.3. Model MTNIP-MILP: MTNIP as a mixed-integer linear program z ¼ min a;h;x;g X uij gij ð30Þ ði;jÞ2A s.t Constraints (17)–(27) and gij P hij  xij ði; jÞ 2 A; gij P 0 ði; jÞ 2 A; ð31Þ ð32Þ If xij = 0 in an optimal solution to MTNIP-MINP, the corresponding term in the objective function (18) is equal to uijhij. If xij = 1 in an optimal solution, then the corresponding term in (18) is 0. Thus, for the linearization to work, it must be true that gij = 0 when xij = 1 and that gij = hij when xij = 0. When xij = 1, constraints (31) are satisfied for 0 6 hij 6 1 and for 0 6 hij 6 0. However, because setting gij to any value greater than 0 unnecessarily increases the objective function value, gij must be zero. When xij = 0, constraint (31) is satisfied for gij P hij. However, due to the minimizing objective function (30), it must be true that gij = hij. This justifies the correctness of the linearization. Next, we argue that forcing constraints (31) to equality does not cause a loss of optimality. To see this, observe that whenever hij = 0, the right side of (31) is either 0 or 1 so that (32) forces gij to be non-negative. Taking gij = 0 will not cause loss of optimality since the coefficient uij of gij in the objective function is positive. Furthermore, taking xij = 0 in this case will also maintain optimality while maintaining feasibility. Hence, (31) can be taken as equality whenever hij = 0. In the remaining case, hij > 0 and xij is either 0 or 1. If xij = 0, then (31) reduces to gij P hij. In this case, taking gij = hij gives an objective value which is at least as good as taking gij > hij. If xij = 1, then the right side of (31) is either zero or negative. In this case, gij is bounded below by zero and optimality is achieved by taking it to be zero. If the right hand side of (31) is negative in this case (i.e., if hij < 1), then we may increase hij to 1 to make the right side of (31) equal to zero, thereby achieving equality in (31). 2.2.4. Model MTNIP-E: Final version of the exact formulation for MTNIP z ¼ min a;g;x s:t X uij gij ; ð33Þ ði;jÞ2A  aik þ ajk þ gij þ xij P 0;  ajk þ aik þ gij þ xij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i; j 2 N  N 0 ; k ¼ 1; . . . ; K; ði; jÞ 2 A; i; j 2 N  N ; 0 N0k ; ð36Þ 0 0 N0k ; ð37Þ 0 ð38Þ ajk þ gij þ xij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; j 2 N  N ; i 2 N   aik þ gij þ xij P 1; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 N0k ; k ¼ 1; . . . ; K; ði; jÞ 2 A; j 2 N0k ; gij þ xij P 1; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 X rij xij 6 R; N0k ; 0 ð35Þ 0 aik þ gij þ xij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; i 2 N  N ; j 2 N   ajk þ gij þ xij P 0; ð34Þ 0 j2N  j2NN; 0 i2NN; N0k and j 2 ð39Þ N0k ; 0 i2N  N0k ; ð40Þ ð41Þ ði;jÞ2A X ij 2 f0; 1g ði; jÞ 2 A; ð42Þ aik ; free k ¼ 1; . . . ; K; i 2 N  N0 ; gij P 0; ði; jÞ 2 A; ð43Þ ð44Þ zMTNIPE gives the maximum flow achievable by the network user after interdiction. When R = 0, MTNIP-E resembles D-MXF (or when nij + xij is replaced by bij). In this case, zMTNIPE ¼ zDMXF ¼ zMXF . As R is increased, the value of zMTNIPE decreases. The decrease in the objective function value is determined depending on which paths between node groups are disconnected by the interdiction of arcs and the flow values on the disconnected paths. The term gij + xij and aik in MTNIP-E can be interpreted as hij and aik in D-MXF (x), respectively, to gain more insight about the solutions to MTNIP-E and its dual, post-interdiction MXF. This is a result of the fact that gij + xij 6 1 in an optimal solution. Because R is limited and the objective value will unnecessarily increase when gij > 0, there is no incentive in setting xij = 1 and gij = 0, simultaneously. Moreover, because 0 6 hij 6 1 and hij is replaced by gij + xij to obtain MTNIP-E, gij + xij 6 1 follows. Thus, using complementary slackness conditions P of optimality for MXF and D-MXF, gij > 0 implies that k¼1;...;K ðyijk þ yjik Þ ¼ uij and yjik > 0 implies that aik þ ajk þ gij þ xij ¼ dijk with dijk being the right-hand-side values of constraints (34)–(40). If a constraint in (34)–(40) corresponding to bijk is not satisfied at equality, yijk ¼ 0. When xij = 1, gij > 0 in MXF. That is, the set of arcs to interdict is chosen from among the set of saturated arcs in MXF. An optimal solution to MTNIP-E gives an optimal interdiction plan x⁄ to the attacker for a specific scenario. By analyzing multiple scenarios with different values/sets of R; N 0k , or other parameters, an interdictor can develop an attack plan. From the point of the network user, x⁄ can be regarded as the smallest set of arcs to be hardened. Thus, the network user can also develop a robust defense plan or system by going through several what-if analyses. _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. 246 The models developed assume that all special node groups can be both source and sink. However, the formulations can easily be adapted to situations in which some special nodes are only source or only sink. Models can easily be extended to handle issues such as interdicting nodes, disallowing interdiction of certain arcs, allowing partial arc interdiction, and using different types of interdiction resources. MTNIP-E is an exact model for MTNIP because it explicitly minimizes the maximum amount of flow among node groups. Although an exact solution is highly desirable, computational studies show that MTNIP-E is difficult to solve. This leads us to develop a new, easy-tosolve approximating model, which is given next. 3. Approximate formulation of MTNIP Multi-partition network-interdiction model (MPNIM) is a binary-integer program. It does not minimize the maximum flow among N 0k ; . . . ; N 0K directly. Instead, it partitions N into K disjoint subsets N1, . . . , NK with N 01 # N 1 ; . . . ; N 0K # N K and interdicts certain arcs connecting the subsets Nk, while observing constraints on interdiction resources. The objective is to minimize the total capacity of the non-interdicted arcs crossing between Nk. Three decision variables are used in MPNIM: (1) xij that takes on the value of 1 if arc (i, j) crosses between two different subsets and is interdicted; 0 otherwise, (2) xij that takes on the value of 1 if node i is assigned to Nk; 0 otherwise, and (3) kij that takes on the value of 1 if arc (i, j) crosses between two different subsets and is not interdicted; 0 otherwise. 3.1. Model MPNIM: Multi-partition network-interdiction model z ¼ min x;k;x s:t: X uij kij ; ð45Þ xik ¼ 1; i 2 N; ð46Þ ði;jÞ2A X k xik  xjk þ kij þ xij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A;  xik þ xjk þ kij þ xij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A; xij þ kij 6 1; ði; jÞ 2 A; X r ij xij 6 R; ð47Þ ð48Þ ð49Þ ð50Þ ði;jÞ2A xik  1; k ¼ 1; . . . ; K; i 2 N0k ; xik0  0; k0 ¼ 1; . . . ; K; i 2 N0k ; k – k0 ; ð51Þ X ij 2 f0; 1g; ði; jÞ 2 A; ð52Þ ð53Þ kij 2 f0; 1g; ði; jÞ 2 A; ð54Þ xik 2 f0; 1g; k ¼ 1; . . . ; K; i 2 N; ð55Þ The objective (45) minimizes the sum of the capacities on non-interdicted arcs crossing between different node subsets. Constraints (46) require each node i to belong to exactly one subset Nk. Constraints (47) and (48) enforce a partitioning of the nodes and determine whether an arc crosses between two subsets: (1) If i, j 2 Nk, then xik  xjk = 0 and xjk  xik = 0, which allows xij = 0 and kij = 0. xij = 1 and/or kij = 1 are also feasible to constraints (47) and (48) in this case, but we may assume that both are 0 because: (a) kij = 0 contributes less to the objective function than does kij = 1, and (b) xij = 0 consumes less resource than does xij = 1. (Alternate optimal solutions with are xij = 1 possible if excess resource exists.) 0 (2) If i 2 Nk and j 2 N 0k0 ; k–k then xij + kij = 1 is required to maintain feasibility. So, either xij = 1, indicating that arc (i, j) is interdicted or kij = 1, indicating that this arc is not interdicted and contributes to the inter-subset capacity after interdiction. MPNIM classifies the arcs in the network into three groups: (a) Arcs that cross between subsets and are interdicted, (b) arcs that cross between subsets and are not interdicted, and (c) arcs that do not cross between subsets. Constraints (49) together with (47) and (48) ensure that each arc is in one of these three groups. Note that constraints (49) are actually implied by the structure of the model. However, the constraints are added explicitly to prevent from violations that may occur in the case of excess interdiction resource without changing the objective function value. Constraint (50) limits the usage of the interdiction resource as before. Constraints (51) set xik = 1 if node i is preassigned to node subset Nk, i.e., i 2 N 0k , and constraints (52) set xik0 ¼ 0 if i 2 N 0k and k – k0 . Constraints (53) and (55) are set restrictions on the decision variables. Proposition 1. MPNIM solves MMCP when R = 0. Proof. Solving MPNIM by setting R = 0 is clearly equivalent to solving the following model. h 3.2. MPNIMC: Multi-way cut model using MPNIM In addition to (45), (46), (51), (52), (54), and (55): xik  xjk þ kij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A;  xik þ xjk þ kij P 0; k ¼ 1; . . . ; K; ði; jÞ 2 A: ð56Þ ð57Þ _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. 247 Table 1 Model statistics and run times for MTNIM-E and MPNIM on GN. ‘‘[ ]’’ indicates that the model could not be solved in 86400 seconds and that shows the best solution found and the integrality gap at termination. 751 471 361 1 281 1 37 1 1 191 461 2 261 243 1 47 3 161 21 3 3 1 363 3 191 18 192 36 471 192 212 443 183 8 3 26 573 2 2 22 47 343 133 2 81 1 21 3 192 243 391 1 172 192 301 2 22 39 2 241 22 2 412 581 2 383 16 392 361 1 192 361 281 361 183 (a) Max-flow solution before interdiction with z * = 376 (R=0). Subset N 2 1 Subset N1 X 2 X X 1 161 X X X 3 X 2 2 1 16 Subset N3 X 3 3 X (b) Max-flow solution after interdiction with z * = 16 (R=9). Fig. 1. Maximum flow solutions for Pr.1 of GN instances before and after interdiction. The shaded nodes with numbers are special nodes while all other nodes are regular nodes. Special nodes with the same number belong to the same node group. The undirected arcs do not carry any flow while directed arcs show the direction of the flow. The values on the directed arcs represent the flow values with superscripts showing the origin of the flow. The arcs with an X are interdicted arcs, which are the same for both MPNIM and MTNIP-E. The nodes enclosed within a rectangle with dashed lines belong to the same node subset imposed by MPNIM. In MPNIMC, the partitioning of the node set is similar to the one in MPNIM. If i, j2Nk, then kij = 0 can be assumed due to the objective 0 function. If i2Nk and j 2 N k0 ; k–k , then kij = 1 for feasibility. The set of arcs with kij = 1 constitutes the minimum multi-way cut due to the minimizing objective function. 248 _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. When there is a solution for either MPNIMC or integral D-MXF, a solution for the other can be obtained. It can be checked that setting bij = 1 in D-MXF for kij = 1 in MPNIMC and setting aik = 1 in D-MXF for xik = 1 with i2N  N0 in MPNIMC gives a feasible solution for the integral D-MXF. Similarly, a solution for MPNIMC can be obtained from a solution of D-MXF. Proposition 2. zMPNIM provides an upper bound on zMTNIPE . Proof. Let Z IDMXF be the optimal objective function value to integral D-MXF. Proposition 1 implies that zMPNIMC ¼ zIDMXF P zDMXF . Moreover, zDMXF P zMTNIPE ; zMPNIMC P zMPNIM , and zMPNIMC P zMTNIPE can be established. MTNIP-E and MPNIM are obtained from D-MXF and P MPNIMC, respectively, in a similar manner and by adding the same set of interdiction constraints X ¼ fx 2 f0; 1gjAj : ði;jÞ2A r ij xij 6 Rg. Specifically, nij + xij in MTNIP-E replaces bij in D-MXF and kij + xij in MPNIM replaces kij in MPNIMC where gij + xij 6 1 and kij + xij 6 1. Thus, the feasible regions of MTNIP-E and MPNIM are the union of the feasible regions of D-MXF and MPNIMC, respectively, with the interdiction set X. It follows that zMPNIM P zMTNIPE . h Computational studies show that zMPNIM ¼ zMTNIPE for some test problems. However, there are many instances for which zMPNIM > zMTNIPE . The results show that, for optimally solved problems, the difference between zMPNIM and zMTNIPE can be as much as 46.2%. However, notice that zMTNIPE gives the post-interdiction flow capacity through the network while zMPNIM gives the post-interdiction multi-cut capacity (due to Proposition 1). We remark that all of the capacities of non-interdicted arcs constituting zMPNIM do not necessarily contribute to the flow capacity. The post-interdiction flow capacity for MPNIM, zMPNIMF , may be less than zMPNIM . Let X MPNIM represent the optimal set of interdicted arcs in MPNIM. Then, zMPNIMF can be evaluated by solving either: (1) MXF after removing X MPNIM from the network, or (2) MPNIM-F obtained by setting xij = 1 for ði; jÞ 2 X MPNIM in MTNIP-E. In our study, we prefer the latter. Proposition 3. zMPNIM P zMPNIMF P zMTNIPE . Proof. It is clear that zMPNIMF P zMTNIPE . zMPNIM P zMPNIMF must also be true because otherwise Proposition 2 is contradicted. h 4. Computational studies We test MTNIP-E and MPNIM using three different types of networks, grid networks (GN), Euclidean-distance networks (EN), and random networks (RN), each with four different sets of data. Table 2 Model statistics and run times for MTNIM-E and MPNIM on RN. ‘‘[ ]’’ indicates that the model could not be solved in 86400 seconds and that shows the best solution found and the integrality gap at termination. _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. 249 Table 3 Model statistics and run times for MTNIM-E and MPNIM on EN. ‘‘[ ]’’ indicates that the model could not be solved in 86400 seconds and that shows the best solution found and the integrality gap at termination. GN are n1  n2 networks where n1 and n2 are the number of nodes in the horizontal and vertical axes, respectively. The number of nodes and arcs change from 28 to 126 and from 63 to 333, respectively. The arc capacities are randomly drawn from the discrete uniform distribution on [13, 99]. EN are complete graphs where the locations of nodes are generated randomly with an uniform distribution in a square on [1, 100]. The arc capacities are assigned as the Euclidean distances between the nodes set to integer units. RN are also complete graphs where the arc capacities are randomly drawn from the discrete uniform distribution on [1, 100]. The number of nodes and arcs in EN and RN change from 30 to 120 and from 435 to 7140, respectively. It is assumed that rij = 1 for all arcs (i, j) in all three types of networks. Different values of R and K are used for all networks. Computational tests are performed on a PC with 3.0 GHz Intel Core 2 Duo processor and 3 GB of RAM by using the solver CPLEX 9.0. The models are run until the optimality is attained or for 24 hour (86400 seconds) at maximum by using default settings of CPLEX, e.g., moving the best bound strategy for branching is used, cuts are allowed (ILOG, 2003). In the tables, run times and z⁄ are given for problems solved to optimality. For problems not solved to optimality, the resulting objective value and the integrality gap jBP  BFj/(1010 + jBPj), where BP is the objective value of the best integer solution and BF is the best remaining objective value of any unexplored node (ILOG, 2003), are given.
The objective values of MTNIP-E and MPNIM are compared by using the statistic Dz ¼ 100%  zMPNIM =zMTNIPE for problems solved to optimality. For problems not solved optimally, Dz and Dz obtained by replacing zMTNIPE in Dz with upper bound zMTNIPE and lower bound zMTNIP-E reached at the end of
allotted time, respectively, are used. To compare the objective values of MPNIM-F and MTNIP-E, DFz ¼ 100%  zMPNIMF =zMTNIPE is used. Table 1 gives results for test problems on GN. MPNIM can optimally solve all of 16 test problems with solution times ranging from 0 to 17.53 seconds. MTNIP-E can optimally solve 13 problems with solution times changing from 0 to 81680 seconds. The remaining three problems not solved by MTNIP-E are solved by MPNIM with the worst solution time being 17.53 seconds. Dz ¼ 100% for 10 of the 13 problems solved optimally by both models, i.e., the objective function values of the models are the same. Dz for the remaining three problems are 103.6%, 102.5%, and 110.4%, respectively. Dz and Dz for problems not solved optimally by MTNIP-E change from 121% to 150% with an average of 136% and from 100% to 129% with an average of 115.67%, respectively. To give a pictorial view of the effect of the interdictor with respect to flow capacity through a network, the maximum flow solutions before and after interdiction for Pr.1 of the GN instances are given in Fig. 1. Fig. 1(a) shows the flow values on the arcs together with the directions and origins when there is no interdictor, i.e., R = 0. The maximum flow value achieved in this case is 376. Fig. 1(b) indicates the interdicted arcs and flow values for R = 9. For this instance, the set of interdicted arcs determined by MPNIM and MTNIP-E are the same. After interdiction, there remain only two flow paths for flow to occur, both of which achieve the maximum flow of 16. Note that one additional unit of resource is needed to cut off the remaining two flow paths. Fig. 1(b) also shows the partitioning of the node set into subsets resulting from the solution of MPNIM. Table 2 gives results for test problems on RN. MPNIM can optimally solve all of 32 test problems with solution times changing from ranging from 0.03 to 0.99 seconds. MTNIP-E can optimally solve only 20 problems with solution times changing from 0.09 to 17605 250 _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. Table 4 Comparison of the objective function values of MPNIM-F and MTNIP-E for problems solved optimally by both MPNIM and MTNIP-E. seconds. 12 problems not solved by MTNIP-E are solved by MPNIM with solution times changing from 0.05 to 0.99 seconds. Dz ¼ 100% for only 1 problem out of 20 solved optimally by both models. Dz for the remaining 19 problems change from 106.69% to 142.17% with an average of 129.86%. Dz and Dz for problems not solved optimally by MTNIP-E change from 113.86% to 147.5% with an average of 130.8% and from 111.91% to 135.75% with an average of 126.74%, respectively. Table 3 gives results for test problems on EN. MPNIM can optimally solve all of 32 test problems with solution times ranging from 0.02 to 3.86 seconds. MTNIP-E can optimally solve only 22 of the problems with solution times changing from 0.22 to 48720.70 seconds. 10 problems not solved by MTNIP-E are solved by MPNIM with the solution times ranging from 0.13 to 3.77 seconds. Dz ¼ 100% for only 2 problems out of 22 solved optimally by both models. Dz for the remaining 20 problems change from 101.41% to 146.19% with an average of 125.93%. Dz and Dz for problems not solved optimally by MTNIP-E change from 108.9% to 146.88% with an average of 127.09% and from 105.64% to 145.97% with an average of 125.45%, respectively. The previous results show that MPNIM is incomparably better than MTNIP-E with respect to solution times and that zMPNIM P zMTNIPE in compliance with Proposition 2. To summarize, Dz ¼ 100% for 13 problems out of 56 optimally-solved problems, for which the average Dz is 120.33%. For the remaining 43 problems, the worst and the average Dz are 146.19% and 126.33%, respectively. The best Dz values are obtained for problems on GN with an average of 101.27%. Average Dz values for RN and EN are 128.44% and 123.57%, respectively. Table 4 gives DFz for test problems solved optimally by both MPNIM and MTNIP-E. For GN, DFz change from 100% to 106.46% with an average of 100.1%. For RN, DFz change from 100% to 106.38% with an average of 102.3%. For EN, DFz range from 100% to 107.9% with an average of 101.97%. Note that the largest DFz is 107.9% while the largest Dz is 146.19%. Notice that DFz is 101.87% for the problem (Pr. Id. 31 on EN) with the worst Dz value of 146.19%. The worst DFz value of 107.9% is obtained for the problem (Pr. Id. 8 on EN) with Dz value of 117.13%. The results show that zMPNIM P zMPNIMF P zMTNIPE in accordance with Proposition 3. Moreover, DFz is significantly smaller than Dz implying that the solution provided by MPNIM may be an adequate approximation to the solution to MTNIP-E in terms of post-interdiction flow capacity. This combined with the fact that MPNIM is incomparably easier to solve shows that MPNIM can be used instead of MTNIP-E especially when there are time constraints. 5. Conclusion This paper defines and studies MTNIP in which a network user attempts to maximize flow in a network among k P 3 pre-specified node groups while an interdictor uses limited resources to interdict network arcs to minimize this maximum flow. The paper proposes an exact (MTNIP-E) and an approximating model (MPNIM) to solve this NP-hard problem and presents computational results carried out on different types of networks to compare both models. MTNIP-E is obtained by formulating MTNIP as bi-level min–max program and then converting it into a mixed integer program where the flow is explicitly minimized. MPNIM is binary-integer program that does not minimize the flow directly. It partitions the node set into disjoint subsets such that each node group is in a different subset and minimizes the sum of _ Akgün et al. / European Journal of Operational Research 211 (2011) 241–251 I. 251 the arc capacities crossing between different subsets. Computational results show that MPNIM can solve all instances in a few seconds while MTNIP-E cannot solve about one third of the problems in 24 hour. The optimal objective function values of both models are equal to each other for some problems while they differ from each other as much as 46.2% in the worst case. However, when the post-interdiction flow capacity incurred by the solution of MPNIM is computed and compared to the objective value of MTNIP-E, the largest difference is only 7.90%. This result implies that MPNIM may be a very good approximation to MTNIP-E. Further research may be on devising methods to improve the solution times of MTNIP-E, e.g., decomposition techniques outlined in Israeli and Wood (2002) can be tried and integer-programming cuts can be developed, improving MPNIM to better approximate MTNIP-E, and extending the models to allow stochastic interdictions. Acknowledgements The authors are grateful to two anonymous referees for providing constructive feedback that has helped improve in major ways the presentation of the material in the paper. References Ahuja, R.K., Magnanti, T.L., Orlin, J.B., 1993. Network Flows. Prentice Hall, New Jersey. Alevras, D., Grötschel, M., Jonas, P., Paul, U., Wessaly, R., 1998. Survivable mobile phone network architectures: Models and solution methods. IEEE Communications Magazine 36, 88–93. Assimakopoulos, N., 1987. A network interdiction model for hospital infection control. Computers in Biology and Medicine 17 (6), 413–422. Anandalingam, G., Apprey, V., 1991. Multi-level programming and conflict resolution. European Journal of Operational Research 51 (2), 233–247. Brown, G., Carlyle, M., Diehl, D., Kline, J., Wood, R.K., 2005. A two-sided optimization for theater-ballistic missile defense. Operations Research 53 (5), 745–763. Brown, G., Carlyle, M., Salmeron, J., Wood, R.K., 2006. Defending critical infrastructure. Interfaces 36 (6), 530–544. Church, R.L., Scaparra, M.P., Middleton, R.S., 2004. Identifying critical infrastructure: The median and covering facility interdiction problems. Annals of the Association of American Geographers 94, 491–502. Cormican, K.J., Morton, D.P., Wood, R.K., 1998. Stochastic network interdiction. Operations Research 46, 184–197. Costa, M., Letocart, L., Roupin, F., 2005. Minimal multicut and maximal integer multiflow: A survey. European Journal of Operational Research 162, 55–69. Dahlhaus, E., Johnson, D.S., Papadimitriou, C.H., Seymour, P.D., Yannakakis, M., 1994. The complexity of multiway cuts. SIAM Journal of Computers 23, 864–894. Dempe, S., 2002. Foundations of Bilevel Programming. Kluwer Academic Publishers, Dordrecht. Desai, J., Sen, S., 2010. A global optimization algorithm for reliable network design. European Journal of Operational Research 200 (1), 1–8. Ford Jr., L.R., Fulkerson, D.R., 1956. Maximal flow through a network. Canadian Journal of Mathematics. Fulkerson, D.R., Harding, G.C., 1977. Maximizing the minimum source-sink path subject to a budget constraint. Mathematical Programming 13, 116–118. Garg, N., Vazirani, V.V., Yannakakis, M., 1996. Approximate max-flow min-(multi)cut theorems and their applications. SIAM Journal of Computers 25 (2), 235–251. Ghare, P.M., Montgomery, D.C., Turner, W.C., 1971. Optimal interdiction policy for a flow network. Naval Research Logistics Quarterly 18, 37–45. Golden, B., 1978. A problem in network interdiction. Naval Research Logistics Quarterly (25), 711–713. Helmbold, R.L., 1971. A counter capacity network interdiction model. Report R-611-PR, March 1971. RAND Corporation, Santa Monica, California. ILOG S.A. ILOG/CPLEX 9.0 Reference Manual, 2003. Israeli, E., 1999. System interdiction and defense, Doctoral dissertation, March 1999, Naval Postgraduate School, Monterey, California. Israeli, E., Wood, R.K., 2002. Shortest path network interdiction. Networks 40, 97–111. Lim, C., Smith, J.C., 2007. Algorithms for discrete and continuous multicommodity flow network interdiction problems. IIE Transactions 39 (1), 15–26. Lubore, S.H., Ratliff, H.D., Sicilia, G.T., 1971. Determining the most vital link in a flow network. Naval Research Logistics Quarterly 17, 261–268. McMasters, A.W., Mustin, T.M., 1970. Optimal interdiction of a supply network. Naval Research Logistics Quarterly 17, 261–268. Myung, Y.S., Kim, H.J., Tcha, D.W., 1999. Design of communication networks with survivability constraints. Management Science 45, 238–252. Ouveysi, I., Wirth, A., 1999. On design of a survivable network architecture for dynamic routing: Optimal solution strategy and efficient heuristic. European Journal of Operational Research 117, 30–44. Salmeron, J., Wood, R.K., Baldick, R., 2004. Analysis of electric grid security under terrorist threat. IEEE Transactions on Power Systems 19, 905–912. Scaparra, M.P., Church, R.L., 2008. A bilevel mixed integer program for critical infrastructure protection planning. Computers and Operations Research 35, 1905–1923. Siman, M., Cruz, J.B., 1973. On the Stackelberg strategy in nonzero-sum games. Journal of Optimization Theory and Applications 11, 533–555. Smith, J.C., Lim, C., Sudargho, F., 2007. Survivable network design under optimal and heuristic interdiction scenarios. Journal of Global Optimization 38, 181–199. Washburn, A., Wood, R.K., 1994. Two-person zero-sum games for network interdiction. Operations Research 42, 243–251. Whiteman, P.S., 1999. Improving single strike effectiveness for network interdiction. Military Operations Research 4, 15–30. Wollmer, R.D., 1964. Removing arcs from a network. Operations Research 12, 807–1076. Wollmer, R.D., 1970a. Algorithm for targeting strikes in a lines-of-communication network. Operations Research 18, 497–515. Wollmer, R.D., 1970b. Interception in a network. Naval Research Logistics Quarterly 17, 207–216. Wood, R.K., 1993. Deterministic network interdiction. Mathematical and Computer Modeling 17, 1–18.