Determining r- and (r,s)-Robustness of Digraphs Using Mixed Integer Linear Programming

TL;DR

This paper presents a novel mixed integer linear programming method to efficiently determine the $r$- and $(r,s)$-robustness of directed graphs, a key property for resilient control algorithms, using only the graph Laplacian.

Contribution

It introduces the first application of mixed integer programming to compute graph robustness, providing an efficient and scalable approach based on the Laplacian matrix.

Findings

The method accurately determines robustness parameters in simulations.

It outperforms prior algorithms in computational efficiency.

The approach is applicable to various directed graphs.

Abstract

There has been an increase in the use of resilient control algorithms based on the graph theoretic properties of $r$- and $(r,s)$-robustness. These algorithms guarantee consensus of normally behaving agents in the presence of a bounded number of arbitrarily misbehaving agents if the values of the integers $r$ and $s$ are sufficiently large. However, determining an arbitrary graph's robustness is a highly nontrivial problem. This paper introduces a novel method for determining the $r$- and $(r,s)$-robustness of digraphs using mixed integer linear programming; to the best of the authors' knowledge it is the first time that mixed integer programming methods have been applied to the robustness determination problem. The approach only requires knowledge of the graph Laplacian matrix, and can be formulated with binary integer variables. Mixed integer programming algorithms such as branch-and-bound are used to iteratively tighten the lower and upper bounds on $r$ and $s$. Simulations are presented which compare the performance of this approach to prior robustness determination algorithms.

Figures

Equations

$$L_{i,j}=\{\begin{array}{ll}|{\mathcal{N}}_i|& \text{if }i=j,\\ -1& \text{if }j\in {\mathcal{N}}_i,\\ 0& \text{if }j\notin {\mathcal{N}}_i.\end{array}$$

$$\displaystyle\mathcal{T}=\big{\{}$$

$$\displaystyle|S_{1}\cap S_{2}|=0\big{\}}.$$

$$\Theta =\{$$

$$(|{\mathcal{X}}_{S_2}^r|=|S_2|)\text{ or }(|{\mathcal{X}}_{S_1}^r|+|{\mathcal{X}}_{S_2}^r|\ge s)\}.$$

$$\mathcal{R}(S)$$

$$r_{\max }(\mathcal{D})=$$

$$|S_1|>0,|S_2|>0,|S_1\cap S_2|=0.$$

$$r_{\max }(\mathcal{D})=$$

$${\sigma }_j(S)$$

$$L_j\mathbf{\sigma }(S)$$

$$\mathcal{R}(S)=\underset{j}{\max }{L}_j\mathbf{\sigma }(S),j\in \{1,\dots ,n\}.$$

$$L_j\mathbf{\sigma }$$

$$=|{\mathcal{N}}_j|{\sigma }_j-\sum _{q\in {\mathcal{N}}_j\cap S}{\sigma }_q-\sum _{q\in {\mathcal{N}}_j\backslash S}{\sigma }_q.$$

$$\mathcal{R}(S)=\underset{j\in S}{\max }|{\mathcal{N}}_j\backslash S|=\underset{j\in S}{\max }(L_j\mathbf{\sigma }(S)).$$

$$\underset{j\in \{1,\dots ,n\}}{\max }{L}_j\mathbf{\sigma }(S)$$

$$=\underset{j\in S}{\max }{L}_j\mathbf{\sigma }(S).$$

$$\mathcal{R}(S)=\underset{j}{\max }(L_j\mathbf{\sigma }(S)),j\in \{1,\dots ,n\}.$$

$$\max \left(\mathcal{R}(S_1),\mathcal{R}(S_2)\right)=$$

$$\max \left(\underset{i\in \{1,\dots ,n\}}{\max }\left(L_i\mathbf{\sigma }(S_1)\right),\underset{j\in \{1,\dots ,n\}}{\max }\left(L_j\mathbf{\sigma }(S_2)\right)\right)$$

$$r_{\max }(\mathcal{D})=$$

$$\displaystyle\underset{(S_{1},S_{2})\in\mathcal{T}}{\min}\max\big{(}\max_{i}\left(L_{i}\bm{\sigma}(S_{1})\right),\max_{j}\left(L_{j}\bm{\sigma}(S_{2})\right)\big{)}.$$

$$\Sigma (S_1,S_2)=(\mathbf{\sigma }(S_1),\mathbf{\sigma }(S_2)),(S_1,S_2)\in \mathcal{T}.$$

$$\displaystyle\mathcal{B}=\bigg{\{}$$

$$\displaystyle 1\leq\bm{1}^{T}\bm{b}^{2}\leq(n-1),\ \bm{b}^{1}+\bm{b}^{2}\preceq\bm{1}\bigg{\}}.$$

$$|S_1|>0$$

$$|S_2|>0$$

$$|S_1|<n$$

$$|S_2|<n$$

$$|S_1\cap S_2|=0⟹\mathbf{\sigma }(S_1)+\mathbf{\sigma }(S_2)⪯1.$$

$$b_j^1=1$$

$$b_j^1=0$$

$$b_j^2=1$$

$$b_j^2=0$$

$$r_{\max }(\mathcal{D})=$$

$${\mathbf{b}}^1+{\mathbf{b}}^2⪯1$$

$$1\le 1^T{\mathbf{b}}^1\le (n-1)$$

$$1\le 1^T{\mathbf{b}}^2\le (n-1)$$

Full text

††thanks: The authors would like to acknowledge the support of the Automotive Research Center (ARC) in accordance with Cooperative Agreement W56HZV-14-2-0001 U.S. Army TARDEC in Warren, MI, and the Award No W911NF-17-1-0526.

, ,

1 Introduction

Consensus on shared information is fundamental to the operation of multi-agent systems. In context of mobile agents, it enables formation control, rendezvous, distributed estimation, and many more objectives. Although a vast number of algorithms for consensus exist, many are unable to tolerate the presence of misbehaving agents subject to attacks or faults. Recent years have seen an increase of attention on resilient algorithms that are able to operate despite such misbehavior. Many of these algorithms have been inspired by [9], which is one of the seminal papers on consensus in the presence of adversaries; [15, 43, 14], which outline discrete- and continuous-time algorithms along with necessary and sufficient conditions for scalar consensus in the presence of Byzantine adversaries; and [39, 38, 32, 34], which outline algorithms for multi-agent vector consensus of asynchronous systems in the presence of Byzantine adversaries. Some of the most recent resilience-based results that draw upon these papers include state estimation [18, 20, 19, 21], rendezvous of mobile agents [23, 22], output synchronization [11], simultaneous arrival of interceptors [17], distributed optimization [31, 29], reliable broadcast [33, 43], clock synchronization [8], randomized quantized consensus [5], self-triggered coordination [26], and multi-hop communication [30].

A large number of results on network resilience are based upon the graph theoretical properties known as $r$-robustness and $(r,s)$-robustness [15, 43]. $r$-robustness and $(r,s)$-robustness are key notions included in the sufficient conditions for convergence of several resilient consensus algorithms including the ARC-P [14], W-MSR [15], SW-MSR [25], and DP-MSR [4] algorithms. Given an upper bound on the global or local number of adversaries in the network, the aforementioned resilient algorithms guarantee convergence of normally behaving agents’ states to a value within the convex hull of initial states if the integers $r$ and $s$ are sufficiently large.

A key challenge in implementing these resilient algorithms is that determining the $r$- and $(r,s)$-robustness of arbitrary digraphs is an NP-hard problem in general [13, 42]. The first algorithmic analysis of determining the values of $r$ and $s$ for arbitrary digraphs was given in [13]. The algorithms proposed in [13] employ an exhaustive search to determine the maximum values of $r$ and $s$ for a given digraph, and have exponential complexity w.r.t. the number of nodes in the network. In [42] it was shown that the decision problem of determining if an arbitrary graph is $r$-robust for a given integer $r$ is coNP-complete in general. Subsequent work has focused on methods to circumvent this difficulty, including graph construction methods which increase the graph size while preserving given values of $r$ and $s$ [15, 6]; lower bounding $r$ with the isoperimetric constant and algebraic connectivity of undirected graphs [27]; and demonstrating the behavior of $r$ as a function of certain graph properties [42, 28, 44, 35, 7, 24]. In particular, it has been shown that the $r$-robustness of some specific classes of graphs can be exactly determined in polynomial time from certain graph parameters. Examples include $k$-circulant graphs [35] and lattice-based formations [7, 24]. Another recent approach has used machine learning to correlate characteristics of certain graphs to the values of $r$ and $s$ [40], but these correlations are inherently stochastic in nature and do not provide explicit guarantees. Despite the impressive results of prior literature, methods to either approximate or determine exactly the $r$- and $(r,s)$-robustness of arbitrary digraphs remain relatively rare. Methods to determine the exact $r$- or $(r,s)$-robustness of arbitrary undirected graphs are also uncommon. Finding more efficient or practical ways of determining the robustness of arbitrary graphs in general, and digraphs in particular, remains an open problem.

In response to this open problem, this paper introduces novel methods for determining the $r$- and $(r,s)$-robustness of digraphs and undirected graphs using mixed integer linear programming (MILP). These methods only require knowledge of the graph Laplacian matrix and are zero-one MILPs, i.e. with all integer variables being binary. To the best of our knowledge, this is the first time the robustness determination problem has been formulated in this way. These results connect the problem of graph robustness determination to the extensive and well-established literature on integer programming and linear programming.

This paper makes the following specific contributions:

We present a method to determine the maximum integer for which a nonempty, nontrivial, simple digraph is $r$-robust using mixed integer linear programming. 2. 2.

We present a method which determines the $(r,s)$-robustness of a digraph using linear programming. Here, the $(r,s)$-robustness of a digraph refers to the maximal $(r,s)$ integer pair according to a lexicographical order for which a given digraph is $(r,s)$-robust, as first described in [13]. Furthermore, we show that our method can also determine the maximum integer $F$ for which a digraph is $(F+1,F+1)$-robust, which is not considered in [13]. 3. 3.

We present two mixed integer linear programs whose optimal values provide lower and upper bounds on the maximum $r$ for which a nonempty, nontrivial, simple digraph is $r$-robust. These two formulations exhibit a lower complexity than the method in the first contribution described above.

The contributions of this paper provide several advantages. First, expressing the robustness determination problem in MILP form allows for approximate lower bounds on a given digraph’s $r$-robustness to be iteratively tightened using algorithms such as branch-and-bound. Lower bounds on the maximum value of $s$ for which a given digraph is $(r,s)$ robust (for a given nonnegative integer $r$) can also be iteratively tightened using the approach in this paper. Prior algorithms are only able to tighten the upper bound on the maximum robustness for a given digraph or undirected graph. Second, this formulation enables commercially available solvers such as Gurobi or MATLAB’s intlinprog to be used to find the maximum robustness of any digraph. Finally, experimental results using this new formulation suggest a reduction in computation time as compared to the centralized algorithm proposed in [13].

Part of this paper was previously submitted as a conference paper [36]. The extensions to the conference version include the following:

• •

The more general case of determining $(r,s)$-robustness is considered.

• •

Two optimization problems with reduced-dimension binary vector variables are given whose optimal values provide a lower bound and an upper bound, respectively, on the maximum value of $r$ for which a graph is $r$-robust.

This paper is organized as follows: notation is presented in Section 2, the problem formulation is given in Section 3, determining the $r$-robustness of digraphs is treated in Section 4, determining the values of $s$ for which a digraph is $(r,s)$-robust for a given $r$ is treated in Section 5, methods to obtain upper and lower bounds on the maximum $r$ for which a digraph is $r$-robust are presented in Section 6, a brief discussion about the advantages of the MILP formulations is given in Section 7, simulations demonstrating our algorithms are presented in Section 8, and a brief conclusion is given in 9.

2 Notation

The sets of real numbers and integers are denoted $\mathbb{R}$ and $\mathbb{Z}$, respectively. The sets of nonnegative real numbers and integers are denoted $\mathbb{R}_{+}$ and $\mathbb{Z}_{+}$, respectively. $\mathbb{R}^{n}$ denotes an $n$-dimensional vector space over the field $\mathbb{R}$, $\mathbb{Z}^{n}$ represents the set of $n$ dimensional vectors with integer entries, and $\{0,1\}^{n}$ represents a binary vector of dimension $n$. Scalars are denoted in normal text (e.g. $x\in\mathbb{R}$) while vectors are denoted in bold (e.g. $\bm{x}\in\mathbb{R}^{n}$). The notation $x_{i}$ denotes the $i$th entry of vector $\bm{x}$.

The inequality symbol $\preceq$ denotes a componentwise inequality between vectors; i.e. for $\bm{x},\bm{y}\in\mathbb{R}^{n}$, $\bm{x}\preceq\bm{y}$ implies $x_{i}\leq y_{i}\ \forall i\in\{1,\ldots,n\}$. A vector of ones is denoted $\bm{1}$, and a vector of zeros is denoted $\bm{0}$, where the length of each vector will be implied by the context. The union, intersection, and set complement operations are denoted by $\cup,\ \cap$, and $\setminus$, respectively. The cardinality of a set $S$ is denoted as $|S|$, and the empty set is denoted $\emptyset$. The infinity norm of a vector is denoted $\left\lVert\cdot\right\rVert_{\infty}$. The notations $C(n,k)=\begin{pmatrix}n\\ k\end{pmatrix}=n!/(k!(n-k)!)$ are both used in this paper to denote the binomial coefficient with $n,k\in\mathbb{Z}_{+}$. Given a set $S$, the power set of $S$ is denoted $\mathcal{P}(S)=\{A:A\subseteq S\}$. Given a function $f:D\rightarrow R$, the image of a set $A\subseteq D$ under $f$ is denoted $f(A)$. Similarly, the preimage of $B\subseteq R$ under $f$ is denoted $f^{-1}(B)$. The logical OR operator, AND operator, and NOT operator are denoted by $\lor,\land,\lnot$, respectively. The lexicographic cone is defined as $K_{lex}=\{0\}\cup\{\bm{x}\in\mathbb{R}^{n}:x_{1}=\ldots=x_{k}=0,\ x_{k+1}>0\}$ for some $0\leq k<n$. The lexicographic ordering on $\mathbb{R}^{n}$ is defined as $\bm{x}\leq_{lex}\bm{y}$ if and only if $\bm{y}-\bm{x}\in K_{lex}$, with $\bm{x},\bm{y}\in\mathbb{R}^{n}$ [3, Ch. 2].

A directed graph (digraph) is denoted as $\mathcal{D}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{1,\ldots,n\}$ is the set of indexed vertices and $\mathcal{E}$ is the edge set. This paper will use the terms vertices, agents, and nodes interchangeably. A directed edge is denoted $(i,j)$, with $i,j\in\mathcal{V}$, meaning that agent $j$ can receive information from agent $i$. The set of in-neighbors for an agent $j$ is denoted $\mathcal{N}_{j}=\{i\in\mathcal{V}:(i,j)\in\mathcal{E}\}$. The minimum in-degree of a digraph $\mathcal{D}$ is denoted $\delta^{in}(\mathcal{D})=\min_{j\in\mathcal{V}}|\mathcal{N}_{j}|$. Occasionally, $\mathcal{G}=(\mathcal{V},\mathcal{E})$ will be used to denote an undirected graph, i.e. a digraph in which $(i,j)\in\mathcal{E}\iff(j,i)\in\mathcal{E}$. The graph Laplacian $L$ for a digraph (or undirected graph) is defined as follows, with $L_{i,j}$ denoting the entry in the $i$th row and $j$th column:

[TABLE]

3 Problem Formulation

The notions of $r$- and $(r,s)$-robustness are graph theoretical properties used to describe the communication topologies of multi-agent networks. Examples of such networks include stations in a power grid, satellites in formation, or a group of mobile robots. In these networks, edges model the ability for one agent $i$ to transmit information to another agent $j$. Prior literature commonly considers simple digraphs, which have no repeated edges or self edges [16, 15, 12, 39, 37]. More specifically, simple digraphs satisfy $(i,i)\notin\mathcal{E}\ \forall i\in\mathcal{V}$, and if the directed edge $(i,j)\in\mathcal{E}$, then it is the only directed edge from $i$ to $j$. Prior work also commonly considers nonempty and nontrivial graphs, where $|V|>1$.

Assumption 1**.**

This paper considers nonempty, nontrivial, simple digraphs.

The property of $r$-robustness is based upon the notion of $r$-reachability. The definitions of $r$-reachability and $r$-robustness are as follows:

Definition 1** ([15]).**

Let $r\in\mathbb{Z}_{+}$ and $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be a digraph. A nonempty subset $S\subset\mathcal{V}$ is $r$-reachable if $\exists i\in S$ such that $|\mathcal{N}_{i}\backslash S|\geq r$.

Definition 2** ([15]).**

Let $r\in\mathbb{Z}_{+}$. A nonempty, nontrivial digraph $\mathcal{D}=(\mathcal{V},\mathcal{E})$ on $n$ nodes $(n\geq 2)$ is $r$-robust if for every pair of nonempty, disjoint subsets of $\mathcal{V}$, at least one of the subsets is $r$-reachable. By convention, the empty graph $(n=0)$ is 0-robust and the trivial graph $(n=1)$ is 1-robust.

The property of $(r,s)$-robustness is based upon the notion of $(r,s)$-reachability. The definitions of $(r,s)$-reachability and $(r,s)$-robustness are as follows:

Definition 3** ([15]).**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be a nonempty, nontrivial, simple digraph on $n\geq 2$ nodes. Let $r,s\in\mathbb{Z}_{+}$, $0\leq s\leq n$. Let $S$ be a nonempty subset of $\mathcal{V}$, and define the set $\mathcal{X}^{r}_{S}=\{j\in S:|\mathcal{N}_{j}\backslash S|\geq r\}$. We say that $S$ is an $(r,s)$-reachable set if there exist $s$ nodes in $S$, each of which has at least $r$ in-neighbors outside of $S$. More explicitly, $S$ is $(r,s)$-reachable if $|\mathcal{X}^{r}_{S}|\geq s$.

Definition 4** ([15]).**

Let $r,s\in\mathbb{Z}_{+}$, $0\leq s\leq n$. Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be a nonempty, nontrivial, simple digraph on $n\geq 2$ nodes. Define $\mathcal{X}^{r}_{S}=\{j\in S:|\mathcal{N}_{j}\backslash S|\geq r\}$, $S\subset\mathcal{V}$. The digraph $\mathcal{D}$ is $(r,s)$-robust if for every pair of nonempty, disjoint subsets $S_{1},S_{2}\subset\mathcal{V}$, at least one of the following conditions holds:

1. A)

$|\mathcal{X}^{r}_{S_{1}}|=|S_{1}|$, 2. B)

$|\mathcal{X}^{r}_{S_{2}}|=|S_{2}|$, 3. C)

$|\mathcal{X}^{r}_{S_{1}}|+|\mathcal{X}^{r}_{S_{2}}|\geq s$.

The properties of $r$- and $(r,s)$-robustness are used to quantify the ability of several resilient consensus algorithms to guarantee convergence of normally behaving agents in the presence of Byzantine and malicious adversaries, collectively referred to in this paper as misbehaving agents [15, 11, 4, 25, 14]. Larger values of $r$ and $s$ generally imply the ability of networks applying these resilient algorithms to tolerate a greater number of misbehaving agents in the network. For a more detailed explanation of the properties of $r$-robustness and $(r,s)$-robustness, the reader is referred to [16, 15, 42].

It should be clear from Definitions 2 and 4 that determining $r$ and $(r,s)$-robustness for a digraph $\mathcal{D}=(\mathcal{V},\mathcal{E})$ by using an exhaustive search method is a combinatorial problem, which involves checking the reachabilities of all nonempty, disjoint subsets of $\mathcal{V}$. For notational purposes, we will denote the set of all possible pairs of nonempty, disjoint subsets of $\mathcal{V}$ as $\mathcal{T}\subset\mathcal{P}(\mathcal{V})\times\mathcal{P}(\mathcal{V})$. More explicitly, $\mathcal{T}$ is defined as

[TABLE]

It was shown in [13] that $|\mathcal{T}|=\sum_{p=2}^{n}\begin{pmatrix}n\\ p\end{pmatrix}(2^{p}-2)$.111Since $(S_{1},S_{2})\in\mathcal{T}\implies(S_{2},S_{1})\in\mathcal{T}$, the total number of unique nonempty, disjoint subsets is $(1/2)|\mathcal{T}|$, denoted as $R(n)$ in [13]. As a simple example, Figure 1 depicts all elements of $\mathcal{T}$ for a graph of 3 agents, i.e. all possible ways to choose two nonempty, disjoint subsets from the graph.

When considering a particular digraph $\mathcal{D}$, there may be multiple values of $r$ for which $\mathcal{D}$ is $r$-robust. Similarly, there may be multiple values of $r$ and $s$ for which $\mathcal{D}$ is $(r,s)$-robust. The following properties of robust graphs demonstrate this characteristic. Note that $r$-robustness is equivalent to $(r,1)$-robustness; i.e. $\mathcal{D}$ is $r$-robust if and only if it is $(r,1)$-robust [13, Property 5.21] [15, Section VII-B].

Property 1** ([10], Prop. 5.13).**

Let $\mathcal{D}$ be an arbitrary, simple digraph on $n$ nodes. Suppose $\mathcal{D}$ is $(r,s)$-robust with $r\in\mathbb{N}$ and $s\in\{0,\ldots,n\}$. Then $\mathcal{D}$ is also $(r^{\prime},s^{\prime})$-robust $\forall r^{\prime}:\ 0\leq r^{\prime}\leq r$ and $\forall s^{\prime}:\ 1\leq s^{\prime}\leq s$.

Property 2** ([10], Prop. 5.20).**

Let $\mathcal{D}$ be an arbitrary, simple digraph on $n$ nodes. Suppose $\mathcal{D}$ is $(r,s)$-robust with $r\in\mathbb{N}$ and $s\in\{1,\ldots,n\}$. Then $\mathcal{D}$ is $(r-1,s+1)$-robust.

As an example, if a digraph $\mathcal{D}_{1}$ is $4$-robust, it is $(4,1)$-robust, and therefore by Property 1 it is also simultaneously 3-robust, 2-robust, and 1-robust. In addition, if a digraph $\mathcal{D}_{2}$ is $(5,4)$-robust, then it is simultaneously $(r^{\prime},s^{\prime})$-robust for all integers $0\leq r^{\prime}\leq 5$ and $0\leq s^{\prime}\leq 4$. Moreover, by Property 2, $\mathcal{D}_{2}$ is also $(4,5)$-robust, $(3,6)$-robust, $(2,7)$-robust, and $(1,8)$-robust. For notational purposes, we denote the set of all values for which a digraph $\mathcal{D}$ is $(r,s)$-robust as $\Theta$, where $\Theta\subset\mathbb{Z}_{+}\times\mathbb{Z}_{+}$. By Definition 4, $\Theta$ is explicitly defined as

[TABLE]

Note that the conditions of (3) are simply an alternate way of expressing the conditions of Definition 4.

To characterize the resilience of graphs however, prior literature has generally been concerned with only a few particular values of $r$ and $s$ for which a given digraph is $r$- or $(r,s)$-robust. For $r$-robustness, the value of interest is the maximum integer $r$ for which the given digraph is $r$-robust.

Definition 5**.**

We denote the maximum integer $r$ for which a given digraph $\mathcal{D}$ is $r$-robust as ${r_{\max}}(\mathcal{D})\in\mathbb{Z}_{+}$.

Several resilient algorithms guarantee convergence of the normal agents when the adversary model is $F$-total or $F$-local in scope,222An $F$-total adversary model implies that there are at most $F\in\mathbb{Z}_{+}$ misbehaving agents in the entire network. An $F$-local adversary model implies that each normal agent has at most $F$ misbehaving agents in its in-neighbor set. and the digraph is $(2F+1)$-robust. The value of ${r_{\max}}(\mathcal{D})$ therefore determines the maximum adversary model under which these algorithms can operate successfully. Furthermore, all other values of $r$ for which a digraph $\mathcal{D}$ is $r$-robust can be determined from ${r_{\max}}(\mathcal{D})$ by using Property 1.

For $(r,s)$-robustness, there are two $(r,s)$ pairs of interest. The authors of [13] order the elements of $\Theta$ using a lexicographical total order, where elements are ranked by $r$ value first and $s$ value second. More specifically, $(r_{1},s_{1})\leq_{lex}(r_{2},s_{2})$ if and only if $\begin{bmatrix}r_{2}-r_{1}\\ s_{2}-s_{1}\end{bmatrix}\in K_{lex}$, where $K_{lex}$ is the lexicographic cone defined in Section 2. Their algorithm $DetermineRobustness$ finds the maximum element of $\Theta$ with respect to this order. For notational clarity, we denote this maximum element as $(r^{*},s^{*})\in\Theta$.

Definition 6**.**

Let $\Theta$ be defined as in (3). The element $(r^{*},s^{*})$ is defined as the maximum element of $\Theta$ under the lexicographical order on $\mathbb{R}^{2}$.

The other $(r,s)$ pair of interest is $(F_{\max}+1,F_{\max}+1)$, where $F_{\max}=\max(\{F\in\mathbb{Z}_{+}:(F+1,F+1)\in\Theta\})$. Several resilient algorithms guarantee convergence of the normally behaving agents when the (malicious [15]) adversary model is $F$-total in scope and the digraph is $(F+1,F+1)$-robust. The value $F_{\max}$ determines the maximum malicious adversary model under which these algorithms can operate successfully. The value of $(F_{\max}+1,F_{\max}+1)$ does not always coincide with the $(r^{*},s^{*})$-robustness of the digraph. A simple counterexample is given in Figure 2, where the $(r^{*},s^{*})$-robustness of the graph is $(2,1)$ but the value of $(F_{\max}+1,F_{\max}+1)$ is equal to $(1,1)$.

The purpose of this paper is to present methods using mixed integer linear programming to determine ${r_{\max}}(\mathcal{D})$, the $(r^{*},s^{*})$-robustness of $\mathcal{D}$, and the $(F_{\max}+1,F_{\max}+1)$-robustness of $\mathcal{D}$ for any nonempty, nontrivial, simple digraph $\mathcal{D}$.

Problem 1**.**

Given an arbitrary nonempty, nontrivial, simple digraph $\mathcal{D}$, determine the value of ${r_{\max}}(\mathcal{D})$.

Problem 2**.**

Given an arbitrary nonempty, nontrivial, simple digraph $\mathcal{D}$, determine the $(r^{*},s^{*})$-robustness of $\mathcal{D}$.

Problem 3**.**

Given an arbitrary nonempty, nontrivial, simple digraph $\mathcal{D}$, determine the $(F_{\max}+1,F_{\max}+1)$-robustness of $\mathcal{D}$.

Finally, the values of $r$ for which a digraph can be $r$-robust lie within the interval $0\leq r\leq\left\lceil n/2\right\rceil$ [10, Property 5.19]. Since $r$-robustness is equivalent to $(r,1)$-robustness, the values of $r$ for which a graph can be $(r,s)$-robust fall within the same interval. The values of $s$ for which a digraph can be $(r,s)$-robust lie within the interval $1\leq s\leq n$.333Footnote 8 in [15] offers an excellent explanation for restricting $s$ to this range by convention. However, we will use an abuse of notation by denoting a graph as $(r,0)$-robust for a given $r\in\mathbb{Z}_{+}$ if the graph is not $(r,1)$-robust.

4 Determining $r$-Robustness using Mixed Integer Linear Programming

In this section we will demonstrate a method for solving Problem 1 using a mixed integer linear program (MILP) formulation. An MILP will be presented whose optimal value is equal to ${r_{\max}}(\mathcal{D})$ for any given nonempty, nontrivial, simple digraph $\mathcal{D}$.

First, an equivalent way of expressing the maximum robustness ${r_{\max}}(\mathcal{D})$ of a digraph $\mathcal{D}$ is derived. This equivalent expression will clarify how ${r_{\max}}(\mathcal{D})$ can be determined by means of an optimization problem. Given an arbitrary, simple digraph $\mathcal{D}=(\mathcal{V},\mathcal{E})$ and a subset $S\subset V$, the reachability function $\mathcal{R}:\mathcal{P}(\mathcal{V})\rightarrow\mathbb{Z}_{+}$ is defined as follows:

[TABLE]

In other words, the function $\mathcal{R}(S)$ returns the maximum $r$ for which $S$ is $r$-reachable. Using this function, the following Lemma presents an optimization formulation which yields ${r_{\max}}(\mathcal{D})$:

Lemma 1**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph with $|\mathcal{V}|=n$. Let ${r_{\max}}(\mathcal{D})$ be defined as in Definition 5. The following holds:

[TABLE]

{pf}

Note that $S_{1},S_{2}\in\mathcal{P}(\mathcal{V})$ and the three constraints of the RHS of (1) imply that the feasible set consists of all subsets $S_{1},S_{2}$ such that $(S_{1},S_{2})\in\mathcal{T}$, as per (2). In addition, $\max\left(\mathcal{R}(S_{1}),\mathcal{R}(S_{2})\right)=m$ implies $\mathcal{R}(S_{1})=m$ or $\mathcal{R}(S_{2})=m$. Let $(S_{1}^{*},S_{2}^{*})$ be a minimizer of (1). Then $\max\left(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*})\right)\leq\max\left(\mathcal{R}(S_{1}),\mathcal{R}(S_{2})\right)$ $\forall(S_{1},S_{2})\in\mathcal{T}$. Therefore $\forall(S_{1},S_{2})\in\mathcal{T},$ either \mathcal{R}(S_{1})\geq\max\big{(}\mathcal{R}(S_{1}^{*}), \mathcal{R}(S_{2}^{*})\big{)} or $\mathcal{R}(S_{2})\geq\max\left(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*})\right)$. This satisfies the definition of $r$-robustness as per Definition 2, therefore $\mathcal{D}$ is at least $\max\left(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*})\right)$-robust. This implies ${r_{\max}}(\mathcal{D})\geq\max\left(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*})\right)$.

We next show that ${r_{\max}}(\mathcal{D})=\max\left(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*})\right)$. We prove by contradiction. Recall from Definition 5 that ${r_{\max}}(\mathcal{D})$ is the maximum integer $r$ for which $\mathcal{D}$ is $r$-robust, which means $\mathcal{D}$ is ${r_{\max}}(\mathcal{D})$-robust by definition. Suppose ${r_{\max}}(\mathcal{D})>\max\left(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*})\right)$. This implies $\mathcal{R}(S_{1}^{*})<{r_{\max}}(\mathcal{D})$ and $\mathcal{R}(S_{2}^{*})<{r_{\max}}(\mathcal{D})$. Since the nonempty, disjoint subsets $(S_{1}^{*},S_{2}^{*})\in\mathcal{T}$ satisfy $\mathcal{R}(S_{1}^{*})<{r_{\max}}(\mathcal{D})$ and $\mathcal{R}(S_{2}^{*})<{r_{\max}}(\mathcal{D})$, by the negation of Definition 2 this implies that $\mathcal{D}$ is not ${r_{\max}}(\mathcal{D})$-robust. However, this contradicts the definition of ${r_{\max}}(\mathcal{D})$ being the largest integer for which $\mathcal{D}$ is $r$-robust (Definition 5). This provides the desired contradiction; therefore ${r_{\max}}(\mathcal{D})=\max\left(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*})\right)$. $\blacksquare$

Remark 1**.**

Using the definition of $\mathcal{T}$ in (2), the constraints on the RHS of (1) can be made implicit [3, section 4.1.3] as follows:

[TABLE]

We demonstrate next that the objective function of (1) can be expressed as a function of the network Laplacian matrix. Recall that $n=|\mathcal{V}|$ and that $\{0,1\}^{n}$ represents a binary vector of dimension $n$. The indicator vector $\bm{\sigma}(\cdot):\mathcal{P}(\mathcal{V})\rightarrow\{0,1\}^{n}$ is defined as follows: for all $S\in\mathcal{P}(\mathcal{V})$,

[TABLE]

In other words the $j$th entry of $\bm{\sigma}(S)$ is 1 if the node with index $j$ is a member of the set $S\in\mathcal{P}(\mathcal{V})$, and zero otherwise. It is straightforward to verify that $\bm{\sigma}:\mathcal{P}(\mathcal{V})\rightarrow\{0,1\}^{n}$ is a bijection. Therefore given $\bm{x}\in\{0,1\}^{n}$, the set $\bm{\sigma}^{-1}(\bm{x})\in\mathcal{P}(\mathcal{V})$ is defined by $\bm{\sigma}^{-1}(\bm{x})=\{j\in\mathcal{V}:x_{j}=1\}$. Finally, observe that for any $S\in\mathcal{P}(\mathcal{V})$, $|S|=\bm{1}^{T}\bm{\sigma}(S)$. The following Lemma demonstrates that for any $S\in\mathcal{P}(\mathcal{V})$, the function $\mathcal{R}(S)$ can be determined as an affine function of the network Laplacian matrix and the indicator vector of $S$:

Lemma 2**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph, let $L$ be the Laplacian matrix of $\mathcal{D}$, and let $S\in\mathcal{P}(\mathcal{V})$. Then the following holds for all $j\in\{1,\ldots,n\}$:

[TABLE]

where $L_{j}$ is the $j$th row of $L$. Furthermore,

[TABLE]

{pf}

The term $\bm{\sigma}(S)$ is shortened to $\bm{\sigma}$ for brevity. Recall that the entry in the $j$th row and $i$th column of $L$ is denoted $L_{j,i}$. The definition of $L$ from (1) implies

[TABLE]

Since by (7), $q\in S$ implies $\sigma_{q}=1$, the term $\sum_{q\in\mathcal{N}_{j}\cap S}\sigma_{q}=|\mathcal{N}_{j}\cap S|$. In addition, since $q\notin S$ implies $\sigma_{q}=0$, the term $\sum_{q\in\mathcal{N}_{j}\backslash S}\sigma_{q}=0$. By this, equation (10) simplifies to $L_{j}\bm{\sigma}=|\mathcal{N}_{j}|\sigma_{j}-|\mathcal{N}_{j}\cap S|$.

The value of the term $|\mathcal{N}_{j}|\sigma_{j}$ depends on whether $j\in S$ or $j\notin S$. If $j\in S$, then $\sigma_{j}=1$, implying $L_{j}\bm{\sigma}=|\mathcal{N}_{j}|-|\mathcal{N}_{j}\cap S|=\left(|\mathcal{N}_{j}\cap S|+|\mathcal{N}_{j}\backslash S|\right)-|\mathcal{N}_{j}\cap S|=|\mathcal{N}_{j}\backslash S|$. If $j\notin S$, then $\sigma_{j}=0$ implying $L_{j}\bm{\sigma}=-|\mathcal{N}_{j}\cap S|$. This proves the result for equation (8).

To prove (9), we first consider nonempty sets $S\in\mathcal{P}(\mathcal{V})\backslash\{\emptyset\}$. By the results above and (4), the maximum reachability of any $S\in\mathcal{P}(\mathcal{V})\backslash\{\emptyset\}$ is

[TABLE]

By its definition, $\mathcal{R}(S)\geq 0$. Observe that if $j\in S$ then $L_{j}\bm{\sigma}(S)=|\mathcal{N}_{j}\backslash S|\geq 0$, implying $\max_{j\in S}L_{j}\bm{\sigma}(S)\geq 0$. Conversely, if an agent $j$ is not in the set $S$, then the function $L_{j}\bm{\sigma}(S)$ takes the nonpositive value $-|\mathcal{N}_{j}\cap S|$. This implies $\max_{j\notin S}L_{j}\bm{\sigma}(S)\leq 0$. By these arguments, we therefore have $\max_{j\notin S}L_{j}\bm{\sigma}(S)\leq 0\leq\max_{j\in S}L_{j}\bm{\sigma}(S)$, which implies

[TABLE]

Therefore by equations (12) and (11), the maximum reachability of $S$ is found by the expression

[TABLE]

Lastly, if $S=\emptyset$, then by (4) we have $\mathcal{R}(S)=0$. In addition, $\bm{\sigma}(S)=0$, implying that $\max_{j}L_{j}\bm{\sigma}(S)=0=\mathcal{R}(S),\ j\in\{1,\ldots,n\}$. $\blacksquare$

Using Lemma 2, it will next be shown that the objective function of (1) can be rewritten as the maximum over a set of affine functions:

Lemma 3**.**

Consider an arbitrary, nonempty, nontrivial, simple digraph $\mathcal{D}=(\mathcal{V},\mathcal{E})$. Let $L$ be the Laplacian matrix of $\mathcal{D}$, and let $L_{i}$ be the $i$th row of $L$. Let $\mathcal{T}$ be defined as in (2). Then for all $(S_{1},S_{2})\in\mathcal{T}$ the following holds:

[TABLE]

{pf}

By Lemma 2, $\mathcal{R}(S_{1})=\max_{i}L_{i}\bm{\sigma}(S_{1})$ and $\mathcal{R}(S_{2})=\max_{j}L_{j}\bm{\sigma}(S_{2})$ for $i,j\in\{1,\ldots,n\}$. The result follows. $\blacksquare$

From Lemma 1, Lemma 3, and Remark 1, we can immediately conclude that ${r_{\max}}(\mathcal{D})$ satisfies

[TABLE]

Note that the terms $\bm{\sigma}(S_{1})$ and $\bm{\sigma}(S_{2})$ are each $n$-dimensional binary vectors. Letting $\bm{b}^{1}=\bm{\sigma}(S_{1})$ and $\bm{b}^{2}=\bm{\sigma}(S_{2})$, the objective function of (4) can be written as \max\big{(}\max_{i}\left(L_{i}\bm{b}^{1}\right),\max_{j}\left(L_{j}\bm{b}^{2}\right)\big{)}. Every pair $(S_{1},S_{2})\in\mathcal{T}$ can be mapped into a pair of binary vectors $(\bm{b}^{1},\bm{b}^{2})$ by the function $\Sigma:\mathcal{T}\rightarrow\{0,1\}^{n}\times\{0,1\}^{n}$, where $\Sigma(S_{1},S_{2})=(\bm{\sigma}(S_{1}),\bm{\sigma}(S_{2}))=(\bm{b}^{1},\bm{b}^{2})$. By determining the image of $\mathcal{T}$ under $\Sigma(\cdot,\cdot)$, the optimal value of (4) can be found by minimizing over pairs of binary vectors $(\bm{b}_{1},\bm{b}_{2})\in\Sigma(\mathcal{T})$ directly. Using binary vector variables instead of set variables $(S_{1},S_{2})$ will allow (4) to be written directly in an MILP form. Towards this end, the following Lemma defines the set $\Sigma(\mathcal{T})$:

Lemma 4**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph, and let $\mathcal{T}$ be defined as in (2). Define the function $\Sigma:\mathcal{T}\rightarrow\{0,1\}^{n}\times\{0,1\}^{n}$ as

[TABLE]

Define the set $\mathcal{B}\subset\{0,1\}^{n}\times\{0,1\}^{n}$ as

[TABLE]

Then both of the following statements hold:

The image of $\mathcal{T}$ under $\Sigma$ is equal to $\mathcal{B}$, i.e. $\Sigma(\mathcal{T})=\mathcal{B}$ 2. 2.

The mapping $\Sigma:\mathcal{T}\rightarrow\mathcal{B}$ is a bijection.

{pf}

We prove 1) by showing first that $\Sigma(\mathcal{T})\subseteq\mathcal{B}$, and then $\mathcal{B}\subseteq\Sigma(\mathcal{T})$. Any $(S_{1},S_{2})\in\mathcal{T}$ satisfies $|S_{1}|>0$, $|S_{2}|>0$, $|S_{1}\cap S_{2}|=0$ as per (2). Observe that

[TABLE]

Because $S_{1},S_{2}\subset\mathcal{V}$ and $|S_{1}\cap S_{2}|=0$, then $|S_{1}|<n$. Otherwise if $|S_{1}|=n$ then either $|S_{2}|=0$ or $|S_{1}\cap S_{2}|\neq 0$, which both contradict the definition of $\mathcal{T}$. Therefore $|S_{1}|<n$, and by similar arguments $|S_{2}|<n$. Observe that

[TABLE]

Finally, $|S_{1}\cap S_{2}|=0$ implies that $j\in S_{1}\implies j\notin S_{2}$ and $j\in S_{2}\implies j\notin S_{1}$ $\forall j\in\{1,\ldots,n\}$. Therefore $\sigma_{j}(S_{1})=1\implies\sigma_{j}(S_{2})=0$ and $\sigma_{j}(S_{2})=1\implies\sigma_{j}(S_{1})=0$. This implies that

[TABLE]

Therefore for all $(S_{1},S_{2})\in\mathcal{T}$, $(\bm{\sigma}(S_{1}),\bm{\sigma}(S_{2}))=$ $\Sigma(S_{1},S_{2})$ satisfies the constraints of the set on the RHS of (16). This implies that $\Sigma(\mathcal{T})\subseteq\mathcal{B}$.

Next, we show $\mathcal{B}\subseteq\Sigma(\mathcal{T})$ by showing that for all $(\bm{b}^{1},\bm{b}^{2})\in\mathcal{B}$, there exists an $(S_{1},S_{2})\in\mathcal{T}$ such that $(\bm{b}^{1},\bm{b}^{2})=\Sigma(S_{1},S_{2})$. Choose any $(\bm{b}^{1},\bm{b}^{2})\in\mathcal{B}$ and define sets $(S_{1},S_{2})$ as follows:

[TABLE]

For the considered sets $(S_{1},S_{2})$, $1\leq\bm{1}^{T}\bm{b}^{1}$ implies $|S_{1}|>0$ and $1\leq\bm{1}^{T}\bm{b}^{2}$ implies $|S_{2}|>0$. In addition since $\bm{b}^{1}+\bm{b}^{2}\preceq\bm{1}$, we have $b^{1}_{j}=1\implies b^{2}_{j}=0$ and $b^{2}_{j}=1\implies b^{1}_{j}=0$. By our choice of $S_{1}$ and $S_{2}$, we have $b_{j}^{1}=1\implies j\in S_{1}$, and from previous arguments $b_{j}^{1}=1\implies b_{j}^{2}=0\implies j\notin S_{2}$. Similar reasoning can be used to show that $b_{j}^{2}=1\implies j\notin S_{1}$. These arguments imply that $|S_{1}\cap S_{2}|=0$. Consequently, $(S_{1},S_{2})$ satisfies all the constraints of $\mathcal{T}$ and is therefore an element of $\mathcal{T}$. Clearly, by (17) we have $\Sigma(S_{1},S_{2})=(\bm{b}^{1},\bm{b}^{2})$, which shows that there exists an $(S_{1},S_{2})\in\mathcal{T}$ such that $(\bm{b}^{1},\bm{b}^{2})=\Sigma(S_{1},S_{2})$. Since this holds for all $(\bm{b}^{1},\bm{b}^{2})\in\mathcal{B}$, this implies $\mathcal{B}\subseteq\Sigma(\mathcal{T})$. Therefore $\Sigma(\mathcal{T})=\mathcal{B}$.

We next prove 2). Since $\Sigma(\mathcal{T})=\mathcal{B}$, the function $\Sigma:\mathcal{T}\rightarrow\mathcal{B}$ is surjective. To show that it is injective, consider any $\Sigma(S_{1},S_{2})\in\mathcal{B}$ and $\Sigma(\bar{S}_{1},\bar{S}_{2})\in\mathcal{B}$ such that $\Sigma(S_{1},S_{2})=\Sigma(\bar{S}_{1},\bar{S}_{2})$. This implies $(\bm{\sigma}(S_{1}),\bm{\sigma}(S_{2}))=(\bm{\sigma}(\bar{S}_{1}),\bm{\sigma}(\bar{S}_{2}))$. Note that $(\bm{\sigma}(S_{1}),\bm{\sigma}(S_{2}))=(\bm{\sigma}(\bar{S}_{1}),\bm{\sigma}(\bar{S}_{2}))$ if and only if $\bm{\sigma}(S_{1})=\bm{\sigma}(\bar{S}_{1})$ and $\bm{\sigma}(S_{2})=\bm{\sigma}(\bar{S}_{2})$. Since the indicator function $\bm{\sigma}:\mathcal{P}(\mathcal{V})\rightarrow\{0,1\}^{n}$ is itself injective, this implies $S_{1}=\bar{S}_{1}$ and $S_{2}=\bar{S}_{2}$, which implies $(S_{1},S_{2})=(\bar{S}_{1},\bar{S}_{2})$. Therefore $\Sigma:\mathcal{T}\rightarrow\mathcal{B}$ is injective. $\blacksquare$

Using Lemma 4 allows us to present the following mixed integer program which solves Problem 1:

Theorem 1**.**

Let $\mathcal{D}$ be an arbitrary nonempty, nontrivial, simple digraph and let $L$ be the Laplacian matrix of $\mathcal{D}$. The maximum $r$-robustness of $\mathcal{D}$, denoted ${r_{\max}}(\mathcal{D})$, is obtained by solving the following minimization problem:

[TABLE]

Furthermore, (1) is equivalent to the following mixed integer linear program:

[TABLE]

{pf}

From Lemmas 1 and 3 we have

[TABLE]

for $i,j\in\{1,\ldots,n\}$. As per Remark 1, the definition of $\mathcal{T}$ can be used to make the constraints implicit:

[TABLE]

for $i,j\in\{1,\ldots,n\}$. Since $\Sigma:\mathcal{T}\rightarrow\mathcal{B}$ is a bijection by Lemma 4, (4) is equivalent to

[TABLE]

Making the constraints of (4) explicit yields (1).

Next, we prove that (19) is equivalent to (1). The variables $\bm{b}^{1}$ and $\bm{b}^{2}$ from (1) are combined into the variable $\bm{b}\in\mathbb{Z}^{2n}$ in (19); i.e. $\bm{b}=\begin{bmatrix}(\bm{b}^{1})^{T}&(\bm{b}^{2})^{T}\end{bmatrix}^{T}$. The first and third constraints of (19) restrict $\bm{b}\in\{0,1\}^{2n}$. Next, it can be demonstrated [3, Chapter 4] that the formulation ${\min_{\bm{x}}}\max_{i}(x_{i})$ is equivalent to $\min_{t,x}t$ subject to $0\leq t$, $\bm{x}\preceq t\bm{1}$. Reformulating the objective of the RHS of (1) in this way yields the objective and first two constraints of (19):

[TABLE]

The fourth, fifth, and sixth constraints of (19) restrict $(\bm{b}^{1},\bm{b}^{2})\in\mathcal{B}$ and are simply a reformulation of the first three constraints in (1). $\blacksquare$

5 Determining $(r,s)$-Robustness using Mixed Integer Linear Programming

In this section we address Problems 2 and 3, which involve determining the $(r^{*},s^{*})$-robustness and $(F_{\max}+1,F_{\max}+1)$-robustness of a given digraph. To determine these values, we will use the following notation:

Definition 7**.**

For a digraph $\mathcal{D}$ and a given $r\in\mathbb{Z}_{+}$, the maximum integer $s$ for which $\mathcal{D}$ is $(r,s)$-robust is denoted as ${s_{\max}}(r)\in\mathbb{Z}_{+}$. If $\mathcal{D}$ is not $(r,s)$-robust for any $1\leq s\leq n$, we will denote ${s_{\max}}(r)=0$.

Using this notation, the $(r^{*},s^{*})$-robustness of a nonempty, nontrivial simple digraph satisfies $r^{*}={r_{\max}}(\mathcal{D})$ and $s^{*}={s_{\max}}({r_{\max}}(\mathcal{D}))$. This can be verified by recalling that $r$-robustness is equivalent to $(r,1)$-robustness [13, Property 5.21], and that $(r^{*},s^{*})$ is the maximum element of $\Theta$ according to the lexicographic ordering defined in Section 2. Since a method for determining ${r_{\max}}(\mathcal{D})$ has already been presented, this section will introduce a method for determining ${s_{\max}}(r)$ for any given $r\in\mathbb{Z}_{+}$. This can then be used to find ${s_{\max}}({r_{\max}}(\mathcal{D}))$ after ${r_{\max}}(\mathcal{D})$ is determined.

Recall that $\lor$ indicates logical OR. An equivalent definition of ${s_{\max}}(r)$ can be given using the following notation:

Definition 8**.**

Let $\Theta$ be the set of all $(r,s)$ values for which a given digraph $\mathcal{D}$ is $(r,s)$-robust, as per (3). Let $r\in\mathbb{Z}_{+}$, and let $\mathcal{X}^{r}_{S}$ be defined as in Definition 4. The set $\Theta_{r}\subset\Theta$ is defined as follows:

[TABLE]

In words, $\Theta_{r}$ is the set of all integers $s$ for which the given digraph $\mathcal{D}$ is $(r,s)$-robust for a given $r\in\mathbb{Z}_{+}$. By this definition, ${s_{\max}}(r)=\max\Theta_{r}$, i.e. ${s_{\max}}(r)$ is simply the maximum element of $\Theta_{r}$.

As per (8), checking directly if an integer $s\in\Theta_{r}$ involves testing a logical disjunction for all possible $(S_{1},S_{2})\in\mathcal{T}$. This quickly becomes impractical for large $n$ since $|\mathcal{T}|$ grows exponentially with $n$. This difficulty can be circumvented, however, by defining the set

[TABLE]

where $\land$ denotes logical AND. The set $\bar{\Theta}_{r}$ contains all integers $\bar{s}$ for which the given digraph is not $(r,\bar{s})$-robust for the given value of $r$.

Definition 9**.**

For a digraph $\mathcal{D}$ and a given $r\in\mathbb{Z}_{+}$, the minimum integer $\bar{s}$ for which $\mathcal{D}$ is not $(r,\bar{s})$-robust is denoted as $\bar{s}_{\min}(r)\in\mathbb{Z}_{+}$.

As a simple example, consider a digraph $\mathcal{D}$ of 7 nodes where ${s_{\max}}(3)=2$. This implies that $\Theta_{3}=\{1,2\}$, i.e. the digraph is $(3,1)$- and $(3,2)$-robust. In this case, the set $\bar{\Theta}_{3}=\{3,4,5,\ldots\}$ since $\mathcal{D}$ is not $(3,3)$-robust, $(3,4)$-robust, etc. Here we have $\bar{s}_{\min}(3)=3$. In general, observe that by definitions 7 and 9 we have

[TABLE]

It is therefore sufficient to find $\bar{s}_{\min}(r)$ in order to determine ${s_{\max}}(r)$. An illustration is given in Figure 3. The methods in this section will solve for $\bar{s}_{\min}(r)$ using a mixed integer linear program. Note that since possible values of ${s_{\max}}(r)$ are limited to $0\leq{s_{\max}}(r)\leq n$ (Definition 7), possible values of $\bar{s}_{\min}(r)$ are limited to $1\leq\bar{s}_{\min}(r)\leq n+1$. We point out that it is easier to test if an integer $\bar{s}\in\bar{\Theta}_{r}$ than to test if an integer $s\in\Theta_{r}$, in the sense that only one element $(S_{1},S_{2})\in\mathcal{T}$ is required to verify that $\bar{s}\in\bar{\Theta}_{r}$ (as per (5)) whereas all $(S_{1},S_{2})\in\mathcal{T}$ must be checked to verify that $s\in\Theta_{r}$ (as per (8)).

The following Lemma is needed for our main result. It shows that given any $r\in\mathbb{Z}_{+}$ and $S\subset\mathcal{V}$, the indicator vector of the set $\mathcal{X}^{r}_{S}$, denoted $\bm{\sigma}(\mathcal{X}^{r}_{S})$, can be expressed using an MILP. Recall that $\mathcal{X}^{r}_{S}$ is the set of agents in $S$ which have $r$ in-neighbors outside of $S$, implying $\sigma_{j}(\mathcal{X}^{r}_{S})=1$ if $L_{j}\bm{\sigma}(S)=|\mathcal{N}_{j}\backslash S|\geq r$.

Lemma 5**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph. Let $L$ be the Laplacian matrix of $\mathcal{D}$, and let $r\in\mathbb{N}$. Consider any subset $S\subset\mathcal{V}$, $|S|>0$ and let $\mathcal{X}^{r}_{S}$ be defined as in Definition 4. Then the following holds:

[TABLE]

{pf}

Recall that the entries of the indicator vector $\bm{\sigma}(\mathcal{X}^{r}_{S})$ are defined as

[TABLE]

Let $\bm{y}^{*}$ be an optimal point of the RHS of (26). To prove that $\bm{y}^{*}=\bm{\sigma}(\mathcal{X}^{r}_{S})$, we demonstrate that $\forall j\in\{1,\ldots,n\}$, $\sigma_{j}(\mathcal{X}^{r}_{S})=1\iff y^{*}_{j}=1$. Observe that this is equivalent to demonstrating $\sigma_{j}(\mathcal{X}^{r}_{S})\neq 1\iff y^{*}_{j}\neq 1\ \forall j$, which is equivalent to demonstrating $\sigma_{j}(\mathcal{X}^{r}_{S})=0\iff y^{*}_{j}=0\ \forall j$. This can be seen by noting $\bm{\sigma}(\mathcal{X}^{r}_{S})\in\{0,1\}^{n}$ which implies $\sigma_{j}(\mathcal{X}^{r}_{S})\neq 1\iff\sigma_{j}(\mathcal{X}^{r}_{S})=0$, and $\bm{y}^{*}\in\{0,1\}^{n}$ which implies $y^{*}_{j}\neq 1\iff y^{*}_{j}=0$. Since proving $\sigma_{j}(\mathcal{X}^{r}_{S})=1\iff y^{*}_{j}=1$ for all $j\in\{1,\ldots,n\}$ is equivalent to proving $\sigma_{j}(\mathcal{X}^{r}_{S})=0\iff y^{*}_{j}=0$ for all $j\in\{1,\ldots,n\}$, and both $\bm{y}^{*}\in\{0,1\}^{n}$ and $\bm{\sigma}(\mathcal{X}^{r}_{S})\in\{0,1\}^{n}$, we therefore have $(\sigma_{j}(\mathcal{X}^{r}_{S})=1\iff y^{*}_{j}=1\ \forall j)$ if and only if $(\bm{y}^{*}=\bm{\sigma}(\mathcal{X}^{r}_{S}))$.

Sufficiency: Consider any $j\in\{1,\ldots,n\}$ such that $\sigma_{j}(\mathcal{X}^{r}_{S})=1$. This implies $j\in\mathcal{X}^{r}_{S}$ and therefore $|\mathcal{N}_{j}\backslash S|\geq r$ by Definition 4. By Lemma 2, $|\mathcal{N}_{j}\backslash S|=L_{j}\bm{\sigma}(S)$, and therefore $L_{j}\bm{\sigma}(S)>(r-1)$. Since $\bm{y}^{*}$ is an optimal point, it is therefore a feasible point. If $y^{*}_{j}=0$, the $j$th row of the first constraint on the RHS of (26) is be violated since $L_{j}\bm{\sigma}(S)-(n)y^{*}_{j}=L_{j}\bm{\sigma}(S)\nleq(r-1)$. Therefore we must have $y^{*}_{j}=1$. Note that $|\mathcal{N}_{j}\backslash S|\leq n\ \forall j\in S$ for any $S\subset\mathcal{V}$.

Necessity: We prove by contradiction. Suppose $y^{*}_{j}=1$ and $\sigma_{j}(\mathcal{X}^{r}_{S})=0$. This implies that $L_{j}\bm{\sigma}(S)=|\mathcal{N}_{j}\backslash S|<r$. Consider the vector $\tilde{\bm{y}}$ where $\tilde{y}_{j}=0$ and $\tilde{y}_{i}=y^{*}_{i}\ \forall i\neq j,\ i\in\{1,\ldots,n\}$. Since $L_{j}\bm{\sigma}(S)=|\mathcal{N}_{j}\backslash S|<r$, then $\tilde{\bm{y}}$ is therefore also a feasible point, and $\bm{1}^{T}\tilde{\bm{y}}<\bm{1}^{T}\bm{y}^{*}$. This contradicts $\bm{y}^{*}$ being an optimal point to (26); therefore we must have $\sigma_{j}(\mathcal{X}^{r}_{S})=1$. $\blacksquare$

The next Theorem presents a mixed integer linear program which determines $\bar{s}_{\min}(r)$ for any fixed $r\in\mathbb{Z}_{+}$.

Theorem 2**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph. Let $L$ be the Laplacian matrix of $\mathcal{D}$. Let $r\in\mathbb{Z}_{+}$, and let $\bar{s}_{\min}(r)$ be the minimum value of $s$ for which $\mathcal{D}$ is not $(r,s)$-robust. Then if $\bar{s}_{\min}(r)<n+1$, the following holds:

[TABLE]

Furthermore, for any $r>0$, $\bar{s}_{\min}(r)=n+1$ if and only if the integer program in (2) is infeasible.

{pf}

Note that the theorem statement assumes $r$ is a fixed integer. First, consider the case where $\bar{s}_{\min}(r)<(n+1)$. The value of $\bar{s}_{\min}$ can be found by solving the problem $\bar{s}_{\min}(r)=\min_{\bar{s}\in\bar{\Theta}_{r}}\bar{s}$. Making the constraints explicit yields

[TABLE]

To put this problem in an MILP form, we show that the terms $|\mathcal{X}^{r}_{S_{1}}|$, $|\mathcal{X}^{r}_{S_{2}}|$, $|S_{1}|$, and $|S_{2}|$ can be represented by functions of binary vectors. This can be done by first observing that for any $S\subseteq\mathcal{V}$, $|S|=\bm{1}^{T}\bm{\sigma}(S)$. Therefore the following relationships hold:

[TABLE]

Equation (29) can therefore be rewritten as

[TABLE]

By Lemma 4, the terms $\bm{\sigma}(S_{1}),\ \bm{\sigma}(S_{2})$ for $(S_{1},S_{2})\in\mathcal{T}$ can be represented by vectors $(\bm{b}^{1},\bm{b}^{2})\in\mathcal{B}$. This yields

[TABLE]

Expanding the last constraint using the definition of $\mathcal{B}$ in (16) yields the sixth, seventh, and eighth constraints in (2) as well as the constraint that $\bm{b}^{1},\bm{b}^{2}\in\{0,1\}^{n}$. In addition, the first constraint of the RHS of (2) limits the search for feasible value of $\bar{s}$ to the range of possible values for $\bar{s}_{\min}(r)$.

The vectors $\bm{y}^{1},\bm{y}^{2}$ are constrained to satisfy $\bm{y}^{1}=\bm{\sigma}(\mathcal{X}^{r}_{S_{1}})$ and $\bm{y}^{2}=\bm{\sigma}(\mathcal{X}^{r}_{S_{2}})$ as follows: by Lemma 4, $\bm{b}^{1}=\bm{\sigma}(S_{1})$ and $\bm{b}^{2}=\bm{\sigma}(S_{2})$ for $(S_{1},S_{2})\in\mathcal{T}$ as per the sixth through ninth constraints. Therefore by Lemma 5,

[TABLE]

The constraints of (32) and (33) are contained in the fifth and last constraints of (2). Since the fourth constraint of (2), $\bm{1}^{T}\bm{y}^{1}+\bm{1}^{T}\bm{y}^{2}\leq(\bar{s}-1)$, simultaneously minimizes $\bm{1}^{T}\bm{y}^{1}$ and $\bm{1}^{T}\bm{y}^{2}$, the fourth, fifth, and last constraints of (2) ensure that $\bm{y}^{1}=\bm{\sigma}(\mathcal{X}^{r}_{S_{1}})$ and $\bm{y}^{2}=\bm{\sigma}(\mathcal{X}^{r}_{S_{2}})$. Therefore $\bm{1}^{T}\bm{y}^{1}=|\mathcal{X}^{r}_{S_{1}}|$ and $\bm{1}^{T}\bm{y}^{2}=|\mathcal{X}^{r}_{S_{2}}|$.

Now, by the above arguments, solving the RHS of (2) yields $\bar{s}_{\min}(r)$ when $\bar{s}_{\min}(r)<(n+1)$. We now prove that for $r>0$, $\bar{s}_{\min}(r)=n+1$ if and only if the RHS of (2) is infeasible. Note that if $r=0$, then it trivially holds by Definition 4 that ${s_{\max}}(0)=n$ and therefore $\bar{s}_{\min}(0)=n+1$.

Sufficiency: $\bar{s}_{\min}(r)=n+1$ implies that ${s_{\max}}(r)=n$. Recall that $\mathcal{D}$ is $(r,{s_{\max}}(r))$-robust by Definition 7, since ${s_{\max}}(r)$ is the largest integer $s$ for which $\mathcal{D}$ is $(r,s)$-robust. By Definition 4, this implies that for all $(S_{1},S_{2})\in\mathcal{T}$, at least one of the following three conditions holds: $|\mathcal{X}^{r}_{S_{1}}|=|S_{1}|$, or $|\mathcal{X}^{r}_{S_{2}}|=|S_{2}|$, or $|\mathcal{X}^{r}_{S_{1}}|+|\mathcal{X}^{r}_{S_{2}}|\geq{s_{\max}}(r)=n$. Given any $(S_{1},S_{2})\in\mathcal{T}$, we consider each condition separately and show that at least one constraint of (2) is violated if the condition holds true:

• •

$|\mathcal{X}^{r}_{S_{1}}|=|S_{1}|$ being true implies that the second constraint of (2) is violated. This can be shown using earlier arguments from this proof. Specifically, we have $\bm{1}^{T}\bm{y}^{1}=|\mathcal{X}^{r}_{S_{1}}|=|S_{1}|=\bm{1}^{T}\bm{b}^{1}>\bm{1}^{T}\bm{b}^{1}-1$. Therefore no feasible point can be constructed from the given set pair $(S_{1},S_{2})$ if $|\mathcal{X}^{r}_{S_{1}}|=|S_{1}|$.

• •

$|\mathcal{X}^{r}_{S_{2}}|=|S_{2}|$ being true implies that the third constraint of (2) is violated. Specifically, we have $\bm{1}^{T}\bm{b}^{2}=|S_{2}|$ and $\bm{1}^{T}\bm{y}^{2}=|\mathcal{X}^{r}_{S_{2}}|$. This implies that $\bm{1}^{T}\bm{y}^{2}=|\mathcal{X}^{r}_{S_{2}}|=|S_{2}|=\bm{1}^{T}\bm{b}^{2}>\bm{1}^{T}\bm{b}^{2}-1$. Therefore no feasible point can be constructed from the given set pair $(S_{1},S_{2})$ if $|\mathcal{X}^{r}_{S_{2}}|=|S_{2}|$.

• •

$|\mathcal{X}^{r}_{S_{1}}|+|\mathcal{X}^{r}_{S_{2}}|\geq n$ being true implies that both $|\mathcal{X}^{r}_{S_{1}}|=|S_{1}|$ and $\mathcal{X}^{r}_{S_{2}}=|S_{2}|$. This follows by observing that $\mathcal{X}^{r}_{S_{1}}\subseteq S_{1}$ and $\mathcal{X}^{r}_{S_{2}}\subseteq S_{2}$, $S_{1}\cap S_{2}=\{\emptyset\}$ by definition of $\mathcal{T}$ in (2), and therefore $\mathcal{X}^{r}_{S_{1}}\cap\mathcal{X}^{r}_{S_{2}}=\{\emptyset\}$. Since $S_{1},S_{2}\subset\mathcal{V}$ and $|\mathcal{V}|=n$, we have $n\leq|\mathcal{X}^{r}_{S_{1}}|+|\mathcal{X}^{r}_{S_{2}}|\leq|S_{1}|+|S_{2}|\leq|\mathcal{V}|=n$. We must therefore have $|\mathcal{X}^{r}_{S_{1}}|=|S_{1}|$ and $|\mathcal{X}^{r}_{S_{2}}|=|S_{2}|$, which from prior arguments both imply that a constraint of (2) is violated. Therefore no feasible point can be constructed from the given set pair $(S_{1},S_{2})$ if $|\mathcal{X}^{r}_{S_{1}}|+|\mathcal{X}^{r}_{S_{2}}|\geq n$

Since for all $(S_{1},S_{2})\in\mathcal{T}$ at least one of these three conditions holds, for all $(S_{1},S_{2})\in\mathcal{T}$ at least one constraint of (2) is violated when ${s_{\max}}(r)=n$, which is equivalent to $\bar{s}_{\min}(r)=n+1$. Therefore $\bar{s}_{\min}(r)=n+1$ implies that (2) is infeasible.

Necessity: We prove the contrapositive, i.e. we prove that $\bar{s}_{\min}\neq n+1$ implies that there exists a feasible point to the RHS of (2). First, no digraph on $n$ nodes is $(r,n+1)$-robust [15, Definition 13], and the contrapositive of Property 1 implies that if a graph is not $(\bar{r},\bar{s})$-robust, then it is also not $(\bar{r}^{\prime},\bar{s}^{\prime})$-robust for all $\bar{r}^{\prime}\geq\bar{r}$, and for all $\bar{s}^{\prime}\geq\bar{s}$. Therefore $\bar{s}_{\min}\neq n+1$ implies $\bar{s}_{\min}\leq n$. Next, $\bar{s}_{\min}\leq n$ implies $n\in\bar{\Theta}_{r}$, which implies that there exists $(S_{1},S_{2})\in\mathcal{T}$ such that $|\mathcal{X}^{r}_{S_{1}}|\leq|S_{1}|-1$ and $|\mathcal{X}^{r}_{S_{2}}|\leq|S_{2}|-1$ and $|\mathcal{X}^{r}_{S_{1}}|+|\mathcal{X}^{r}_{S_{2}}|\leq n-1$, as per (5). Letting $\bar{s}=n$, $\bm{b}^{1}=\bm{\sigma}(S_{1})$, $\bm{b}^{2}=\bm{\sigma}(S_{2})$, $\bm{y}^{1}=\bm{\sigma}(\mathcal{X}^{r}_{S_{1}})$, and $\bm{y}^{2}=\bm{\sigma}(\mathcal{X}^{r}_{S_{2}})$ yields a feasible point to (2). $\blacksquare$

The MILPs in Theorems 1 and 2 can be used to determine the $(r^{*},s^{*})$-robustness of any digraph satisfying Assumption 1, thereby solving Problem 2. Recall from the beginning of Section 5 that $r^{*}={r_{\max}}(\mathcal{D})$ and $s^{*}={s_{\max}}({r_{\max}}(\mathcal{D}))$. Theorem 1 can first be used to determine the value of ${r_{\max}}(\mathcal{D})=r^{*}$. Using ${r_{\max}}(\mathcal{D})$, Theorem 2 can then be used to find the value of ${s_{\max}}({r_{\max}}(\mathcal{D}))=s^{*}$.

More generally however, the MILP formulation in Theorem 2 allows for ${s_{\max}}(r)$ to be determined for any $r\in\mathbb{Z}_{+}$. Since $(r,1)$-robustness is equivalent to $r$-robustness, the MILP in Theorem 2 can also be used to determine whether a digraph $\mathcal{D}$ is $r$ robust for a given $r\in\mathbb{Z}_{+}$. If $\bar{s}_{\min}(r)\geq 2$, then ${s_{\max}}(r)\geq 1$ which implies that $\mathcal{D}$ is $(r,1)$-robust. On the other hand, $\bar{s}_{\min}(r)=1$ implies that $1\in\bar{\Theta}_{r}$ and therefore $\mathcal{D}$ is not $(r,1)$-robust (and not $r$-robust).

Finally, to solve Problem 3 Theorem 2 can be used to determine the $(F_{\max}+1,F_{\max}+1)$-robustness of a nonempty, nontrivial, simple digraph. Recall that $F_{\max}=\max(\{F\in\mathbb{Z}_{+}:(F+1,F+1)\in\Theta\})$. The value of $F_{\max}$ is determined by Algorithm 1, presented below.

In essence, Algorithm 1 finds the largest values of $r^{\prime}$ and $s^{\prime}$ such that $r^{\prime}=s^{\prime}$ and $(r^{\prime},s^{\prime})\in\Theta$. It begins by setting $r^{\prime}\leftarrow{r_{\max}}(\mathcal{D})$, and finding ${s_{\max}}(r^{\prime})$ using Theorem 2. If ${s_{\max}}(r^{\prime})\geq r^{\prime}$, then by Proposition 1 the digraph $\mathcal{D}$ is $(r^{\prime},s)$-robust for $s=r^{\prime}$ and therefore $(r^{\prime},r^{\prime})$-robust. This implies $r^{\prime}=F_{\max}+1$. However, if ${s_{\max}}(r^{\prime})<r^{\prime}$ then $r^{\prime}$ is decremented, ${s_{\max}}(r^{\prime})$ recalculated, and the process is repeated until the algorithm terminates with the highest integer $r^{\prime}$ such that $\mathcal{D}$ is $(r^{\prime},r^{\prime})$-robust, yielding $F_{\max}=(r^{\prime}-1)$.

6 Reducing Worst-Case Performance of Determining $\bf{r_{\max}}(\mathcal{D})$

When solving a MILP with zero-one integer variables using a branch-and-bound technique, the worst-case number of subproblems to be solved is equal to $2^{q}$, where $q$ is the dimension of the zero-one integer vector variable. In Section 4, the MILP in Theorem 1 which solves for ${r_{\max}}(\mathcal{D})$ has a binary vector variable with dimension $2n$. In this section, we present two MILPs whose optimal values provide upper and lower bounds on the value of ${r_{\max}}(\mathcal{D})$. Each MILP has a binary vector variable with dimension of only $n$, which implies a reduced worst-case performance as compared to the MILP in Theorem 1.

6.1 A Lower Bound on Maximum $r$-Robustness

In [6], a technique is presented for lower bounding ${r_{\max}}(\mathcal{D})$ of undirected graphs by searching for the minimum reachability of subsets $S\subset\mathcal{V}$ such that $|S|\leq\left\lfloor n/2\right\rfloor$. We extend this result to digraphs in the next Lemma.

Lemma 6**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph. Let $\Psi=\{S\subset\mathcal{V}:1\leq|S|\leq\left\lfloor n/2\right\rfloor\}$. Then the following holds:

[TABLE]

{pf}

By Lemma 1 and Remark 1, proving (34) is equivalent to proving

[TABLE]

Denote $S^{*}=\arg\min_{S\in\Psi}\mathcal{R}(S)$ and $(S_{1}^{*},S_{2}^{*})=\arg\min_{(S_{1},S_{2})\in\mathcal{T}}\max(\mathcal{R}(S_{1}),\mathcal{R}(S_{2}))$. We prove by contradiction. Suppose $\mathcal{R}(S^{*})>\max(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*}))$. Since $S_{1}^{*}$ and $S_{2}^{*}$ are nonempty, $|S_{1}^{*}|\geq 1$ and $|S_{2}^{*}|\geq 1$. Since they are disjoint, we must have either $|S_{1}^{*}|\leq\left\lfloor n/2\right\rfloor$, or $|S_{2}^{*}|\leq\left\lfloor n/2\right\rfloor$, or both $|S_{1}^{*}|$ and $|S_{2}^{*}|$ less than or equal to $\left\lfloor n/2\right\rfloor$. Therefore either $S_{1}^{*}\in\Psi$ or $S_{2}^{*}\in\Psi$. This implies that either $\mathcal{R}(S_{1}^{*})\geq\mathcal{R}(S^{*})$ (if $S_{1}^{*}\in\Psi$) or $\mathcal{R}(S_{2}^{*})\geq\mathcal{R}(S^{*})$ (if $S_{2}^{*}\in\Psi$), since $S^{*}$ is an optimal point. But this contradicts the assumption that $\mathcal{R}(S^{*})>\max(\mathcal{R}(S_{1}^{*}),\mathcal{R}(S_{2}^{*}))$. Therefore we must have

[TABLE]

which concludes the proof. $\blacksquare$

Using this result, a lower bound on ${r_{\max}}(\mathcal{D})$ can be obtained by the following optimization problem:

Theorem 3**.**

Let $\mathcal{D}$ be an arbitrary nonempty, nontrivial, simple digraph and let $L$ be the Laplacian matrix of $\mathcal{D}$. A lower bound on the maximum integer for which $\mathcal{D}$ is $r$-robust, denoted ${r_{\max}}(\mathcal{D})$, is found as follows:

[TABLE]

Furthermore, (37) is equivalent to the following mixed integer linear program:

[TABLE]

{pf}

To prove the result we show that the RHS of (37) is equivalent to the RHS of (34). By Lemma 2, $\mathcal{R}(S)=\max_{i}L_{i}\bm{\sigma}(S)$. Therefore (34) is equivalent to

[TABLE]

Next, we demonstrate that the set

[TABLE]

satisfies $\mathcal{B}_{\Psi}=\bm{\sigma}(\Psi)$, where $\bm{\sigma}(\Psi)$ is the image of $\Psi$ under $\bm{\sigma}:\mathcal{P}(\mathcal{V})\rightarrow\{0,1\}^{n}$. Since $1\leq|S|\leq\left\lfloor n/2\right\rfloor$ for all $S\in\Psi$, then by (7) we have $1\leq\bm{1}^{T}\bm{\sigma}(S)\leq\left\lfloor n/2\right\rfloor$ for all $S\in\Psi$. Also, $\bm{\sigma}(S)\in\{0,1\}^{n}$, and therefore $\bm{\sigma}(S)\in\mathcal{B}_{\Psi}\ \forall S\in\Psi$, implying that $\bm{\sigma}(\Psi)\subseteq\mathcal{B}$. Next, for any $\bm{b}\in\mathcal{B}_{\Psi}$, choose the set $S=\bm{\sigma}^{-1}(\bm{b})$ (recall from LABEL:subsec:reformulate that $\bm{\sigma}^{-1}:\{0,1\}^{n}\rightarrow\mathcal{P}(\mathcal{V})$). Then clearly $\bm{\sigma}(S)=\bm{\sigma}(\bm{\sigma}^{-1}(\bm{b}))=\bm{b}$, and therefore $\mathcal{B}_{\Psi}\subseteq\bm{\sigma}(\Psi)$. Therefore $\mathcal{B}_{\Psi}=\bm{\sigma}(\Psi)$.

The function $\bm{\sigma}:\Psi\rightarrow\mathcal{B}_{\Psi}$ is therefore surjective. Since $\bm{\sigma}:\mathcal{P}(\mathcal{V})\rightarrow\{0,1\}^{n}$ is injective, $\Psi\subset\mathcal{P}(\mathcal{V})$, and $\mathcal{B}_{\Psi}\subset\{0,1\}^{n}$, then $\bm{\sigma}:\Psi\rightarrow\mathcal{B}_{\Psi}$ is also injective and therefore a bijection. This implies that (39) is equivalent to

[TABLE]

Making the constraints of (41) explicit yields (37). More specifically, since $\mathcal{B}_{\Psi}=\{\bm{b}\in\{0,1\}^{n}:1\leq\bm{1}^{T}\bm{b}\leq\left\lfloor n/2\right\rfloor\}$ by (40), equation (41) can be rewritten with explicit constraints on $\bm{b}$ as follows:

[TABLE]

Equation (42) is the same as (37).

We next prove that (38) is equivalent to (37). As per the proof of Theorem 1, the objective and first two constraints of (38) are a reformulation of the objective of the RHS of (37). The first and third constraint ensure $b$ to be in $\{0,1\}^{n}$, and the fourth constraint ensures $\bm{b}\in\mathcal{B}_{\Psi}$. $\blacksquare$

6.2 An Upper Bound on Maximum r-Robustness

This section will present an MILP whose solution provides an upper bound on the value of ${r_{\max}}(\mathcal{D})$, and whose binary vector variable has a dimension of $n$. This will be accomplished by searching a subset $\mathcal{T}^{\prime}\subset\mathcal{T}$ which is defined as

[TABLE]

In other words, $\mathcal{T}^{\prime}$ is the set of all possible partitionings of $\mathcal{V}$ into $S_{1}$ and $S_{2}$. Considering only elements of $\mathcal{T}^{\prime}$ yields certain properties that allow us to calculate an upper bound on ${r_{\max}}(\mathcal{D})$ using an MILP with only an $n$-dimensional binary vector.

Observe that $|\mathcal{T}^{\prime}|=2^{n}-2$, since neither $\mathcal{T}$ nor $\mathcal{T}^{\prime}$ include the cases where $S_{1}=\{\emptyset\}$ or $S_{2}=\{\emptyset\}$. Similar to the methods discussed earlier, the partitioning of $\mathcal{V}$ into $S_{1}$ and $S_{2}$ can be represented by the indicator vectors $\bm{\sigma}(S_{1})$ and $\bm{\sigma}(S_{2})$, respectively. Note that since $S_{1}\cup S_{2}=\mathcal{V}$ for all $(S_{1},S_{2})\in\mathcal{T}^{\prime}$, it can be shown that ${\bm{\sigma}(S_{1})+\bm{\sigma}(S_{2})=\bm{1}}\ \forall(S_{1},S_{2})\in\mathcal{T}^{\prime}$. These properties allow the following Lemma to be proven:

Lemma 7**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph. Let $L$ be the Laplacian matrix of $\mathcal{D}$ and let $L_{j}$ be the $j$th row of $L$. Let $\mathcal{T}^{\prime}$ be defined as in (43). Then for all $(S_{1},S_{2})\in\mathcal{T}^{\prime}$, the following holds:

[TABLE]

{pf}

Lemma 2 implies that $L_{j}\bm{\sigma}(S_{1})=|\mathcal{N}_{j}\backslash S_{1}|$ if $j\in S_{1}$, and $L_{j}\bm{\sigma}(S_{2})=|\mathcal{N}_{j}\backslash S_{2}|$ if $j\in S_{2}$. Since $(S_{1},S_{2})\in\mathcal{T}^{\prime}\implies\bm{\sigma}(S_{1})+\bm{\sigma}(S_{2})=\bm{1}$, we have

[TABLE]

This relation holds because, by the definition of $L$, $\bm{1}$ is always in the null space of $L$. Therefore for $j\in S_{2}$ we have $L_{j}\bm{\sigma}(S_{1})=-L_{j}\bm{\sigma}(S_{2})=-|\mathcal{N}_{j}\backslash S_{2}|$, and for $j\in S_{1}$ we have $L_{j}\bm{\sigma}(S_{2})=-L_{j}\bm{\sigma}(S_{1})=-|\mathcal{N}_{j}\backslash S_{1}|$. $\blacksquare$

An interesting result of Lemma 7 is that for any subsets $(S_{1},S_{2})\in\mathcal{T}^{\prime}$, the maximum reachability of the two subsets can be recovered using the infinity norm.

Lemma 8**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph and let $L$ be the Laplacian matrix of $\mathcal{D}$. For all $(S_{1},S_{2})\in\mathcal{T}^{\prime}$, the following holds:

[TABLE]

{pf}

Denote the nodes in $S_{1}$ as $\{i_{1},\ldots,i_{p}\}$ and the nodes in $S_{2}$ as $\{j_{1},\ldots,j_{(n-p)}\}$ with $p\in\mathbb{Z},\ 1\leq p\leq(n-1)$. Note that since $S_{1}\cup S_{2}=\mathcal{V}$, we have $\{i_{1},\ldots,i_{p}\}\cup\{j_{1},\ldots,j_{n-p}\}=\{1,\ldots,n\}$.

The right hand side of equation (46) can be expressed as

[TABLE]

Similarly, using Lemma 7 yields

[TABLE]

Finally, observe that

[TABLE]

which completes the proof. $\blacksquare$

A mixed integer linear program yielding an upper bound on the value of ${r_{\max}}(\mathcal{D})$ is therefore given by the following Theorem:

Theorem 4**.**

Let $\mathcal{D}=(\mathcal{V},\mathcal{E})$ be an arbitrary nonempty, nontrivial, simple digraph. Let $L$ be the Laplacian matrix of $\mathcal{D}$. The maximum integer for which $\mathcal{D}$ is $r$-robust, denoted ${r_{\max}}(\mathcal{D})$, is upper bounded as follows:

[TABLE]

Furthermore, (50) is equivalent to the following mixed integer linear program:

[TABLE]

{pf}

Consider the optimization problem

[TABLE]

Since $\mathcal{T}^{\prime}\subset\mathcal{T}$, the optimal value of (52) is a valid upper bound on the value of ${r_{\max}}(\mathcal{D})$ as per Remark 1. From (52) and Lemma 8 we obtain

[TABLE]

Since $S_{1},S_{2}$ are nonempty and $S_{1}\cup S_{2}=\mathcal{V}$ for all $(S_{1},S_{2})\in\mathcal{T}^{\prime}$, the set of all possible $S_{1}$ subsets within elements of $\mathcal{T}^{\prime}$ is $(\mathcal{P}(\mathcal{V})\backslash\{\emptyset,\mathcal{V}\})$. Similarly, the set of all possible $S_{2}$ subsets within elements of $\mathcal{T}^{\prime}$ is also $(\mathcal{P}(\mathcal{V})\backslash\{\emptyset,\mathcal{V}\})$. For brevity, denote $\mathcal{P}_{\emptyset,\mathcal{V}}=\mathcal{P}(\mathcal{V})\backslash\{\emptyset,\mathcal{V}\}$.

Next, we demonstrate that the set $\mathcal{B}^{\prime}=\{\bm{b}\in\{0,1\}^{n}:1\leq\bm{1}^{T}\bm{b}\leq(n-1)\}$ satisfies $\bm{\sigma}(\mathcal{P}_{\emptyset,\mathcal{V}})=\mathcal{B}^{\prime}$. Since $1\leq|S|\leq(n-1)$ for all $S\in\mathcal{P}_{\emptyset,\mathcal{V}}$, then by (7) we have $1\leq\bm{1}^{T}\bm{\sigma}(S)\leq(n-1)$ for all $S\in\mathcal{P}_{\emptyset,\mathcal{V}}$. Also, $\bm{\sigma}(S)\in\{0,1\}^{n}$, and therefore $\bm{\sigma}(S)\in\mathcal{B}^{\prime}\ \forall S\in\mathcal{P}_{\emptyset,\mathcal{V}}$, implying that $\bm{\sigma}(\mathcal{P}_{\emptyset,\mathcal{V}})\subseteq\mathcal{B}^{\prime}$. Next, for any $\bm{b}\in\mathcal{B}^{\prime}$, choose the set $S=\bm{\sigma}^{-1}(\bm{b})$. Then clearly $\bm{\sigma}(S)=\bm{\sigma}(\bm{\sigma}^{-1}(\bm{b}))=\bm{b}$, and therefore $\mathcal{B}^{\prime}\subseteq\bm{\sigma}(\mathcal{P}_{\emptyset,\mathcal{V}})$. Therefore $\mathcal{B}^{\prime}=\bm{\sigma}(\mathcal{P}_{\emptyset,\mathcal{V}})$.

The function $\bm{\sigma}:\mathcal{P}_{\emptyset,\mathcal{V}}\rightarrow\mathcal{B}^{\prime}$ is therefore surjective. Since $\bm{\sigma}:\mathcal{P}(\mathcal{V})\rightarrow\{0,1\}^{n}$ is injective, $\mathcal{P}_{\emptyset,\mathcal{V}}\subset\mathcal{P}(\mathcal{V})$, and $\mathcal{B}^{\prime}\subset\{0,1\}^{n}$, then $\bm{\sigma}:\mathcal{P}_{\emptyset,\mathcal{V}}\rightarrow\mathcal{B}^{\prime}$ is also injective. Therefore $\bm{\sigma}:\mathcal{P}_{\emptyset,\mathcal{V}}\rightarrow\mathcal{B}^{\prime}$ is a bijection, implying that (53) is equivalent to

[TABLE]

Making the constraints of (54) explicit yields (50). We next prove that (51) is equivalent to (50). It can be shown [3, Chapter 4] that $\min_{\bm{x}}\left\lVert\bm{x}\right\rVert_{\infty}$ is equivalent to

[TABLE]

Likewise, the objective and first two constraints of the RHS of (51) are a reformulation of the objective of the RHS of (50). The first and third constraint restrict $\bm{b}\in\{0,1\}^{n}$, and the fourth constraint restricts $\bm{b}$ to be an element of $\mathcal{B}^{\prime}$. These arguments imply that the RHS of (51) is equivalent to the RHS of (50). $\blacksquare$

7 Discussion

MILP problems are $NP$-hard problems to solve in general. As such, the formulations presented in this paper do not reduce the theoretical complexity of the robustness determination problem. However, it has been pointed out that algorithmic advances and improvement in computer hardware have led to a speedup factor of 800 billion for mixed integer optimization problems during the last 25 years [1]. The results of this paper allow for the robustness determination problem to benefit from ongoing and future improvements in the active areas of optimization and integer programming.

In addition, one crucial advantage of the MILP formulations is the ability to iteratively tighten a global lower bound on the optimal value over time by using a branch-and-bound algorithm. The reader is referred to [41] for a concise overview of how such a lower bound can be calculated. In context of robustness determination, lower bounds on ${r_{\max}}(\mathcal{D})$ and ${s_{\max}}(r)$ are generally more useful than upper bounds since they can be used to calculate lower bounds on the maximum adversary model that the network can tolerate. The ability to use branch-and-bound algorithms for solving the robustness determination problem offers the flexibility of terminating the search for ${r_{\max}}(\mathcal{D})$ and/or ${s_{\max}}(r)$ when sufficiently high lower bounds have been determined. In this manner, approximations of these values can be found when it is too computationally expensive to solve for them exactly. The investigation of additional methods to find approximate solutions to the MILP formulations in this paper is left for future work.

8 Simulations

This section presents simulations which demonstrate the performance of the MILP formulations as compared to a robustness determination algorithm from prior literature called $DetermineRobustness$ [13]. Computations for these simulations are performed in MATLAB 2018b on a desktop computer with 8 Intel core i7-7820X CPUs (3.60 GHz) capable of handling 16 total threads. All MILP problems are solved using MATLAB’s intlinprog function.

Four types of random graphs are considered in the simulations: Erdős-Rényi random graphs, random digraphs, $k$-out random graphs [2], and $k$-in random graphs. The Erdős-Rényi random graphs in these simulations consist of $n$ agents, with each possible undirected edge present independently with probability $p\in[0,1]$ and absent with probability $1-p$. Three values of $p$ are considered: 0.3, 0.5, and 0.8. The random digraphs considered consist of $n$ nodes with each possible directed edge present independently with probability $p$ and absent with probability $1-p$. The values of $p$ considered are again 0.3, 0.5, and 0.8. The $k$-out random graphs consist of $n$ nodes. For each vertex $i\in\mathcal{V}$, $k\in\mathbb{N}$ other nodes are chosen with all having equal probability and all choices being independent. For each node $j$ of the $k$ chosen nodes, a directed edge $(i,j)$ is formed. $k$-in random graphs are formed in the same manner as $k$-out random graphs with the exception that the direction of the directed edges are reversed; i.e. edges $(j,i)$ are formed. The values of $k$ considered are $\{3,4,5\}$.

Two sets of simulations are considered. The first set compares two algorithms which determine the pair $(r^{*},s^{*})$ for a digraph: the $DetermineRobustness$ algorithm from [13] and Algorithm 3, $(r,s)$-Rob. MILP, which is an MILP formulation using results from Theorems 1 and 2. Details about the implementation of these algorithms can be found in the Appendix, section A. The algorithms are tested on the four types of graphs described above with values of $n$ ranging from $7$ to $15$. In addition, the MILP formulation is tested on digraphs with values of $n$ ranging from $17$ to $25$. The $DetermineRobustness$ algorithm is not tested on values of $n$ above 15 since the convergence rate trend is clear from the existing data, and the projected convergence times are prohibitive for large $n$. 100 graphs per graph type and combination of $n$ and $p$ (or $n$ and $k$ depending on the respective graph type), are randomly generated, and the algorithms are run on each graph. Overall, 10,800 total graphs are analyzed with $DetermineRobustness$ and 16,800 total graphs are analyzed with Algorithm 3. The time for each algorithm to determine the pair $(r^{*},s^{*})$ is averaged for each combination of $n$ and $p$ (for Erdős-Rényi random graphs and random digraphs), and for each combination of $n$ and $k$ (for $k$-out and $k$-in random graphs). The interpolated circles represent the average convergence time in seconds over 100 trials for each value of $n$, while the vertical lines represent the spread between maximum convergence time and minimum convergence time over trials for the respective value of $n$. Note the logarithmic scale of the y-axis.

To facilitate the large number of graphs being tested, a time limit of $10^{3}$ seconds (roughly 17 minutes) is imposed on Algorithm 3 (the MILP formulation). However, out of the 16,800 graphs tested by Algorithm 3, there are only 62 instances where the algorithm did not converge to optimality before this time limit. Instances where the time limit was violated are given the maximum time of $10^{3}$ seconds and included in the data. The graphs where optimality was not reached by the time limit all have between 21 and 25 nodes, are either Erdős-Rényi random graphs or random digraphs, and have edge formation probabilities of $p=0.8$.

Several patterns in the data warrant discussion. It is clear that in some cases, the minimum time of $DetermineRobustness$ is less than that of $(r,s)$-Rob. MILP. $DetermineRobustness$ terminates if two subsets $S_{1}$ and $S_{2}$ which are both 0-reachable are encountered, since this implies that the graph is at most $(0,n)$-robust. This can result in fast termination if such subsets are encountered early in the search. Second, there are instances where the maximum time for the $(r,s)$-Rob. MILP is much higher than that of the $DetermineRobustness$ algorithm (e.g. for $k$-out random digraphs with $k=4$). It is not immediately clear why this is the case; future work will investigate graph characteristics which affect the convergence time of the MILP formulations. Finally, for small values of $n$ (e.g. $n\in\{7,8\}$) the average time for $DetermineRobustness$ is lower than the average time for the $(r,s)$-Rob. MILP. This likely reflects that it may be quicker to simply test all unique nonempty, disjoint subsets in these cases (966 for $n=7$, 3025 for $n=8$) than to incur computational overhead associated with solving the MILP formulations. We point out that, with a few exceptions, the difference in this case is small: the average convergence time for both algorithms is generally under $10^{-1}$ seconds for $n\in\{7,8\}$.

The second simulation set compares the performance of four algorithms which determine only the value of ${r_{\max}}(\mathcal{D})$ for digraphs. These include Algorithm 4, a modified version of $DetermineRobustness$ which determines ${r_{\max}}(\mathcal{D})$, the MILP formulation from Theorem 1 (denoted r-Rob. MILP), the lower bound MILP formulation from Theorem 3 (denoted r-Rob. Lower Bnd), and the upper bound MILP formulation from Theorem 4 (denoted r-Rob. Upper Bnd). These algorithms are tested on the four types of graphs described above with values of $n$ ranging from $7$ to $15$. Additionally, the MILP formulations are tested on digraphs with values of $n$ ranging from $17$ to $25$. Again, 100 graphs per graph type and combination of $n$ and $p$ (or $n$ and $k$, depending on the respective graph type) are randomly generated, and the algorithms are run on each graph. Overall, 10,800 graphs are analyzed with Algorithm 4 and 16,800 graphs are analyzed by each of the three MILP formulations. The time for each algorithm to determine ${r_{\max}}(\mathcal{D})$ is averaged for each combination of $n$, $p$ or $k$, and graph type. The average, minimum, and maximum times per combination are plotted in Figure 5. A time limit of $10^{3}$ seconds is again imposed on all three of the MILP formulations, but out of the 16,800 graphs tested there are no instances where this time limit was violated.

Some of the same patterns as in the first set of simulations (with the $DetermineRobustness$ and $(r,s)$-Rob. MILP algorithms) are evident in the second set of simulations. The Mod. Det. Rob. algorithm also terminates if a pair of subsets $S_{1}$ and $S_{2}$ are found which are both 0-reachable, which is likely the reason for the small minimum computation time of this algorithm for several of the graphs. Mod. Det. Rob. generally has a lower average computational time for $n\in\{7,8\}$, again likely due to the speed of checking the relatively low number of unique nonempty, disjoint subset pairs as compared to solving the MILPs. It is not clear why the $r$-Rob. Upper Bnd MILP exhibits high average and maximum computational times for the $k$-out random digraphs. Future work will further analyze graph characteristics which negatively affect the convergence time of the MILP formulations.

9 Conclusion

This paper presents a novel approach on determining the $r$- and $(r,s)$-robustness of digraphs using mixed integer linear programming (MILP). The advantages of the MILP formulations and branch-and-bound algorithms over prior algorithms are discussed, and the performance of the MILP methods to the $DetermineRobustness$ algorithm is compared.

Much work remains to be done in the area of robustness determination. The results in this paper merely open the door for the extensive literature on mixed integer programming to be applied to the robustness determination problem. Future work will focus on applying more advanced integer programming techniques to the formulations in this paper to yield faster solution times. In particular, the Laplacian matrix exhibits a high degree of structure and plays a central role in the MILP formulations presented in this paper. Future efforts will explore ways to exploit this structure to determine the robustness of digraphs more efficiently.

Appendix A Description of Algorithm implementations

This section gives additional details about the implementations of the algorithms tested in the Simulation. Algorithm 2 provides the details about the implementation of the $DetermineRobustness$ algorithm implemented in the simulations.

One modification was made to the $DetermineRobustness$ algorithm to ensure accuracy of results. In the first line of the original $DetermineRobustness$ algorithm, $r$ is initialized with $\min\left(\delta^{\text{in}}(\mathcal{D}),\left\lceil n/2\right\rceil\right)$. This yields incorrect results however for directed spanning trees where the in-degree of the root node is zero. Consider the digraph depicted in Figure 6. Here, $\delta^{\text{in}}(\mathcal{D})=0$, since the left agent has no in-neighbors. This implies that the original $DetermineRobustness$ algorithm would initialize $r\leftarrow 0$ and return $(0,n)$ as the values of $({r_{\max}}(\mathcal{D}),{s_{\max}}({r_{\max}}(\mathcal{D})))$. However, Figure 6 shows that for all nonempty, disjoint subsets $S_{1}$ and $S_{2}$, at least one is $1$-reachable. The depicted graph is therefore $(1,1)$-robust with ${r_{\max}}(\mathcal{D})=1$ and ${s_{\max}}({r_{\max}}(\mathcal{D}))=1$. In fact, initializing $r\leftarrow\delta^{\text{in}}(\mathcal{D})$ will always yield this error for any directed spanning tree where the in-degree of the root node is 0. This happens because any digraph is 1-robust if and only if it contains a rooted out-branching [15, Lemma 7], yet $DetermineRobustness$ initializes $r\leftarrow\delta^{\text{in}}(\mathcal{D})=0$ which results in termination at line 23.

In Algorithm 2, $r$ is instead initialized with $\min\left(\max\left(\delta^{\text{in}}(\mathcal{D}),1\right),\left\lceil\frac{n}{2}\right\rceil\right)$. This initializes $r\leftarrow 1$ if $\delta^{\text{in}}(\mathcal{D})=0$, since it is still possible for the digraph to be $1$-robust.

To compare the MILP methodologies of this paper with the $DetermineRobustness$ algorithm, Algorithm 3 was used which determines the values of $({r_{\max}}(\mathcal{D}),{s_{\max}}({r_{\max}}(\mathcal{D})))$.

The first part of Algorithm 3 uses the formulation in Theorem 1 to determine ${r_{\max}}(\mathcal{D})$. If ${r_{\max}}(\mathcal{D})=0$, then $s\leftarrow n$ and the algorithm returns $(r,s)$. If ${r_{\max}}(\mathcal{D})>0$, then the algorithm determines in line 5 whether the minimum in-degree of the graph $\delta^{\text{in}}(\mathcal{D})\geq\left\lfloor n/2\right\rfloor+r-1$. By Property 5.23 in [10], if this is satisfied then the digraph is $(r,s)$-robust for all $1\leq s\leq n$. Since for $r>0$ the MILP in Theorem 2 is infeasible if and only if $s=n$, this test attempts to help the MILP solver avoid a fruitless search for a feasible solution if $s$ is indeed equal to $n$. This test works for graphs with large minimum in-degrees (e.g. complete graphs), but since it is a sufficient condition only it may not always detect when ${s_{\max}}(r)=n$. Determining a more rigorous test to determine when ${s_{\max}}(r)=n$ is left for future work. Finally, if the test in line 5 fails then the MILP formulation in Theorem 2 is performed to determine $\bar{s}_{\min}$. Since ${s_{\max}}(r)=\bar{s}_{\min}(r)-1$, the value of ${s_{\max}}(r)$ is stored in $s$ and the algorithm returns $(r,s)$.

The MILP algorithms in Section 4 consider $r$-robustness, which is equivalent to $(r,1)$-robustness [13, Property 5.21], [15, Section VII-B]. Since they effectively do not consider values of $s$ greater than $1$, it is unfair to compare them directly with the $DetermineRobustness$ algorithm. Algorithm 4 is a modified version of $DetermineRobustness$ which only considers $(r,1)$-robustness. This is accomplished by initializing $s\leftarrow 1$ in lines 2 and 20 instead of $s\leftarrow n$. Algorithm 4 is labeled Mod. Det. Rob. in the simulation legends.