Optimal Transport-based Coverage Control for Swarm Robot Systems: Generalization of the Voronoi Tessellation-based Method

Swarm robot systems, which consist of many cooperating mobile robots, have attracted attention for their environmental adaptability and fault tolerance advantages. One of the most important tasks for such systems is coverage control, in which robots autonomously deploy to approximate a given spatial distribution. In this study, we formulate a coverage control paradigm using the concept of optimal transport and propose a novel control technique, which we have termed the optimal transport-based coverage control (OTCC) method. The proposed OTCC, derived via the gradient flow of the cost function in the Kantorovich dual problem, is shown to covers a widely used existing control method as a special case. We also perform a Lyapunov stability analysis of the controlled system, and provide numerical calculations to show that the OTCC reproduces target distributions with better performance than the existing control method.


I. INTRODUCTION
Swarm robot systems, in which many mobile robots work cooperatively to perform given tasks, are expected to have strong environmental adaptability and high fault tolerance in comparison with single-robot systems [1][2][3]. One of the most fundamental and important challenges of such systems is coverage control, in which robots move and reposition themselves autonomously so that their placement approaches a predetermined spatial distribution. The application ranges from the optimal placement of sensor networks to efficient rescue of human life in the event of a disaster [4]. From the early 2000s to the present, various methods have been proposed, such as potential-function-based control [5], probability-based control [6,7], and broadcast-based control [8,9], in order to solve the coverage tasks. A detailed review of these methods is provided in [10].
Among them, the Voronoi tessellation-based coverage control (VTCC) method proposed by Cortes et al. [11] is a seminal work and widely used. In this method, the coverage area is divided into subspaces referred to as Voronoi regions, each of which is assigned to a robot, and the robot is moved toward the center of gravity of its assigned Voronoi region. The cost function defined for the entire robot swarm is shown to decrease over time, meaning that eventually the robots are appropriately scattered throughout the coverage area. The VTCC has been commonly used in practice for its mathematical guarantee of stability, as well as the simplicity and the scalability of the algorithm [12].
By regarding a robot swarm as an abstract group of points in Euclidean space, the coverage control is interpreted as the problem of transporting a given discrete distribution to approximate a target continuous distribution. This is commonly referred to as the optimal transport problem, and its mathematical properties and numerical solutions have been widely investigated [13][14][15]. Recent areas of interest concerning the optimal transport problem extend to applications such as machine learning [16][17][18], image processing [19,20], and natural language processing [21].
In this study, we formulate coverage control as an optimal transport problem in order to propose a novel control technique, which we call the optimal transport-based coverage control (OTCC) method. Multi-agent control methods using optimal transport have been proposed in [22,23]. However, the relation between the control laws proposed in these references and the VTCC has not been investigated. Our control method differs from existing methods in that it considers gradient flows for the cost function of the optimal transport problem, which allows us to compare the structure and performance of the proposed OTCC with that of the VTCC. The contributions of the present coverage control formulation is summarized as follows: • The cost function for the VTCC is shown to be a special case of the cost function for the Kantorovich dual problem.
• The new control law is derived as the gradient flows of the cost function for the Kantorovich dual problem.
• A sufficient condition for the Lyapunov stability is provided for the controlled system, followed by a more specific condition in one-dimensional case. * daisuke-inoue@mosk.tytlabs.co.jp arXiv:2011.08337v1 [math.OC] 16 Nov 2020 • Numerical analysis is conducted to show that the proposed method reaches a closer state to the global optimum than can be achieved via the VTCC.
Notation: Let R, R + , and N be a set of real, non-negative real, and positive integer numbers, respectively. The Euclidean norm of x ∈ R n is represented as x . We call U z ⊂ R n is a neighborhood of z ∈ R n if there exists an open ball B r := {x | x − z < r} for some r > 0 and B r ⊂ U z holds. For a real-valued function G : R n × R m → R, we denote the partial derivative of G(p, q) with respect to p as ∇ p G and with respect to q as ∇ q G. The higher-order derivatives follow the convention ∇ 2 p G = ∂ 2 G ∂p 2 , ∇ 2 q G = ∂ 2 G ∂q 2 , and so on.

II. REVIEW OF OPTIMAL TRANSPORT
This section provides a brief overview of optimal transport in order to assist in formulating the optimal transportbased coverage control. We begin by considering two continuous density functions ρ 0 : R d → R + and ρ T : Definition 1 (Kantorovich problem). The problem of finding a simultaneous probability density function p(x, y) that minimizes the following cost function C K (p) is called the Kantorovich problem: where p(x, y) satisfies the following conditions: We denote the solution of the Kantorovich problem as W (ρ 0 , ρ T ) := inf p C K (p) and call it the Wasserstein metric. In fact, W is a function that measures the distance between the two distributions, and it is known that W actually satisfies the axiom of distance [14].
Definition 2 (Kantorovich dual problem). The problem of finding integrable functions φ : R d → R and ψ : R d → R that maximize the following cost functions C KD (φ, ψ) is called the Kantorovich dual problem: where φ and ψ satisfy the following condition: Strong duality is known to hold for the Kantorovich problem and its dual problem.
Proposition 1 ([13, Theorem 1.3]). For the Kantorovich problem solution p * (x, y) and the dual problem solutions φ * (x) and ψ * (y), the following equation holds: Thus, the problem of finding the simultaneous distribution of p in (1) is replaced by the problem of finding the functions φ and ψ in (4). In fact, it is known that we only need to find one of the two functions.
Proposition 2 ([13, Remark 1.12]). In the Kantorovich dual problem, the following equality holds for the pair of functions φ * and ψ * that maximize the cost function C KD (φ, ψ).
In the next section we focus on the optimal transport problem where the distribution is restricted to a particular family.

III. PROPOSED CONTROL METHOD
This section provides the optimal transport-based coverage control with the aid of the idea in the previous section. We begin by considering n ∈ N mobile robots located on R d space. Let the position of the i-th robot at time t ∈ R + be x i (t) ∈ R d and let its dynamics be given as follows: where u i (t) ∈ R d is the input that determines the speed of the robot. Next, we define a distribution formed by robots as where we define δ(·) as a Dirac's delta function. We focus on the problem of minimizing the distance (measured by W ) between the target distribution ρ T and the distribution ρ 0 (x) = ρ(x, t) at each time: where we define x := [x 1 , . . . , x n ] as a position vector of robots. Our goal is to design inputs u i (t) (∀i ∈ {1, . . . , n}) that achieves (10) at each time for the system of (8).
The results provided in the previous section yield the following expression for the problem of (10).
Proposition 3. The following equation holds: where we define a real vector φ = [φ 1 , . . . , φ n ] ∈ R n , and we define a function F : In addition, the set V φ i (x) is called Laguerre regions: Proof. We show that the function W is transformed into (12) by using the Kantorovich duality. The first term in (4) is calculated as Using (7), the second term in (4) is rearranged as Finally, (6) ensures that (11) holds.
Thus, the optimization problem to be solved by the robot swarm is rearranged as the following min-max problem: Accordingly, we propose the following controller as a solution to the problem of (16). (8), the optimal transport-based coverage control (OTCC) is defined as

Definition 3. For the system in
where k ∈ R + and k ∈ R + are design parameters that provide the feedback gain, andã i (t) andb i (t) are defined as We define and consider the following time derivative of the function F in (12) along the trajectories x(t) and φ(t) of the system controlled via the OTCC:Ḟ The following proposition then justifies the conclusion that the OTCC provides a solution to the problem in (16).
Proposition 4. The following hold for any t ∈ R + : Proof. We calculate the gradients of the function F as where we used Reynolds' transport theorem [24] to differentiate the functions including variables in the integration domain. The time derivative of F along the trajectories x and φ are evaluated aṡ which ends the proof.
Remark 1. We show that the Voronoi tessellation-based coverage control (VTCC) method [11] is regarded as a special case of the proposed OTCC. In [11], the following cost function is introduced: where the set V i (x) is called the Voronoi region: In their study, they proposed the following control law to minimize the cost function of (28): The OTCC in (17) and (18) agree with the VTCC of (30) when the variable φ is set to φ(t) ≡ 0. Under this condition, the cost function for both methods ( (12) and (28)) coincide. Therefore, the VTCC is one of the gradient flows that realize the transport from ρ(x, t) to ρ T (x). However, the VTCC is not optimal from the aspect of optimal transport because the cost function is not maximized for φ at each time. In contrast, the proposed OTCC overcomes this problem and is expected to provide better control performance.
Remark 2. We discuss the difference between the Voronoi region in (29) and the Laguerre region in (13). In both sets, the boundary is perpendicular to the line between neighboring robots i and j. While the boundary in the Voronoi region is a bisector, it is not a bisector in the Laguerre region if φ is not zero. Specifically, when φ i is larger than φ j , the boundary moves towards the robot j, and consequently robot i acquires more area.

Remark 3.
Under suitable conditions, the proposed algorithm is scalable in the sense that each robot only needs the information of neighboring robots. Indeed, the controllers in (17) and (18) are computed using x j and φ j of the robots ∅} is the set of adjacent robots in the Laguerre sense. Hence, we focus on clarifying the condition that each robot is capable of obtaining the information of the neighboring robots. Suppose that each robot i has a measurement range R i and obtains the information of robots j satisfying x j − x i ≤ R i . Then, the proposed algorithm is feasible if R i satisfies x j − x i ≤ R i for any i and j ∈ N φ i . Such a range R i does not increase in a common situation where many robots are distributed on the workspace and |φ i − φ j | (for any i, j) is sufficiently smaller than the size of the workspace.

IV. STABILITY ANALYSIS
In this section, we prove the Lyapunov stability of the system's equilibrium when using our proposed OTCC. Suppose that the dynamical system represented by (8) are controlled with the OTCC in (17) and (18). The system equilibrium (x * , φ * ) is then characterized as a point that satisfies the following conditions: where we assume that x * i = x * j holds for all i = j.
where H(x, φ) ∈ R nd×nd is defined as Here, (H(x, φ)) ij denotes (i, j)-block matrix element of H(x, φ), and I d ∈ R d×d denotes d-dimensional identity matrix.
We define ∂V φ ij := V φ i (x) ∩ V φ j (x) as the Laguerre region boundary of i and j, and the integral in (35) is 0 when ∂V φ ij is an empty set.
We prove Theorem 1 by utilizing following proposition.

Proposition 5 ([25]
). For a continuous and differentiable function G : R n × R n → R, consider a system with the following dynamics:ṗ If G is convex-concave around the saddle point (p * , q * ), then the system of (36) is Lyapunov stable at (p * , q * ).
Here, the convex-concave function and the saddle point are defined below.
Definition 4 (Convex-concave function). A function g : R n → R is convex aroundz ∈ R n if a neighborhood Uz exists and the inequality of λg(x) + (1 − λ)g(y) ≥ g(λx + (1 − λ)y) holds for any x, y ∈ Uz and λ ∈ [0, 1]. The function g is concave aroundz if the inverse inequality holds. Furthermore, the function G : R n × R n → R is convex-concave around (p,q) if p → G(p,q) is convex aroundp and q → G(p, q) is concave aroundq.
Definition 5 (Saddle point). For a continuous and differentiable function G : R n ×R n → R, a point (p * , q * ) ∈ R n ×R n is a saddle point if there exists neighbors U p * and U q * , and the following holds for any p ∈ U p * and q ∈ U q * : Proposition 5 ensures that it is sufficient to show that the function F in (12) is convex-concave around the equilibrium (x * , φ * ) of (33), and that (x * , φ * ) is the saddle point of F . The former is guaranteed with the aid of the following two propositions: Proof. We denote the above function as F x * (φ). In order to use the second-order sufficient condition for the concave function [26], we evaluate ∇ 2 φ F x * (φ), the second-order derivative of F x * (φ) for φ. Using Reynolds' transport theorem, it is shown that Equation (38) is a weighted graph Laplacian multiplied by −1. Thus, ψ ∇ 2 φ F x * (φ)ψ ≤ 0 holds for any φ, ψ ∈ R n , which means that ∇ 2 φ F x * (φ) is negative semidefinite at any point φ ∈ R n . Proposition 7. If the condition of (34) holds, the function Proof. We denote the above function as F φ * (x). Using Reynolds' transport theorem, the second-order derivative of (34) holds. The equilibrium is directly shown to be the saddle point using the convex-concavity. Proposition 8. If the condition of (34) holds, (x * , φ * ) is the saddle point of the function F of (12).
Proof. We show that (37) holds. By using Taylor's formula, the following holds for the equilibrium point (x * , φ * ) and any φ in the neighborhood of φ * : whereφ is a point between φ * and φ. The inequality comes from the definition of the equilibrium point in (33) and the concavity of F x * (φ * ). The other inequality in (37) is also shown by considering Taylor's formula in the neighborhood of x * .
Finally, Theorem 1 is proved by using above propositions. Proof of Theorem 1: Propositions 6 and 7 show that the function F of (12) is convex-concave around (x * , φ * ) when (34) is satisfied. In addition, Proposition 8 ensures that (x * , φ * ) is the saddle-point of F . As a result, Proposition 5 guarantees that (x * , φ * ) is Lyapunov stable.
By limiting the space dimension to d = 1, we reduce the condition in Theorem 1 to the following explicit form.
Theorem 2. Suppose that d = 1 holds for the space dimension. Then, (x * , φ * ) is Lyapunov stable provided that the following inequality holds for any i ∈ {1, . . . , n}: where the function h i : R n × R n → R is defined as In (41), we define The proof of Theorem 2 is obvious through the following proposition.
Proposition 9. Suppose that d = 1 holds. If the inequality of (40) holds for any i ∈ {1, . . . , n}, the function Proof. Using d = 1, we assume x i < x i for i < i without a loss of generality. Then, ∂V φ ij exists as (35) is then calculated as follows: where we denote k := i − 1 and h := i + 1. For j = i, By using Gershgorin theorem [27], all eigenvalues of ∇ 2 x F φ * (x) lie within a disk with a center of (∇ 2 x F φ * (x)) ii and a radius of j (∇ 2 x F φ * (x)) ij , so that the following is obtained as a sufficient condition for all eigenvalues to be positive: Therefore, in order for ∇ 2 x F φ * (x) to be positive definite at the equilibrium x * , it is sufficient that (40) holds, where we use the fact thatã i = 1/n holds at the equilibrium. From the assumption that x * i = x * j (∀i = j), we see that both sides of (44) are continuous functions. Thus, if (44) holds at the equilibrium x * , it remains in its neighborhood, which shows that the function F is convex around x * .

V. NUMERICAL EXPERIMENTS
In this section, numerical analysis is conducted to verify the performance of the proposed OTCC. We first consider a one-dimensional case with 40 robots uniformly distributed over the interval [−10, −5]. The target density distribution ρ T is set as ρ T (x) = N (0, 3), where N (µ, σ 2 ) denotes the density function for the normal distribution with a mean µ and a variance σ 2 . We use k = 0.5 for the feedback gain in the VTCC, and k = 0.5 and k = 1.0 × 10 −4 as the feedback gains in the proposed OTCC. The integrations in (19) and (20) are carried out numerically by restricting the robot workspace to the interval of [−10, 10] and discretizing it into small cells. Each robot then determines the ownership of each cell based on the definition of the Laguerre regions and performs numerical integration over the area belonging to the robot. The two control methods are performed for the above system and the trajectory of the robots' position is shown in Fig. 1. In both methods, the robots that were centered around x = −7.5 at time t = 0 move over time to present the desired distribution ρ T centered at x = 0. However, in the VTCC, many robots remain in the region of x < 0, and the mean value of the distribution does not approach to x = 0. In contrast, the OTCC reproduces the shape of the target distribution ρ T more closely over time. For more quantitative analysis, we examine the value of the cost function in (28) in the steady state. The value of the cost function in the VTCC is 1.36 × 10 −3 , whereas the value in the proposed OTCC is 1.08 × 10 −3 , which suggests that the latter provides better control. To check whether the conditions of Theorem 2 are satisfied, we next examine the value of the left-hand side of (40) in the steady state. The maximum value for all robots is 5.69 × 10 −1 (< 1), which indicates that the stationary point is Lyapunov stable.
Next, we consider a two-dimensional case, with 25 robots uniformly distributed over the interval of [−8, −2] 2 and with an interval of [−10, 10] 2 for robot workspace. We use ρ T (x) = 1 2 {N (µ 1 , Σ) + N (µ 2 , Σ)} as the target density distribution, while µ 1 = [−5, −5] , µ 2 = [5,5] , and Σ = diag(4, 4) are used for the mean and covariance. We set k = 0.5 for the feedback gain in the VTCC, and k = 0.5 and k = 5.0 × 10 −2 for the feedback gains in the proposed OTCC. Fig. 2 shows the visualized trajectory of the positions of the robots using two control laws. The black circles represent the positions of the robots, and the areas painted in different colors represent the Voronoi/Laguerre region to which each robot belongs. In both methods, the robots move in a way that reproduces the desired distribution ρ T . In the VTCC, only 6 of the 25 robots move to the upper-right distribution, while the remaining robots stay in the lower-left distribution. In contrast, in the OTCC, 12 of the 25 robots move to the upper-right distribution, thus suggesting that the shape of the target distribution ρ T is better reproduced. We calculate the value of the cost function of (28) at the steady state. The value in the VTCC is 9.70 × 10 −1 , whereas the value in the proposed OTCC is 8.84 × 10 −1 , thus indicating the distribution is better reproduced in the proposed method.
Remark 4. We discuss the reasons of higher performance of the OTCC. This is because that the equilibrium condition in the OTCC is stricter than that in the VTCC, which makes the robot less likely to be trapped in the stationary point. The equilibrium condition in the VTCC is that the robot is placed at the center of gravity of each region (x i = b i ), and once this condition is satisfied, each robot does not move thereafter. In the OTCC, however, there is an additional condition that the weighted area of each region is equal (ã i = 1/n in (33)). Thus, even if the former condition is satisfied, robots continue to move unless the latter condition is satisfied, which is why robots are less likely to be trapped at the stationary point.

VI. CONCLUSION
In this paper, we propose the optimal transport-based coverage control (OTCC) method as an improvement to the Voronoi tessellation-based coverage control (VTCC) method. Our proposed method, which is derived via the Kantorovich dual problem, is consistent with the VTCC with setting φ(t) ≡ 0. This correspondence successfully reconsiders the VTCC in the optimal transport framework. We also derive the conditions of Lyapunov stability for the controlled system in Theorem 1 and 2, showing that once the robot swarm reaches the equilibrium point, they remain in the neighborhood. Numerical calculations clearly show that our proposed method is more capable of escaping the local optimum point and achieving better control than the VTCC. As one of our next research topics, we plan to