A Linear Input Dependence Model for Interdependent Networks

In this section we describe the randomized rounding approximation algorithms which will be the subject of the remainder of Section 3, as well as explaining why a randomized rounding approach was chosen. In the BIDM the binary linking variable $y_{ij}^{kl}$ for each interdependent pair $(ij,kl)\in I$ either forces the parent arc to be saturated (when $y_{ij}^{kl}=1$) or the child arc to be unused (when $y_{ij}^{kl}=0$). A common approximation technique for such a program with binary variables would be to solve the linear relaxation (restricting $y_{ij}^{kl}$ to $[0,1]$ instead of $\{0,1\}$) and then deterministically round the relaxed binary variables to 1 if above some threshold or 0 if below some threshold, fixing the rounded variables and then finally solving the resulting pure LP. Unfortunately there is no deterministic threshold that could be applied in this case that would guarantee the existence of a feasible solution. As an alternative we consider the related idea of randomized rounding algorithms.

Let $\mathbf{x^{\boldsymbol{*}}}$ be the optimal flow for the LP relaxation of the BIDM. For each interdependence $(ij,kl)\in I$ we set $y_{ij}^{kl}$ to 1 with probability $P_{ij}^{kl}$ and 0 otherwise, where $P_{ij}^{kl}$ describes the relative saturation of either the parent or the child. If the resulting binary variables define an infeasible LP then the random binary variables can be re-realized, repeating as necessary until a feasible solution is found. As it would be reasonable to base the value of $y_{ij}^{kl}$ on the saturation of either the parent or the child, we define three basic types of randomized rounding algorithm: randomized rounding based on child flow (RR-Child) in which $P_{ij}^{kl}:=x_{kl}^{*}/u_{kl}$, randomized rounding based on parent flow (RR-Parent) in which $P_{ij}^{kl}:=x_{ij}^{*}/u_{ij}$, and the control test fair randomized rounding (RR-Fair) in which $P_{ij}^{kl}:=0.5$.

There is a danger that such a randomized rounding algorithm may not terminate regardless of how many times the random variables are realized. This can occur in certain cases where $P_{ij}^{kl}$ is allowed to equal exactly 0 or 1, as shown in the examples in Appendix B. To prevent this we define modified versions of the RR-Child and RR-Parent schemes with alternate definitions for $P_{ij}^{kl}$. For any $\epsilon\in[0,0.5]$, let $\textsc{RR-Child}(\epsilon)$ be equivalent to RR-Child described above, but with its value of $P_{ij}^{kl}$ restricted to the interval $[\epsilon,1-\epsilon]$ (that is, equal to $\max\{\min\{x_{kl}^{*}/u_{kl},1-\epsilon\},\epsilon\}$) in order to avoid probabilities too close to 0 or 1. For the purposes of our study we will explore the special cases of $\textsc{RR-Child}(0.00)$, which is equivalent to RR-Child as described above, as well as $\textsc{RR-Child}(0.01)$ and $\textsc{RR-Child}(0.05)$, which restrict $P_{ij}^{kl}$ to the intervals $[0.01,0.99]$ and $[0.05,0.95]$, respectively. We define $\textsc{RR-Parent}(\epsilon)$ similarly to $\textsc{RR-Child}(\epsilon)$, with $P_{ij}^{kl}=\max\{\min\{x_{ij}^{*}/u_{ij},1-\epsilon\},\epsilon\}$.

Finally we note some simple theoretical bounds. Clearly the LP relaxation of the BIDM provides a lower bound for its optimal value, and any feasible solution to the BIDM obtained from a randomized rounding scheme provides an upper bound. Unfortunately in general these bounds may not be very tight, and there exist example networks for which they are arbitrarily far from the true optimum (see Appendix C). For this reason theoretical bounds are not very useful in evaluating the LP relaxation and randomized rounding approximations, and so in order to learn more about the typical errors in these approximations we turn to empirical results obtained through computational trials.

In order to constitute reasonable approximation methods we should expect that the LP relaxation and the randomized rounding schemes produce typically similar objective values to those of the BIDM solution, and that the randomized rounding schemes typically find a feasible solution in a small number of iterations. Moreover, in order to justify the utility of the LP relaxation in finding a feasible solution to the BIDM we should expect all RR-Child and RR-Parent schemes to typically produce higher-quality solutions than the RR-Fair scheme, as the RR-Child and RR-Parent schemes require an LP relaxation pre-solve while the RR-Fair scheme does not. It is also reasonable to expect that all approximation methods perform better at lower interdependence densities, as fewer interdependencies leads to fewer relaxed constraints. These expectations were substantiated through use of the computational experiments explained below.

The approximation methods described above were empirically tested on a large number of random interdependent networks. Each network was solved once as a MILP by treating its interdependencies as binary, once as an LP by treating its interdependencies as linear, and then once again as a MILP using each of the randomized rounding schemes (based on the LP solution). For each solution method the objective value was recorded, and for the randomized rounding schemes the number of realizations required to first achieve a feasible MILP solution was recorded. In order to evaluate the effect of the interdependence density as well as the network’s size, arc density, and interdependence structure on these results, several network parameters were varied between problem sets.

We used a modified version of NETGEN [15], a random MCNF problem generator, to create our test cases, adding a procedure to generate pairs of interdependent arcs alongside the standard MCNF problem parameters and used these to define instances of the BIDM, as well as a preliminary check to ensure that all problem instances were feasible. We tested two main types of problem: those in which parent arcs corresponded to demand nodes, and those in which parent arcs were distributed randomly throughout the network. We will refer to these problems as Structured-Type and Unstructured-Type, respectively. Structured-Type problems represent a more realistic interdependence structure that might be seen in infrastructure networks, while Unstructured-Type problems assume no particular interdependence structure in order to evaluate the models’ behaviors in a more general setting. For Structured-Type problems, we applied a transformation to convert parent demand nodes into equivalent parent arcs, relaxing the demand values at these nodes as well as the supply values. The full source code for the computational trial generation can be viewed online [28].

For each problem type, we generated test networks on 256 and 512 nodes. For each test network the number of arcs was either 4 or 8 times the number of nodes. For Structured-Type problems we randomly selected either 2%, 5%, 10%, or 15% of the demand nodes as parents, with their corresponding children being chosen uniformly at random from the arc set. For Unstructured-Type problems we randomly selected either 1%, 2%, 5%, or 10% of the arcs within the network as parents and children.

For each combination of parameters, 60 test cases were generated. In all cases, 20% of the nodes were sources and 20% were sinks (in Structured-Type problems, this refers to the original number of sink nodes, before some were converted into transshipment nodes during the parent arc transformation). Each arc’s cost was chosen uniformly at random from the interval $[1,100]$, with its capacity chosen from $[100,500]$ and with 100% of the NETGEN skeleton arcs having maximum cost. The total supply was 10,000 per 256 nodes. The total number of test cases generated over all parameter combinations listed above was 1960.

Each test began by solving the binary input dependence MILP for the current network using CPLEX. Next the LP relaxation was solved using CPLEX, and the flow values from the optimal solution were recorded. Finally each of the seven randomized rounding schemes described above was executed in turn, using the flow values from the LP solution to determine the probability of setting each binary variable to 0 or 1. During each attempt a random binary vector was generated according to the rules of the current scheme, the resulting constraints were added to the MILP formulation, and the resulting problem was evaluated with CPLEX. If infeasible, a new random binary vector was generated from the same probability distribution and the process was repeated. If no feasible solution was found after an arbitrarily-chosen cutoff of 1000 attempts then the attempt was labeled as a failure and the next scheme was tested.

In this section we summarize the results of the computational trials described above. Full data tables can be found in Appendix E, while the raw data can be viewed online [13]. Across all groups of trials the RR-Child and RR-Parent schemes performed so similarly that only the RR-Child and RR-Fair results will be discussed here.

The computational results of the previous section largely confirm our expectations for the approximation algorithms as explained in Section 3.1. The LP relaxation and the resulting randomized rounding solutions typically produce reasonable bounds for the MILP solution. The mean relative error of the LP relaxation had a mean of less than 1.7% for all trial groups. For Structured-Type problems, with one exception the $\textsc{RR-Child}(0.00)$ scheme was able to find a feasible solution at a relative error of less than 6.5%, while $\textsc{RR-Child}(0.01)$ achieved nearly the same accuracy with a better success rate. For Unstructured-Type problems the $\textsc{RR-Child}(0.00)$ scheme achieved the smallest relative errors and the $\textsc{RR-Child}(0.01)$ scheme achieved the smallest failure rate, but both performed well especially for the smallest trials. In all trial sets, all RR-Child and RR-Parent schemes achieved significantly smaller errors and higher success rates than the RR-Fair scheme, which further justifies the utility in solving an instance of the LIDM while approximating the solution to the BIDM. However, these results also show that the error and failure rates of the approximation methods increase with the interdependence density. While this was expected, it also implies that the use of the randomized rounding algorithms should be limited to networks with relatively low interdependence density. That being said, in spite of the upward trend, the relative error in the LP relaxation was so small across all trial sets that it is still a useful approximation method even for networks with relatively large interdependence densities, which motivates the development of an efficient solution algorithm for the LIDM.

The remainder of the paper will be dedicated to developing and proving the correctness of a generalized version of network simplex capable of solving an arbitrary instance of the MCNFLI. This problem is technically a special case of the MCNF problem with side constraints studied in Mathies and Mevert 1998 [18], but we can exploit the structure of our side constraints to define a simpler solution algorithm. Similar algorithms have been studied in Calvete 2003 [6] for the general equal flow problem and in Bahçeci and Feyziog̃lu 2012 [4] for the proportional flow problem.

In this section we define some terminology that will be used throughout the remainder of the paper. For brevity we will use the terms “arc”, “flow variable”, and “column” interchangeably to refer to an arc, an arc’s associated flow variable, and a flow variable’s associated column of $\mathbf{A}$. Similarly we will use the terms “node”, “flow conservation constraint”, and “row” interchangeably to refer to a node, a node’s associated flow conservation constraint, and a constraint’s associated row of $\mathbf{A}$.

During each iteration of simplex, let $B$ denote the set of basic variables. For such $B$, let $\mathbf{B}$ be the submatrix of $\mathbf{A}$ containing only basic columns. Let $L$ and $U$ be the sets of nonbasic variables at their lower or upper bound, respectively. Let a variable be called interdependent if it is involved in any interdependence as either a parent arc, child arc, or slack variable, and independent otherwise. For each interdependence $(ij,kl)\in I$, we say that flow variables $x_{ij}$ and $x_{kl}$ and slack variable $s_{ij}^{kl}$ are all linked to each other. We call arcs $ij$ and $kl$ partners of each other. For any interdependent arc $ij$, let $\hat{a}_{ij}$ be defined as the coefficient of variable $x_{ij}$ in its corresponding linking constraint (5), which is $-\alpha_{ij}^{kl}$ if $ij$ is a parent and 1 if it is a child.

In this section we prove the main result of the paper, the characterization of the basis of the LIDM-B, which will require first establishing some preliminary observations regarding the problem’s structure. A basis of the LIDM-B consists of a set $B$ of basic variables whose corresponding columns of $\mathbf{A}$ are full-rank, plus a partition of all nonbasic variables into sets $L$ and $U$. We begin by observing the rank of $\mathbf{A}$.

The proof of this follows from the fact that, as an $m\times n$ incidence matrix, submatrix $\mathbf{\hat{A}}$ has rank $m-1$ [5, p. 280] while the $p\times p$ identity matrix $\mathbf{I}$ in the lower right has rank $p$, and the block structure of $\mathbf{A}$ makes it clear that the last $p$ rows are linearly independent when taken together with any full-rank subset of the first $m$ rows. This tells us that our basis must always contain $m+p-1$ basic variables.

Our next goal is to find a graph-based characterization of this basis analogous to the spanning tree basis of the standard MCNF. We assume without loss of generality that $G$ contains at least one spanning tree made up entirely of independent arcs. As in the MCNF, the subgraph of basic arcs must contain only a single component.

The proof of this relies on the fact that the basis matrix $\mathbf{B}$ can be rearranged into a block form with one block for each component of the basic subgraph. Each block is an incidence matrix and thus rank-deficient by 1, meaning that adding all rows corresponding to one component results in a zero row. Since $\mathbf{B}$ is rank-deficient by only 1, only one zero row should be possible, which implies that only one component is possible.

Unlike the standard MCNF it is possible for the MCNFLI basis to contain cycles of arcs, however no such cycle can consist entirely of independent arcs.

The proof of this is almost identical to that of the MCNF problem [5, p. 283]. Observations 2 and 3, combined, imply that the subgraph of basic independent arcs must form a spanning forest of $G$. Let a spanning forest with $k$ components be referred to as a $k$-spanning forest.

For the remainder of this paper we will let $r$ represent the number of basic interdependent variables, which could include a combination of interdependent arcs and slack variables. Note that it is always the case that $r\geq p$ or else Observation 3 would be violated. This allows us to define the structure of the basic independent arcs.

This follows from Observations 2 and 3. Then the collection of basic variables consists of a spanning forest of independent arcs, plus possibly some interdependent arcs linking the components of the forest, plus possibly some slack variables. For the remainder of the paper we will use $F=\{T_{1},\dots,T_{r-p+1}\}$ to refer to the spanning $(r-p+1)$-forest of $G$ formed by the basic independent arcs, where $T_{h},h=1,\dots,r-p+1$ are the subgraphs of $G$ that form the components of the forest.

Define matrix $\mathbf{D}\in\mathbb{R}^{r\times r}$ as

where $r_{a}$ is the number of basic interdependent arcs, $r_{s}$ is the number of basic slack variables, $\mathbf{\hat{I}}\in\mathbb{R}^{p\times r_{s}}$ and $\mathbf{Q}\in\mathbb{R}^{p\times r_{a}}$ are the submatrices of $\mathbf{\hat{I}}$ and $\mathbf{\hat{Q}}$ containing only basic columns, and $\delta_{ij}^{h}$ is defined as

This matrix proves to be the key for determining the rank of $\mathbf{B}$.

Adopting the terminology of [2, 4, 6], we refer to an $(r-p+1)$-spanning forest of $G$ as a good $(r-p+1)$-forest with respect to a set $S$ of $r$ non-independent variables if $\mathbf{D}$ is full-rank. Then the preceding results allow us to finally characterize the basis of the LIDM-B.

The proof of Theorem 1 follows directly from Lemma 2 and the definition of a basis. It should also be noted that the basis characterization presented in Theorem 1 can be modified to depend on only a submatrix of $\mathbf{D}$, as shown in Appendix A. In either case the spanning forest element of the the LIDM-B basis makes it structurally similar to that of the standard MCNF basis. This enables us to develop a generalized network simplex algorithm for its solution, as will be explained in the next section.

The vector of reduced costs $\mathbf{c^{\boldsymbol{\pi}}}\in\mathbb{R}^{n+p}$ is defined by $(\mathbf{c^{\boldsymbol{\pi}}})^{\prime}:=\mathbf{c}^{\prime}-\boldsymbol{\pi}^{\prime}\mathbf{A}$, where $\boldsymbol{\pi}\in\mathbb{R}^{m+p}$ is the potential vector [5, p. 84]. There is one reduced cost $c_{ij}^{\pi}$ for each flow variable $x_{ij}$ and one $(c^{\pi})_{ij}^{kl}$ for each slack variable $s_{ij}^{kl}$. There is one potential $\pi_{i}$ for reach node $i\in V$ and one $\pi_{ij}^{kl}$ for each interdependence $(ij,kl)\in I$. Given our definitions for $\mathbf{c^{\boldsymbol{\pi}}}$ and $\mathbf{A}$, we have the relationships

The reduced cost of a basic variable is zero, allowing us to use equations (39)–(42) to solve for the potentials by adopting a two-step strategy based on that of Calvete 2003 [6] and Bahçeci and Feyziog̃lu 2012 [4]. First we make an initial guess $\boldsymbol{\pi}$ for the potential vector. To do this, within each tree $T_{h}$ of $F$, arbitrarily select a root and arbitrarily assign it a potential. Assign a potential of zero to all interdependencies of basic slack variables and arbitrarily assign a potential to all remaining interdependencies. Then equations (39) can be used to solve for all remaining node potentials within each tree in order from root to leaf, since $c_{ij}^{\pi}=0$ for all $ij\in T_{h}$.

Second, we test whether the guessed potentials were correct by testing whether they produce reduced costs of zero for all basic variables when substituted into (39)–(42). If so, then $\boldsymbol{\pi}$ is the correct potential vector and we can move on. If not, then we must calculate a “corrected potential” vector $\mathbf{\tilde{\boldsymbol{\pi}}}$ defined by

The vector $\boldsymbol{\sigma}$ of correction terms is defined as the solution of the linear system $\mathbf{D}^{\prime}\boldsymbol{\sigma}=\mathbf{c^{\boldsymbol{\pi}}}$, which includes the reduced cost vector $\mathbf{c^{\boldsymbol{\pi}}}$ resulting from the initial guess $\boldsymbol{\pi}$.

At the end of the two-step process the resulting corrected potential vector $\mathbf{\tilde{\boldsymbol{\pi}}}$ can be applied to equations (39)–(42) to compute the reduced cost of any nonbasic variable. As noted at the end of Section 4.2, there is a way to characterize the basis using only a submatrix of $\mathbf{D}$. This reduced version of $\mathbf{D}$ also gives rise to a reduced version of the above algorithm explained in Appendix A. In addition, note that the full version of the above algorithm only needs to be used during the initial iteration. Between consecutive iterations of simplex, within most components of $F$ the previous potential vector will still satisfy (39) after the change of basis. As a result the previous potential vector $\boldsymbol{\pi}$ can be used to calculate the correction vector $\boldsymbol{\sigma}$ without the need to recalculate all node potentials within each component of $F$.

To explain, if a component $T_{h}$ of $F$ remains unchanged between iterations (meaning that no arc within $T_{h}$ leaves the basis and no independent arc incident to a node in $T_{h}$ enters the basis), or if $T_{h}$ loses an arc and splits into two separate components, then all previous node potentials within $T_{h}$ still satisfy (39). Only if an independent arc enters the basis, causing two components of $F$ to merge, can one of these equations become violated. In this event equations (39) can be satisfied by increasing all node potentials within one of the merged components by an appropriate constant. If independent arc $ij$ with endpoints in separate components $T_{h}$ and $T_{g}$ of $F$, respectively, enters the basis, then either $(\pi_{i}-c_{ij})-\pi_{j}$ can be added to all node potentials in $T_{g}$ or $(\pi_{j}+c_{ij})-\pi_{i}$ can be added to all node potentials in $T_{h}$. All basic slack potentials can then be set to 0 as before, after which the correction vector can be calculated and used it to calculate the corrected potential vector.

In any case, if all reduced costs of variables in $L$ are nonnegative and all reduced costs of variables in $U$ are nonpositive, then the current basic solution is optimal and the algorithm may terminate. Otherwise we choose some variable in $L$ with a negative reduced cost or in $U$ with a positive reduced cost to enter the basis and conduct change of basis.

After deciding on a variable to enter the basis we conduct pivoting and update the solution vector. Before describing this process it will be helpful to first discuss how to calculate the values of the basic variables given a particular basis. This may be a necessary step during the algorithm’s first iteration starting from an initial basis, and it will also introduce some of the notation and results required to define the change of basis procedure in Section 5.2.2.