Paper deep dive
Density-Driven Optimal Control: Convergence Guarantees for Stochastic LTI Multi-Agent Systems
Kooktae Lee
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 95%
Last extracted: 4/10/2026, 2:01:24 AM
Summary
The paper introduces Stochastic Density-Driven Optimal Control (D2OC), a decentralized Lagrangian framework for non-uniform area coverage in stochastic LTI multi-agent systems. By formulating a stochastic MPC-like problem that minimizes the Wasserstein distance as a running cost, the method ensures convergence to a target density while accounting for process and measurement noise. The approach provides formal convergence guarantees via reachability analysis and outperforms existing heuristic methods in optimality and consistency.
Entities (4)
Relation Signals (2)
Stochastic D2OC → controls → Multi-Agent Systems
confidence 95% · we propose Stochastic Density-Driven Optimal Control (D2OC)... for non-uniform coverage in stochastic linear time-invariant (LTI) multi-agent systems.
Stochastic D2OC → minimizes → Wasserstein Distance
confidence 95% · By formulating a stochastic MPC-like problem that minimizes the Wasserstein distance as a running cost
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper addresses the decentralized non-uniform area coverage problem for multi-agent systems, a critical task in missions with high spatial priority and resource constraints. While existing density-based methods often rely on computationally heavy Eulerian PDE solvers or heuristic planning, we propose Stochastic Density-Driven Optimal Control (D$^2$OC). This is a rigorous Lagrangian framework that bridges the gap between individual agent dynamics and collective distribution matching. By formulating a stochastic MPC-like problem that minimizes the Wasserstein distance as a running cost, our approach ensures that the time-averaged empirical distribution converges to a non-parametric target density under stochastic LTI dynamics. A key contribution is the formal convergence guarantee established via reachability analysis, providing a bounded tracking error even in the presence of process and measurement noise. Numerical results verify that Stochastic D$^2$OC achieves robust, decentralized coverage while outperforming previous heuristic methods in optimality and consistency.
Tags
Links
- Source: https://arxiv.org/abs/2604.08495v1
- Canonical: https://arxiv.org/abs/2604.08495v1
Trouble viewing inline? Open PDF directly →
Full Text
46,922 characters extracted from source content.
Expand or collapse full text
Density-Driven Optimal Control: Convergence Guarantees for Stochastic LTI Multi-Agent Systems Kooktae Lee kooktae.lee@nmt.edu Department of Mechanical Engineering, New Mexico Institute of Mining and Technology, Socorro, NM 87801, USA Abstract This paper addresses the decentralized non-uniform area coverage problem for multi-agent systems, a critical task in missions with high spatial priority and resource constraints. While existing density-based methods often rely on computationally heavy Eulerian PDE solvers or heuristic planning, we propose Stochastic Density-Driven Optimal Control (D2OC). This is a rigorous Lagrangian framework that bridges the gap between individual agent dynamics and collective distribution matching. By formulating a stochastic MPC-like problem that minimizes the Wasserstein distance as a running cost, our approach ensures that the time-averaged empirical distribution converges to a non-parametric target density under stochastic LTI dynamics. A key contribution is the formal convergence guarantee established via reachability analysis, providing a bounded tracking error even in the presence of process and measurement noise. Numerical results verify that Stochastic D2OC achieves robust, decentralized coverage while outperforming previous heuristic methods in optimality and consistency. keywords: Multi-Agent Systems; Stochastic Control; Optimal Transport; Wasserstein Distance; Area Coverage. †thanks: This work was supported by NSF CAREER Grant CMMI-DCSD-2145810. 1 Introduction The multi-agent area coverage problem has recently garnered significant attention due to its broad range of applications, including search and rescue, environmental monitoring, infrastructure inspection, smart farming, and planetary exploration. The central challenge lies in guiding a team of agents to maximize coverage performance within a given domain. While traditional strategies like the lawnmower path have been widely used for uniform coverage, these methods are often inefficient in large-scale environments and may become infeasible when resources such as operation time, agent number, and communication range are constrained. Consequently, non-uniform area coverage has emerged as a more practical alternative, especially under limited resources and mission-specific priorities. Recent approaches to non-uniform area coverage aim to match agent behavior with a reference distribution encoding spatial priorities. Spectral Multiscale Coverage (SMC) [10, 11] leverages Fourier-based metrics to align long-term visitation frequencies based on ergodic control. However, ergodicity is only theoretically achieved as t→∞t→∞, making it ill-suited for missions with strict finite-time constraints. To address this, density-based optimal control frameworks have been explored. For instance, [1] and [13] study optimal transport over dynamical systems using Eulerian-based gradient flows or dynamic programming in probability spaces. While theoretically rigorous, these methods often rely on centralized computations of a global density flow and are primarily restricted to linear dynamics or simplified state spaces, which are computationally prohibitive for high-dimensional decentralized swarms. More recently, [12] employs mean-field Schrödinger bridges with Gaussian Mixture Models (GMM) to steer agent populations. However, such approaches are limited by the parametric assumptions of GMMs and the high complexity of solving coupled PDEs, making them less flexible for non-parametric target densities. Density-Driven Control (D2C) has emerged as a practical Lagrangian-based alternative for matching empirical agent distributions with a target density [6, 9]. By adopting a Lagrangian perspective, where the swarm is represented as a collection of discrete, identifiable particles rather than a continuous field, D2C effectively circumvents the curse of dimensionality and the heavy computational burden associated with solving high-dimensional PDEs in Eulerian density control. Despite its computational efficiency and scalability, the reliance on heuristic planning often limits its performance and optimality, particularly in complex environments with high uncertainty. Specifically, existing D2C frameworks frequently overlook the dynamic coupling between individual agent trajectories and the collective distribution matching objective under stochastic perturbations. This gap necessitates a unified decentralized framework capable of simultaneously accounting for strict physical operational constraints and the inherent stochasticity arising from process and measurement noises. To bridge this gap, we propose the Stochastic D2OC, a rigorous optimal control framework for non-uniform coverage in stochastic linear time-invariant (LTI) multi-agent systems. Unlike previous heuristic methods, our framework reformulates distribution matching within the Optimal Transport (OT) context by minimizing the Wasserstein distance as a running cost in a stochastic MPC-like formulation. This leads to an optimal control law that ensures the time-averaged spatial distribution converges to a reference density in a provably consistent manner, even under severe noise. Furthermore, by characterizing the reachable set of the agent dynamics, we provide a formal convergence analysis ensuring the empirical distribution remains within a bounded neighborhood of the target. Numerical simulations under aggressive noise and perturbations demonstrate that Stochastic D2OC achieves robust, decentralized matching while strictly satisfying individual dynamics and physical constraints. 2 Preliminaries Notation: Let ℝR and ℤZ be the sets of real and integer numbers, with ℤ>0Z_>0 and ℤ≥0Z_≥ 0 denoting positive and non-negative integers, respectively. For a matrix A∈ℝm×nA ^m× n, its transpose is A⊤A . The Euclidean norm is ‖\|x\|, and the trace is tr(A)tr(A). The identity and zero matrices of size n are nI_n and n0_n. A Gaussian distribution with mean μ and covariance Σ is (μ,Σ)N(μ, ), and ‖U‖R:=U⊤RU\|U\|_R:= U RU denotes the weighted norm for R≻0R 0. The operators diag(⋅)diag(·) and blkdiag(Ah)h=r+H−1blkdiag(A_h)_h=r^r+H-1 construct diagonal and block-diagonal matrices, respectively. Finally, ⊗ and ⊙ represent the Kronecker and Hadamard products. Consider a multi-agent system, where each agent is governed by discrete-time stochastic LTI dynamics, evolving over a discrete-time index k∈ℤ≥0k _≥ 0, as follows: xik+1 x_i^k+1 =Aixik+Biuik+wik,wik∼(0,Σi,w), =A_ix_i^k+B_iu_i^k+w_i^k, w_i^k (0, _i,w), (1) yik y_i^k =Cixik+vik,vik∼(0,Σi,v), =C_ix_i^k+v_i^k, v_i^k (0, _i,v), where xik∈ℝnx_i^k ^n is the state, uik∈ℝmu_i^k ^m is the control input, and yik∈ℝdy_i^k ^d is the output at time k. The matrices Ai∈ℝn×nA_i ^n× n, Bi∈ℝn×mB_i ^n× m, and Ci∈ℝd×nC_i ^d× n define the system dynamics, control influence, and observation model, respectively. The process noise wikw_i^k and measurement noise vikv_i^k are independent, zero-mean Gaussian random variables with covariances Σi,w _i,w and Σi,v _i,v, capturing model uncertainty and sensing errors. To ensure well-posedness of the control problems considered in this paper, we impose the following standard assumption on each agent’s dynamics. Assumption 1. For each agent i, the pair (Ai,Bi)(A_i,B_i) is completely controllable. Assumption 2. For each agent i, AiA_i is at least marginally stable, i.e., all eigenvalues lie in the closed unit disk, and any eigenvalue on the unit circle is non-defective (i.e., the algebraic and geometric multiplicities are equal). These assumptions ensure that each agent is fully controllable and that its open-loop dynamics are at least marginally stable, which, together with bounded noise covariance, forms the basis for the convergence analysis. 2.1 Wasserstein Distance and the Kantorovich Optimal Transport Problem To develop the D2OC strategy, we utilize optimal transport theory [14]. The p-Wasserstein distance between two discrete measures ρ and ν on a metric space (,d)(X,d) is defined via the Kantorovich problem: p(ρ,ν)=(minπij≥0∑i=1M∑j=1Nπijd(yi,qj)p)12,W_p(ρ,ν)= ( _ _ij≥ 0 _i=1^M _j=1^N _ij\,d(y_i,q_j)^p ) 12, (2) subject to ∑j=1Nπij=αi _j=1^N _ij= _i, ∑i=1Mπij=βj _i=1^M _ij= _j, and ∑i,jπij=1 _i,j _ij=1, where πij _ij is the transport plan from agent points yi∈ℝdy_i ^d to fixed sample points qj∈ℝdq_j ^d with masses αi,βj _i, _j. We set p=2p=2 and use Euclidean distance for d(⋅,⋅)d(·,·). To facilitate D2OC, we distinguish between evolving agent trajectories yiky_i^k and fixed reference locations qjq_j. Each of the nan_a agents generates MiM_i points over a finite horizon, totaling M=∑MiM=Σ M_i points. Each agent i maintains its own capacity αik=1/Mi _i^k=1/M_i and locally estimates sample capacities βi,jk _i,j^k (initially 1/N1/N). These capacities βi,jk _i,j^k are updated locally to reflect coverage history, i.e., βi,jk _i,j^k decreases as agent i visits or passes near qjq_j, effectively tracking explored regions. To evaluate global density-driven coverage, we define the time-averaged empirical distribution ρkρ^k and the reference distribution ν as: ρk:=1(k+1)na∑t=0k∑i=1naδyit,ν:=1N∑j=1Nδqj,ρ^k:= 1(k+1)n_a _t=0^k _i=1^n_a _y_i^t, ν:= 1N _j=1^N _q_j, (3) where δy _y is the Dirac measure at y. The objective of D2OC is to coordinate agents to minimize the global Wasserstein distance 2(ρk,ν)W_2(ρ^k,ν) within a finite operational time, achieving effective and robust decentralized non-uniform coverage. 2.2 Three-stage D2OC Overview Although the core methodology presented here differs fundamentally from prior works [6, 7, 9], which lack guarantees on optimality or convergence, the overall structure for realizing D2OC follows a similar high-level logic, summarized in the following three stages: 1. Local Target Sample Selection and Optimal Control: At each time step k, each agent selects some sample points among the reference samples qj\q_j\ based on proximity and remaining weight (favoring less-visited points). Each agent then computes a control input using the selected local target samples, to minimize the local Wasserstein distance, subject to dynamic and input constraints. 2. Weight Update: Next time step after moving, the agent solves a local Wasserstein distance to update sample weights, reducing those recently covered and encouraging coverage of under-explored areas. 3. Weight Sharing: Agents exchange sample weights within communication range, synchronizing coverage estimates via a min-weight consensus that selects the smallest weight for each sample point. This ensures agreement on maximal coverage progress, enhancing coordination and reducing redundant visits. The first two stages are executed independently by each agent without communication, while the final stage enables collaborative coverage via information exchange among neighbors. This iterative cycle guides agents to match the reference density over time within dynamic and communication constraints. Although this paper primarily focuses on the first stage and theoretical convergence analysis of D2OC, further details about three-stage approach can be found in [6, 7, 9]. 3 Optimal Control Minimizing Local Wasserstein Distance This section presents the optimal control law for the first stage of the three-stage D2OC framework. Let ik⊆1,…,NS_i^k \1,…,N\ be the index set of sample points assigned to agent i, and ik=qj∈ℝd∣j∈ikQ_i^k=\q_j ^d j _i^k\ the corresponding local targets with capacity weights βi,jk≥0 _i,j^k≥ 0. The selection of local targets is detailed later based on reachability analysis. Definition 1. (Output Relative Degree in Stochastic Discrete-Time LTI Systems) Consider the stochastic discrete-time LTI system (1). The output relative degree r∈ℤ>0r _>0 is the smallest positive integer such that CiAir−1Bi≠0C_iA_i^r-1B_i≠ 0, and CiAiℓ−1Bi=0,∀ℓ=1,…,r−1.C_iA_i -1B_i=0,\,∀ =1,…,r-1. This means the control input affects the output starting at step k+rk+r. Thus, the expected squared local Wasserstein distance for agent i is constructed from time k+rk+r with the prediction horizon H∈ℤ>0H _>0 as follows: [∑h=rH+r−1(ik+h)2]:=∑h=rH+r−1∑j∈ik+hπjk+h[‖yik+h−qj‖2], [ _h=r^H+r-1(W_i^k+h)^2 ]:= _h=r^H+r-1 _j _i^k+h _j^k+h\,E [\|y_i^k+h-q_j\|^2 ], where πjk+h _j^k+h is the transport weight from agent i to local sample qjq_j given by the local optimal transport plan at time k+hk+h. This captures stochastic system dynamics and expected transport cost from predicted agent positions to assigned targets. Leveraging this property, we establish the following result. Relationship to [8]: The following quadratic reformulations and optimal control laws (see Proposition 1 and Theorem 1) maintain notation consistency with [8] while providing a fundamental stochastic generalization. Unlike the deterministic study in [8], explicitly incorporating noise statistics via the expectation operator is a critical prerequisite for the rigorous convergence analysis in Section 5. Proposition 1. Let ik+hS_i^k+h denote the index set for the local sample points for agent i at time k and the optimal transport plan πjk+h≥0 _j^k+h≥ 0 be given for all j∈ik+hj _i^k+h and for h=r,…,H+r−1h=r,…,H+r-1. Define the weighted barycenter at time k+hk+h as q¯ik+h:=1∑j∈ik+hπjk+h∑j∈ik+hπjk+hqj. q_i^k+h:= 1 _j _i^k+h _j^k+h _j _i^k+h _j^k+hq_j. Then, the following equality holds: [∑h=rH+r−1(ik+h)2]=[‖ik|r:H(Yik|r:H−Q¯ik|r:H)‖2] [ _h=r^H+r-1 (W_i^k+h )^2 ]=E [ \| _i^k|r:H(Y_i^k|r:H- Q_i^k|r:H) \|^2 ] +const., +const., where Yik|r:H:=[(yik+r)⊤⋯(yik+H+r−1)⊤]⊤∈ℝdH, \,Y_i^k|r:H:= bmatrix(y_i^k+r) &·s&(y_i^k+H+r-1) bmatrix ^dH, Q¯ik|r:H:=[(q¯ik+r)⊤⋯(q¯ik+H+r−1)⊤]⊤∈ℝdH, Q_i^k|r:H:= bmatrix( q_i^k+r) &·s&( q_i^k+H+r-1) bmatrix ^dH, ik|r:H:=blkdiag(∑j∈i(k+h)πj(k+h)d)h=r+H−1. _i^k|r:H:=blkdiag\! ( _j _i(k+h)\! _j(k+h)\,I_d )_h=r^r+H-1. (4) Proof. By expanding the quadratic term, [∑h=rH+r−1(ik+h)2]=[∑h=rH+r−1∑j∈ik+hπjk+h‖yik+h−qj‖2] [ _h=r^H+r-1 (W_i^k+h )^2 ]=E [ _h=r^H+r-1 _j _i^k+h _j^k+h\|y_i^k+h-q_j\|^2 ] =[∑h=rH+r−1(∑jπjk+h)‖yik+h−q¯ik+h‖2+Ch] =E [ _h=r^H+r-1 ( _j _j^k+h )\|y_i^k+h- q_i^k+h\|^2+C_h ] =[∑h=rH+r−1(∑jπjk+h)‖yik+h−q¯ik+h‖2]+const. =E [ _h=r^H+r-1 ( _j _j^k+h )\|y_i^k+h- q_i^k+h\|^2 ]+const. where each Ch:=∑jπjk+h‖qj−q¯ik+h‖2C_h:= _j _j^k+h\|q_j- q_i^k+h\|^2 is constant w.r.t. Uik|HU_i^k|H. Stacking the terms into vectors and using the diagonal matrix ik|r:H _i^k|r:H yields the compact quadratic form. ∎ Remark 1 (Local Transport Plan). The transport plan πjk+h _j^k+h is an internal variable computed locally by agent i to determine the optimal mass distribution to target samples qj\q_j\. Given a constant mass αi=1/Mi _i=1/M_i allocated at each step, the agent independently solves a matching problem to minimize the local Wasserstein distance. The existence and uniqueness of this optimal πjk+h _j^k+h are guaranteed by the optimal matching strategy established in [7, Prop. 1]. The control objective is to minimize the expected squared Wasserstein cost along with a weighted input penalty: minUik|HJ(Uik|H):=[∑h=rH+r−1(ik+h)2]+‖Uik|H‖R2 _U_i^k|HJ(U_i^k|H):=E [ _h=r^H+r-1(W_i^k+h)^2 ]+\|U_i^k|H\|^2_R (5) where Uik|H∈ℝmHU_i^k|H ^mH is the stacked input and R≻0R 0 is a positive definite weighting matrix. To avoid the intractability of global Wasserstein minimization, we focus on the local distance ik+hW_i^k+h. Based on Proposition 1, the cost in (5) is reformulated into a strictly convex quadratic form using the following matrix definitions: Θi _i :=[CiAir−1Bidm⋯dmCiAirBiCiAir−1Bi⋯dm⋮⋱⋮CiAir+H−2Bi⋯CiAir−1Bi], := bmatrixC_iA_i^r-1B_i&0_dm&·s&0_dm\\ C_iA_i^rB_i&C_iA_i^r-1B_i&·s&0_dm\\ & & & \\ C_iA_i^r+H-2B_i&·s&·s&C_iA_i^r-1B_i bmatrix, (6) Φi _i :=[(CiAir)⊤,…,(CiAir+H−1)⊤]⊤. :=[(C_iA_i^r) ,…,(C_iA_i^r+H-1) ] . (7) Specifically, for a reference barycenter Q¯ik|r:H Q_i^k|r:H, the cost J is equivalently expressed as J=(Uik|H)⊤HiUik|H+2fi⊤Uik|H,where Hi=(iΘi)⊤(iΘi)+R,fi=(iΘi)⊤i(Φi[xik]−Q¯ik|r:H), split&J=(U_i^k|H) H_iU_i^k|H+2f_i U_i^k|H,\\ &where H_i=( _i _i) ( _i _i)+R,\\ & f_i=( _i _i) _i( _iE[x_i^k]- Q_i^k|r:H), split (8) with i:=ik|r:H _i:= _i^k|r:H. This quadratic reformulation leads to the following optimal control law. Theorem 1 (Strictly Convex Optimal Control Law). Consider the stochastic discrete-time LTI system (1) for agent i. For the cost function J(Uik|H)J(U_i^k|H) defined in (5), the unique unconstrained optimal control sequence Uik|H⋆U_i^k|H that minimizes the local Wasserstein-based cost is given by Uik|H⋆=−Hi−1fi,U_i^k|H =-H_i^-1f_i, (9) where HiH_i and fif_i are as defined in (8). If box constraints on the control input umin≤uik+τ≤umax,τ=0,…,H−1,u_ ≤ u_i^k+τ≤ u_ , τ=0,…,H-1, are imposed, then the constrained optimization becomes the strictly convex quadratic program: minUik|H _U_i^k|H (Uik|H)⊤HiUik|H+2fi⊤Uik|H (U_i^k|H) H_iU_i^k|H+2f_i U_i^k|H (10) subject to umin(H)≤Uik|H≤umax(H)u_ ^(H)≤ U_i^k|H≤ u_ ^(H), where umin(H):=H⊗umin,umax(H):=H⊗umaxu_ ^(H):=1_H u_ ,\,u_ ^(H):=1_H u_ with H∈ℝH1_H ^H is the all-ones vector. The unique optimal solution (Uik|H)∗(U_i^k|H)^* satisfies the Karush-Kuhn-Tucker (KKT) conditions: Hi(Uik|H)∗+fi+λ+−λ−=0,λ+,λ−≥0, H_i(U_i^k|H)^*+f_i+λ^+-λ^-=0, λ^+,λ^-≥ 0, λ+⊙((Uik|H)∗−umax(H))=0,λ−⊙(umin(H)−(Uik|H)∗)=0. λ^+ ((U_i^k|H)^*-u_ ^(H))=0,\,λ^- (u_ ^(H)-(U_i^k|H)^*)=0. Proof. From Proposition 1 with [Yik|r:H]=ΘiUik|H+Φi[xik]E[Y_i^k|r:H]= _iU_i^k|H+ _iE[x_i^k], the Wasserstein cost term can be expanded as ∑h=rH+r−1(ik+h)2=tr(ik|r:HΣi,Y(ik|r:H)⊤) _h=r^H+r-1(W_i^k+h)^2=tr( _i^k|r:H _i,Y( _i^k|r:H) ) +‖ik|r:H(ΘiUik|H+Φixik−Q¯ik|r:H)‖2+const, +E\| _i^k|r:H( _iU_i^k|H+ _ix_i^k- Q_i^k|r:H)\|^2+const, where Σi,Y:=diag([CiΣi,wCi⊤+Σi,v,…,CiΣi,wCi⊤+Σi,v])∈ℝdH×dH _i,Y:=diag([C_i _i,wC_i + _i,v,…,C_i _i,wC_i + _i,v]) ^dH× dH is the block-diagonal output noise covariance over the prediction window. Since the trace and constant terms are independent of Uik|HU_i^k|H, it is omitted from the optimization. Including the weighted input penalty term ‖Uik|H‖R2\|U_i^k|H\|_R^2, the cost becomes J(Uik|H)=(Uik|H)⊤HiUik|H+2fi⊤Uik|H,J(U_i^k|H)=(U_i^k|H) H_iU_i^k|H+2f_i U_i^k|H, where HiH_i and fif_i are defined in (8). By Assumption 1, the pair (Ai,Bi)(A_i,B_i) is controllable, which implies that the matrix Θi _i has full column rank or sufficient rank properties ensuring that (ik|r:HΘi)⊤(ik|r:HΘi)( _i^k|r:H _i) ( _i^k|r:H _i) is positive semidefinite. Adding R≻0R 0 then guarantees that Hi≻0.H_i 0. Consequently, the cost function is strictly convex and admits a unique global minimizer obtained from the first-order condition ∇Uik|HJ=2HiUik|H+2fi=0, _U_i^k|HJ=2H_iU_i^k|H+2f_i=0, which yields the unconstrained solution in (9). When box constraints are included, the problem becomes a strictly convex quadratic program. The KKT conditions are both necessary and sufficient for optimality and yield the unique constrained solution (Uik|H)∗(U_i^k|H)^*. ∎ Corollary 1. (Optimal Solution with Per-Step Ball Input Constraints ‖uik+τ‖≤umax\|u_i^k+τ\|≤ u_ ) Consider the optimization problem (5) over a finite horizon H, subject to per-step Euclidean norm constraints ‖uik+τ‖≤umax,∀τ=0,…,H−1,\|u_i^k+τ\|≤ u_ , ∀τ=0,…,H-1, where the system evolves according to (1). Then, the optimal solution (Uik|H)∗(U_i^k|H)^* is obtained by projecting each block uik+τ,unconu_i^k+τ,uncon of the unconstrained solution (Uik|H)uncon(U_i^k|H)^uncon onto the Euclidean ball of radius umaxu_ : uik+τ,∗=min(1,umax‖uik+τ,uncon‖)uik+τ,uncon,u_i^k+τ,*= (1, u_ \|u_i^k+τ,uncon\| )u_i^k+τ,uncon, where uik+τ,unconu_i^k+τ,uncon is the (k+τ)(k+τ)-th m-dimensional block of the unconstrained solution defined in (9). This result leverages the separable nature of the per-step Euclidean norm constraints, enabling a closed-form projection of the unconstrained solution without the need for iterative solvers. The proposed density-driven optimal control law is implemented via a model predictive control (MPC) scheme. At each time step, an optimal control sequence over a finite horizon H is computed, with only the first input applied before re-optimizing at the next step using updated agent states and local sample information. Remark 2. To track the reference distribution, each agent updates its local target samples and corresponding barycenter at every step. However, frequent updates can induce oscillations due to shifting targets. This is mitigated by the receding-horizon structure, which tempers responsiveness, and by an input regularization term that smooths motion while preserving adaptability to local changes. Since the effectiveness of D2OC depends critically on the choice of local target samples and their barycenter, the next section presents a detailed analysis based on the mean reachable set. 4 Reachability-Aware Target Sample Selection and Adjustment For systems with output relative degree r, control inputs do not affect the output until time k+rk+r. Thus, mean reachability analysis over h=r,…,H+r−1h=r,…,H+r-1 is crucial for selecting local target points and computing their barycenter. This section presents a framework for constructing mean reachable output sets, selecting target samples accordingly, and computing a reachable barycenter that reflects future agent capabilities. An adjustment procedure is also provided for cases where the barycenter lies outside the reachable set. 4.1 h-Step Reachable Set Analysis Given the current state xikx_i^k at time k, and assuming an output relative degree r, we consider the stochastic h-step reachable set for h∈r,…,H+r−1h∈\r,…,H+r-1\ at confidence level α, defined in [3] as ℛik+h(α) _i^k+h(α) :=x∈ℝn∣∃uik,…,uik+h−1⊆, =\\,x ^n ∃\u_i^k,…,u_i^k+h-1\ , (11) (x−μik+h)⊤Σi,w,h−1(x−μik+h)≤χn,α2, (x- _i^k+h) _i,w,h^-1(x- _i^k+h)≤χ^2_n,α\,\, where μik+h:=[xik+h]=Aihμik+∑τ=0h−1Aih−1−τBiuik+τ _i^k+h:=E[x_i^k+h]=A_i^h _i^k+ _τ=0^h-1A_i^h-1-τB_iu_i^k+τ is the mean h-step state prediction under nominal dynamics, Σi,w,h _i,w,h is the aggregated process noise covariance over h steps, and χn,α2χ^2_n,α denotes the chi-squared quantile at confidence level α in dimension n. The control set ⊂ℝmhU ^mh is assumed convex and compact. The corresponding mean reachable state set at time k+hk+h is defined as ℳik+h:=μik+h=Aihμik+∑τ=0h−1Aih−1−τBiuik+τ∣uik+τ∈.M_i^k+h:= \ _i^k+h=A_i^h _i^k+ _τ=0^h-1A_i^h-1-τB_iu_i^k+τ u_i^k+τ \. (12) Under Assumption 1, this set spans a rich subset of the state space consistent with input constraints, enabling meaningful target selection and barycenter computation. Thus, ℛik+h(α)R_i^k+h(α) can be viewed as a union of confidence ellipsoids centered at points in ℳik+hM_i^k+h, capturing control variability and stochastic disturbances over h steps. 4.2 Local Sample Selection and MPC Implementation Due to the output relative degree r, the reachable output set remains a singleton for the first r−1r-1 steps, determined solely by the open-loop propagation [yik+h]=CiAihμikE[y_i^k+h]=C_iA_i^h _i^k for all h<rh<r. As a result, the reachability and sample selection become meaningful only from time step k+rk+r onward. To resolve the indeterminacy between the future agent trajectory and the transport plan, let AihμikA_i^h _i^k denote the fixed open-loop nominal prediction of the state at time k+hk+h, h=r,…,H+r−1h=r,…,H+r-1, obtained by propagating the current state mean μik _i^k forward h steps with zero input. By treating these nominal predictions as deterministic reference points, we propose selecting local target points ik+h:=qjj∈ik+hQ_i^k+h:=\q_j\_j _i^k+h centered around this nominal prediction, so that their barycenter q¯ik+h q_i^k+h lies close to or within the h-step mean reachable output set Ciℳik+hC_iM_i^k+h. To this end, we adopt the following local optimal transport-based local sample selection strategy: minqj,πjk+h∑j∈ik+hπjk+h‖Aihμik−qj‖2 _q_j,\, _j^k+h _j _i^k+h _j^k+h\|A_i^h _i^k-q_j\|^2 subject to ∑j∈ik+hπjk+h=αik+h:=1Mi _j _i^k+h _j^k+h= _i^k+h:= 1M_i, where αik+h _i^k+h is the constant mass allocated to agent i at each step, and πjk+h _j^k+h is the transport plan. Since AihμikA_i^h _i^k is fixed, an analytic solution exists that allocates mass to nearest neighbors qjq_j sequentially, respecting sample capacities βi,jk _i,j^k, until the transported mass reaches 1Mi 1M_i as in [7]. According to Proposition 1, the Wasserstein cost is minimized when the expected outputs coincide with the barycenter q¯ik+h q_i^k+h for h∈r,…,H+r−1h∈\r,…,H+r-1\. This requires the barycenter to lie within the mean reachable output set Ciℳik+hC_iM_i^k+h, which may not always hold. In the following, we thus discuss how to handle cases where q¯ik+h q_i^k+h lies outside Ciℳik+hC_iM_i^k+h to ensure feasibility. 4.3 Reachable Barycenter Approximation under MPC 4.3.1 Projection onto Feasible Output Set If the selected barycenter q¯ik+h q_i^k+h lies outside the mean reachable output set Ciℳik+hC_iM_i^k+h, we project it onto the closest feasible output: q~ik+h:=argminz∈Ciℳik+h‖z−q¯ik+h‖2. q_i^k+h:= *argmin_z∈ C_iM_i^k+h\|z- q_i^k+h\|^2. (13) This projection ensures that there exists a feasible control sequence that reaches q~ik+h q_i^k+h, acknowledging MPC’s receding horizon limitation. 4.3.2 Sample Selection with Soft Constraint Alternatively, we incorporate feasibility into sample selection via a soft constraint. Denoting q¯ik+h q_i^k+h as the empirical mean of selected target points qj\q_j\, we solve: minuik+τ,qj _\u_i^k+τ\,\q_j\\, ∑j∈ik+hπjk+h‖qj−q¯ik+h‖2+λ‖q¯ik+h−q^ik+h‖2 _j _i^k+h _j^k+h\|q_j- q_i^k+h\|^2+λ \| q_i^k+h- q_i^k+h \|^2 s.t.qj .t. q_j ∈Ciℳik+h,q¯ik+h:=∑j∈ik+hπjk+hqj∑j∈ik+hπjk+h, ∈ C_iM_i^k+h,\,\, q_i^k+h:= _j _i^k+h _j^k+hq_j _j _i^k+h _j^k+h, q^ik+h q_i^k+h :=Ci(Aihxik+∑τ=0h−1Aih−1−τBiuik+τ), :=C_i (A_i^hx_i^k+ _τ=0^h-1A_i^h-1-τB_iu_i^k+τ ), where λ>0λ>0 controls the trade-off between proximity to the desired barycenter and feasibility. Integrating the reachability-aware target selection and adjustment procedures discussed in this section, the complete three-stage D2D^2OC framework is summarized in Algorithm 1. Algorithm 1 Decentralized D2D^2OC with Reachability-Aware MPC 1:Initialize: Each agent i sets μi0 _i^0 and initial coverage weights βi,j0 _i,j^0. 2:for each time step k=0,1,2,…k=0,1,2,… do 3: Step 1: Reachability-Aware Target Selection and MPC (Sec. 3 & 4) 4: for h=rh=r to H+r−1H+r-1 do ⊳ Receding horizon H-step analysis 5: Compute nominal prediction AihμikA_i^h _i^k and identify local samples ik+hQ_i^k+h. 6: Obtain barycenter q¯ik+h q_i^k+h via OT; apply projection (13) or soft constraint. 7: end for 8: Compute the optimal control uiku_i^k by solving (5) 9: Apply uiku_i^k and update state mean μik+1 _i^k+1. 10: Step 2: Weight Update 11: Solve local Wasserstein distance to update sample weights βi,jk+1 _i,j^k+1. 12: Reduce weights of covered samples to prioritize under-explored areas. 13: Step 3: Weight Sharing and Consensus 14: Exchange βi,jk+1 _i,j^k+1 with neighbors within communication range. 15: Synchronize via min-weight consensus: βi,j←min(βi,j,βneighbor,j) _i,j← ( _i,j, _neighbor,j). 16:end for 5 Convergence Analysis Now the convergence behavior of the proposed control strategy is analyzed under this section. The main result shows that the global output distribution converges toward the target in the mean squared Wasserstein sense with a certain bound. Theorem 2. (Convergence of D2OC with MPC and Decentralized Barycenter Updates) Consider a multi-agent system with discrete-time stochastic linear dynamics (1) and decentralized, local communication. Let the output relative degree be r, and define the reachable output barycenter set at time k+hk+h, for h=r,…,H+r−1h=r,…,H+r-1, as Ciℳik+hC_iM_i^k+h, where ℳik+hM_i^k+h is given by (12), μik=[xik] _i^k=E[x_i^k], and U is a compact input set. At each k, an MPC solves the finite-horizon problem (5) over horizon H to minimize deviation between predicted outputs and target barycenters q¯ik+hh=rH+r−1\ q_i^k+h\_h=r^H+r-1, applying only the first input uiku_i^k before re-optimization. Each agent determines these barycenters from local sample points with weights βi,jk\ _i,j^k\, updated by (i) local coverage and (i) decentralized communication through pointwise minimum of weights between agents within communication range. Define the empirical multi-agent output distribution ρkρ^k and let ν be the target distribution as shown in (3). Then, under decentralized D2OC with local weight sharing, [22(ρk+1,ν)]≤(1−ck+1)[22(ρk,ν)]+C(k+1)2, [W_2^2(ρ^k+1,ν)]≤ (1- ck+1 )E[W_2^2(ρ^k,ν)]+ C(k+1)^2, (14) for constants c,C>0c,\,C>0, where c depends on controllability and convergence rate, and C bounds projection and noise effects. Consequently, limk→∞[22(ρk,ν)]≤ϵh, _k→∞E[W_2^2(ρ^k,ν)]≤ _h, for some ϵh>0 _h>0 determined by communication frequency, noise, and projection errors. Proof. At each step k, the MPC solves the finite-horizon problem (5), which is strictly convex under box input constraints (Theorem 1), applying only the first input before re-optimizing. Deviations arise from: (i) stochastic noise, and (i) dynamic local weight updates. Let the mean reachable output set at k+hk+h be Ciℳik+hC_iM_i^k+h, and q¯ik+h q_i^k+h be the target barycenter selected from target samples. Under decentralized communication, agents update sample weights βi,jk _i,j^k by taking the minimum of local and received weights within communication range. This conservative update aligns barycenters across neighbors, reducing redundant coverage and enhancing distributional alignment with ν. Then, the error is defined for h=r,…,H+r−1h=r,…,H+r-1 by eik+h:=yik+h−q¯ik+h=Cixik+h−q¯ik+h.e_i^k+h:=y_i^k+h- q_i^k+h=C_ix_i^k+h- q_i^k+h. By triangle inequality, we have ‖eik+h‖≤‖Cixik+h−q~ik+h‖+‖q~ik+h−q¯ik+h‖,\|e_i^k+h\|≤\|C_ix_i^k+h- q_i^k+h\|+\| q_i^k+h- q_i^k+h\|, where q~ik+h q_i^k+h is the projection of q¯ik+h q_i^k+h onto the mean reachable output set Ciℳik+hC_iM_i^k+h. (i) Projection error: If q¯ik+h∈Ciℳik+h q_i^k+h∈ C_iM_i^k+h, then ‖q~ik+h−q¯ik+h‖=0\| q_i^k+h- q_i^k+h\|=0; otherwise, it is bounded by some ϵε. (i) Tracking error: Writing xik+h=x^ik+h+ξik+hx_i^k+h= x_i^k+h+ _i^k+h with nominal x^ik+h x_i^k+h and noise ξik+h _i^k+h, we have ‖Cixik+h−q~ik+h‖≤‖Cix^ik+h−q~ik+h‖+‖Ciξik+h‖.\|C_ix_i^k+h- q_i^k+h\|≤\|C_i x_i^k+h- q_i^k+h\|+\|C_i _i^k+h\|. By construction, q~ik+h∈Ciℳik+h q_i^k+h∈ C_iM_i^k+h, so the nominal output Cix^ik+hC_i x_i^k+h can be steered arbitrarily close to q~ik+h q_i^k+h under the given dynamics. For the stochastic part, [‖Ciξik+h‖2]≤λmax(CiΣξCi⊤)E[\|C_i _i^k+h\|^2]≤ _ (C_i _ξC_i ), where Σξ=∑ℓ=0h−1Aih−1−ℓΣw(Aih−1−ℓ)⊤ _ξ= _ =0^h-1A_i^h-1- _w(A_i^h-1- ) , and λmax(⋅) _ (·) denotes the largest eigenvalue of a symmetric positive semi-definite matrix. With finite H, Σξ _ξ is bounded under Assumption 2, ensuring uniform boundedness of [‖eik+h‖]E[\|e_i^k+h\|] with a bound depending on H. (i) Wasserstein descent: The bounded tracking error established in (i) and (i) ensures a dissipative drift. Under limited communication, the decentralized min-weight consensus mitigates overestimation of coverage weights. By reducing variance in ρkρ^k, this mechanism promotes consistent alignment with the target distribution ν. To derive the recurrence relation in (14), we utilize the update law ρk+1=k+1ρk+1k+1δyik+1ρ^k+1= kk+1ρ^k+ 1k+1 _y_i^k+1. By leveraging the L-smoothness of the squared Wasserstein distance in the sense of Fréchet derivatives [4, 5], [22(ρk+1,ν)]E[W_2^2(ρ^k+1,ν)] is expanded as: 22(ρk+1,ν) _2^2(ρ^k+1,ν) ≤22(ρk,ν)+2k+1⟨∇22(ρk,ν),δyik+1−ρk⟩ _2^2(ρ^k,ν)+ 2k+1 _2^2(ρ^k,ν), _y_i^k+1-ρ^k +L(k+1)2‖δyik+1−ρk‖2 + L(k+1)^2\| _y_i^k+1-ρ^k\|^2 (15) By taking the total expectation on both sides of (15) and applying the dissipative drift property [⟨∇22,δyik+1−ρk⟩∣ℱk]≤−c~22(ρk,ν)+ϵprojE[ _2^2, _y_i^k+1-ρ^k _k]≤- cW_2^2(ρ^k,ν)+ _proj, the recurrence relation in (14) is directly obtained. Here, ℱkF_k is the filtration up to time k, c~>0 c>0 is the contraction gain, and ϵproj≥0 _proj≥ 0 denotes the residual error bound from control projection and discretization. By defining c:=2c~>1c:=2 c>1 and a constant C>0C>0 that bounds the second-order expansion and projection terms, applying Chung’s Lemma [2] directly yields an O(1/k)O(1/k) convergence rate, specifically [22(ρk,ν)]≤Cc−1k−1+o(k−1)E[W_2^2(ρ^k,ν)]≤ Cc-1k^-1+o(k^-1). This ensures that limk→∞[22(ρk,ν)]≤ϵh, _k→∞E[W_2^2(ρ^k,ν)]≤ _h, where the residual error ϵh≥0 _h≥ 0 is determined by the communication frequency, the intensity of process/measurement noises, and projection errors associated with physical constraints. ∎ This convergence in expected squared Wasserstein distance implies that the mean squared error between the empirical output distribution and the reference distribution vanishes asymptotically up to ϵh _h. 6 Simulation Results To validate the proposed method, simulation results are presented in this section. Three linearized quadrotors are considered for the multi-agent platform. The state of agent i at time step k is defined by xik:=[xik,x˙ik,yik,y˙ik,zik,z˙ik,ϕik,ϕ˙ik,θik,θ˙ik,ψik,ψ˙ik]⊤,x_i^k:=[x_i^k, x_i^k,y_i^k, y_i^k,z_i^k, z_i^k, _i^k, φ_i^k, _i^k, θ_i^k, _i^k, ψ_i^k] , where xikx_i^k, yiky_i^k, and zikz_i^k denote position coordinates, and ϕik _i^k, θik _i^k, and ψik _i^k represent roll, pitch, and yaw angles, respectively. Dotted symbols indicate the corresponding translational or angular velocities in discrete time. The four control inputs for agent i are ui,1k=τi,ϕku_i,1^k= _i,φ^k, ui,2k=τi,θku_i,2^k= _i,θ^k, ui,3k=τi,ψku_i,3^k= _i,ψ^k, and ui,4k=Tiku_i,4^k=T_i^k, corresponding to torques for roll, pitch, yaw, and total thrust. Each input is box-constrained to reflect motor torque and thrust limits. and the communication range for decentralized weight-sharing is set to 5. The discretized quadrotor model exhibits an output relative degree of r=4r=4, satisfying CiAiℓ−1Bi=0C_iA_i -1B_i=0 for ℓ=1,2,3 =1,2,3. This implies that control inputs at time k begin to affect the system output only at time k+4k+4. To evaluate the robustness of the proposed D2OC under severe uncertainty, substantial stochasticity is incorporated with process noise wik∼(0,0.2×12)w_i^k (0,0.2×I_12) and measurement noise vik∼(0,0.5×3)v_i^k (0,0.5×I_3). Furthermore, the initial states are randomized as xi0=x¯i0+δix_i^0= x_i^0+ _i with δi∼(0,4×12) _i (0,4×I_12) to represent initial spatial uncertainty in every simulation run. With the prediction horizon H=1H=1, the optimal control input is obtained using quadprog in MATLAB by solving the quadratic program defined in (8). When the target barycenter is outside the reachable set, it is projected to the closest feasible output as in (13). (a) Conventional D2C (b) Proposed D2OC Figure 1: Comparison of 3D quadrotor trajectories under (a) Conventional D2C and (b) Proposed Stochastic D2OC. Solid lines represent agent paths, while ‘×’ and ‘∘ ’ markers denote initial and final positions, respectively. The comparative performance between the conventional D2C approach [7, 9] and the proposed method is illustrated in Fig. 1. In the conventional approach (Fig. 1(a)), a goal point is greedily selected from sample points closest to the agent and tracked via a standard PID-type controller. This method fails to accurately reconstruct the torus-shaped distribution (gray dots), with trajectories drifting toward the outer boundaries due to the high sensitivity of the greedy logic to injected noise. In contrast, the proposed Stochastic D2OC in Fig. 1(b) demonstrates superior performance, precisely matching the reference distribution despite severe noise and initial state uncertainties. (a) Communication Events (b) Expected W22W_2^2 distance Figure 2: Performance analysis of the proposed Stochastic D2OC. (a) Communication events between agents. (b) Empirical mean (solid line) and standard deviation (shaded region) of W22W_2^2 over 100 simulation runs. To further evaluate the performance in decentralized and stochastic environments, Fig. 2 provides a detailed analysis. Fig. 2(a) shows intermittent communication events, confirming that D2OC achieves time-averaged distribution matching even with occasional weight-sharing. The statistical results in Fig. 2(b) show that the empirical mean (solid blue line) asymptotically converges toward a bounded region ϵh _h, consistent with Theorem 2. Notably, the narrow standard deviation (shaded region) demonstrates the high robustness of the proposed method, ensuring consistent performance despite the presence of large process/measurement noises and injected initial state perturbations across 100 simulation runs. 7 Conclusion This paper introduced Density-Driven Optimal Control, a novel framework for solving the multi-agent non-uniform area coverage problem under realistic operational constraints. Unlike traditional uniform coverage strategies that are inefficient or infeasible in large-scale or resource-limited settings, D2OC leverages optimal transport theory to derive control laws that steer agents toward matching a task-specific reference distribution. The framework accounts for stochastic LTI system dynamics, agent constraints, and limited communication, providing a unified and principled solution. Theoretical guarantees on convergence were established via reachability analysis, and simulation results demonstrated the effectiveness and practicality of the proposed approach. References [1] Y. Chen, T. T. Georgiou, and M. Pavon (2016) Optimal transport over a linear dynamical system. IEEE Transactions on Automatic Control 62 (5), p. 2137–2152. Cited by: §1. [2] K. L. Chung (1954) On a stochastic approximation method. The Annals of Mathematical Statistics, p. 463–483. Cited by: §5. [3] M. Fiacchini and T. Alamo (2021) Probabilistic reachable and invariant sets for linear systems with correlated disturbance. Automatica 132, p. 109808. Cited by: §4.1. [4] M. Gelbrich (1990) On a formula for the l2 wasserstein metric between measures on euclidean and hilbert spaces. Mathematische Nachrichten 147 (1), p. 185–203. Cited by: §5. [5] A. Gupta and W. B. Haskell (2021) Convergence of recursive stochastic algorithms using wasserstein divergence. SIAM Journal on Mathematics of Data Science 3 (4), p. 1141–1167. Cited by: §5. [6] R. H. Kabir and K. Lee (2021) Efficient, decentralized, and collaborative multi-robot exploration using optimal transport theory. In 2021 American Control Conference (ACC), p. 4203–4208. Cited by: §1, §2.2, §2.2. [7] R. H. Kabir and K. Lee (2021) Wildlife monitoring using a multi-uav system with optimal transport theory. Applied Sciences 11 (9), p. 4070. Cited by: §2.2, §2.2, §4.2, §6, Remark 1. [8] K. Lee and E. Brook (2025) Connectivity-preserving multi-agent area coverage via density-driven optimal control (d2oc). IEEE Control Systems Letters 9, p. 2723–2728. Cited by: §3, §3. [9] K. Lee and R. Hasan Kabir (2022) Density-aware decentralised multi-agent exploration with energy constraint based on optimal transport theory. International Journal of Systems Science 53 (4), p. 851–869. Cited by: §1, §2.2, §2.2, §6. [10] G. Mathew and I. Mezic (2009) Spectral multiscale coverage: a uniform coverage algorithm for mobile sensor networks. In Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, p. 7872–7877. Cited by: §1. [11] G. Mathew and I. Mezić (2011) Metrics for ergodicity and design of ergodic dynamics for multi-agent systems. Physica D: Nonlinear Phenomena 240 (4-5), p. 432–442. Cited by: §1. [12] G. Rapakoulias, A. R. Pedram, and P. Tsiotras (2025) Steering large agent populations using mean-field schrödinger bridges with gaussian mixture models. IEEE Control Systems Letters. Cited by: §1. [13] A. Terpin, N. Lanzetti, and F. Dörfler (2024) Dynamic programming in probability spaces via optimal transport. SIAM Journal on Control and Optimization 62 (2), p. 1183–1206. Cited by: §1. [14] C. Villani (2008) Optimal transport: old and new. Vol. 338, Springer Science & Business Media. Cited by: §2.1.