---
abstract: 'We investigate methods to solve the maximum-reliability stochastic network interdiction problem (SNIP).  In this problem, a defender interdicts arcs on a directed graph to minimize an attacker’s probability of undetected traversal through the network.  The attacker’s origin and destination are unknown to the defender and assumed to be random.  SNIP can be formulated as a stochastic mixed-integer program via a deterministic equivalent formulation (DEF).  As the size of this DEF makes it impractical for solving large instances, current approaches to solving SNIP rely on modifications of Benders decomposition.  We present two new approaches to solve SNIP.  First, we introduce a new DEF that is significantly more compact than the standard DEF.  Second, we propose a new path-based formulation of SNIP.  The number of constraints required to define this formulation grows exponentially with the size of the network, but the model can be solved via delayed constraint generation.  We present valid inequalities for this path-based formulation which are dependent on the structure of the interdicted arc probabilities.  We propose a branch-and-cut (BC) algorithm to solve this new SNIP formulation.  Computational results demonstrate that directly solving the more compact SNIP formulation and this BC algorithm both provide an improvement over a state-of-the-art implementation of Benders decomposition for this problem.' 
author:
- 'Eli Towle[^1]'
- 'James Luedtke[^2]'
bibliography:
- 'refs.bib'
date: 'Department of Industrial and Systems Engineering, University of Wisconsin – Madison'
title: 'New Solution Approaches for the Maximum-Reliability Stochastic Network Interdiction Problem'
---

This work was supported by National Science Foundation grants CMMI-1130266 and SES-1422768. 

Introduction
============

Network interdiction problems feature a defender and an attacker with opposing goals.  The defender first modifies a network using a finite budget in an attempt to diminish the attacker’s ability to perform a task.  These modifications could include increasing arc traversal costs, reducing arc capacities, or removing arcs altogether.  The attacker then optimizes his or her objective with respect to the newly modified network.  Network interdiction problems have been studied in the context of nuclear weapons interdiction [@morton2007; @pan2003], disruption of nuclear weapons projects [@brown2009], interdiction of drug smuggling routes [@wood1993], and general critical infrastructure defense [@brown2009-2]. 

In the maximum-reliability stochastic network interdiction problem (SNIP), formulated by Pan [et al.]{}[@pan2003], an attacker seeks to maximize the probability of avoiding detection while traveling through a directed network from an origin to a destination.  There is a chance the attacker will be detected when traversing each arc.  A defender may expend resources to place sensors on certain arcs, thereby increasing the probability of detecting an attacker on those arcs.  The origin and destination of the attacker are unknown to the defender, and are assumed to be random.  The defender’s goal is to minimize the expected value of the attacker’s maximum-reliability path through the network over all origin–destination scenarios.  This can be interpreted as minimizing the overall probability of the attacker successfully attacking a critical infrastructure target or smuggling contraband to a target location undetected. 

Morton [et al.]{}[@morton2007] present a deterministic equivalent formulation (DEF) of SNIP derived from the dual of the attacker’s optimization problem.  Pan and Morton [@pan2008] improve the DEF by incorporating the values of uninterdicted maximum-reliability paths into the model constraints.  They also derive to strengthen the linear programming (LP) relaxation of the Benders master problem before initiating the Benders algorithm.  Bodur [et al.]{}[@bodur2016] investigate the strength of integrality-based cuts in conjunction with Benders cuts for stochastic integer programs with continuous recourse, and use test instances of the SNIP problem to demonstrate the value of integrality-based cuts within Benders decomposition. 

We propose two new approaches to solve SNIP.  First, we present a more compact DEF.  This formulation combines constraints of the DEF for scenarios ending at the same destination, significantly reducing the number of constraints and variables in the DEF.  Second, we present a new formulation of SNIP that includes constraints for every origin–destination path.  Although the formulation size grows exponentially with the number of arcs in the network, the model can be solved with delayed constrained generation.  A path’s reliability function is a supermodular set function representable with a strictly convex and decreasing function.  We propose a branch-and-cut algorithm that exploits this structure for each path through the network to derive valid inequalities.  We consider three cases for the probabilities of traversing interdicted arcs.  Two of these cases use inequalities derived by Ahmed and Atamt[ü]{}rk [@ahmed2011] for the problem of maximizing a submodular utility function having similar structure. 

The maximum-reliability network interdiction problem is closely related to the shortest-path variant, in which the network defender attempts to maximize the length of the attacker’s shortest path from origin to destination.  When the probabilities of traversing interdicted arcs in a maximum-reliability problem are all strictly positive , the problem is equivalent, via a logarithmic transformation, to a shortest-path network interdiction problem [@ahuja1993; @morton2010].  This transformation, however, is not valid when there exists an arc that reduces the attacker’s probability of successful traversal to zero when it is interdicted by the defender.  A deterministic version of the shortest-path network interdiction problem was initially explored by Fulkerson and Harding [@fulkerson1977], and later by Golden [@golden1978]. 

An alternative interdiction problem mostly unrelated to our work is the maximum-flow network interdiction problem, introduced by Wollmer [@wollmer1964].  In this problem, a defender changes arc capacities in order to minimize the attacker’s maximum flow.  Deterministic and stochastic versions of this problem have been studied by many authors, e.g., [@cormican1998; @hemmecke2003; @jeff2008; @wood1993]. 

In Section \[sec:existingresults\], we review relevant existing approaches to solving SNIP.  We propose a compact DEF of SNIP in Section \[sec:compact\].  In Section \[sec:pathbased\], we propose a path-based formulation for SNIP and describe valid inequalities for the relevant mixed-integer set.  We describe our computational experiments and results in Section \[sec:computational\]. 

Problem statement and existing results {#sec:existingresults}
======================================

Let $N$ and $A$ denote the set of nodes and the set of arcs of a directed network.  Let $D \subseteq A$ denote the set of arcs available for interdiction.  In the first stage of SNIP, the defender may choose to install sensors on a subset of the interdictable arcs.  The cost to install a sensor on arc $a \in D$ is ${c}_a > 0$. The defender is constrained by installation budget $b > 0$. The probability of the attacker’s undetected traversal through arc $a \in A$ without a sensor installed is $r_a \in (0,1]$. When a sensor is installed on arc $a \in D$, this probability reduces to $q_a \in [0, r_a)$. Let ${\Omega}$ be the finite set of origin–destination pairs, where ${\omega}= (s, t) \in {\Omega}$ occurs with probability $p_{\omega}> 0$ ($\sum_{{\omega}\in {\Omega}} p_{\omega}= 1$).  We assume a path exists from $s$ to $t$ for all $(s, t) \in {\Omega}$. If not, the defender can discard that scenario from consideration.  In the second stage of SNIP, the attacker’s origin and destination are realized.  The attacker traverses the network from the origin to the corresponding destination via the maximum-reliability path given the defender’s set of interdicted arcs. 

For $s,t \in N$, a simple $s$-$t$ path $P$ is a set of arcs, say $\added{\{}(i_0,i_1),(i_1,i_2),\ldots,$ $\replaced{(i_{|P|-1},i_{|P|})}{(i_{k-1},i_{k})}\added{\}}$, where $i_0 = s$, $\replaced{i_{|P|}}{i_{k}} = t$, and the nodes $i_0,i_1,i_2,\ldots,i_{|P|} \in N$ are all distinct.  Let ${\mathcal{P}}_{st}$ be the set of all simple paths from $s \in N$ to $t \in N$. Let $S \subseteq D$ be the set of interdicted arcs, as chosen by the defender.  The function $$\begin{aligned}
    h_P(S) \coloneqq \left( \prod_{a \in P} r_a \right) \left( \prod_{a \in P \cap S} \frac{q_a}{r_a} \right)\end{aligned}$$ calculates the probability of the attacker traversing the set of arcs $P$ undetected given the interdicted arcs $S$.

In the stochastic network interdiction problem (SNIP), the defender selects to interdict to minimize the expected value of the attacker’s maximum-reliability path, subject to the budget restriction on the selected arcs.  SNIP can be formulated as $$\begin{aligned}
{3}
    \label{eq:path}
    \begin{aligned}
        \min_{S} && \ \smashoperator{\sum\limits_{{\omega}= (s, t) \in {\Omega}}} p_{{\omega}} & \max\{h_P(S) \colon P \in {\mathcal{P}}_{st}\} && \\
        \textrm{s.t.} &&  \sum\limits_{a \in S} {c}_a & \leq b. &&
    \end{aligned}\end{aligned}$$

Deterministic equivalent formulation {#subsec:extensive}
------------------------------------

A DEF of SNIP can be obtained by using first-stage binary variables $x_{a}$ to represent whether the defender installs a sensor on arc $a \in D$. That is, $x_a = 1$ if the defender elects to install a sensor on arc $a \in D$, and $x_a = 0$ if no sensor is installed on that arc.  Second-stage continuous variables $\pi_i^{\omega}$ represent the maximum probability of the attacker traveling undetected from $i \in N$ to $t$ in scenario ${\omega}= (s, t) \in {\Omega}$.

The formulation uses second-stage constraints to calculate the $\pi$ variables for each scenario ${\omega}\in {\Omega}$. The $\pi$ variables are calculated using an LP formulation of the dynamic programming (DP) optimality conditions for calculating a maximum-reliability path.  Using the convention $q_a \coloneqq r_a$ and $x_a \coloneqq 0$ for all $a \in A \setminus D$, each $\pi_i^{\omega}$ is calculated as $$\begin{aligned}
\label{eq:dp}
    \pi_i^{\omega}& = \max_{a = (i,j) \in A} \left\lbrace \pi^{\omega}_j\left[r_a + (q_a - r_a)x_a \right] \right\rbrace, \quad i \in N,\ {\omega}\in {\Omega}.\end{aligned}$$ The maximum probability of the attacker reaching $t$ undetected from node $i \in N$ is the maximum over all forward adjacent nodes $j \in N$ of $\replaced{\pi_j^{{\omega}}}{\pi_j^{t}}$, multiplied by $r_a$ if the arc is not interdicted, or $q_a$ if it is interdicted. 

The optimality conditions imply the set of inequalities $$\begin{aligned}
\label{eq:dpineqs}
    \pi_i^{\omega}& \geq \pi^{\omega}_j\left[r_a + (q_a - r_a)x_a \right], \quad a = (i,j) \in A,\ {\omega}\in {\Omega}.\end{aligned}$$ The inequalities are nonlinear.  Pan and Morton [@pan2008] derive the following DEF of SNIP that uses a linear reformulation of these inequalities:

\[eq:ext\] $$\begin{aligned}
{3}
        \min_{x, \pi}\enspace\ && \quad \smashoperator{\sum\limits_{{\omega}= (s, t) \in {\Omega}}} p_{\omega}\pi_{s}^{\omega}&&& \label{eq:objfun} \\
        \textrm{s.t.}\enspace\ && \smashoperator{\sum\limits_{a \in D}} {c}_a x_a & \leq b && \label{eq:extcon1} \\
        && \pi_i^{\omega}- r_a\pi_j^{\omega}& \geq 0, && a=(i,j) \in A \setminus D,\ {\omega}\in {\Omega}\label{eq:pi1} \\
        && \pi_i^{\omega}- r_a\pi_j^{\omega}& \geq -(r_a - q_a) {u}_j^{\omega}x_a, \qquad && a=(i,j) \in D,\ {\omega}\in {\Omega}\label{eq:pi2} \\
        && \pi_i^{\omega}- q_a\pi_j^{\omega}& \geq 0, && a=(i,j) \in D,\ {\omega}\in {\Omega}\label{eq:pi3} \\
        && \pi_t^{\omega}& = 1, && {\omega}= (s, t) \in {\Omega}\label{eq:pi4} \\
        && x_a & \in \{0,1\}, && a \in D. \label{eq:extcon2}
    \end{aligned}$$

The objective function minimizes the expected value of the attacker’s maximum-reliability path.  For each scenario ${\omega}= (s, t) \in {\Omega}$, the parameter ${u}_j^{{\omega}}$ represents the value of the maximum-reliability path from $j \in N$ to $t$ when no sensors are installed, and hence is an upper bound on $\pi_j^{{\omega}}$. These parameters are calculated in a model preprocessing step.  The linear constraints – formulate the nonlinear DP inequalities .  Constraints enforce $\pi_i^{\omega}\geq r_a
\pi_j^{\omega}$ for all $a = (i, j) \in A \setminus D,\ {\omega}\in {\Omega}$. For ${\omega}\in {\Omega}$ and $a = (i, j) \in D$, if $x_a =
0$, then becomes $\pi_i^{\omega}\geq r_a \pi_j^{\omega}$, which dominates .  On the other hand, if $x_a
= 1$, then implies $\pi_i^{\omega}- r_a \pi_j^{\omega}\geq -(r_a - q_a) {u}_j^{\omega}$, dominates .  Thus, in either case – are equivalent to .  Since the variables $\pi_s^{{\omega}}$, ${\omega}=(s,t) \in {\Omega}$, have positive coefficients in the objective, the equations will be satisfied optimal solution. 

Benders decomposition {#subsec:benders}
---------------------

Directly solving the DEF with a mixed-integer programming solver may be too time-consuming due to its large size.  Benders decomposition can be used to decompose large problems like SNIP.  After introducing the SNIP formulation, Pan and Morton [@pan2008] outline a Benders decomposition algorithm for the SNIP DEF .  For a fixed vector $x \in [0, 1]^{\added{|D|}}$ satisfying $\sum_{a \in D} c_a x_a \leq b$, the $\pi_j^{\omega}$ variables in can be obtained by solving

\[eq:extfixed\] $$\begin{aligned}
{3}
        {E}^{st}(x) \coloneqq \min_{\pi}\enspace && \pi_{s} \qquad\quad &&& \label{eq:fixedobjfun} \\
        \textrm{s.t.}\enspace && \pi_i - r_a\pi_j & \geq 0, && a=(i,j) \in A \setminus D \label{eq:fixedpi1} \\
        && \pi_i - r_a\pi_j & \geq -(r_a - q_a) {u}_j^{\omega}x_a, \qquad && a=(i,j) \in D \label{eq:fixedpi2} \\
        && \pi_i - q_a\pi_j & \geq 0, && a=(i,j) \in D \label{eq:fixedpi3} \\
        && \pi_t & = 1 && \label{eq:fixedpi4}
    \end{aligned}$$

for each scenario ${\omega}= (s, t) \in {\Omega}$. The dual of is $$\begin{aligned}
{3}
    \label{eq:dualsubproblem}
    \begin{aligned}
        \max_{y, z}\ && y_t - \smashoperator{\sum\limits_{a = (i, j) \in D}} (r_a - q_a) {u}_j^{\omega}y_{ij} x_a \qquad \qquad\ &&& \\
        \textrm{s.t.}\ && \smashoperator{\sum\limits_{(s, j) \in A}} ( y_{s j} + z_{s j} ) & = 1 && \\
        && \smashoperator{\sum\limits_{(i, j) \in A}} (y_{ij} + z_{ij}) - \smashoperator{\sum\limits_{a=(j, i) \in A}} (r_a y_{ji} + q_a z_{ji}) & = 0,\quad && i \in N \setminus \{s, t\} \\
        && y_{t} - \smashoperator{\sum\limits_{a=(j, t) \in A}} (r_a y_{j t} + q_a z_{j t}) & = 0 && \\
        && y_{ij}, z_{ij} & \geq 0, && (i, j) \in A \\
        && y_{t} & \geq 0. 
    \end{aligned}\end{aligned}$$ Dual variables $y_{ij}$ correspond to DEF constraints and , $z_{ij}$ is the dual variable for constraints , and $y_{t}$ is the dual variable for constraint .  We fix $z_{ij} \coloneqq 0$ for all $(i, j) \in A \setminus D$, as constraint only applies to arcs $D$.

The Benders master problem is as follows:

\[eq:benders\] $$\begin{aligned}
{3}
        \min_{x, \theta}\enspace && \sum\limits_{{\omega}\in {\Omega}} p_{\omega}\theta^{\omega}&&& \\
        \textrm{s.t.}\enspace && \smashoperator{\sum\limits_{a \in D}} {c}_a x_a & \leq b && \\
        && \theta^{\omega}\geq \bar{y}_{t}^{{\omega}} &- \smashoperator{\sum\limits_{a = (i, j) \in D}} (r_a - q_a) {u}_j^{\omega}\bar{y}_{ij}^{{\omega}}x_a,\quad  && (\bar{y}^{\omega}, \bar{z}^{\omega}) \in K^{{\omega}},\ {\omega}= (s, t) \in {\Omega}\label{eq:optimalitycut} \\
        && \theta^{\omega}& \geq 0, && {\omega}\in {\Omega}\\
        && x_a & \in \{0, 1\}, && a \in D. \label{eq:bendersintegrality}
    \end{aligned}$$

The Benders algorithm begins by solving the master problem with no constraints of the form to obtain a candidate solution $(\bar{x}, \bar{\theta})$. At iteration $k$ of the Benders cutting plane algorithm, a dual subproblem is solved for each scenario using the candidate solution $(\bar{x}, \bar{\theta})$ to find a dual-feasible extreme point $(\bar{y}^{{\omega}}, \bar{z}^{{\omega}})$. If a cut in the form of constraint cuts off the candidate solution $\bar{\theta}^{\omega}$, this cut is added to the Benders master problem by including the point in the set $K^{\omega}$. The updated master problem is solved to generate a new candidate solution $(\bar{x}, \bar{\theta})$. This process is repeated until none of the scenario cuts constructed at an iteration of the algorithm cut off the candidate solution obtained by solving the master problem at the previous iteration. 

Benders decomposition can also be used in a branch-and-cut algorithm.  The integrality constraint is relaxed and enforced within the branch-and-bound tree.  At each integer-feasible solution $(\bar{x}, \bar{\theta})$ obtained at a node in the master problem branch-and-bound tree, a dual subproblem is solved for each scenario.  Violated cuts are added to the LP formulation of that node and it is re-solved.  When no violated inequalities are found for an integer-feasible solution, the upper bound may be updated and the node pruned. 

Pan and Morton [@pan2008] enhance the Benders branch-and-cut algorithm by using optimal solutions to scenario subproblems from previous iterations to create step inequalities for the relaxed Benders master problem.  These additional inequalities are added to the formulation as cuts to tighten the LP relaxation of the master problem.  This leads to a smaller branch-and-bound tree, which in turn reduces the time spent solving mixed-integer programs. 

The general-purpose Benders approach tested by Bodur [et al.]{}[@bodur2016] is very effective in solving SNIP.  A key implementation detail of this approach is to first solve the relaxed Benders master problem, obtained by removing the integrality restriction on $x$. Once no more Benders inequalities can be added to the LP relaxation of the master problem, the integrality restriction on $x$ is restored, and the model, including the identified Benders cuts, is passed to a mixed-integer programming solver to begin the Benders branch-and-cut algorithm.  The solver adds its own general-purpose integrality-based cuts to the formulation using the initial set of Benders cuts.  This results in a stronger LP relaxation bound, and reduces the size of the branch-and-bound tree.  Bodur [et al.]{}[@bodur2016] employed a smaller gap tolerance as a stopping criterion ($10^{-3}$) for their Benders branch-and-cut implementation than Pan and Morton [@pan2008] used in their experiments ($10^{-2}$).  We implemented the Benders branch-and-cut algorithm as in [@bodur2016], and tested it with a relative optimality gap tolerance of $10^{-2}$ to compare the results to those published in [@pan2008].  With these matching tolerances, the Benders algorithm ran $3$–$4$ times faster than Pan and Morton’s algorithm.  Although the difference in hardware used in these experiments makes it impossible to directly compare the performance of these methods, we conclude that the implementation of Benders branch-and-cut described in [@bodur2016] is currently among the most efficient ways to solve instances of SNIP, and hence we compare against this approach in our numerical experiments. 

Compact deterministic equivalent formulation {#sec:compact}
============================================

For scenarios sharing a destination node, the DP optimality conditions for SNIP are identical.  Hence, the DEF contains redundancies by repeating the LP formulation of these optimality conditions, constraints –, for each destination node.  Motivated by this observation, we present a compact DEF of SNIP that groups together second-stage constraints for scenarios with a common destination.  Let $T$ be the set of unique destination nodes: $T = \{ t \in N \colon (s,t) \in {\Omega}$ for some $s \in N \}$. Second-stage variables $\pi_j^t$ now represent the probability of the attacker successfully traveling from node $j \in N$ to destination $t \in T$ in any scenario having destination $t$. We obtain the following new DEF for SNIP:

\[eq:com\] $$\begin{aligned}
{3}
        \min_{x, \pi}\ && \quad\ \smashoperator{\sum\limits_{{\omega}= (s, t) \in {\Omega}}} \ p_{\omega}\pi_{s}^{t} &&& \label{eq:comobj} \\
        \textrm{s.t.} &&  \smashoperator{\sum\limits_{a \in D}} {c}_a x_a & \leq b && \label{eq:combudget} \\
        && \pi_i^t - r_a \pi_j^t & \geq 0, && a = (i,j) \in A \setminus D,\ t \in T \label{eq:com1} \\
        && \pi_i^t - r_a \pi_j^t & \geq -(r_a - q_a) {u}_j^t x_a, \qquad && a = (i,j) \in D,\ t \in T \label{eq:com2} \\
        && \pi_i^t - q_a \pi_j^t & \geq 0, && a = (i,j) \in D,\ t \in T \label{eq:com3} \\
        && \pi_t^t & = 1, && t\in T \label{eq:com4} \\
        && x_a &\in \{0,1\}, && a \in D. \label{eq:comcon1}
    \end{aligned}$$

Parameter ${u}_j^t$ is the value of the attacker’s maximum-reliability path from $j \in N$ to $t \in T$ when no sensors are installed.  By definition, ${u}_j^t = {u}_j^{\omega}$ for all ${\omega}= (s, t) \in {\Omega}$ and $j \in N$. The objective function weights each $\pi_s^t$ variable by its respective scenario probability.  The summation of these terms represents the overall probability of a successful attack, which the defender seeks to minimize. 

\[prop:1\]

$$\begin{aligned}
{2}
        \added{\min_{s \in \mathcal{S}(t)} \bar{\pi}_i^{st}} & \added{\geq \min_{s \in \mathcal{S}(t)} r_a \bar{\pi}_j^{st},} && \added{a = (i,j) \in A \setminus D,\ t \in T} \\
        \added{\min_{s \in \mathcal{S}(t)} \bar{\pi}_i^{st}} & \added{\geq \min_{s \in \mathcal{S}(t)} [r_a \bar{\pi}_j^{st} - (r_a - q_a) {u}_j^t \bar{x}_a], \qquad} && \added{a = (i,j) \in D,\ t \in T} \\
        \added{\min_{s \in \mathcal{S}(t)} \bar{\pi}_i^{st}} & \added{\geq \min_{s \in \mathcal{S}(t)} q_a \bar{\pi}_j^{st},} && \added{a = (i,j) \in D,\ t \in T} \\
        \added{\min_{s \in \mathcal{S}(t)} \bar{\pi}_t^{st}} & \added{= 1,} && \added{t \in T.}
    \end{aligned}$$ $$\begin{aligned}
        \added{\min_{s \in \mathcal{S}(t)} [r_a \bar{\pi}_j^{st} - (r_a - q_a) {u}_j^t \bar{x}_a]} &\added{= r_a \min_{s \in \mathcal{S}(t)} \bar{\pi}_j^{st} - (r_a - q_a) {u}_j^t \bar{x}_a,}
    \end{aligned}$$ $(\bar{x}, \hat{\pi})$ $$\begin{aligned}
        \added{\smashoperator{\sum_{{\omega}= (s, t) \in {\Omega}}} \ p_{{\omega}} \hat{\pi}_{s}^{t} = \enspace \smashoperator{ \sum_{{\omega}= (s, t) \in {\Omega}}} \ p_{{\omega}} \min_{s' \in \mathcal{S}(t)} \bar{\pi}_{s}^{s' t} \leq \enspace \smashoperator{ \sum_{{\omega}= (s, t) \in {\Omega}}} \ p_{{\omega}} \bar{\pi}_{s}^{st}.} 
    \end{aligned}$$ Now, let $(\hat{x}, \hat{\pi})$ be a solution to the LP relaxation of .  Let $\bar{\pi}_i^{{\omega}} \coloneqq \hat{\pi}_i^t$ for all $i \in N$ and ${\omega}= (s, t) \in {\Omega}$. $(\replaced{\hat{x}}{\bar{x}}, \bar{\pi})$ is a solution to with objective value $\sum_{{\omega}\added{= (s,t)} \in {\Omega}} p_{{\omega}} \bar{\pi}_{s}^{{\omega}} = \sum_{{\omega}\added{= (s,t)} \in {\Omega}} p_{{\omega}} \hat{\pi}_{s}^{t}$.

A Benders algorithm can be applied to the compact formulation by introducing variables $\theta^t$ to represent the probability-weighted sum of maximum-reliability path values for scenarios ending at node $t \in T$. We experimented with a Benders decomposition on the DEF using the Benders algorithm of Bodur [et al.]{}[@bodur2016].  We found that the Benders algorithm performed much worse on this formulation than on the DEF .  The reason for this poor performance is that the root relaxation obtained after the mixed-integer programming solver added its general-purpose cuts was much weaker when using a Benders reformulation of , as compared to the Benders reformulation .  This is consistent with the results of [@bodur2016], which indicates that the integrality-based cuts derived in the formulation can be stronger than those derived in a “projected” Benders master problem that only uses $\theta^t$ variables. 

Path-based formulation and cuts {#sec:pathbased}
===============================

We now derive a new mixed-integer linear programming formulation of SNIP based on the path-based formulation .  With the introduction of binary variables $x_a$ to denote whether network arc $a \in D$ is interdicted, we map the set function $h_P$ to a vector function $\bar{h}_P$. In particular, for $S \subseteq D$ we let ${\rchi^S}\in \{0, 1\}^{\added{|D|}}$ be the characteristic vector of the set $S$: ${\rchi^S}_a = 1$ if $a \in S$, and $0$ otherwise.  Then we define $\bar{h}_P: \{0,1\}^{\added{|D|}} \rightarrow \mathbb{R}_{\added{+}}$ as $$\begin{aligned}
    \bar{h}_P(x) \coloneqq \left[ \prod_{a \in P} r_a \right] \left[ \prod_{a \in P \cap D} \left( \frac{q_a}{r_a}
    \right)^{x_a} \right]\end{aligned}$$ so that $h_P(S) = \bar{h}_P({\rchi^S})$ for $S \subseteq D$.

Then model can be formulated as $$\begin{aligned}
{3}
    \label{eq:path2}
    \begin{aligned}
        \min_{x,{\pi}} && \smashoperator{\sum\limits_{{\omega}= (s,t) \in {\Omega}}}\ p_{\omega}& {\pi}_s^t && \\
        \textrm{s.t.} &&  \sum\limits_{a \in D} {c}_a x_a & \leq b && \\
        && {\pi}_s^t & \geq \max\{\bar{h}_P(x) \colon P \in {\mathcal{P}}_{st}\}, && (s, t) \in {\Omega}\\
        && x_a & \in \{0, 1\}, && a \in D.
    \end{aligned}\end{aligned}$$

We use the special structure of $h_P$ for each $P \in {\mathcal{P}}_{st}$ to build a linear formulation of model . 

$h \colon 2^{\mathcal{N}}\rightarrow {\mathbb{R}}$ is a *supermodular* set function over a ground set ${\mathcal{N}}$ if $h(S_1 \cup \{{a}\}) - h(S_1) \leq h(S_2 \cup \{{a}\}) - h(S_2)$ for all $S_1 \subseteq S_2 \subseteq {\mathcal{N}}$ and ${a}\in {\mathcal{N}}\setminus S_2$.

$h_P(\cdot)$ is a supermodular set function. 

Let $S_1, S_2 \subseteq D$ with $S_1 \subseteq S_2$. Let ${a'}\in D \setminus S_2$. The result is immediate if ${a'}\notin P$. Therefore, assume ${a'}\in P$ to obtain $$\begin{aligned}
{2}
        h_P(S_1 \cup \{{a'}\}) - h_P(S_1) &= \left( \prod_{a \in P} r_a \right) \left( \prod_{a \in P \cap S_1} \frac{q_a}{r_a} \right) \left(\frac{q_{{a'}}}{r_{{a'}}} - 1 \right) && \\
        &\leq \left( \prod_{a \in P} r_a \right) \left( \prod_{a \in P \cap S_2} \frac{q_a}{r_a} \right) \left(\frac{q_{{a'}}}{r_{{a'}}} - 1 \right) && \\
        &= h_P(S_2 \cup \{{a'}\}) - h_P(S_2). && 
    \end{aligned}$$

The maximum of a set of supermodular functions is not in general a supermodular function, so $\max_{P \in {\mathcal{P}}_{st}}
h_P(S)$ is not necessarily supermodular.  To exploit the supermodular structure for each individual path, we consider an equivalent formulation of that contains an inequality for every scenario $(s, t) \in {\Omega}$ and every path from $s$ to $t$:

\[eq:path3\] $$\begin{aligned}
{3}
        \min_{x,{\pi}} && \quad\enspace\ \smashoperator{\sum\limits_{{\omega}= (s,t) \in {\Omega}}}\ p_{\omega}&{\pi}_s^t && \\
        \textrm{s.t.} &&  \sum\limits_{a \in D} {c}_a x_a & \leq b && \label{eq:path3budget} \\
        && {\pi}_s^t & \geq \bar{h}_P(x), \qquad && P \in {\mathcal{P}}_{st},\ (s, t) \in {\Omega}\label{eq:path3path} \\
        && x_a & \in \{0, 1\}, && a \in D.
    \end{aligned}$$

Each second-stage constraint includes a supermodular function on the right-hand side. 

Consider the mixed-integer feasible region of for a fixed scenario $(s, t) \in {\Omega}$ and path $P \in {\mathcal{P}}_{st}$: $$\begin{aligned}
    {H_P}\coloneqq \left\lbrace (x,{\pi}) \in \{0,1\}^{\added{|D|}} \times {\mathbb{R}}\colon {\pi}\geq \bar{h}_P(x) \right\rbrace \label{eq:FP}.\end{aligned}$$ We are interested in a linear formulation of ${H_P}$. Let ${\rho_P^{a}}(S) \coloneqq h_P(S \cup \{a\}) - h_P(S)$ be the marginal difference function of the set function $h_P$ with respect to arc $a \in D$. Nemhauser [et al.]{}[@nemhauser1978] provide an exponential family of linear inequalities that can be used to define ${H_P}$. Applied to ${H_P}$, these inequalities are $$\begin{aligned}
    {\pi}&\geq h_P(S) - \smashoperator{\sum\limits_{a \in S}} {\rho_P^{a}}(D \setminus \{a\})(1 - x_a) + \smashoperator{\sum\limits_{a \in D \setminus S}} {\rho_P^{a}}(S) x_a,\phantom{\emptyset} S \subseteq D \label{eq:ineq1} \\
    {\pi}&\geq h_P(S) - \smashoperator{\sum\limits_{a \in S}} {\rho_P^{a}}(S \setminus \{a\})(1 - x_a) + \smashoperator{\sum\limits_{a \in D \setminus S}} {\rho_P^{a}}(\emptyset) x_a,\phantom{D} S \subseteq D \label{eq:ineq2}.\end{aligned}$$ Only one of the sets of inequalities and is required to define ${H_P}$ [@nemhauser1988].  Using these inequalities, the feasible region of ${H_P}$ can be formulated as $$\begin{aligned}
    {H_P}\coloneqq \left\lbrace (x,{\pi}) \in \{0,1\}^{\added{|D|}} \times {\mathbb{R}}\colon \eqref{eq:ineq1} \textrm{ or } \eqref{eq:ineq2} \right\rbrace.\end{aligned}$$ The number of inequalities required to define model in this manner grows exponentially with the number of arcs in a path, and the number of $s$–$t$ paths grows exponentially with the size of the network.  Enumerating all $s$–$t$ paths for scenario $(s, t) \in {\Omega}$ is impractical.  Nevertheless, lends itself to a delayed constraint generation algorithm.  Instead of adding inequalities for all possible paths, we add violated valid inequalities as needed (see Section \[subsec:bncsummary\]). 

To demonstrate the potential of this formulation, we next show that if we could obtain the convex hull of each set ${H_P}$, a relaxation that is at least as strong as the LP relaxation of the DEF (and ).  implies that the path-based formulation has the potential to yield a better LP relaxation if we can identify strong valid inequalities for $\operatorname{conv}({H_P})$.

\[theorem:1\] Consider the following LP relaxation of : $$\begin{aligned}
{3}
        \label{eq:path6}
        \begin{aligned}
            \min_{\added{x,} {\pi}}\ && \smashoperator{\sum\limits_{{\omega}= (s,t) \in {\Omega}}} p_{\omega}\replaced{{\pi}_s^t}{{\pi}_s^{\omega}} &&& \\
            \emph{s.t.}\ && \sum\limits_{a \in D} {c}_a x_a & \leq b && \\
            && (x, \replaced{{\pi}_s^t}{{\pi}_s^{\omega}}) & \in \operatorname{conv}({H_P}), &&P \in {\mathcal{P}}_{st},\ {\omega}= (s, t) \in {\Omega}\\
            && x_a & \in [0, 1], && a \in D.
        \end{aligned}
    \end{aligned}$$ For all $(\bar{x},\bar{{\pi}})$ feasible to , there exists $\hat{\pi}$ such that $(\bar{x},\hat{\pi})$ is feasible to and has objective value in not greater than the objective value of $(\bar{x},\bar{{\pi}})$ in . 

The above theorem is a consequence of Lemma \[lemma:1\].  For a fixed $(s, t) \in {\Omega}$ and $x \in [0, 1]^{\added{|D|}}$ satisfying $\sum_{a \in D} {c}_a x_a \leq b$, define $$\begin{aligned}
{3}
    {O}^{st}(x) \coloneqq \min_{\replaced{{\pi}_s}{{\pi}_s^t}}\ && \replaced{{\pi}_s}{{\pi}_s^t} \qquad &&& \\
    \textrm{s.t.}\ && (x, \replaced{{\pi}_s}{{\pi}_s^t}) & \in \operatorname{conv}({H_P}),\ \enspace && P \in {\mathcal{P}}_{st}.\end{aligned}$$ Recall that ${E}^{st}(x)$ is defined in . 

\[lemma:1\] Let ${\omega}= (s, t) \in {\Omega}$ and $\replaced{\bar{x}}{x} \in [0, 1]^{\added{|D|}}$ satisfy $\sum_{a \in D} {c}_a \replaced{\bar{x}_a}{x_a} \leq b$. Then, $${E}^{st}(\replaced{\bar{x}}{x}) \leq {O}^{st}(\replaced{\bar{x}}{x}).$$

The vector of all ones in ${\mathbb{R}}^{|N|}$ is feasible to .  A feasible solution to can be constructed as follows.  Let $P = \{\added{(i_0,i_1),} (i_1,i_2), \deleted{(i_2,i_3),} \ldots, \replaced{(i_{|P|-1},i_{|P|})}{(i_{n-1},i_n)} \}
    \in {\mathcal{P}}_{st}$ be a simple $s$–$t$ path, where $\replaced{i_0}{i_1} = s$ and $\replaced{i_{|P|}}{i_n} = t$. Let $\replaced{\bar{y}_{i_0,i_1}}{\bar{y}_{i_1,i_2}} =
    1$. For $k = \replaced{1}{2}, \ldots, \replaced{|P|-1}{n}$, let $\bar{y}_{i_k, i_{k+1}} = r_{i_{k-1}, i_k} \bar{y}_{i_{k-1}, i_k}$. Let $\bar{y}_{ij} = 0$ for all $(i, j) \in A \setminus P$, and $\bar{z}_{ij} =
    0$ for all $(i,j) \in A$.  Then $(\bar{y}, \bar{z})$ is feasible to .  Because the feasible regions of and are nonempty , both problems have an optimal solution.  Let $\pi^*$ be an optimal solution to , with corresponding optimal dual solution $(y^*, z^*)$.

Assume first there exists an $s$–$t$ path $P^*$ such that one of the constraints – is satisfied at equality for all $a \in P^*$. Thus, $$\begin{aligned}
            \pi^*_i = \max \left\lbrace r_a \pi^*_j - (r_a - q_a)u_j^{{\omega}} \replaced{\bar{x}_a}{x_a},\ q_a \pi^*_j \right\rbrace
    \end{aligned}$$ for all $a=(i, j) \in P^*$. $$\added{G^{LP}(P) \coloneqq \left\lbrace \mystrut(44,44) \right. } 
        \begin{alignedat}{2}
            \added{(x, \pi) \in [0,1]^{|D|}} &\added{\times {\mathbb{R}}^{|N|} \colon} && \\
            \added{\pi_i - r_a\pi_j} & \added{\geq 0,} && \added{a=(i,j) \in P \setminus D} \\
            \added{\pi_i - r_a\pi_j} & \added{\geq -(r_a - q_a) {u}_j^{\omega}x_a, \quad} && \added{a=(i,j) \in P \cap D} \\
            \added{\pi_i - q_a\pi_j} & \added{\geq 0,} && \added{a=(i,j) \in P \cap D} \\
            \added{\pi_t} & \added{= 1} &&
        \end{alignedat}
        \added{\left.  \mystrut(44,44) \right\rbrace. } $$ $$\begin{aligned}
        \added{{E}^{st}(\bar{x})} &\added{= \min_{\pi} \{\pi_s \colon (\bar{x}, \pi) \in G^{LP}(P^*) \} } \\
            &= \min_{\pi_s} \{ \pi_s \colon (\bar{x}, \pi_s) \in \operatorname{proj}_{(x,\pi_s)} \big( G^{LP}(P^*) \big) \}. 
    \end{aligned}$$ $$\begin{aligned}
        &\added{\operatorname{proj}_{(x,\pi_s)} G^{IP}(P^*) = {H_{P^*}}} \\
        \added{\implies} &\added{\operatorname{proj}_{(x,\pi_s)} G^{LP}(P^*) \supseteq \operatorname{conv}({H_{P^*}})} \\
        \added{\implies} &{E}^{st}(\replaced{\bar{x}}{x}) \leq \min\{{\pi}_s \colon (\replaced{\bar{x}}{x}, {\pi}_s) \in \operatorname{conv}({H_{P^*}})\} \leq {O}^{st}(\replaced{\bar{x}}{x}). 
    \end{aligned}$$

Finally, assume there is no $s$–$t$ path $P^* \in {\mathcal{P}}_{st}$ such that one of the constraints – is satisfied at equality for all $a \in P^*$. Consider any $s$–$t$ path $P \in {\mathcal{P}}_{st}$. By assumption, there exists an arc $(i, j) \in P$ such that constraints – are not binding.  By complementary slackness, the corresponding dual optimal $y^*_{ij}$ and $z^*_{ij}$ are $0$ the flow of path through the network is $0$.

From Lemma \[lemma:1\], we have $$\smashoperator{\sum_{{\omega}= (s,t) \in {\Omega}}} p_{\omega}\bar{{\pi}}_s^{{\omega}} \geq \enspace \smashoperator{\sum_{{\omega}= (s,t) \in {\Omega}}} p_{\omega}{O}^{st}(\bar{x}) \geq
        \enspace \smashoperator{\sum_{{\omega}= (s,t) \in {\Omega}}} p_{\omega}{E}^{st}(\bar{x}).$$ This implies there exists $\hat{\pi} \in {\mathbb{R}}^{\added{|N| \times |{\Omega}|}}$ such that $(\bar{x},\hat{\pi})$ is feasible to with objective value $\sum_{{\omega}= (s,t) \in {\Omega}} p_{\omega}\hat{\pi}_s^{\omega}\leq \sum_{{\omega}= (s,t) \in {\Omega}} p_{\omega}\bar{{\pi}}_s^{\omega}$.

Computational experiments by Ahmed and [@ahmed2011] show that the inequalities and may provide a poor approximation of the set $\operatorname{conv}({H_P})$. Thus, in the following sections, we describe additional inequalities that can used to approximate $\operatorname{conv}({H_P})$ for all $P \in {\mathcal{P}}_{st}$. We consider three different classes of inequalities, based on the traversal probabilities of interdicted arcs. 

Inequalities for the $q > 0$ case {#subsec:q>0}
---------------------------------

In this section, we assume $q_a > 0$ for all $a \in D$. For a ground set ${\mathcal{N}}$, Ahmed and [@ahmed2011] study valid inequalities for the mixed-integer set $$\begin{aligned}
    {F}= \left\lbrace (x, w) \in \{0,1\}^{\added{|{\mathcal{N}}|}} \times {\mathbb{R}}\colon w \leq f(\textstyle\sum_{i \in {\mathcal{N}}}{\alpha}_i x_i + {\beta}) \right\rbrace, \label{eq:barF}\end{aligned}$$ where ${\alpha}\in {\mathbb{R}}^{\added{|{\mathcal{N}}|}}_+$, ${\beta}\in {\mathbb{R}}$, and $f$ is a strictly concave, increasing, differentiable function.  They derive valid inequalities for the set ${F}$ and prove that they dominate the submodular inequality equivalents of and applied to ${F}$. These improved inequalities are shown empirically to yield significantly better relaxations than and . 

We translate the set ${H_P}$ to match the structure of ${F}$. Let ${\mathcal{N}}\coloneqq D$, $w \coloneqq -{\pi}$, $f(u) \coloneqq -\exp(-u)$, ${\beta}\coloneqq -\sum_{a \in P} \log(r_a)$. For $a \in D$, let $$\begin{aligned}
    {\alpha}_a \coloneqq \begin{cases} \log(r_a) - \log(q_a) & \textrm{if } a \in P \cap D, \\ 0 & \textrm{otherwise.}  \end{cases}\end{aligned}$$ Because $\log(r_a) > \log(q_a)$ for all $a \in D$, we have ${\alpha}\in {\mathbb{R}}^{\added{|D|}}_{+}$. With these definitions, ${H_P}$ is expressed in the form of ${F}$.

We now describe how to calculate valid inequalities for ${H_P}$ using the results of Ahmed and Atamt[ü]{}rk [@ahmed2011].  Without loss of generality, let $D \setminus S \coloneqq \{1, 2, \ldots, m\}$ be indexed such that ${\alpha}_1 \geq {\alpha}_2 \geq \ldots \geq {\alpha}_m$, and let ${A}_k = \sum_{a = 1}^k {\alpha}_a$ for $k \in D
\setminus S$, with ${A}_0 = 0$. Also define ${\alpha}(S) = \sum_{a \in S} {\alpha}_a$. The subadditive lifting inequality $$\begin{aligned}
        {\pi}\geq h_P(S) - \smashoperator{\sum_{a \in S}} {\phi}(-{\alpha}_a)(1 - x_a) + \smashoperator{\sum_{a \in D \setminus S}} {\rho_P^{a}}(S) x_a \label{eq:lifted1}\end{aligned}$$ is valid for ${H_P}$ for any set of interdicted arcs $S \subseteq D$, where ${\phi}$ is calculated as follows [@ahmed2011].  Consider the function $\zeta \colon {\mathbb{R}}_- \rightarrow {\mathbb{R}}$, calculated according to Algorithm \[alg:1\]. 

\[alg:1\] $k \leftarrow 0$ $\zeta({\eta}) \leftarrow -\exp(-{\alpha}(S) - {A}_k - {\beta}- {\eta}) + \sum_{a = 1}^k {\rho_P^{a}}(S) + \exp(-{\alpha}(S) - {\beta})$

For a provided value of ${\eta}$, Algorithm \[alg:1\] returns a specific $k$. With this $k$, let ${\phi}\colon {\mathbb{R}}_- \rightarrow {\mathbb{R}}$ be defined as $$\begin{aligned}
    {\phi}({\eta}) \coloneqq \begin{cases} \zeta(\mu_k - {A}_{k-1}) + {\rho_P^{k}}(S) \frac{b_k({\eta})}{{\alpha}_k} & \textrm{if }\mu_k - {A}_k \leq {\eta}\leq \mu_k - {A}_{k-1}, \\ \zeta({\eta}) & \textrm{otherwise,} \end{cases}\end{aligned}$$ where $\mu_k = -\log(-{\rho_P^{k}}(S)/{\alpha}_k) - {\alpha}(S) - {\beta}$ and $b_k({\eta}) = \mu_k - {A}_{k-1} - {\eta}$. Inequality dominates the general supermodular inequality [@ahmed2011]. 

We now construct Ahmed and Atamt[ü]{}rk’s inequalities that dominate the supermodular inequalities .  Let $S = \{1, 2, \ldots, n\}$ be indexed such that ${\alpha}_1 \geq {\alpha}_2 \geq \ldots \geq {\alpha}_n$. Also, let ${A}_k = \sum_{a = 1}^k {\alpha}_a$ for $k \in S$, with ${A}_0 = 0$. We define the function $\xi \colon {\mathbb{R}}_+
\rightarrow {\mathbb{R}}$ to be calculated according to Algorithm \[alg:2\]. 

\[alg:2\] $k \leftarrow 0$ $\xi({\eta}) \leftarrow -\exp(-{\alpha}(S) + {A}_k - {\beta}- {\eta}) - \sum_{a = 1}^k {\rho_P^{a}}(S \setminus \{a\}) + \exp(-{\alpha}(S) - {\beta})$

let $$\begin{aligned}
    {\psi}({\eta}) \coloneqq \begin{cases} \xi({A}_k - \nu_k) + {\rho_P^{k}}(S \setminus \{k\}) \frac{b_k({\eta})}{{\alpha}_k} & \textrm{if }{A}_{k-1} - \nu_k \leq {\eta}\leq {A}_k - \nu_k, \\ \xi({\eta}) & \textrm{otherwise,} \end{cases}\end{aligned}$$ where $\nu_k = {\alpha}(S) + \log(-{\rho_P^{k}}(S \setminus \{k\}) / {\alpha}_k) - {\beta}$ and $b_k({\eta}) = {A}_k - \nu_k - {\eta}$ for every $k \in \{1, 2, \ldots, r\}$. For any $S \subseteq D$, the inequality $$\begin{aligned}
    {\pi}\geq h_P(S) - \smashoperator{\sum_{a \in S}} {\rho_P^{a}}(S \setminus \{a\}) (1 - x_a) - \smashoperator{\sum_{a \in D \setminus S}} {\psi}({\alpha}_a)x_a \label{eq:lifted2}\end{aligned}$$ is valid for ${H_P}$. For all $a \in D \setminus P$, the coefficients in inequalities and are $0$.

We next discuss separation of the inequalities and for use in a delayed constraint generation framework.  Given a point $(\bar{x}, \bar{{\pi}}) \in [0, 1]^{\added{|D|}} \times {\mathbb{R}}$, we want to determine if there is a set $S$ such that $(\bar{x},\bar{{\pi}})$ violates either or .  For an integral vector $\bar{x}$, we construct the set $S = \{a \in D \colon \bar{x}_a = 1\}$. If $\bar{{\pi}} < h_P(S)$, then inequalities and are violated and can be added to the formulation of ${H_P}$ to cut off $(\bar{x}, \bar{{\pi}})$.

If $\bar{x}$ is not integral, we use Ahmed and Atamt[ü]{}rk’s scheme to heuristically construct a set $S$ using $\bar{x}$. This set is used in inequalities and to cut off $(\bar{x},
\bar{{\pi}})$, if possible.  We first consider constructing a set $S$ to apply inequality .  We first solve the following nonlinear program: $$\begin{aligned}
    \max_{z \in [0,1]^{\added{|D|}}}\ \frac{-\bar{{\pi}} + 1}{ -\exp(-\textstyle\sum_{a \in D}{\alpha}_a z_a - {\beta}) + 1} + \smashoperator{\sum\limits_{a \in D}} \frac{{\rho_P^{a}}(\emptyset)}{h_P(\emptyset) + 1} \bar{x}_a (1 - z_a).  \label{eq:sep1}\end{aligned}$$ Given a solution $\bar{z}$, we greedily round the fractional components in $\bar{z}$ that result in the least reduction in objective value when rounded to either $0$ or $1$. Given the rounded solution, say $\bar{z}'$, we set $S = \{a \in D
\colon \bar{z}'_a = 1\}$. If ${\pi}< h_P(S) - \sum_{a \in S} {\phi}(-{\alpha}_a)(1 - \bar{x}_a) + \sum_{a \in D
\setminus S} {\rho_P^{a}}(S) \bar{x}_a$, valid inequality can be added to the relaxation to cut off $(\bar{x}, \bar{{\pi}})$.

We solve a related nonlinear program to find a set $S$ to apply inequality : $$\begin{aligned}
    \max_{z \in [0,1]^{\added{|D|}}}\ \frac{-\bar{{\pi}} + 1}{ -\exp(-\textstyle\sum_{a \in D}{\alpha}_a z_a - {\beta}) + 1} - \smashoperator{\sum\limits_{a \in D}} \frac{{\rho_P^{a}}(\emptyset)}{h_P(\{a\}) + 1} (1 - \bar{x}_a) z_a.  \label{eq:sep2}\end{aligned}$$ We again greedily round the solution of to obtain the characteristic vector of $S$. If ${\pi}< h_P(S) - \sum_{a \in S} {\rho_P^{a}}(S \setminus \{a\}) (1 - \bar{x}_a) - \sum_{a \in D \setminus S} {\psi}({\alpha}_a) \bar{x}_a$, inequality cuts off $(\bar{x}, \bar{{\pi}})$.

Yu and Ahmed [@yu2017] consider deriving stronger valid inequalities by imposing a knapsack constraint on the set ${F}$ given in equation .  In particular, $$\begin{aligned}
    {F}= \left\lbrace (x, w) \in \{0,1\}^{\added{|{\mathcal{N}}|}} \times {\mathbb{R}}\colon w \leq f(\textstyle\sum_{i \in {\mathcal{N}}}{\alpha}_i x_i + {\beta}),\ \sum_{i \in {\mathcal{N}}} x_i \leq k \right\rbrace,\end{aligned}$$ where $k > 0$. Although our test instances include this structure, we do not explore applying their results here.  In particular, when interdicting arcs along an [*individual*]{} path, the knapsack constraint is only relevant if the number of arcs in the path under consideration is larger than the defender’s budget. 

Inequalities for $q = 0$ case {#subsec:q=0}
-----------------------------

When $q_a = 0$ for all $a \in D$, $\log(q_a)$ is no longer defined, and the results of Section \[subsec:q&gt;0\] do not apply.  We consider a different class of inequalities for this special case of $q$.

If all interdicted arc probabilities are $0$, $\bar{h}_P(x)$ simply reduces to $$\begin{aligned}
    \bar{h}_P(x) = \left[\, \prod_{a \in P} r_a \right] \left[ \enspace\ \smashoperator{\prod_{a \in P \cap D}}\, (1 - x_a) \right].\end{aligned}$$ If any arc $a \in P \cap D$ is interdicted, $\bar{h}_P(x)$ equals $0$. By setting $\hat{{\pi}} \coloneqq {\pi}/ \prod_{a \in P} r_a$, our relevant mixed-integer linear set is $$\begin{aligned}
    {H_P}&= \left\lbrace (x,\hat{{\pi}}) \in \{0,1\}^{\added{|D|}} \times {\mathbb{R}}\colon \hat{{\pi}} \geq \textstyle\prod_{a \in P \cap D} (1 - x_a) \right\rbrace.\end{aligned}$$ We are again interested in valid inequalities for ${H_P}$. The inequalities $$\begin{aligned}
    \hat{{\pi}} &\geq 1 - \smashoperator{\sum_{a \in P \cap D}} x_a \label{eq:french1}\end{aligned}$$ and $\hat{{\pi}} \geq 0$ are valid for ${H_P}$ [@fortet1960] and define its convex hull [@alkhayyal1983]. 

Inequalities for the mixed case {#subsec:qmixed}
-------------------------------

We now consider the case where $q_a > 0$ for some, but not all, $a \in D$. Let ${P_+}\coloneqq \{a \in P \cap D \colon q_a > 0\}$ and ${P_0}\coloneqq \{a \in P \cap D \colon q_a = 0\}$. We write $\bar{h}_P(x)$ as $$\begin{aligned}
    \bar{h}_P(x) &= r(P \setminus D) \bar{h}_{{P_+}}(x) \bar{h}_{{P_0}}(x).  \label{hprewrite}\end{aligned}$$ We introduce two new variables, ${{\pi}_+}\in \left[0, r({P_+}) \right]$ and ${{\pi}_0}\in [0, r({P_0})]$ to represent $\bar{h}_{{P_+}}(x)$ and $\bar{h}_{{P_0}}(x)$, respectively, and arrive at the formulation: $$\begin{aligned}
    {\pi}&\geq r(P \setminus D) {{\pi}_+}{{\pi}_0}\label{eq:theta1} \\
    {{\pi}_+}&\geq \bar{h}_{{P_+}}(x) \label{eq:thetaplus1} \\
    {{\pi}_0}&\geq \bar{h}_{{P_0}}(x).  \label{eq:thetazero1}\end{aligned}$$ Our approach is to use results from Sections \[subsec:q&gt;0\] and \[subsec:q=0\] to derive valid inequalities for and , respectively, and relax the nonconvex constraint using the McCormick inequalities [@mccormick1976], which in this case reduce to $$\begin{aligned}
    {\pi}& \geq r(P \setminus D) \left[ r({P_0}){{\pi}_+}- r({P_+})(r({P_0}) - {{\pi}_0}) \right] \label{eq:q0ineq1}\end{aligned}$$ and ${\pi}\geq 0$.

Let ${\gamma}^S \in {\mathbb{R}}^{\added{|{P_+}|}}$ and constant ${\delta}^S \in {\mathbb{R}}$ be coefficients from one inequality of the form (with $P_+$ taking the place of $P$) for some $S \subseteq \replaced{D}{P_+}$. An inequality of the form may also be used.  We relax to the constraint $$\begin{aligned}
    {{\pi}_+}& \geq \smashoperator{\sum_{a \in {P_+}}} {\gamma}^S_a x_a + {\delta}^S. \label{eq:fme1}\end{aligned}$$ Inequality is relaxed as ${{\pi}_0}\geq 0$, and, applying , $$\begin{aligned}
    {{\pi}_0}& \geq r({P_0}) \left( 1 - \smashoperator{\sum_{a \in {P_0}}} x_a \right) .   \label{eq:fme2}\end{aligned}$$ Finally, we project variables ${{\pi}_+}$ and ${{\pi}_0}$ out of inequality using Fourier-Motzkin elimination with constraints – and ${{\pi}_0}\geq 0$. This results in the two inequalities $$\begin{aligned}
    {\pi}& \geq r(P \setminus D) r({P_0}) \left[\ \smashoperator{\sum_{a \in {P_+}}} {\gamma}^S_a x_a +
    {\delta}^S - r({P_+})\smashoperator{\sum_{a \in {P_0}}} x_a \right] \label{eq:q01} \\
    {\pi}& \geq r(P \setminus D) r({P_0}) \left[\ \smashoperator{\sum_{a \in {P_+}}} {\gamma}^S_a x_a +
    {\delta}^S - r({P_+}) \right].  \label{eq:q02}\end{aligned}$$ Inequality is dominated by ${\pi}\geq 0$, because $\sum_{a \in {P_+}} {\gamma}^S_a x_a +
{\delta}^S
\leq {{\pi}_+}\leq r(P_+)$. Therefore, inequality is the only inequality we obtain from this procedure for the given inequality defined by $S$.

We next argue that for any integer solution $(\bar{x},\bar{{\pi}}) \in \{0, 1\}^{\added{|D|}} \times {\mathbb{R}}_+$, if it is not feasible (i.e., $(\bar{x},\bar{{\pi}}) \notin {H_P}$), then we can efficiently find an inequality of the form that this solution violates.  Indeed, first observe that if $\bar{x}_a = 1$ for any $a \in P_0$ then $\bar{h}_P(\bar{x}) = 0 \leq \bar{{\pi}}$, by assumption, and hence the solution is feasible.  So, we may assume $\bar{x}_a = 0$ for all $a \in P_0$. We set $S =
\{ a \in P_+ : \bar{x}_a = 1\}$, which yields $\sum_{a \in {P_+}} {\gamma}^S_a \bar{x}_a + {\delta}^S = h_{{P_+}}(S) = \bar{h}_{{P_+}}(\bar{x})$. Inequality thus yields $${\pi}\geq   r(P \setminus D) r({P_0})  \left[ \bar{h}_{{P_+}}(\bar{x}) - 0 \right] = \bar{h}_P(\bar{x})$$ by since $\bar{h}_{{P_0}}(\bar{x}) = r({P_0})$.

Path-based branch-and-cut algorithm {#subsec:bncsummary}
-----------------------------------

We next describe how the inequalities derived in Sections \[subsec:q&gt;0\]–\[subsec:qmixed\] can be used within a branch-and-cut algorithm to solve the path-based formulation .  Similar to the use of Benders decomposition in a branch-and-cut algorithm described in Section \[subsec:benders\], the algorithm is based on solving a master problem via branch-and-cut in which the constraints are relaxed and approximated with cuts. 

Let $(\bar{x}, \bar{{\pi}})$ be a solution obtained in the branch-and-cut algorithm with integral $\bar{x}$. Since are relaxed, we must check if this solution satisfies these constraints, and if not, finding a cut that this solution violates.  We assign arc traversal probabilities ${\sigma}_a$ for $a \in D$ according to the formula $$\begin{aligned}
    {\sigma}_a &= r_a^{1 - \bar{x}_a} q_a^{\bar{x}_a}.  \label{reldef1}\end{aligned}$$ For each ${\omega}= (s,t) \in {\Omega}$, we find a maximum-reliability path $\bar{P} \in {\mathcal{P}}_{st}$. If $\bar{{\pi}}_s^t \geq
\bar{h}_{\bar{P}}(\bar{x})$ then this solution is feasible to for this ${\omega}$, since, by construction, $\bar{h}_{\bar{P}}(\bar{x}) = \max_{P \in {\mathcal{P}}_{st}} \bar{h}_P(\bar{x})$. Otherwise, depending on if $q_a =
0$ for any or all $a \in \bar{P}$, we derive a cut from one of Sections \[subsec:q&gt;0\]–\[subsec:qmixed\], using this path $\bar{P}$ and scenario ${\omega}$. By construction, for the solution $\bar{x}$, this cut enforces that $\pi_s^t
\geq \bar{h}_{\bar{P}}(\bar{x})$, and hence (i) cuts off the current infeasible solution $(\bar{x}, \bar{{\pi}})$, and (ii) implies that any solution $(x,{\pi})$ of the updated master problem with $x =\bar{x}$ will satisfy for this ${\omega}$. This ensures that the branch-and-cut algorithm is finite (since at most finitely many such cuts are needed at finitely many integer solutions) and correct (since any infeasible solution obtained in the algorithm will be cut off). 

One may also attempt to generate cuts at solutions $(\bar{x}, \bar{{\pi}})$ in which $\bar{x}$ is not necessarily integer, in order to improve the LP relaxation.  To find a path for a scenario ${\omega}= (s,t)$ in this case, we again use formula to define arc reliabilities, and then find the most reliable $s$-$t$ path.  With this path, we again apply the methods in Sections \[subsec:q&gt;0\]–\[subsec:qmixed\] to attempt to derive a violated cut.  When $\bar{x}$ is not integral, this approach is not guaranteed to find a violated cut, even if one exists.  Thus, to increase the chances a cut is found, we may also consider identifying another path by assigning arc costs as follows: $$\begin{aligned}
    {\sigma}_a &= (1 - \bar{x}_a)r_a +  \bar{x}_a q_a,\end{aligned}$$ and then finding a maximum-reliability path using these arc reliabilities.  Our preliminary tests did not reveal any obvious benefit of using one particular arc-reliability calculation method. 

In our implementation, before starting the branch-and-cut algorithm we solve a relaxation of the master problem in which constraints are dropped, and the integrality constraints are relaxed.  After solving this LP relaxation, we attempt to identify violated cuts for the relaxation solution, and if found, we add them to the master relaxation and repeat this process until no more violated cuts are found.  We then begin the branch-and-cut process with all cuts found when solving the LP relaxation included in the mixed-integer programming formulation.  This allows the solver to use these cuts to generate new cuts in the master problem relaxation.  At any point in the branch-and-bound process where a solution $(\bar{x},\bar{{\pi}})$ with $\bar{x}$ integral is obtained, we identify if there is any violated cut (as discussed above) and if so add it to the formulation via the lazy constraint callback function. 

Computational experiments {#sec:computational}
=========================

We test the different solution methods on SNIP instances from Pan and Morton [@pan2008], which consist of a network of $783$ nodes and $2586$ arcs, $320$ of which can be interdicted by the defender.  Each instance considers the same $456$ scenarios.  We consider three cases for the $q$ vector outlined by Pan and Morton [@pan2008].  In particular, we consider instances with $q \coloneqq 0.5r$, $q \coloneqq 0.1r$, and $q
\coloneqq 0$. Pan and Morton also considered $q_a$ values independently sampled from a $U(0,0.5)$ distribution.  Since the DEF for most of these instances could be solved within $30$ seconds with little or no branching, we exclude these in our experiments.  The cost of interdiction is ${c}_a = 1$ for all $a \in D$. Seven different budget levels were tested for each network and $q$ level, and there are five test instances for each budget level and $q$ level. 

We consider the following algorithms in our tests. 

-   Solve DEF with a MIP solver

-   Solve compact DEF with a MIP solver

-   Benders branch-and-cut algorithm on DEF ; refer to

-   Branch-and-cut algorithm on path-based decomposition

DEF and C-DEF are solved using the commercial MIP solver IBM ILOG CPLEX 12.6.3.  BEN is the Benders branch-and-cut algorithm as implemented in Bodur [et al.]{}[@bodur2016 Section 4].  The branch-and-cut algorithms, BEN and PATH, were implemented within the CPLEX solver using lazy constraint callbacks.  Ipopt 3.12.1 was used to solve nonlinear models and in PATH.  The computational tests were run on a 12-core machine with two 2.66GHz Intel Xeon X5650 processors and 128GB RAM, and one-hour time limit was imposed.  We allowed four threads for DEF, while all other algorithms were limited to one thread.  We used the CPLEX default relative gap tolerance of $10^{-4}$ to terminate the branch-and-bound process.  The algorithms were written in Julia 0.4.5 using the Julia for Mathematical Optimization framework [@lubin2015]. 

In our implementation of the BEN and PATH algorithms, when generating cuts at the root LP relaxation, we used the following procedure to focus the computational effort on scenarios that yield cuts.  After each iteration, we make a list of scenarios for which violated cuts were found in that iteration.  Then, in the next iteration, we first only attempt to identify cuts for scenarios in this list.  If successful, the list is further reduced to the subset of scenarios that yielded a violated cut.  If we fail to find any violated cuts in the current list, we re-initialize the list with all scenarios and attempt to identify violated cuts for each scenario. 

  -------- ----- ----------- ----------- ----------- -----------
                                                     
  $q$      $b$           DEF       C-DEF         BEN        PATH
  $0.5r$   30      (2)1652.1        23.9       298.5        54.1
           40      (2)2678.1        35.6       337.7        32.9
           50      (2)2242.4        31.7       355.8        42.8
           60      (2)1891.7        34.0       454.1        70.4
           70      (1)1730.0        27.4       487.4        63.9
           80          826.0        14.4       386.7        27.0
           90          460.3         6.9       385.0        26.8
  $0.1r$   30            (5)       117.5       438.7        55.3
           40            (5)       466.0       812.3       301.2
           50            (5)       327.9       787.2       265.3
           60            (5)       638.6       793.7       194.1
           70            (5)       851.6       983.0       335.0
           80            (5)   (1)1257.4   (1)1106.3    (1)891.1
           90            (5)   (2)2329.4   (3)1269.4   (2)1497.5
  $0$      30            (5)        96.1       487.0        61.7
           40            (5)       196.4       556.8       177.0
           50            (5)       332.3       663.6       250.8
           60            (5)       369.0      1141.8       374.8
           70            (5)       341.0       739.5       354.4
           80            (5)       262.9       535.1       151.5
           90            (5)       540.0       765.7       635.2
  -------- ----- ----------- ----------- ----------- -----------

  : Computational results[]{data-label="tab:results"}

  -------- ----- -------- -------- -- -------- -------- --------
                                                        
    $q$     $b$       BEN     PATH       C-DEF      BEN     PATH
   $0.5r$   30     10.64%   10.64%       2.38%    6.11%    6.03%
            40     11.34%   11.35%       2.80%    6.28%    6.20%
            50     11.22%   11.23%       2.40%    5.60%    5.62%
            60     10.54%   10.55%       2.07%    4.76%    4.80%
            70      8.88%    8.89%       1.83%    4.14%    4.34%
            80      6.25%    6.25%       0.97%    3.13%    3.19%
            90      3.92%    3.91%       0.43%    1.92%    2.02%
   $0.1r$   30     22.47%   22.49%       6.72%    6.11%    6.67%
            40     26.22%   26.18%       9.35%    8.73%    8.21%
            50     27.54%   27.48%       9.38%    9.01%    9.28%
            60     28.16%   28.08%       9.93%    8.75%    9.35%
            70     28.92%   28.83%       9.67%    9.30%    9.70%
            80     30.88%   30.90%      10.85%   11.57%   11.96%
            90     33.07%   32.98%      12.45%   15.18%   15.46%
    $0$     30     25.15%   25.30%       7.07%    6.72%    6.50%
            40     28.45%   28.55%       9.62%    8.07%    9.25%
            50     30.54%   30.78%      11.02%   10.03%   11.51%
            60     32.07%   32.45%      12.37%   12.48%   14.20%
            70     32.60%   33.30%      11.21%   11.82%   14.92%
            80     33.28%   33.99%      10.43%   11.44%   15.87%
            90     36.17%   36.97%      12.89%   15.45%   20.94%
  -------- ----- -------- -------- -- -------- -------- --------

  : Average root gaps before and after CPLEX cuts.  The LP relaxation of BEN is the same as DEF and C-DEF.[]{data-label="tab:lprelax"}

  ------ ----- ----------- ------------- ------------- ------------- -- ---------- ---------
                                                                                   
   $q$    $b$          DEF         C-DEF           BEN          PATH           BEN      PATH
   0.5r   30     (2)1080.7        1406.6       39800.8       51512.2         266.5      18.9
          40     (2)1982.7        2470.2       21909.8        7697.8         313.6      21.2
          50     (2)1996.3        2516.4         24112         17161         324.8      21.1
          60     (2)1723.7        3270.6       35094.2       35860.2         404.5        22
          70       (1)1464        2242.4       24387.2       34456.2         456.9      22.1
          80           727        1235.4        5989.4        3509.4         378.7      21.2
          90           556           372        4930.2        7512.2         379.1      17.3
   0.1r   30           (5)        3870.8        6576.2       10156.8         417.2      26.2
          40           (5)       17798.2       87973.4      108169.6         583.9      29.9
          50           (5)        8950.4       89062.4         70182         507.6      29.7
          60           (5)       22457.8       95993.8       53885.8         528.7      30.9
          70           (5)       28358.4      176635.4       83490.6         448.8      32.4
          80           (5)    (1)45951.8   (1)163500.3     (1)242734      (1)497.9   (1)33.9
          90           (5)   (2)101616.7   (3)170902.5   (2)327946.3      (3)474.8   (2)37.5
    0     30           (5)          2205        5613.4        9114.8           462      27.2
          40           (5)        3821.6       13803.4       29569.6         488.2      28.1
          50           (5)        6037.8       30095.8       47806.4         511.4      28.2
          60           (5)        5144.6        127214         65936         546.7      27.9
          70           (5)          5379         45454       50358.6         477.9        28
          80           (5)          2667       18581.6       11890.6         406.7      28.1
          90           (5)        7289.4       52156.2       64416.2         424.5      29.5
  ------ ----- ----------- ------------- ------------- ------------- -- ---------- ---------

  : Average branch-and-bound nodes and cut generation time[]{data-label="tab:bnb"}

Table \[tab:results\] reports the average computational time, in seconds, over the five instances for each combination of the vector $q$ and budget level $b$. The average times include only instances solved within the time limit – the number of unsolved instances are enclosed in parentheses.  We find that PATH and C-DEF are consistently faster than both DEF and BEN.  PATH and C-DEF have comparable performance, with C-DEF being somewhat faster on the $q=0.5r$ instances (the easiest set) and PATH being somewhat faster on the $q=0.1r$ instances (the hardest set). 

Table \[tab:lprelax\] reports the root gap after adding LP cuts for the decomposition algorithms BEN and PATH, relative to the optimal solution value, before and after the addition of cuts from CPLEX.  The root gap is calculated as $(z^* - z^{LP}) / z^*$, where $z^*$ is the instance’s optimal value and $z^{LP}$ is the objective value obtained after running the initial cutting plane loop on the relaxed master problem.  The reported value is the arithmetic mean over the five instances for each combination of interdicted traversal probabilities $q$ and budget level $b$. DEF, C-DEF, and BEN all have the same LP relaxation value, so we show this gap only for BEN.  The columns under the heading “LP relaxation after CPLEX cuts” refer to the gap that is obtained at the root node, after CPLEX finishes its process of adding general-purpose cuts.  We exclude DEF from these results since in many instances the root node has not completed processing in the time limit.  These results show that the initial LP relaxation value obtained using the cuts in the PATH method is very similar to that obtained using the BEN (or the DEF or C-DEF formulations).  We also find that the general purpose cuts in CPLEX significantly improve the relaxation value for all formulations, with the greatest improvement coming in the C-DEF formulation.  The gaps obtained after the addition of CPLEX cuts are quite similar between PATH and BEN, with the exception of the $q=0$ instances, where the CPLEX cuts seem to be more effective for the BEN formulation. 

Table \[tab:bnb\] displays the average number of branch-and-bound nodes of each algorithm and the time spent generating cuts in BEN and PATH.  This includes time spent generating cuts during the initial cutting plane loop as well as within the branch-and-cut process.  The mean is over instances that solved within the time limit.  Consistent with the results on root relaxation gaps, we find that C-DEF requires the fewest number of branch-and-bound nodes, and that PATH and BEN methods require a similar number of branch-and-bound nodes.  On the other hand, significantly less time is spent generating cuts in the PATH method than in the BEN method, which explains the better performance of PATH.  In particular, generating cuts in PATH requires solving a maximum-reliability path problem for each destination node, as opposed to the Benders approach which requires solving a linear program for each scenario. 

Conclusion
==========

We have proposed two new methods for solving maximum-reliability SNIP.  The first method is to solve a more compact DEF, and the second is based on a reformulation derived from considering the reliability for each path separately.  Both of these methods outperform state-of-the-art methods on a class of SNIP instances from the literature.  An interesting direction for future research is to investigate whether ideas similar to the path-based formulation might be useful for solving different variants of SNIP, such as the maximum-flow network interdiction problem. 

[^1]: etowle@wisc.edu

[^2]: jim.luedtke@wisc.edu
