Randomized Constraints Consensus for Distributed Robust Mixed-Integer Programming

TL;DR

This paper introduces a randomized distributed algorithm for solving uncertain mixed-integer convex programs across a network, ensuring high-confidence optimality despite unreliable communication and local uncertainties.

Contribution

It presents a novel asynchronous, randomized distributed method for robust mixed-integer programming with convergence guarantees and practical implementation on multi-core systems.

Findings

Algorithm achieves finite-time consensus on optimal solutions.

High-confidence solutions are feasible and nearly optimal.

Effective in uncertain, distributed environments like sensor networks.

Abstract

In this paper, we consider a network of processors aiming at cooperatively solving mixed-integer convex programs subject to uncertainty. Each node only knows a common cost function and its local uncertain constraint set. We propose a randomized, distributed algorithm working under asynchronous, unreliable and directed communication. The algorithm is based on a local computation and communication paradigm. At each communication round, nodes perform two updates: (i) a verification in which they check---in a randomized fashion---the robust feasibility of a candidate optimal point, and (ii) an optimization step in which they exchange their candidate basis (the minimal set of constraints defining a solution) with neighbors and locally solve an optimization problem. As main result, we show that processors can stop the algorithm after a finite number of communication rounds (either because verification has been successful for a sufficient number of rounds or because a given threshold has been reached), so that candidate optimal solutions are consensual. The common solution is proven to be---with high confidence---feasible and hence optimal for the entire set of uncertainty except a subset having an arbitrary small probability measure. We show the effectiveness of the proposed distributed algorithm using two examples: a random, uncertain mixed-integer linear program and a distributed localization in wireless sensor networks. The distributed algorithm is implemented on a multi-core platform in which the nodes communicate asynchronously.

Tables

Table 1. Table 1: A verage—over all nodes—number of times a basis is transmitted to the neighbors, average number of times verification is performed ( k i subscript 𝑘 𝑖 k_{i} at the convergence) and empirical violation of the computed solution ( 𝐱 sol subscript 𝐱 sol {\mathbf{x}}_{\text{sol}} ) over 10 , 000 10 000 10,000 random samples for different number of nodes and neighbors in each node. The simulation is performed 100 100 100 times for each row and average results are reported.

# Nodes	# Neighbors	Graph	# Constraints	# Transmissions	# Verifications	Empirical
$n$	in each node (degree)	diameter	in each node	(averaged)	(averaged)	violation
$10$	$3$	$4$	$100$	$29.86$	$15.6$	$8.18\times {10}^{-4}$
$50$	$4$	$4$	$100$	$37.58$	$10.73$	$2\times {10}^{-4}$
$100$	$6$	$4$	$100$	$41.73$	$9.79$	$3\times {10}^{-4}$
$200$	$8$	$4$	$100$	$43.94$	$8.92$	$2.3\times {10}^{-4}$

Table 2. Table 2: A verage—over all nodes—number of times a basis is transmitted to the neighbors, average number of times verification is performed ( k i subscript 𝑘 𝑖 k_{i} at the convergence) and empirical violation of the computed solution ( 𝐱 sol subscript 𝐱 sol {\mathbf{x}}_{\text{sol}} ) over 10 , 000 10 000 10,000 random samples for different number of nodes and neighbors in each node. We allow 10 10 10 violation samples to be returned by the verification step. The simulation is performed 100 100 100 times for each row and average results are reported.

# Nodes	# Neighbors	Graph	# Constraints	# Transmissions	# Verifications	Empirical
$n$	in each node (degree)	diameter	in each node	(averaged)	(averaged)	violation
$10$	$3$	$4$	$100$	$16.1$	$8.1$	$10\times {10}^{-3}$
$50$	$4$	$4$	$100$	$24.37$	$7.05$	$9\times {10}^{-3}$
$100$	$6$	$4$	$100$	$26.25$	$6.3$	$4\times {10}^{-3}$

Table 3. Table 3: The average—over all nodes—number of times a basis is transmitted to the neighbors, average number of times verification is performed ( k i subscript 𝑘 𝑖 k_{i} at the convergence) and empirical violation of the computed solution ( 𝐱 sol subscript 𝐱 sol {\mathbf{x}}_{\text{sol}} ) over 10 , 000 10 000 10,000 random samples for different number of known sensors.

# Nodes	# Transmissions	# Verifications	Empirical
$n$	(averaged)	(averaged)	violation
$10$	$7.77$	$8.77$	$1\times {10}^{-4}$
$50$	$12.48$	$7.36$	$0$
$100$	$9.37$	$6.52$	$0$
$200$	$7.95$	$7.42$	$0$

Equations

$$\min$$

$$\mathbf{x}\in \mathcal{F},$$

$$\mathbf{x}\in S.$$

$$\forall i_1,\dots ,i_h\in$$

$$K_1\cap \cdots \cap K_m\ne \varnothing$$

$$\min$$

$$\mathbf{x}\in \bigcap _{i=1}^n{\mathcal{F}}^i(q),\forall q\in \mathbb{Q},$$

$$\mathbf{x}\in S,$$

$${\mathbf{q}}_{k_i}^i\doteq \{q_{k_i,i}^{(1)},\dots ,q_{k_i,i}^{(M_{k_i}^i)}\}\in {\mathbb{Q}}^{M_{k_i}^i},$$

$$M_{k_{i}}^{i}=\Big{\lceil}\frac{2.3+1.1\ln k_{i}+\ln\frac{1}{\delta_{i}}}{\ln\frac{1}{1-\varepsilon_{i}}}\Big{\rceil}$$

$$\displaystyle\mathbb{S}:=\big{\{}\mathbf{q}:\;\mathbf{q}=$$

$$\displaystyle\mathbf{q}_{i}=\{\mathbf{q}_{k_{i}}^{i}\}_{\color[rgb]{0,0,0}{k_{i}=1,\ldots,\infty}},\;\;i=1,\dots,n\big{\}}.$$

$${\mathbb{S}}_{\text{sol}}:=\{\mathbf{q}\in \mathbb{S}:\text{Algorithm refalg: distributed algorithm terminates}\},$$

$$\displaystyle\mathbb{P}_{\infty}\!\bigg{\{}\mathbf{q}\in\mathbb{S}_{\textnormal{sol}}:$$

$$\displaystyle\leq\epsilon\bigg{\}}\!\geq\!1-\delta.$$

$$\displaystyle{\mathbb{P}}_{\infty}\!\bigg{\{}\!{{\rm\bf q}}\in\mathbb{S}_{\textnormal{sol}}:$$

$$\displaystyle\leq\varepsilon\bigg{\}}\geq 1\!-\!\delta$$

$$\min$$

$$\mathbf{x}\in \bigcap _{i=1}^n\bigcap _{\ell =1}^{M_i^{\text{scen}}}{\mathcal{F}}^i(q^{(\ell )}),$$

$$\mathbf{x}\in S,$$

$${\delta }_i\ge \sum _{\ell =0}^{h(S)-2}\left(\genfrac{}{}{0ex}{}{M_i^{\text{scen}}}{\ell }\right){\epsilon }_i^{\ell }(1-{\epsilon }_i{)}^{M_i^{\text{scen}}-\ell }.$$

$$\displaystyle{\mathbb{P}}^{M^{\textnormal{scen}}}\bigg{\{}{{\rm\bf q}}\in\mathbb{Q}^{M^{\textnormal{scen}}}:$$

$$\displaystyle\leq\varepsilon\bigg{\}}\geq 1-\delta,$$

$$\displaystyle\!\!\!\!\!\!\!{\mathbb{P}}^{M}\bigg{\{}{{\rm\bf q}}\in\mathbb{Q}^{M}:$$

$$\displaystyle\leq\varepsilon\bigg{\}}\geq 1-\delta,$$

$$\min$$

$$\mathbf{x}\in \bigcap _{i=1}^n\bigcap _{\ell =1}^{M_{k_i}^i}{\mathcal{F}}^i(q^{(\ell )}),$$

$$\mathbf{x}\in S,$$

$$\min$$

$$(A_i^0+A_i^q{)}^T\mathbf{x}\le b_i,i=1,\dots ,n$$

$$\mathbf{x}\in {\mathbb{Z}}^2\times {\mathbb{R}}^3,$$

$$M^{\text{scen}}_{i}\geq\frac{1.582}{\varepsilon}\bigg{(}\ln\frac{1}{\delta_{i}}+h(s)-2\bigg{)}.$$

$$k_i\ge \exp \left(\frac{0.58\ln \frac{1}{{\delta }_i}+1.58(h(S)-2)-2.3}{1.1}\right).$$

$$\mathcal{F}^{i}_{R}\doteq\Bigg{\{}{\mathbf{x}}\in\mathbb{R}^{2}:\left[\begin{array}[]{c c}rI_{2}&({\mathbf{x}}-p_{i})\\ ({\mathbf{x}}-p_{i})^{T}&r\end{array}\right]\succeq 0\Bigg{\}},$$

$${\mathbf{x}}_x^l$$

$$\text{subject to }\mathbf{x}$$

$$\Vert p_i-{\bar{p}}_i{\Vert }_2\le \rho ,$$

$$J({\mathcal{B}}^i(t_0))\le J({\mathcal{B}}^p(t_0+nL))\le J({\mathcal{B}}^i(t_0+2nL)).$$

$${\text{Viol}}_i$$

Full text

, , ,

††thanks: This work is supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme, grant agreement No 638992 - OPT4SMART, and a grant from the Singapore National Research Foundation (NRF) under the ASPIRE project, grant No NCR-NCR001-040. A preliminary conference version of this paper appeared as [15]. This paper addresses a more general optimization set-up, includes a more in-depth analysis and new numerical computations. Corresponding author Mohammadreza Chamanbaz, Tel. +6584019212.

1 Introduction

Robust optimization plays an important role in several areas such as estimation and control and has been widely investigated. Its rich literature dates back to the $1950$s, see [7] and references therein. Very recently, there has been a renewed interest in this topic in a parallel and/or distributed frameworks. In the early reference [29], a synchronous distributed random projection algorithm with almost sure convergence was proposed for the case where each node knows a local cost function and an uncertain constraint. In [30] , the same authors extended the algorithm to a gossip-based protocol (based on pairwise node selection). To prove almost sure convergence, the algorithms in [29, 30] require suitable assumptions on the random set, network topology and algorithmic weights. Authors in [34] proposed a distributed approach based on barrier methods. The algorithm needs to start from a strictly feasible point. In [14], a parallel/distributed scheme is considered for solving an uncertain optimization problem by means of the scenario approach [11, 12]. The scheme consists of extracting a number of samples from the uncertain set and assigning them to nodes in a network. Each node is assigned a portion of the extracted samples. Then, a variant of the constraints consensus algorithm introduced in [32] is used to solve the deterministic optimization problem. A similar parallel framework based on the scenario approach is considered in [36] to solve convex optimization problems with a common uncertain constraint. In this setup, the sampled optimization problem is solved in a distributed way by using a primal-dual subgradient (resp. random projection) algorithm over an undirected (resp. directed) graph. In [9], a cutting plane consensus algorithm is introduced for solving convex optimization problems where constraints are distributed through the network processors and all processors have a common cost function. In the case where constraints are uncertain, a worst-case approach based on a pessimizing oracle is used. A distributed scheme based on the scenario approach is introduced in [31]. Each node initially extracts random samples from its local uncertain constraint set and a distributed proximal minimization algorithm is designed to solve the sampled optimization problem. Xie and coauthors proposed in [36] a distributed algorithm for solving convex programs over time-varying and unbalanced digraphs. To this end, an epigraph form of the problem is considered to convert the objective function into a linear form and a two-step recursive algorithm is used which asymptotically converges to a solution. One of the assumptions in this setup is that each node needs to know its identifier. A distributed algorithm—based on an iterative scheme that combines dual decomposition and proximal minimization—is presented in [19] which takes into account privacy of information.

All aforementioned papers deal with usual (continuous) optimization problems. A centralized method—based on the scenario approach—is presented in [10] for solving robust mixed-integer convex optimization problems. Few recent works address the parallel or distributed solution of deterministic mixed-integer programs. A decentralized (parallel)—but not distributed—approach, based on dual decomposition, is proposed in [35, 21] to approximately solve mixed-integer linear programs (MILPs) with guaranteed suboptimality bounds. A central unit is used in the two parallel methods above to update and then broadcast back to agents the dual variables.

In [20] a distributed version of the algorithm in [21] is proposed. In [13] a distributed algorithm based on primal decomposition is proposed for the same MILP set-up. A Lagrange relaxation approach combined with proximal bundle has been used in [27] to solve a demand response problem involving mixed-integer constraints. In [23], the problem of heterogeneous multi-vehicle routing is formulated as a MILP and a decentralized algorithm, based on a gossip protocol, is proposed to find a feasible suboptimal solution almost surely. A problem of cooperative distributed robust trajectory optimization is addressed in [28]. The problem is formulated as a MILP and solved in a receding horizon framework where sequentially each vehicle solves its subproblem—optimizing its own control input—and perturbs the decision of other subproblems to improve the global cost. Finally , a distributed algorithm based on generation of cutting planes and exchange of active constraints is proposed in [33] to solve deterministic MILPs.

The main contribution of this paper —which extends the conference paper [15]—is the design of a fully distributed algorithm to solve uncertain mixed-integer programs in a network with directed, asynchronous, and possibly unreliable, communication. The problem under investigation is a mixed-integer program in which the constraint set is the intersection of local uncertain constraints, each one known only by a single node. Starting from a deterministic constraint exchange idea introduced in [32], the algorithm proposed in this paper introduces a randomized, sequential approach in which each node locally: (i) performs a probabilistic verification step (based on a “local” sampling of its uncertain constraint set), and (ii) solves a deterministic optimization problem with a limited number of constraints. If suitable termination conditions are satisfied, we are able to prove that the nodes reach consensus on a common solution , which is probabilistically feasible and optimal with high confidence. We also propose a stopping criterion ensuring that the distributed algorithm can in fact be stopped after a finite number of communication rounds. Although tailored to mixed-integer programs, the proposed algorithm can solve a more general class of problems known as S-optimization problems [16] (which include continuous, integer and mixed-integer optimization problems). As compared to the literature reviewed above, the proposed algorithm has three main advantages. First, no assumptions are needed on the probabilistic nature of the local constraint sets. Second, each node can sample locally its own uncertain set. Thus, no central unit is needed to extract samples and no common constraint set needs to be known by the nodes. Third and final, nodes do not need to perform the whole sampling at the beginning and subsequently solve the (deterministic) optimization problem. Online extracted samples are used only for verification, which is computationally inexpensive. The optimization is always performed on a number of constraints that remains constant for each node and depends only on the dimension of the decision variable and on the number of neighboring nodes. Furthermore, the distributed algorithm can be immediately implemented using existing off-the-shelf optimization tools, and no ad-hoc implementation of specific update rules is needed.

The paper is organized as follows. In Section 2, we provide some preliminaries and we formulate the distributed robust mixed-integer convex program. Section 3 presents our distributed, randomized algorithm for finding a solution—with probabilistic robustness—to robust distributed mixed-integer problems. The probabilistic convergence properties of the distributed algorithm are investigated in Section 3.2. Finally, extensive numerical simulations are performed in Section 4 to show the effectiveness of the proposed methodology.

2 Preliminaries & Problem Set-Up

In this section , we introduce some preliminary notions of combinatorial optimization that will serve us to set-up the Robust Mixed-Integer Convex Program addressed in the paper and the methods used for the proposed algorithm.

2.1 Preliminaries

Given a constraint set ${\mathcal{F}}={\mathcal{F}}^{1}\cap\dots\cap{\mathcal{F}}^{N}\subset\mathbb{R}^{d}$ and a vector $c\in\mathbb{R}^{d}$, we will denote with $({\mathcal{F}},c)$ the following optimization problem over $S\subset\mathbb{R}^{d}$

[TABLE]

Moreover, let us denote $J({\mathcal{F}})$ the optimal cost of the above problem. For a large class of problems, often known as abstract programs or LP-type problems, a solution to $({\mathcal{F}},c)$ can be identified by considering a subset of the constraints ${\mathcal{F}}^{1},\dots,{\mathcal{F}}^{N}$. This concept is characterized by the notion of basis, which is, informally, a “minimal” set of constraints that defines a solution. The concept of basis is supported by Helly-type theorems, initially introduced by Eduard Helly in [25], see also [2]. Before formally defining a basis, we define the Helly number and its variation $S$-Helly number [3, 16].

Definition 2.1** (Helly number).**

Given a nonempty family ${\mathcal{K}}$ of sets, Helly number $h=h({\mathcal{K}})\in\mathbb{N}$ of ${\mathcal{K}}$ is defined as the smallest number satisfying the following

[TABLE]

for all $m\in\mathbb{N}$ and $K_{1},\dots,K_{m}\in{\mathcal{K}}$. If no such $h$ exists then $h({\mathcal{K}})=\infty$.

For the classical Helly’s theorem, ${\mathcal{K}}$ is a finite family of convex subsets of $\mathbb{R}^{d}$. One of the important extensions of the Helly number is the $S$-Helly number defined for $S\subset\mathbb{R}^{d}$. $\square$

Definition 2.2** ($S$-Helly number).**

Given a set $S\subset\mathbb{R}^{d}$, let $S\cap{\mathcal{K}}$ be the family of sets $S\cap K$ such that $K\subset\mathbb{R}^{d}$ is convex. The $S$-Helly number is defined as $h(S)=h(S\cap{\mathcal{K}})$. $\square$

It is worth pointing out the main difference between the definitions given above. The Helly number states that if the intersection of every $h$ sets is nonempty, then there is a point in common between all the sets. Differently, the $S$-Helly number states that if the intersection of any $h$ convex sets contains a point in $S$, then all the sets have a point in common in $S$.

Definition 2.3** (Basis).**

Given a constraint set ${\mathcal{F}}={\mathcal{F}}^{1}\cap\cdots\cap{\mathcal{F}}^{N}$, a basis ${\mathcal{B}}$ of $({\mathcal{F}},c)$ is a minimal collection of constraints from ${\mathcal{F}}^{1},\dots,{\mathcal{F}}^{N}$, such that the optimal cost of the problem $({\mathcal{F}},c)$ is identical to the one of $({\mathcal{B}},c)$, i.e., $J({\mathcal{F}})=J({\mathcal{B}})$. $\square$

We point out that in this paper we are slightly abusing the notation since we are denoting by ${\mathcal{B}}$ both the collection of constraints (when referring to the basis) and their intersection (when denoting the constraint set of the optimization problem).

The size of a largest basis of problem $({\mathcal{F}},c)$ is called its combinatorial dimension. The following result connects the combinatorial dimension with the $S$-Helly number. This result is presented in [10, Theorem 2], [16, Theorem 1.3] and [3, Theorem 3.11].

Theorem 2.4**.**

Let ${\mathcal{F}}={\mathcal{F}}^{1}\cap\!\cdots\!\cap{\mathcal{F}}^{N}$, then the combinatorial dimension of problem $({\mathcal{F}},c)$ is $h(S)-1$. $\square$

It is worth noticing that if in problem $({\mathcal{F}},c)$, the sets ${\mathcal{F}}^{1},\dots,{\mathcal{F}}^{N}$ are convex, well known classes of problems arise depending on $S$. If $S={\mathbb{R}}^{d}$, problem $({\mathcal{F}},c)$ is a continuous convex optimization problem. If $S={\mathbb{Z}}^{d}$, then $({\mathcal{F}},c)$ becomes an integer optimization problem. Choosing $S={\mathbb{Z}}^{d_{Z}}\times{\mathbb{R}}^{d_{R}}$ with $d=d_{Z}+d_{R}$ leads to a mixed-integer convex problem. The S-Helly number for $S={\mathbb{Z}}^{d_{Z}}\times{\mathbb{R}}^{d_{R}}$ is given below, see [5].

Theorem 2.5** (Mixed-Integer Helly Theorem).**

The Helly number $h({\mathbb{Z}}^{d_{Z}}\times{\mathbb{R}}^{d_{R}})$ equals $(d_{R}+1)2^{d_{Z}}$. $\square$

This implies that the combinatorial dimension of $({\mathcal{F}},c)$ with $S={\mathbb{Z}}^{d_{Z}}\times{\mathbb{R}}^{d_{R}}$ is $(d_{R}+1)2^{d_{Z}}-1$.

2.2 Problem Set-Up

We consider a network of $n$ processors with limited computation and/or communication capabilities aiming at cooperatively solving the following Robust Mixed-Integer Convex Program (RMICP)

[TABLE]

where ${\mathbf{x}}\in S$ is the vector of decision variables, $q\in{\mathbb{Q}}$ is the vector of uncertain parameters acting on the system, while $\mathcal{F}^{i}(q)=\{{\mathbf{x}}\in S:f^{i}({\mathbf{x}},q)\leq 0\}$ is the constraint set known only by agent $i$, with $f^{i}({\mathbf{x}},q):\mathbb{R}^{d}\times{\mathbb{Q}}\rightarrow{\mathbb{R}}$ its related constraint function, which is assumed to be convex for any fixed value of $q\in\mathbb{Q}$. The objective function is considered to be linear. This assumption is without loss of generality. In fact, a nonlinear convex objective function can be transformed into the epigraph form by introducing an extra decision variable.

We make the following assumption regarding the solution of any deterministic subproblem of (1), i.e., a problem in which only a finite number of constraints ${\mathcal{F}}^{i}(q)$ (for a given $q$) are considered.

Assumption 2.6** (Non-degeneracy).**

The minimum point of any subproblem of (1) with at least $h(S)-1$ constraints is unique and there exist only $h(S)-1$ constraints intersecting at the minimum point. $\square$

Assumption 2.6 is not restrictive. In fact, to ensure uniqueness of the optimal point we could use a strictly convex objective function, a lexicographic ordering, or any universal tie-breaking rule, see [2, Observation 8.1] for further details.

Processor $i$ has only knowledge of the constraint set ${\mathcal{F}}^{i}(q)$ and of the objective direction $c$, which is common among all nodes. The goal is to design a fully distributed algorithm—consisting of purely local computation and communication with neighbors—used by each agent in order to cooperatively solve problem (1). We want to stress that there is no central node having access to all constraints.

We let the nodes communicate according to a time-dependent, directed communication graph $\mathcal{G}(t)=\{\mathcal{V},\mathcal{E}(t)\}$ where $t\in{\mathbb{N}}$ is a universal time which does not need to be known by nodes, $\mathcal{V}=\{1,\ldots,n\}$ is the set of agent identifiers and $(i,j)\in\mathcal{E}(t)$ indicates that $i$ sends information to $j$ at time $t$. The time-varying set of incoming ( resp. outgoing) neighbors of node $i$ at time $t$, $\mathcal{N}_{\text{in}}(i,t)$ ($\mathcal{N}_{\text{out}}(i,t)$), is defined as the set of nodes from ( resp. to) which agent $i$ receives ( resp. transmits) information at time $t$. A directed static graph is said to be strongly connected if there exists a directed path (of consecutive edges) between any pair of nodes in the graph. For time-varying graphs we use the notion of uniform joint strong connectivity formally defined next.

Assumption 2.7** (Uniform Joint Strong Connectivity).**

There exists an integer $L\geq 1$ such that the graph \bigg{(}\mathcal{V},\bigcup_{\tau=t}^{t+L-1}\mathcal{E}(\tau)\bigg{)} is strongly connected for all $t\geq 0$. $\square$

There is no assumption on how uncertainty $q$ enters problem (1), thus making its solution particularly challenging. In fact, if the uncertainty set $\mathbb{Q}$ is an uncountable set, problem (1) is a semi-infinite optimization problem involving an infinite number of constraints. In general, there are two main paradigms to solve an uncertain optimization problem of the form (1). The first approach is a deterministic worst-case paradigm in which the constraints are enforced to hold for all possible uncertain parameters in the set $\mathbb{Q}$. This approach is computationally intractable for cases where uncertainty does not appear in a “simple” form, e.g. affine, multi-affine, or convex. Moreover, since some uncertainty scenarios are very unlikely to happen, the deterministic paradigm may be overly conservative. The second approach— the one pursued in this paper—is a probabilistic approach where uncertain parameters are considered to be random variables and the constraints are enforced to hold for the entire set of uncertainty except a subset having an arbitrary small probability measure.

3 Randomized Constraints Consensus

In this section, we present a distributed, randomized algorithm for solving problem (1) in a probabilistic sense. Note that, since the uncertain constraint sets in (1) are uncountable, it is in general very difficult to verify if a candidate solution is feasible for the entire set of uncertainty or not. We instead use a randomized approach based on Monte Carlo simulation s to check probabilistic feasibility.

3.1 Algorithm description

The distributed algorithm we propose has a probabilistic nature and consists of two main steps: verification and optimization. The main idea is the following. A node has a candidate basis and candidate solution point. First, it verifies if the candidate solution point belongs to its local uncertain set with high probability. Then, it collects bases from neighbors and solves a convex mixed-integer problem with its basis and its neighbors’ bases as constraint set. If the verification step is not successful, the first violating constraint is also considered in the optimization.

Formally, we assume that $q$ is a random variable and a probability measure ${\mathbb{P}}$ over the Borel $\sigma-$algebra of $\mathbb{Q}$ is given. We denote by $k_{i}$ a local counter keeping track of the number of times the verification step is performed by agent $i$.

In the verification step, each agent $i$ generates a multi-sample ${{\rm\bf q}}_{k_{i}}^{i}$ with cardinality $M_{k_{i}}^{i}$ from the set of uncertainty

[TABLE]

according to the measure ${\mathbb{P}}$, where $\mathbb{Q}^{M_{k_{i}}^{i}}\doteq\mathbb{Q}\times\mathbb{Q}\times\cdots\times\mathbb{Q}$ ($M_{k_{i}}^{i}$ times). Node $i$ checks feasibility of the candidate solution ${\mathbf{x}}^{i}(t)$ only at the extracted samples (by simply checking the sign of $f^{i}({\mathbf{x}}^{i}(t),q_{k_{i},i}^{(\ell)})$ for each extracted sample $q_{k_{i},i}^{(\ell)}$, $\ell\in\{1,\ldots,M_{k_{i}}^{i}\}$) . If a violation happens, the first violating sample is used as a violation certificate111In fact, more than one violated sample—if they exist—can be returned by the verification step. See Remark 3.2 for more details. .

In the optimization step, agent $i$ transmits its current basis to all outgoing neighbors and receives bases from incoming ones. Then, it solves a convex mixed-integer problem whose constraint set is composed of: i) a constraint generated at the violation certificate (if it exists), ii) its current basis, and iii) the collection of bases from all incoming neighbors. We define a primitive $[{\mathbf{x}}^{*},{\mathcal{B}}]={\texttt{Solve\textsubscript{MIP}}}({\mathcal{F}},c)$ which solves the deterministic mixed-integer convex problem defined by the pair $({\mathcal{F}},c)$ and returns back the optimal point ${\mathbf{x}}^{*}$ and the corresponding basis ${\mathcal{B}}$. Node $i$ repeats these two steps until a termination condition is satisfied, namely if the candidate basis has not changed for $2nL+1$ times, with $L$ defined in Assumption 2.7. The distributed algorithm is formally presented in Algorithm 1.

At this point, it is worth highlighting some key interesting features of the proposed algorithm. First, the verification step, even if it may run possibly a large number of times, consists of simple, inexpensive inequality checks which are not computationally demanding. Moreover, we remark that if at some $t$ the candidate solution has not changed, that is ${\mathbf{x}}^{i}(t)={\mathbf{x}}^{i}(t-1)$, then ${\mathbf{x}}^{i}(t)$ has successfully satisfied a verification step and, thus, the algorithm does not perform it again. Second, the number of constraints involved in the optimization problem being solved locally at each node is fixed. Hence, the complexity of the problem does not change with time. Third, the amount of data a processor needs to transmit does not depend on the number of agents, but only on the dimension of the space. Indeed, processor $i$ transmits only $h(S)-1$ constraints at each iteration and for mixed-integer convex programs this number only depends on the dimension of the space, see Theorem 2.5. Fourth and final, the proposed distributed randomized algorithm is completely asynchronous and works under unreliable communication. Indeed, $t$ is a universal time that does not need to be known by the processors and the graph can be time-varying. Thus, if nodes run the computation at different speeds, this can be modeled by having no incoming and outgoing edges in that time interval. Similarly, if transmission fails at a given iteration, this is equivalent to assum ing that the associated edge is not present in the graph at that iteration.

Remark 3.1** (Local optimization and initial constraint).**

In the deterministic constraints consensus algorithm presented in [32], at each iteration of the algorithm, each node needs to include in the local optimization problem its original constraint set. Here, we can drop this requirement because of the uncertain nature of the problem handled in the verification step. $\square$

Remark 3.2** (Multiple Violation Certificates).**

In the verification step of Algorithm 1, we set $q^{\textnormal{{viol}}}_{t,i}$ as the first sample for which ${\mathbf{x}}^{i}(t)\notin\mathcal{F}^{i}(q^{\textnormal{{viol}}}_{t,i})$. However, if more than one violated sample is found, then, in step $O_{2}$ of Algorithm 1, we may include $r>1$ violated constraints in the SolveMIP primitive, thus solving a larger optimization problem. By doing this, one expects a faster convergence and hence less communication burden. In fact, as we will show in the numerical computations, the number of violated samples $r$ constitutes a trade-off between the number of times each node communicates with neighbors and the complexity of the local optimization problem it needs to solve. $\square$

3.2 Algorithm Convergence and Solution Probabilistic Guarantees

Here , we analyze the convergence properties of the distributed algorithm and investigate the probabilistic properties of the solution computed by the algorithm.

We start by defining the set of all possible successive independent random extractions of all the sequences \{\mathbf{q}_{k_{i}}^{i}\}_{\color[rgb]{0,0,0}{k_{i}=1,\ldots,\infty}}, $i=1,\dots,n$ as follows:

[TABLE]

Moreover, we define the set

[TABLE]

which is the domain of the random variables $\mathbf{x}_{\text{sol}}^{i}$. The probability measure on $\mathbf{x}_{\text{sol}}^{i}$ is therefore defined on this space, and it is denoted $\mathbb{P}_{\infty}$. Now , we are ready to present the first main result of the paper.

Theorem 3.3**.**

Let Assumptions 2.6 and 2.7 hold. Given the probabilistic levels $\varepsilon_{i}>0$ and $\delta_{i}>0$, $i=1,\ldots,n$, let $\varepsilon=\sum_{i=1}^{n}\varepsilon_{i}$ and $\delta=\sum_{i=1}^{n}\delta_{i}$. Then, the following statements hold:

Along the evolution of Algorithm 1, the cost $J({\mathcal{B}}^{i}(t))$ at each node $i\in\{1,\ldots,n\}$ is monotonically non-decreasing, i.e., $J({\mathcal{B}}^{i}(t+1))\geq J({\mathcal{B}}^{i}(t))$, and converges to a common value asymptotically. That is, $\lim_{t\rightarrow\infty}J({\mathcal{B}}^{i}(t))=\bar{J}$ for all $i\in\{1,\ldots,n\}$. 2. 2.

If the candidate solution of node $i$, ${\mathbf{x}}^{i}(t)$, has not changed for $2Ln+1$ communication rounds, all nodes have a common candidate solution, i.e., ${\mathbf{x}}_{\textnormal{sol}}^{i}={\mathbf{x}}_{\textnormal{sol}}$ for all $i=1,\dots,n$. 3. 3.

The following inequality holds for ${\mathbf{x}}_{\textnormal{sol}}$:

[TABLE] 4. 4.

Let ${\mathcal{B}}_{\textnormal{sol}}$ be the basis corresponding to ${\mathbf{x}}_{\textnormal{sol}}$. The following inequality holds for ${\mathcal{B}}_{\textnormal{sol}}$

[TABLE]

*where ${\mathcal{F}}(q)\doteq\bigcap_{i=1}^{n}{\mathcal{F}}^{i}(q)$. *$\square$

The proof is given in Appendix A.

We briefly discuss the results of this theorem. Statement 1 shows that agents are increasing their local cost and asymptotically achieve the same cost. This result gives insights on the following statement. Statement 2 provides a condition under which processors can stop the distributed algorithm with the guarantee of having a common solution, whose probabilistic properties are shown in the next statements. In the next subsection, we provide a variation of the algorithm that ensures to stop the algorithm in finite time. We point out, though, that numerical computations show that the algorithm in the original version always satisfies the condition of Statement 2. It is also worth mentioning that for a fixed communication graph the bound $2nL+1$ in Statement 2 to stop the algorithm can be replaced by $2D+1$ where $D$ is the graph diameter. Statement 3 says that the probability of violating the (common) solution ${\mathbf{x}}_{\text{sol}}$ with a new sample (constraint) is smaller than $\varepsilon$ and this statement holds with probability at least $1-\delta$. Similarly, based on Statement 4, the probability that the cost $J({\mathcal{B}}_{\text{sol}})$ associated to ${\mathbf{x}}_{\text{sol}}$ increases, i.e., the probability of ${\mathbf{x}}_{\text{sol}}$ not being optimal, if a new sample (constraint) is extracted, is smaller than $\varepsilon$ and this statement holds with probability at least $1-\delta$.

Remark 3.4** (S-Optimization).**

The class of problems that Algorithm 1 can solve is not limited to mixed-integer problems. The algorithm is immediately applicable to any problem having a finite $S$-Helly number. This class of problems have recently been introduced—under the name of S-optimization problems—in [16]. It is worth mentioning that some classes of non-convex problems have finite $S$-Helly number—see [3] and references therein—and hence can be solved using Algorithm 1. $\square$

3.3 Modified Algorithm with Finite-Time Convergence

Guarantee

In this subsection, we derive a stopping criterion for Algorithm 1 so that it is guaranteed to be halted in finite time. Indeed, there is no guarantee for the condition in Statement 2 of Theorem 3.3 to be satisfied (and thus for Algorithm 1 to be halted in finite-time), even though we found this to be the case for all simulations we run (see Section 4 for more details).

The stopping criterion is derived by borrowing tools from the scenario approach presented in [11, 12]. In the scenario approach, the idea is to extract a finite number of samples from the constraint set of the semi-infinite optimization problem (1) and solve the obtained (deterministic) optimization problem. The developed theory provides bounds on the number of samples to guarantee a desired probabilistic robustness. The main idea of the stopping criterion is the following. Once each $M_{k_{i}}^{i}$ (the cardinality of the multi-sample used in the verification step of Algorithm 1) exceeds a (common) threshold associated to a proper scenario bound, agents stop generating new samples in the verification step and Algorithm 1 proceeds with the multi-sample extracted at the last iteration. Thereby, Algorithm 1 solves (in a distributed way) a scenario problem and in finite time finds a solution feasible for the extracted samples across all the nodes. Each node has a flag which is communicated locally among all neighbors. The flag is initialized to zero and is set to one once $M_{k_{i}}^{i}$ is greater than or equal to the scenario bound. We remark that using, for example, the algorithm presented in [22] nodes can identify—in a fully distributed manner—if all flags across the network are set to one and thus stop generating new samples.

In [31] , scenario bounds are proposed for a distributed robust convex optimization framework in which agents use a different set of uncertainty scenarios in their local optimization programs. In the next Lemma we provide analogous bounds for our problem setup. To this aim, let us consider the following optimization problem

[TABLE]

where $M_{i}^{\text{scen}}$ is the smallest integer satisfying

[TABLE]

Lemma 3.5**.**

*Given the probabilistic levels $\varepsilon_{i},\delta_{i}>0,\,i=1,\ldots,n$, let $\varepsilon=\sum_{i=1}^{n}\varepsilon_{i}$, $\delta=\sum_{i=1}^{n}\delta_{i}$ and $M^{\textnormal{scen}}=\sum_{i=1}^{n}M_{i}^{\textnormal{scen}}$, where $M_{i}^{\textnormal{scen}}$ is the smallest integer satisfying (4). If agents formulate the sampled optimization problem (3) and all of them agree on an optimal solution ${\mathbf{x}}^{*}$, *

[TABLE]

where ${\mathbb{Q}}^{M^{\textnormal{scen}}}={\mathbb{Q}}\times{\mathbb{Q}}\times\ldots\times{\mathbb{Q}}$ ($M^{\textnormal{scen}}$ times) and ${\mathbb{P}}^{M^{\textnormal{scen}}}$ is the product probability measure on ${\mathbb{Q}}^{M^{\textnormal{scen}}}$. $\square$

We omit the proof since it follows arguments similar to those in in [31, Proposition 1]. Specifically, by using the result in [10, Theorem 3] , one can show that the same arguments in [31, Proposition 1] can be followed by properly changing the cardinality of the support constraint set. In particular, in the proof of [31, Proposition 1] the (Helly) number $d+1$ of the (continuous) convex program, with $d$ being the dimension of the decision variable, can be replaced by the $S$-Helly number $(d_{R}+1)2^{d_{Z}}$ (as from Theorem 2.5) of problem (1).

The following Proposition summarizes the convergence properties and the solution guarantees of the modified algorithm (Algorithm 1 with the stopping criterion). We report only the feasibility result corresponding to Statement 3 of Theorem 3.3, but also the equivalent of Statement 4 can be proven.

Proposition 3.6**.**

Given the probabilistic levels $\varepsilon_{i}>0$ and $\delta_{i}>0$, $i=1,\ldots,n$, let $\varepsilon=\sum_{i=1}^{n}\varepsilon_{i}$ and $\delta=\sum_{i=1}^{n}\delta_{i}$. If nodes stop performing steps $V_{2}$ and $V_{4}$ of Algorithm 1 when $M_{k_{i}}^{i}\geq M_{i}^{\textnormal{scen}}$ for all $i=1,\ldots,n$, then the following statements hold:

nodes agree on a solution ${\mathbf{x}}_{\textnormal{sol}}$ in ** finite time; 2. 2.

the solution ${\mathbf{x}}_{\textnormal{sol}}$ satisfies

[TABLE]

where $M=\sum_{i}M_{k_{i}}^{i}$.

Proof 1**.**

First, note that if the Algorithm does not stop due to Condition 2 of Theorem 3.3, then all nodes stop generating new samples in finite time. When all nodes have stopped generating samples they start solving the following deterministic optimization problem

[TABLE]

*where $M_{k_{i}}^{i}\geq M_{i}^{\textnormal{scen}}$ is the cardinality of the verification multi-sample at the stopping point. To prove that nodes agree in finite time on an optimal solution of (6) , we resort to the following two arguments. First, as from Statement 1 of Theorem 3.3, the cost is monotonically non-decreasing. Second, the number of constraints in the network is fixed, thus the number of possible bases is finite. The proof of the first statement, thus, follows from arguments analogous as those in [32]. Indeed, when nodes stop generating samples, the Randomized Constraints Consensus algorithm becomes deterministic and turns out to be a variant of the Constraints Consensus algorithm proposed in [32]. Due to Assumption 2.6, nodes agree on a ${\mathbf{x}}_{\textnormal{sol}}$, which is the unique optimal solution of problem (6). Thus, the second statement follows by noting that, from Lemma 3.5, ${\mathbf{x}}_{\textnormal{sol}}$ satisfies (5) since the number of random constraints at each node $M_{k_{i}}^{i}$ is grater than or equal to the scenario bound $M_{i}^{\textnormal{scen}}$. *

4 Numerical Simulation s

We test the effectiveness of the distributed algorithm presented in Section 3 through extensive numerical simulations. To this end, we consider two different problems: (1) randomly generated mixed-integer linear programs (MILP) with uncertain parameters, and (2) distributed convex position estimation in wireless sensor networks. These two numerical simulations are discussed in the next subsections.

Numerical computations are run on a Linux -based high performance computing cluster222 http://idc.sutd.edu.sg/titan/ with $256$ CPUs. Each node uses only one CPU of the cluster and executes Algorithm 1 in an independent Matlab environment. The communication is modeled by sharing files between different Matlab environments over a fixed digraph.

4.1 Uncertain Mixed-Integer Linear Programming

We randomly generate robust MILPs with the following structure

[TABLE]

where $A_{i}^{0}$ is a fixed (nominal) matrix and $A_{i}^{q}$ is an interval matrix—a matrix whose entries are bounded in given intervals—defining the uncertainty in the optimization problem. We follow the methodology presented in [18] in order to generate the pair $A_{i}^{0},~{}b_{i}$ and the objective direction $c$. In particular, elements of $A_{i}^{0}$ are drawn from a standard Gaussian distribution , i.e. zero mean and unit standard deviation. The entries of $A_{i}^{q}$ are bounded in $[-\rho,\rho]$ for some given $\rho$. The $\ell$-th element of $b_{i}$ is defined by b_{i,\ell}=\gamma\,(\sum_{m=1}^{d}\big{(}A^{0}_{\ell m}\big{)}^{2})^{1/2}. This procedure generates feasible linear programs. However, introducing discrete variables in a feasible linear program may lead to an infeasible MILP. To ensure feasibility of the MILP, we increase the volume of the feasible region using the parameter $\gamma>1$. In the set of simulations reported here, we set $\gamma=20$. The objective direction $c$—which is the same for all nodes—is drawn from the standard Gaussian distribution. The communication graph $\mathcal{G}$ is a (connected) random $k$-nearest neighbor graph (with $k$ being the fixed degree or number of neighbors) [24]. Over this fixed graph, agents implement the distributed algorithm according to the asynchronous and unreliable communication protocol (based on Matlab environments) described above. In particular, only one node at a time is able to read/write from/to each node file containing its basis. Hence, if more than one node is willing to read/write from/to a file, only one would be allowed to do so. This gives raise to an asynchronous and unreliable communication protocol.

We use the intlinprog function of Mosek [4] to solve optimization problems appearing at each iteration of the distributed algorithm. As for the basis computation , it is worth mentioning that in (continuous) convex optimization problems, bases can be found by identifying a minimal set of active constraints. For mixed-integer problems, however, finding the bases can be computationally expensive. A computationally manageable procedure for finding a bases set—which might not be of minimum cardinality—is to check all constraints one by one to see whether or not they belong to a basis. To this end, we drop one constraint at a time and solve the resulting optimization problem. If the optimal objective of this optimization problem is smaller than the objective value of the original problem, we conclude that the dropped constraint is part of the basis. However, we note that, this brute-force approach is performed on a limited number of constraints—the constraint formed at the violation certificate, constraints in the node basis, and constraints in the neighboring bases.

For the computations reported in Table 1, we varied the number of nodes and neighbors per node, while using only those graphs having a diameter equal to $4$. The number of constraints per node is set to $100$. We also consider all elements of $A_{q}$ to be bounded in $[-0.2,0.2]$. The underlying probability distribution of uncertainty appearing in $A_{q}$ is selected to be uniform, see, e.g., [6]. The probabilistic accuracy and confidence levels of each agent are $\varepsilon_{i}=0.1/n$ and $\delta_{i}=10^{-9}/n$ respectively, with $n$ being the number of nodes (first column of Table 1). We report the average—over all nodes—number of times each node updates its basis and transmits it to the outgoing neighbors. It is assumed that each node keeps the latest information received from neighbors and hence, if the basis is not updated, there is no need to re-transmit it to the neighbors. This also accounts for the asynchronicity of the distributed algorithm. We also report the average—over all nodes—number of times node $i$ performs the verification step, i.e., the average —over all the nodes—of the $k_{i}$s at convergence. This allows us to show that with a relatively small number of “design” samples used in step $O_{2}$ of Algorithm 1, nodes compute a solution with a high degree of robustness. In order to examine the robustness of the obtained solution, we run an a posteriori analysis based on Monte Carlo simulations. To this end, we check the feasibility of the obtained solution for $10,000$ random samples extracted from the uncertain sets of the nodes. The empirical violation (last column of Table 1) is measured by dividing the number of samples that violate the solution by $10,000$. Since Algorithm 1 has a stochastic nature, for each row of Table 1 we generate randomly $100$ problems, solve them using Algorithm 1, run an a posteriori analysis and then report the average values.

In Figure 1, we report the distance of objective value $J({\mathcal{B}}^{i}(t))$ and candidate solutions ${\mathbf{x}}^{i}(t),\forall i\in\{1,\ldots,n\}$ from $J({\mathcal{B}}_{\text{sol}})$ and ${\mathbf{x}}_{\text{sol}}$ respectively along the distributed algorithm execution for a problem instance corresponding to the last row of Table 1.

In the following computations , we test the effect of using more violation certificates in the optimization step—see Remark 3.2. Specifically, we use in the optimization step up to $10$ violating certificates obtained from the verification step. As expect ed, this decreases the number of communications required for convergence at the expense of the local optimization complexity. This result is reported in Table 2.

In all $700$ simulations reported in Tables 1 and 2, Algorithm 1 converges to a solution with desired probabilistic robustness in finite time. We further remark that in none of the simulations we ran, the stopping criterion derived in subsection 3.3 was met. In fact , a closer investigation of the stopping condition presented in Lemma 3.6—namely $M_{k_{i}}^{i}\geq M_{i}^{\text{scen}}\text{ for all }i\in\{1,\ldots,n\}$—reveals that this condition happens only for extremely large values of the verification counter $k_{i}$. In [1, Theorem 4], an analytical sub-optimal solution for the sample size $M_{i}^{\text{scen}}$ in the inequality (4) is provided

[TABLE]

Solving the inequality $M_{k_{i}}^{i}\geq M^{\text{scen}}_{i}$ for $k_{i}$ and noting that \ln\big{(}\frac{1}{1-\varepsilon_{i}}\big{)}\simeq\varepsilon_{i} for small $\varepsilon_{i}$, we obtain

[TABLE]

The Helly number associated to the MILP envisaged here is $16$. Considering, for example, $\varepsilon_{i}=0.01$ and $\delta_{i}=10^{-10}$, the verification counter ensuring that $M_{k_{i}}^{i}\geq M^{\text{scen}}_{i}$ becomes $k_{i}\geq 1.2\times 10^{13}$. By looking at the sixth column of Tables 1 and 2 , which indicates the value of verification counter at convergence, one can see that the distributed algorithm converges to a solution with desired probabilistic properties in a much smaller number of verification steps.

4.2 Distributed Localization in Wireless Sensor Networks

The second numerical example is a problem of distributed position estimation in wireless sensor networks. A centralized version of this problem—with no uncertainty—is formulated in [17].

Consider a two-dimensional333Extension to a three-dimensional space is straightforward. space containing a number of heterogeneous wireless sensors which are randomly placed over a given environment. The sensors are of two classes. The first class are wireless sensors with known positions (SwKP) which are capable of positioning themselves up to a given accuracy, i.e., with some uncertainty in their position. They are equipped with processors with limited computation/communication capabilities and play the role of computational nodes in the distributed optimization framework. These sensors can communicate with each other based on a metric distance. That is, two sensors which are close enough can establish a bidirectional communication link. The second class is a wireless sensor with unknown position (SwUP) which has no positioning capabilities. This sensor is only equipped with a short-range transmitter having an isotropic communication range.

We further consider a heterogeneous setup in which some of the SwKPs are equipped with a laser transceiver providing an estimate of the relative angle with respect to the SwUP, which is within the range of the laser transceiver. This “angular constraint” can be represented by the intersection of three half-spaces $\mathcal{F}^{i}_{A}\doteq\{{\mathbf{x}}\in\mathbb{R}^{2}:a_{k}{\mathbf{x}}-b_{k}\leq 0,k=1,2,3\}$, where ${\mathbf{x}}\in\mathbb{R}^{2}$ is the position of SwUP and parameters $a_{k}\in\mathbb{R}^{1\times 2}\text{ and }b_{k}\in\mathbb{R}$ define the three half-spaces (two to bound the angle and one corresponding to the transceiver range) as presented in Fig. 2. SwKPs, which are not equipped with laser transceiver, can place a “radial constraint” for the SwUP if it is within the communication range. If the signal from the SwUP can be sensed by a SwKP, a proximity constraint exists between them, that is, the SwUP is placed in a circle with radius $r$ centered at the position of SwKP with $r$ being the communication range. This constraint can be formulated as $\|{\mathbf{x}}-p_{i}\|_{2}\leq r$, where $p_{j}\in\mathbb{R}^{2}$ is the position of $i$th SwKP. Using the Schur-complement [8], the feasible set of radial constraint can be represented as

[TABLE]

where $I_{2}$ is the $2\times 2$ identity matrix and $\succeq$ denotes positive semidefiniteness of the matrix .

The objective of the network ed system is to “estimate”—in a distributed fashion—the position of the SwUP. More precisely, SwKPs try to find the smallest box $\{{\mathbf{x}}\in\mathbb{R}^{2}:[{\mathbf{x}}_{x}^{l},{\mathbf{x}}_{y}^{l}]\leq{\mathbf{x}}\leq[{\mathbf{x}}_{x}^{u},{\mathbf{x}}_{y}^{u}]\}$ that is guaranteed to contain the SwUP, i.e., the red dotted box in Fig. 2. As shown in Fig. 2, the position of the SwUP is often sensed and hence constrained by several SwKP. This is used to shrink the feasibility set of the position of the SwUP. We remark that in a distributed setup, constraints regarding the position of the SwUP are distributed through different computational nodes and there is no agent having full knowledge of all constraints. The bounding box can be computed by solving four optimization problems with linear objectives. For instance ${\mathbf{x}}_{x}^{l}$ can be obtained by solving the following optimization problem

[TABLE]

where $n$ is the number of SwKPs capable of constraining the position of the SwUP. We point out that for SwKPs that are not equipped with a laser transceiver $\mathcal{F}^{i}_{A}=\mathbb{R}^{2}$ holds.

SwKPs are capable of localizing themselves. Nevertheless, their localization is not accurate and there is an uncertainty associated with the position vector reported by each of these nodes. This uncertainty can be modeled as a norm-bounded vector

[TABLE]

where $\bar{p}_{i}$ is the nominal position reported by sensor $i$, $p_{i}$ is its actual position and $\rho$ is the radius of the uncertain set. In other words, the uncertainty is defined by an $\ell_{2}$ -ball centered at the nominal point $\bar{p}_{i}$ with radius $\rho$.

For simulation purposes, we consider a $10\times 10$ square environment containing $n$ SwKPs (computational nodes) whose purpose is to estimate the position of the SwUP. We remark that the extension to the case in which there is more than one SwUPs is straightforward.

Two SwKPs can establish a two-way communication link if their distance is less than $7$ units of distance. The communication range of SwUP is also considered to be $7$ units. Among the $n$ SwKPs, half of them are able to place angular constraints—using laser transceiver—on their relative angle to the SwUP within the laser range. The angle representing the accuracy of the laser transceiver—the angle $\alpha$ in Fig. 2—is $20$ degrees. Finally, the uncertainty radius $\rho$ in (8) is $0.1$ and the distribution of the uncertainty is selected to be uniform. The probabilistic accuracy and confidence parameters are selected to be $\varepsilon_{i}=0.1/n$ and $\delta_{i}=10^{-9}/n$ for all $i=1,\ldots,n$, leading to $\varepsilon=0.1$ and $\delta=10^{-9}$.

The distributed robust localization problem is solved for different values of the number of SwKPs $n$. The result is reported in Table 3. Similar to Tables 1 and 2, we report the average—over all $n$ nodes—number of times nodes transmit their basis and run the verifications step. We also run an a posteriori analysis with Monte Carlo simulation s using $10,000$ random samples to check the robustness of the computed solution. The empirical violation reported by the Monte Carlo simulation s for each formulated scenario is reported in the last column of Table 3. We remark that the number of uncertain parameters is $2n$, thus rendering the problem nontrivial to solve. For instance, considering the last row of Table 3, where $n=200$, there are $400$ uncertain parameters in the problem.

To conclude, we would like to point out an interesting feature of Algorithm 1. In the results presented in Tables 1, 2, and 3, one can observe that when the number of computational nodes increases while keeping the graph diameter constant, the number of transmissions required for convergence does not change significantly. This suggests that the number of transmissions required for convergence is independent from the number of computational nodes but rather depends on the graph diameter.

5 Conclusions

In this paper, we proposed a randomized distributed algorithm for solving robust mixed-integer convex programs in which constraints are scattered across a network of processors communicating by means of a directed time-varying graph. The distributed algorithm has a sequential nature consisting of two main steps: verification and optimization. Each processor iteratively verifies a candidate solution through a Monte Carlo simulations, and solves a local MICP whose constraint set includes its current basis, the collection of bases from neighbors and possibly, a constraint—provided by the Monte Carlo algorithm—violating the candidate solution. The two steps, i.e. verification and optimization, are repeated until a local stopping criterion is met and all nodes converge to a common solution. We analyzed the convergence properties of the proposed algorithm.

Appendix A Proof of Theorem 3.3

Proof of Statement 1. The proof of the first statement follows arguments very similar to those in [9] and is thus omitted.

Proof of Statement 2. The proof of this statement heavily relies on results in [26], but we report it here for completeness. Since the graph is uniformly jointly strongly connected, for any pair of nodes $u$ and $v$ and for any $t_{0}>0$, there exists a time-dependent path from $u$ to $v$ [26]—a sequence of nodes $\ell_{1},\ldots,\ell_{k}$ and a sequence of time instances $t_{1},\ldots,t_{k+1}$ with $t_{0}\leq t_{1}<\ldots<t_{k+1}$, such that the directed edges $\{(u,\ell_{1}),(\ell_{1},\ell_{2}),\ldots,(\ell_{k},v)\}$ belong to the directed graph at time instances $\{t_{1},\ldots,t_{k+1}\}$, respectively. We now prove that the path from $u$ to $v$ is of length at most $nL$. We recall that $n$ is the number of nodes and $L$ is defined in Assumption 2.7. To this end, following [26], define a time-varying set $S_{t}$ and initialize it with the node $u$ itself, i.e. $S_{t_{0}}={u}$. Next, given $S_{t}$ construct $S_{t+1}$ by adding all nodes $\ell_{t+1}$ such that $\ell_{t+1}\in\mathcal{N}_{\text{out}}(\ell_{t},t+1)$ with $\ell_{t}$ being all the nodes in the set $S_{t}$. Now, notice that if $S_{t}$ does not contain all nodes, i.e., $S_{t}\neq\mathcal{V}$ then after $L$ time instants at least one new node must be added to the set $S_{t}$, otherwise taking the union graph over the $L$ time instants, the set $S_{t}$ and that node would not be connected by any path, thus violating Assumption 2.7. Therefore, in at most $nL$ time instants $S_{t}$ will contain all nodes including $v$. Meaning that the length of time-varying path connecting $u$ and $v$ is at most $nL$. Consider nodes $i$ and $p$. If $\ell_{1}\in\mathcal{N}_{\text{out}}(i,t_{0})$, then $J({\mathcal{B}}^{i}(t_{0}))\leq J({\mathcal{B}}^{\ell_{1}}(t_{0}+1))$ as the constraint set of node $\ell_{1}$ at time $t_{0}+1$ is a superset of the constraint set of node $i$ at time $t_{0}$. Iterating this argument, we obtain $J({\mathcal{B}}^{i}(t_{0}))\leq J({\mathcal{B}}^{p}(t_{0}+nL))$. Again since the graph is uniformly jointly strongly connected, there will be a time varying path of length at most $nL$ from node $p$ to node $i$. Therefore,

[TABLE]

If $J({\mathcal{B}}^{i}(t_{0}))=J({\mathcal{B}}^{i}(t_{0}+2nL))$ and considering the point that node $p$ can be any node of the graph, then all nodes have the same cost. That is, $J({\mathcal{B}}^{1}(t))=\ldots=J({\mathcal{B}}^{n}(t))$. This combined with Assumption 2.6 proves Statement 2.

Proof of Statement 3. We start by defining the following violation sets:

[TABLE]

We observe that

[TABLE]

hence

[TABLE]

therefore

[TABLE]

At this point, we note that, applying simple properties coming from the probability theory, we have

[TABLE]

As a consequence, if we prove that

[TABLE]

then the thesis will follow.

Now, let $i=1,\dots,n$ and $k{{}_{i}}\in\mathbb{N}$. We define the events:

[TABLE]

We immediately note that (9) is equivalent to:

[TABLE]

Thus we will prove (10) instead of (9).

First, we observe that if $\mathbf{q}\in\text{ExitBad}_{i}$, then there exists $k_{i}$ such that $\mathbf{q}\in\text{Feas}_{k_{i}}^{i}$. Consequently, the following relation holds:

[TABLE]

Therefore, we have

[TABLE]

and we can bound the terms of the last series as follows:

[TABLE]

A few comments are given in order to explain the above set of equalities and inequalities: equation (12) follows from the chain rule; inequality (13) is due to the fact that $\mathbb{P}_{\infty}\leq 1$; in equation (14) we have written explicitly the events $\text{Feas}_{k_{i}}^{i}$ and $\text{ExitBad}_{i}$; equation (15) exploits the property \mathbb{P}\big{\{}\mathbb{Q}-A\big{\}}=1-\mathbb{P}\big{\{}A\big{\}}; inequality (16) is due to the fact that in order to declare $\mathbf{x}_{\text{sol}}$ feasible, the algorithm needs to perform $M^{i}_{k_{i}}$ independent successful trials each having a probability of success smaller than or equal to $(1-\epsilon_{i})$.

At this point, combining (16) and (11), we obtain

[TABLE]

and the last series can be made arbitrarily small by an appropriate choice of $M^{i}_{k_{i}}$. In particular, by choosing

[TABLE]

where $\xi(\alpha):=\sum_{s=1}^{\infty}\frac{1}{s^{\alpha}}$ is the Riemann Zeta function evaluated at $\alpha>1$, we have

[TABLE]

Hence, the choice of sample size is obtained by solving (17) for $M^{i}_{k_{i}}$:

[TABLE]

Note that the optimal value of $\alpha$ minimizing $M^{i}_{k_{i}}$ is computed by trial and error to be $\alpha=1.1$; by means of this choice of $\alpha$, we obtain the number $M^{i}_{k_{i}}$ specified in the algorithm, and we can thus conclude that (10) is true, and so the thesis follows.

Proof of Statement 4. Define the following two sets

[TABLE]

where ${\mathcal{F}}(q)\doteq\bigcap_{i=1}^{n}{\mathcal{F}}^{i}(q)$. We first prove that $V=W$.

Let $q_{v}\in V$, that is, ${\mathbf{x}}_{\text{sol}}\notin\bigcap_{i=1}^{n}\mathcal{F}^{i}(q_{v})$, then $J({\mathcal{B}}_{\text{sol}}\cap{\mathcal{F}}(q_{v}))\geq J({\mathcal{B}}_{\text{sol}})$ with ${\mathcal{F}}(q_{v})\doteq\bigcap_{i=1}^{n}{\mathcal{F}}^{i}(q_{v})$. Also, due to Assumption 2.6, any subproblem of (1) has a unique minimum point and hence, $J({\mathcal{B}}_{\text{sol}}\cap{\mathcal{F}}(q_{v}))\neq J({\mathcal{B}}_{\text{sol}})$. Hence, $q_{v}\in W$ and subsequently $V\subseteq W$.

Now let $q_{v}\in W$ that is J\big{(}{\mathcal{B}}_{\text{sol}}\cap{\mathcal{F}}(q_{v})\big{)}>J({\mathcal{B}}_{\text{sol}}) and suppose by contradiction that ${\mathbf{x}}_{\text{sol}}\in\mathcal{F}(q_{v})$, i.e. $q_{v}\notin V$. Considering the fact that ${\mathcal{B}}_{\text{sol}}$ is the basis corresponding to ${\mathbf{x}}_{\text{sol}}$, we conclude that ${\mathbf{x}}_{\text{sol}}\in\mathcal{B}_{\text{sol}}\cap\mathcal{F}(q_{v})$ implying that $J({\mathcal{B}}_{\text{sol}}\cap{\mathcal{F}}(q_{v}))\leq J({\mathcal{B}}_{\text{sol}})$. This contradicts the fact that J\big{(}{\mathcal{B}}_{\text{sol}}\cap{\mathcal{F}}(q_{v})\big{)}>J({\mathcal{B}}_{\text{sol}}). Hence, $q_{v}\in V$ and subsequently, $W\subseteq V$. Since $V\subseteq W$ and $W\subseteq V$, therefore $V=W$.

This argument combined with the result of Statement 4 proves the fifth statement of the theorem. That is, the probability that the solution ${\mathbf{x}}_{\text{sol}}$ is no longer optimal for a new sample equals the probability that the solution is violated by the new sample.