Paper deep dive
Distributed Harmonic-Covariance Feedback for Hole Closure in Planar Robot Swarms
Jitian Liu, Subhrajit Bhattacharya
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 10/12/2026, 5:37:58 AM
Summary
This paper proposes a distributed harmonic-covariance feedback controller for planar robot swarms to achieve hole closure in their sensing coverage. The method uses harmonic edge signals from the regularized Hodge Laplacian of the Čech complex to generate attractive forces between agents. The controller is realized distributively using a Chebyshev polynomial approximation of the harmonic projector. Theoretical analysis guarantees unique global trajectories and provides offline computable rational lower bounds for cumulative hole area, proving that holes close in a time-averaged sense. Numerical simulations demonstrate effective hole closure within 9.9 seconds.
Entities (8)
Relation Signals (6)
Jitian Liu → affiliatedwith → Lehigh University
confidence 99% · Jitian Liu and Subhrajit Bhattacharya are with the Department of Mechanical Engineering and Mechanics, Lehigh University
Subhrajit Bhattacharya → affiliatedwith → Lehigh University
confidence 99% · Jitian Liu and Subhrajit Bhattacharya are with the Department of Mechanical Engineering and Mechanics, Lehigh University
Harmonic-Covariance Feedback → appliedto → Robot Swarm
confidence 95% · Distributed Harmonic-Covariance Feedback for Hole Closure in Planar Robot Swarms
Harmonic-Covariance Feedback → uses → Hodge Laplacian
confidence 95% · For distributed evaluation, we use a fixed-degree polynomial of the regularized Hodge Laplacian.
Harmonic-Covariance Feedback → utilizes → Čech Complex
confidence 95% · We consider hole closure in a planar multi-agent system whose sensing disks form a time-varying Čech complex.
Harmonic-Covariance Feedback → approximatedby → Chebyshev Polynomial
confidence 90% · We choose a Chebyshev expansion... The symmetric matrix Hd=pd(Lw) approximates the exponential.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We consider hole closure in a planar multi-agent system whose sensing disks form a time-varying Čech complex. We propose a harmonic-covariance feedback that uses harmonic edge signals to locate holes and converts those signals into agent forces. For distributed evaluation, we use a fixed-degree polynomial of the regularized Hodge Laplacian. This polynomial is evaluated through local operator applications. The closed loop has a unique global trajectory from every safe initial state. The model data and initial state are supplied with procedures for approximation to any requested accuracy with explicit error bounds. We assume that hole area is measurable and bounded, holes evolve continuously without splitting or merging, and only finitely many temporary new holes appear. Under these conditions, offline finite calculations produce rational lower bounds for the cumulative hole area on the time horizon without querying trajectory data. These bounds converge to that cumulative area, whether finite or infinite. If it is finite, holes close in the time-averaged sense: the average normalized hole area over every fixed-length time window tends to zero as time tends to infinity. In a 15-second numerical demo, all holes close by 9.9s and remain closed through the end.
Tags
Links
- Source: https://arxiv.org/abs/2610.12637v1
- Canonical: https://arxiv.org/abs/2610.12637v1
Trouble viewing inline? Open PDF directly →
Full Text
61,174 characters extracted from source content.
Expand or collapse full text
Distributed Harmonic-Covariance Feedback for Hole Closure in Planar Robot Swarms Jitian Liu Subhrajit Bhattacharya †thanks: $ˆ1$Jitian Liu and Subhrajit Bhattacharya are with the Department of Mechanical Engineering and Mechanics, Lehigh University, Bethlehem, PA, USA. Emails: jila24, sub216@lehigh.edu. Abstract We consider hole closure in a planar multi-agent system whose sensing disks form a time-varying Čech complex. We propose a harmonic-covariance feedback that uses harmonic edge signals to locate holes and converts those signals into agent forces. For distributed evaluation, we use a fixed-degree polynomial of the regularized Hodge Laplacian. This polynomial is evaluated through local operator applications. The closed loop has a unique global trajectory from every safe initial state. The model data and initial state are supplied with procedures for approximation to any requested accuracy with explicit error bounds. We assume that hole area is measurable and bounded, holes evolve continuously without splitting or merging, and only finitely many temporary new holes appear. Under these conditions, offline finite calculations produce rational lower bounds for the cumulative hole area on the time horizon without querying trajectory data. These bounds converge to that cumulative area, whether finite or infinite. If it is finite, holes close in the time-averaged sense: the average normalized hole area over every fixed-length time window tends to zero as time tends to infinity. In a 15-second numerical demo, all holes close by 9.9s and remain closed through the end. Index Terms: multi-agent systems, distributed control, algebraic topology, simplicial complexes, coverage control. I Introduction In nature, among groups of animals such as bird flocks and fish schools, local interaction rules can generate collective motions [1, 2]. The same principle motivates groups of mobile robots known as swarms. For multi-robot systems, sensor coverage is a basic problem: each robot senses a disk around its position, and together the robots should cover their operating region. As the robots move, the covered set can develop holes: regions enclosed by the sensing disks that receive no coverage at all. These holes are exactly the holes of the Čech complex built from the current robot positions [3, 4]. The holes of a moving cover are detected by persistent homology [3, 4, 5] and verified without localization or with switched higher-order Laplacians [6, 7]. Related work on repair and coverage includes simplicial coverage repair for robot teams [8], gradient descent on weighted combinatorial Laplacians [9], homology-based repair of agricultural coverage [10], conformal mapping of domains with holes onto multi-hole disks [11], and diffusive coverage of bounded domains [12]. Related tools include cellular sheaves, which generalize weighted Laplacians for distributed coordination [13, 14]; simplicial filters for processing on complexes [15, 16]; and occupation-measure relaxations, which express dynamical constraints and integral costs through measures and lead to linear optimization problems over measures [17, 18]. The coverage-repair methods above do not provide a criterion for hole closure based solely on the model and initial state. We derive such a criterion for our distributed harmonic-covariance controller. The contributions of this paper are as follows: • A feedback controller that converts the swarm Čech complex’s harmonic edge flows, which localize the holes, into locally Lipschitz force laws on the agents (such that the controlled trajectory is unique), together with a distributed realization of the controller (Sections I and IV). • Under this paper’s controlled dynamics and assumptions, we construct an offline sequence of rational lower bounds cTT\c_T\_T. Each bound is produced by explicit finite calculations without querying trajectory data. The sequence converges to the cumulative hole area over the infinite time horizon. If its limit is finite, the hole area averaged over any fixed time window tends to zero as the window moves forward (Section IV). • A numerical demo in which two holes close within 10 seconds under the regularized harmonic feedback controller (Section V). Fig. 1: Left: the union of sensing disks of the 30-agent annulus; right: its strict Čech complex. The paper is organized as follows. Section I sets up the task and model. Section I designs the controller and its distributed realization. Section IV proves uniqueness of the state space trajectory under the regularized controller, and constructs cTT\c_T\_T from the swarm’s initial state. Section V reports the numerical simulation. Appendix A records details to approximate the normalized hole areas in the swarm and the vector field yielded by the control law, whereas Appendix B supplies proofs to main theoretical results in this paper. I Problem Formulation and Sensing Model We consider N mobile sensing agents. Agent i has position qi(t)∈ℝ2q_i(t) ^2 and velocity vi(t)∈ℝ2v_i(t) ^2; we stack these variables as q,v∈ℝ2Nq,v ^2N and write the state as z=(q,v)z=(q,v). Each agent possesses a sensing disk centered at itself, and the control objective is to remove the bounded regions left uncovered by the union of those disks. Section IV-A defines these regions, which we call holes, and their cumulative cost. I-A Sensing, the Čech complex, and the Hodge Laplacian Each agent senses the open disk of radius rs>0r_s>0 centered at its position. The Čech 2-complex records how these disks overlap. For a configuration q, (q)=(q,rs)C(q)=C(q;r_s) has vertices [N]=1,…,N[N]=\1,…,N\, an edge i,j\i,j\ whenever ‖qi−qj‖<2rs\|q_i-q_j\|<2r_s, and a triangle i,j,k\i,j,k\ whenever the intersection of the three disks has strictly positive measure, equivalently whenever the minimal enclosing circle of the three points has radius strictly smaller than rsr_s [5] (see Fig. 1). At a fixed configuration, call the edges and triangles in (q)C(q) active. Orient each edge by i<ji<j and each triangle consistently with its boundary edges. These choices define incidence (boundary) matrices B1∈ℝN×EB_1 ^N× E and B2∈ℝE×FB_2 ^E× F, where E and F count the active edges and triangles. Because B1B_1 and B2B_2 represent consecutive boundary maps, they satisfy B1B2=0B_1B_2=0. The one-dimensional Hodge Laplacian is L1=B1⊤B1+B2B2⊤∈ℝE×EL_1=B_1 B_1+B_2B_2 ^E× E, and dimkerL1=β1(q) L_1= _1(q), the first Betti number of (q)C(q), equals the number of holes of the cover [16]. As the disks move, their overlap pattern changes at edge and triangle thresholds. To make the controller vary smoothly across these changes, fix 0<εE<2rs0< _E<2r_s and define the scalar cutoff η(x)=exp(−1/x),x>0,0,x≤0.η(x)= cases (-1/x ),&x>0,\\ 0,&x≤ 0. cases (1) For e=i,je=\i,j\, the edge weight is ae=aij=η(sij)η(sij)+η(1−sij),sij=2rs−‖qi−qj‖εE.a_e=a_ij= η(s_ij)η(s_ij)+η(1-s_ij),\;\;s_ij= 2r_s-\|q_i-q_j\| _E. (2) The weight aija_ij in (2) lies in [0,1][0,1] and equals one for ‖qi−qj‖≤2rs−εE\|q_i-q_j\|≤ 2r_s- _E and zero for ‖qi−qj‖≥2rs\|q_i-q_j\|≥ 2r_s. Between these transition endpoints, it decreases smoothly as the agents separate. It is locally Lipschitz, with all distance derivatives vanishing at both thresholds. Thus crossing the edge threshold creates no jump in the weighted operators. To keep matrix dimensions fixed as simplices enter and leave the active complex, use an ambient roster containing every possible edge and triangle. Its sizes are E¯=(N2) E= N2 and F¯=(N3) F= N3, with oriented incidence matrices B¯1∈ℝN×E¯ B_1 ^N× E and B¯2∈ℝE¯×F¯ B_2 E× F. For each ambient triangle f, let Rf(q)≔minc∈ℝ2maxℓ∈f‖qℓ−c‖R_f(q) _c ^2 _ ∈ f\|q_ -c\| be the radius of the minimal enclosing circle of its three vertices. Fix 0<εF<rs0< _F<r_s and, using η from (1), define bf(q)=η(sf)η(sf)+η(1−sf),sf=rs−Rf(q)εF.b_f(q)= η(s_f)η(s_f)+η(1-s_f),\;\;s_f= r_s-R_f(q) _F. The weight bfb_f is globally Lipschitz in q. Both aea_e in (2) and bfb_f can be evaluated to any prescribed accuracy with known error bounds. The active restrictions of B¯1,B¯2 B_1, B_2 are B1,B2B_1,B_2. To preserve the corresponding identity B¯1B¯2=0 B_1 B_2=0 under smooth weights, define D1=B¯1diag(ae)∈ℝN×E¯D_1= B_1diag(a_e) ^N× E and define D2∈ℝE¯×F¯D_2 E× F entrywise by D2[e,f]=B¯2[e,f]bfae∏e′⊂f,e′≠eae′2.D_2[e,f]= B_2[e,f]\,b_f\,a_e _e ⊂ f,\;e ≠ ea_e ^2. For a triangle f and a node i, entry (i,f)(i,f) is (D1D2)[i,f] (D_1D_2)[i,f] =∑e⊂fD1[i,e]D2[e,f] = _e⊂ fD_1[i,e]\,D_2[e,f] =bf∑e⊂fB¯1[i,e]B¯2[e,f]ae2∏e′⊂f,e′≠eae′2 =b_f _e⊂ f B_1[i,e] B_2[e,f]a_e^2 _e ⊂ f,\;e ≠ ea_e ^2 =bf(∏e′⊂fae′2)(B¯1B¯2)[i,f]=0. =b_f ( _e ⊂ fa_e ^2 )( B_1 B_2)[i,f]=0. The regularized Hodge Laplacian is Lw=D1⊤D1+D2D2⊤∈ℝE¯×E¯.L_w=D_1 D_1+D_2D_2 E× E. A harmonic edge form is a vector in the kernel of the active-edge restriction of LwL_w. When the weights are binary and match the current Čech complex, this restriction of LwL_w is L1L_1, so its kernel has dimension β1(q) _1(q). Figure 2 shows a numerical orthonormal pair from the weighted active-edge kernel of a connected two-hole complex. Fig. 2: Two orthonormal harmonic forms of the weighted active-edge restriction of LwL_w for the same connected two-hole complex. Triangles are shown lightly; colored edges give signed coefficients on the common scale, while gray denotes the zero level. Orientations are fixed by i<ji<j. I-B Agents and dynamics The agents have second-order dynamics with three acceleration terms: repulsion prevents hard collisions, damping reduces relative velocity along sensed edges [19, 20, 21], and the controller supplies harmonic-covariance feedback. For a damping gain α≥0α≥ 0, define L0reg(q)=B¯1diag(ae(q)2)B¯1⊤,D(q)=αL0reg(q)⊗I2.L_0^reg(q)= B_1diag (a_e(q)^2 ) B_1 ,\;\;D(q)=α L_0^reg(q) I_2. The matrix D(q)D(q) is symmetric positive semidefinite. An absent edge has ae(q)=0a_e(q)=0 and contributes nothing to the damping. Fix a hard distance a>0a>0, a repulsion activation distance b>ab>a, a gain g>0g>0, and a smoothing parameter τU>0 _U>0, and define the pairwise repulsive potential U(q)=∑i<jg2ρij4τU2+ρij2,ρij=(b−‖qi−qj‖)+‖qi−qj‖−a.U(q)= _i<j g2\, _ij^4 _U^2+ _ij^2,\;\; _ij= (b-\|q_i-q_j\|)_+\|q_i-q_j\|-a. (3) The ratio ρij _ij combines finite-range activation, vanishing at and beyond b, with a barrier as the distance decreases to a. The factor ρij4/(τU2+ρij2) _ij^4/( _U^2+ _ij^2) turns the force on softly at b; the potential retains the hard-boundary divergence used in Section IV-C. The agents obey q˙=v,v˙=−∇U(q)−D(q)v+u(q), q=v,\;\; v=-∇ U(q)-D(q)v+u(q), where u(q)∈ℝ2Nu(q) ^2N is the control input designed in Section I. I Controller Design I-A The harmonic-covariance force Harmonic information is stored as scalar values on edges, while each agent needs a two-dimensional force at a node. We first convert edge values into node vectors and use these vectors to define force laws on each agent. The ideal law uses the exact harmonic projector of the binary Čech complex; Section I-B replaces it with a fixed-degree polynomial of LwL_w from Section I-A. Here “exact” and “ideal” refer only to the projector. Let u(q)=FH(q)u(q)=F^H(q) be the ideal harmonic-covariance force. At a fixed configuration with E active edges, let P be the orthogonal projector onto kerL1 L_1 defined in Section I-A and choose an orthonormal basis ω(k)\ω^(k)\ for that space. A harmonic edge vector ω assigns a scalar ωe _e to every edge. Orient e=(p,q)e=(p,q) by p<qp<q and set de=q−qpd_e=q_q-q_p. At endpoint i, let t^e t_e be the unit vector pointing along the edge away from i and read ωe _e with that local orientation (i.e, for e=(p,q)e=(p,q), ω(p,q)=−ω(q,p) _(p,q)=- _(q,p)). Fix λR>0 _R>0 and define the lift vector yi(ω)y_i(ω) of ω for agent i as the unique minimizer of the following least-squares problem yi(ω)=argminw∈ℝ2λR‖w‖2+∑e∋i(aew⊤t^e−ωe‖de‖)2,y_i(ω)= _w ^2\; _R\|w\|^2+ _e i (a_e\,w t_e- _e\|d_e\| )^2, (4) where the factor ae∈[0,1]a_e∈[0,1] is the edge weight in (2), and λR‖w‖2 _R\|w\|^2 makes the problem strictly convex. Hence (4) has a unique solution yi=Gi−1Tiy_i=G_i^-1T_i with Gi=λRI2+∑e∋iae2t^et^e⊤,Ti=∑e∋iaet^eωe‖de‖.G_i= _RI_2+ _e ia_e^2\, t_e t_e ,\;\;T_i= _e ia_e\, t_e\, _e\|d_e\|. (5) In global edge coordinates, let Ri∈ℝ2×ER_i ^2× E with (Ri):,e=aeGi−1de/‖de‖2,i∈e,0,i∉e.(R_i)_:,e= casesa_eG_i^-1d_e/ \|d_e \|^2,&i∈ e,\\ 0,&i∉ e cases. It then follows that yi(ω(k))=Riω(k)y_i(ω^(k))=R_iω^(k). Because a single harmonic vector may concentrate around one hole and miss others, we aggregate the lifts of the full orthonormal basis in the node “covariance” Ci=∑kyi(ω(k))yi(ω(k))⊤=RiPRi⊤⪰0,C_i= _ky_i(ω^(k))\,y_i(ω^(k)) =R_iPR_i 0, (6) where ∑kω(k)ω(k)=P _kω^(k)ω^(k) T=P, P∈ℝE×EP ^E× E; thus Ci∈ℝ2×2C_i ^2× 2. Division by ‖de‖\|d_e\| puts the edge values on a common distance scale, while the outer-product sum removes the sign and basis choices. Reading the endpoint covariances along an edge gives a nonnegative scalar score, which sets the strength of an attractive pair force on that edge. For e=i,je=\i,j\, gain κ>0κ>0, and common force limit cH>0c_H>0, define se s_e =(1/2)t^e⊤(Ci+Cj)t^e≥0, =(1/2) t_e (C_i+C_j ) t_e≥ 0, (7) Fie F_i^e =κaese(qj−qi),Fje=−Fie,Firaw=∑e∋iFie, =κ a_es_e(q_j-q_i),\;\;F_j^e=-F_i^e,\;\;F_i^raw= _e iF_i^e, λ λ =cHmaxcH,maxℓ‖Fℓraw‖,FiH=λFiraw. = c_H \c_H, _ \|F_ ^raw\| \,\;\;F_i^H=λ F_i^raw. (8) Equation (7) gives the edge score, while (8) assembles the incident pair forces and scales every node by the same factor λ. Since Ci,Cj⪰0C_i,C_j 0, the score is nonnegative and independent of edge orientation. The weight aea_e attenuates the force near loss of support (see Fig. 3). Each pair is equal and opposite and collinear with its edge, so the harmonic force has zero net force and zero torque. The common scaling by λ in (8) preserves these properties and enforces ‖FiH‖≤cH\|F_i^H\|≤ c_H at every node. Edge forces FieF_i^e Zoomed-in detail Fig. 3: Edge forces of the harmonic-covariance law. Lemma 1. For every q¯∈ℝ2 q ^2, ∑i(qi−q¯)⊤FiH≤0 _i(q_i- q) F_i^H≤ 0. What Lemma 1 implies, in plain language, is that the net harmonic force tends to pull all agents together. Yet the lemma guarantees neither monotonic decrease of the hole areas nor their eventual closure. I-B Distributed realization of harmonic-covariance forces Direct computation of P in (6) requires eigenvectors in the kernel of the full matrix L1L_1, which requires centralized computation. For a fixed binary L1L_1, P=limτ→∞e−τL1P= _τ→∞e^-τ L_1. The exponential leaves the zero eigenspace unchanged and attenuates every other subspace. This limit also applies to LwL_w and one may use e−τLwe^-τ L_w as an approximate projector onto kerLw L_w. Evaluating e−τLwe^-τ L_w directly, however, still requires centralized computation. We therefore replace it by a polynomial of degree d, which requires d local applications of LwL_w [22, 15, 23]. The factor λ in (8) is computed separately. We choose a Chebyshev expansion because it treats the full spectral interval uniformly. A Taylor expansion is centered at a single point and requires a separate remainder estimate to control its error over the full interval. The Chebyshev interval requires a state-independent spectral upper bound, i.e., some B≥λmax(Lw)B≥ _ (L_w). Since the entries of D1D_1 and D2D_2 have magnitude at most one, for N≥2N≥ 2 the row-sum bounds give ‖D1⊤D1‖∞≤2N−2 \|D_1 D_1 \|_∞≤ 2N-2 and ‖D2D2⊤‖∞≤3(N−2) \|D_2D_2 \|_∞≤ 3 (N-2 ). Consequently, λmax(Lw)≤‖Lw‖∞≤‖D1⊤D1‖∞+‖D2D2⊤‖∞≤5N−8≕B _ (L_w)≤\|L_w\|_∞≤\|D_1 D_1\|_∞+\|D_2D_2 \|_∞≤ 5N-8 B. Since LwL_w is positive semidefinite, its spectrum lies in [0,B][0,B], and the affine map 2λ/B−12λ/B-1 sends this interval to [−1,1][-1,1]. Accordingly, define L~w=2Lw/B−I L_w=2L_w/B-I. Put σ=τB/2σ=τ B/2. Let TjT_j be the Chebyshev polynomial defined by Tj(cosθ)=cos(jθ)T_j( θ)= (jθ), and let IjI_j be the modified Bessel function of the first kind. Define cj=e−σI0(σ),j=0,2e−σ(−1)jIj(σ),j≥1,c_j= casese^-σI_0 (σ ),&j=0,\\ 2e^-σ (-1 )^jI_j (σ ),&j≥ 1, cases In the Jacobi–Anger expansion of e−τλe^-τλ [24], its degree-d truncation is pd(λ)=∑j=0dcjTj(2λB−1).p_d(λ)= _j=0^dc_j\,T_j ( 2λB-1 ). The symmetric matrix Hd=pd(Lw)H_d=p_d(L_w) approximates the exponential. Its parameters τ,B,d,cjτ,B,d,c_j are fixed and can be evaluated offline, so no online spectral estimate is needed. We aim to span kerLw L_w and generate edge signals by approximately projecting random initial vectors onto kerLw L_w. We initialize Z∈ℝE¯×KZ E× K once at random with K>1K>1 columns on the fixed ambient roster of Section I-A, and set Y=e−τLwZ≈HdZ≕YdY=e^-τ L_wZ≈ H_dZ Y_d. Computationally, YdY_d may be evaluated by the following recursive procedure: X0=Z,X1 X_0=Z,\;\;X_1 =L~wZ(d≥1), = L_wZ\;\;(d≥ 1), Xr+1=2L~wXr X_r+1=2 L_wX_r −Xr−1, 1≤r<d, -X_r-1,\;\;1≤ r<d, Yd=∑r=0dcrXr Y_d= _r=0^dc_rX_r =pd(Lw)Z=HdZ. =p_d(L_w)Z=H_dZ. All these computations may be performed locally. To convert the columns of YdY_d back to node quantities, extend the active-edge map RiR_i to ℝ2×E¯R^2× E by assigning zero columns to absent edges. The lift of YdY_d’s k-th column is yi,k=Ri(Yd):,ky_i,k=R_i(Y_d)_:,k, and their outer products give the probe covariance C^i=∑kyi,kyi,k⊤=RiHdZZ⊤Hd⊤Ri⊤⪰0. C_i= _ky_i,ky_i,k =R_iH_dZ H_d R_i 0. (9) Here yi,ky_i,k is analogous to yi(ω(k))y_i(ω^(k)) in Section I-A. Since C^i⪰0 C_i 0, substituting it into (7) and (8) preserves nonnegative edge scores. The rules in (8) are otherwise unchanged, except (Ci+Cj)(C_i+C_j) is replaced by (C^i+C^j)( C_i+ C_j). The factor λ uses the maximum raw force norm over all nodes, so the distributed realization assumes a connected communication graph over which this maximum can propagate. Damping is likewise local because L0reg(q)L_0^reg(q) couples each agent only to its current one-hop neighbors. IV Well-Posedness and Bounds on Cumulative Hole Areas After the regularized controller is defined, it is natural to raise the following questions: 1. Given a safe initial state, does the controller yield one unique trajectory? 2. Given an initial state, without simulation, can we determine whether the holes will eventually close under the controlled dynamics? Section IV-C gives the first question a positive answer, and Section IV-D answers the second question with an offline construction of a rational sequence cTT\c_T\_T converging to the all-time cumulative hole area, with each term obtained by finite calculation. IV-A Holes and the hole-area cost All sensing disks are open, and their union at configuration q is A(q)=⋃i∈[N]B(qi,rs)A(q)= _i∈[N]B(q_i,r_s). The holes are the bounded connected components of ℝ2∖A(q)R^2 A(q). Let H(q)H(q) be their union and let |⋅||·| denote two-dimensional Lebesgue measure. A terminal hole count does not record how much uncovered area persisted before the terminal time. We therefore use the normalized instantaneous area and its all-time cumulative cost: h(q)=|H(q)|(8Nrs)2;C=∫0∞h(q(t))t∈[0,∞).h(q)= |H(q)|(8Nr_s)^2;\;C= _0^∞h(q(t))\,dt∈[0,∞). (10) If our proposed controller is successful, our objective will be to show that C is finite. Every hole lies in the convex hull of one sensing-graph component. A component with kjk_j vertices has center diameter at most 2rs(kj−1)2r_s(k_j-1) and fits in a box of area at most 4rs2(kj−1)24r_s^2(k_j-1)^2; the normalization by (8Nrs)2(8Nr_s)^2 therefore gives 0≤h(q)≤HN0≤ h(q)≤ H_N, where HN≔(N−1)2/(16N2)<1/16H_N (N-1)^2/(16N^2)<1/16. Section I-A records the correspondence between the holes and β1(q) _1(q). IV-B Assumptions The admissible motion class satisfies the following standing assumptions: A1 ‖qi(t)−qj(t)‖>a\|q_i(t)-q_j(t)\|>a for all i<ji<j and all t≥0t≥ 0. A2 The normalized hole area h(q(t))h(q(t)) is well defined, measurable, and bounded for all t≥0t≥ 0. A3 Hole components grow and shrink continuously and do not split or merge. A4 At most finitely many holes appear after time zero, and every such hole is spurious (present only over a bounded time interval). Without splitting or merging, each hole keeps a well-defined identity over its lifetime, and since every hole appearing after time zero is spurious, their cumulative total area is finite, and therefore will not affect the finiteness of C. The measurability and boundedness assumptions make the all-time integral (10) well defined. IV-C Well-posedness For our controlled dynamics, the natural state space is the safe open set Ω=(q,v):mini<j‖qi−qj‖>a = \(q,v): _i<j\|q_i-q_j\|>a \. Proposition 2. For every safe initial state z0∈Ωz_0∈ , the controlled dynamics of Section I-B with the regularized law of Section I-B has a unique solution z:[0,+∞)→Ωz:[0,+∞)→ . That means no hard collision occurs in finite time. The proof has two parts. We first show that the vector field is locally Lipschitz. Indeed, the edge weights and entries of LwL_w have this property on Ω , and the fixed-degree polynomial filter preserves it. In (8), the maximum raw-force norm has the same regularity, while the denominator is bounded below by cH>0c_H>0. On the safe open set, the denominators in the repulsive potential are bounded away from zero locally, so U is locally C1,1C^1,1; the damping term D(q)vD(q)v is a finite sum using the same weights. Hence the full vector field z↦(v,FH(q)−∇U(q)−D(q)v)z (v,F^H(q)-∇ U(q)-D(q)v ) is locally Lipschitz. Picard-Lindelöf therefore gives a unique solution for a sufficiently short time after t=0t=0. We next show that this solution can be continued to any time. Fix any finite horizon T>0T>0. Define the energy by E(q,v)=∥v∥2/2+U(q)E(q,v)= v ^2/2+U(q). Its damping dissipation satisfies v⊤D(q)v=α∑i,jaij(q)2∥vi−vj∥2≥0v D(q)v=α _\i,j\a_ij(q)^2 v_i-v_j ^2≥ 0. Moreover, (8) gives, via ‖FiH‖≤cH\|F_i^H\|≤ c_H, that v⊤FH=∑ivi⊤FiH≤cH∑i∥vi∥≤12∥v∥2+NcH22,v F^H= _iv_i F_i^H≤ c_H _i v_i ≤ 12 v ^2+ Nc_H^22, where the last inequality comes from (cH−‖vi‖)2≥0(c_H-\|v_i\|)^2≥ 0. Thus, with CH=NcH2/2C_H=Nc_H^2/2, dEdt=v⊤FH−v⊤D(q)v≤12∥v∥2+CH≤E+CH. dEdt=v F^H-v D(q)v≤ 12 v ^2+C_H≤ E+C_H. Consequently, for t≤Tt≤ T, E(t)≤et(E(0)+CH)−CH≤eT(E(0)+CH)−CHE(t)≤ e^t (E(0)+C_H )-C_H ≤ e^T (E(0)+C_H )-C_H. Hence E, and therefore ∥v∥ v , is bounded up to time T, and q will not escape to infinity before T. The finiteness of E also enforces that U(q)U(q) does not blow up, i.e., one always has ‖qi−qj‖>a\|q_i-q_j\|>a. Hence the state stays in a compact subset of Ω before any finite endpoint. The vector field is bounded on that compact set, so the state has a limit in Ω at the endpoint. Starting from this limiting state, the same local existence result continues the solution. This gives a unique solution on every finite interval [0,T][0,T]. Since T is arbitrary, letting T→∞T→∞ gives existence for all t≥0t≥ 0. ∎ IV-D Construction of cumulative hole area lower bounds We now fix the feedback law and examine whether its cumulative hole area can be determined from the model and initial state. Fix the safe initial state z0=(q0,v0)z_0=(q_0,v_0) and write the controlled dynamics as z˙ z =F(z),F(q,v)=(v,f(q,v)), =F(z),\;\;F(q,v)= (v,f(q,v) ), (11) f(q,v) f(q,v) =−∇U(q)−D(q)v+FH(q). =-∇ U(q)-D(q)v+F^H(q). For a positive integer horizon T, the physical quantity of interest is the actual hole-area cost from Section IV-A, JT=∫0Th(q(t))t.J_T= _0^Th (q(t) )\,dt. IV-D1 Outline of bound construction Directly using h is difficult because it need not vary continuously with the agent positions. A small displacement can connect a hole to the exterior, so that region is no longer counted by h, even though it remains uncovered. We therefore approach h from below by continuous functions. For each positive integer T, it is possible to construct qT:Ω→[0,HN]q_T: →[0,H_N] independent of v such that 0≤qT(z)≤qT+1(z)≤h(q),qT(z)↑h(q).0≤ q_T(z)≤ q_T+1(z)≤ h(q),\;\;q_T(z) h(q). (12) Here ↑ denotes non-decreasing convergence. Appendix A constructs this graph-based family qTq_T. Define AT=∫0TqT(z(t))t≤JT≤CA_T= _0^Tq_T (z(t) )\,dt≤ J_T≤ C (13) such that AT↑CA_T C. For fixed T, it is possible to compute rational bounds LT,m≤AT≤UT,m.L_T,m≤ A_T≤ U_T,m. (14) The index m will tighten a fixed-T finite calculation, whereas increasing T both extends the time horizon and improves the lower approximation qTq_T of h. Consequently, the stopping tolerance controls the upper and lower bounds of ATA_T. The details proceed in three steps. Section IV-D2 constructs integral constraints that recover the cumulative behavior of the unique trajectory. Section IV-D3 replaces the controlled field in those constraints by finite rational data. Section IV-D4 then computes LT,mL_T,m and UT,mU_T,m in (14). IV-D2 Recovering the trajectory through integral constraints The construction begins by enclosing every possible state that the true trajectory can visit before time T. The energy, displacement, and repulsion estimates of Section IV-C give rational bounds QT,VTQ_T,V_T and rT>ar_T>a, depending only on F,z0,TF,z_0,T, such that the trajectory lies in a compact subset KT=(q,v)∈Ω:∥q∥≤QT,∥v∥≤VT,mini<j∥qi−qj∥≥rT⊂Ω.K_T= \\,(q,v)∈ : aligned & q ≤ Q_T,\;\; v ≤ V_T,\\ & _i<j q_i-q_j ≥ r_T aligned\, \⊂ . (15) Write XT=[0,T]×KTX_T=[0,T]× K_T. We now add constraints to rule out candidate trajectories that are inconsistent with the controlled dynamics. We use scalar test functions to impose the controlled dynamics as integral constraints and narrow the feasible set. By Lemma A and the uniqueness established in Proposition 2, the full family of constraints recovers the cumulative integrals of the true trajectory. Choose a scalar rational function ψ(t,q,v)ψ(t,q,v) with rational coefficients, C1C^1 on a neighborhood of XTX_T, and satisfying ψ(T,⋅)=0ψ(T,·)=0, independently of the unknown trajectory. By the chain rule, dψ(t,z(t))/dt=∂tψ+v⊤∇qψ+f⊤∇vψdψ (t,z(t) )/dt= _tψ+v _qψ+f _vψ, with the right-hand side evaluated at (t,z(t)) (t,z(t) ). The fundamental theorem of calculus then yields the integral constraints where the right hand side only involves the initial state: ∫0T(∂tψ+F⊤∇zψ)(t,z(t))t _0^T ( _tψ+F _zψ ) (t,z(t) )\,dt (16) =ψ(T,z(T))⏟=0−ψ(0,z0)=−ψ(0,z0). = ψ (T,z(T) )_=0-ψ(0,z_0)=-ψ(0,z_0). To obtain a sequence of finite constraint sets, we enumerate all polynomials with rational coefficients as p1,p2,…p_1,p_2,… and define ψℓ(t,z)=pℓ(t,z)−pℓ(T,z) _ (t,z)=p_ (t,z)-p_ (T,z), so ψℓ(T,z(T))=0 _ (T,z(T))=0. For each index ℓ , define aℓ=∂tψℓ+F⊤∇zψℓ a_ = _t _ +F _z _ (17) ⟹∫0Taℓ(t,z(t))dt=− _0^Ta_ (t,z(t) )\,dt=- ψℓ(0,z0)≕−βℓ. _ (0,z_0) - _ . IV-D3 Approximating the controlled vector field The integral identities in (17) constrain cumulative behavior, but they are not yet finite rational data: aℓa_ contains the generally non-polynomial field F, and βℓ _ may not be rational. We make both quantities finite rational data by approximating the field and the initial-value term. The integer k below records the requested field accuracy. Proposition 3. Suppose uniform rational approximations of the model data, with explicit error bounds, can be constructed on compact sets to any prescribed accuracy. For every integer k≥1k≥ 1, a finite construction from the model data alone produces a continuous piecewise-affine field FkF_k with rational coefficients on a neighborhood of KTK_T, satisfying ∥Fk−F∥∞,KT≤2−k F_k-F _∞,K_T≤ 2^-k. Appendix A-A proves Proposition 3 by explicit construction. For this FkF_k, choose a rational Γℓ≥∥∇zψℓ∥∞,XT _ ≥ _z _ _∞,X_T and a rational βℓk _ k satisfying |βℓk−βℓ|≤2−k | _ k- _ |≤ 2^-k. Set aℓk=∂tψℓ+Fk⊤∇zψℓ,‖aℓk−aℓ‖∞,XT≤Γℓ2−k.a_ k= _t _ +F_k _z _ ,\;\; a_ k-a_ _∞,X_T≤ _ 2^-k. (18) At stage m, compute aℓka_ k and βℓk _ k for every pair 1≤ℓ,k≤m1≤ ,k≤ m. We use these quantities to construct constraints for the finite optimization problem below to compute the lower and upper bounds of ATA_T. IV-D4 Finite optimization of the area bounds To formulate this optimization problem, we approximate qTq_T by a continuous rational piecewise polynomial q^T,m q_T,m with ‖q^T,m−qT‖∞,KT≤2−m. q_T,m-q_T _∞,K_T≤ 2^-m. (19) Appendix A-C gives this construction. Stage m now has m2m^2 constraints and one approximate-area objective. For the true trajectory, the exact identity ∫0Taℓ(t,z(t))t=−βℓ _0^Ta_ (t,z(t))\,dt=- _ gives, for each 1≤ℓ,k≤m1≤ ,k≤ m, ∫0T _0^T aℓk(t,z)dt+βℓk=∫0T(aℓk−aℓ)(t,z)dt+βℓk−βℓ a_ k(t,z)dt+ _ k= _0^T (a_ k-a_ )(t,z)dt+ _ k- _ (20) ⟹|∫0Taℓk(t,z)dt+βℓk|≤(TΓℓ+1)2−k. | _0^Ta_ k(t,z)dt+ _ k |≤ (T _ +1 )2^-k. Lemma B in Appendix B shows that there exists a collection of points xr=(tr,zr)∈XTx_r=(t_r,z_r)∈ X_T and weights wr≥0w_r≥ 0, with total weight T, that can simultaneously represent the m2m^2 constraints and the one objective value, using at most M=m2+2M=m^2+2 points. For every 1≤ℓ,k≤m1≤ ,k≤ m, this gives ∫0Taℓk(t,z(t))t _0^Ta_ k(t,z(t))\,dt =∑r=1Mwraℓk(xr), = _r=1^Mw_ra_ k(x_r), (21) ∫0Tq^T,m(z(t))t _0^T q_T,m(z(t))\,dt =∑r=1Mwrq^T,m(zr), = _r=1^Mw_r q_T,m(z_r), where r indexes the common representing points. The desired optimization problems then follow from the constraints and approximated cost q^T,m q_T,m, for each 1≤ℓ,k≤m1≤ ,k≤ m: ST,m(x,w)=∑r=1Mwrq^T,m(zr), S_T,m(x,w)= _r=1^Mw_r q_T,m(z_r), (22) j−T,m=minx,wST,m(x,w),j+T,m=maxx,wST,m(x,w), j^-_T,m= _x,wS_T,m(x,w),\;\;j^+_T,m= _x,wS_T,m(x,w), s.t.xr=(tr,zr)∈XT,wr≥0,∑r=1Mwr=T, .t.\;x_r=(t_r,z_r)∈ X_T,\;\;w_r≥ 0,\;\; _r=1^Mw_r=T, |∑r=1Mwraℓk(xr)+βℓk|≤(TΓℓ+1)2−k. | _r=1^Mw_ra_ k(x_r)+ _ k |≤ (T _ +1 )2^-k. The last line of (22) directly comes from (20) and (21). Minimizing and maximizing ST,mS_T,m subject to the constraints in (22) therefore gives the bounds jT,m−T2−m≤AT≤jT,m++T2−m.j^-_T,m-T2^-m≤ A_T≤ j^+_T,m+T2^-m. (23) The term T2−mT2^-m accounts for replacing qTq_T by q^T,m q_T,m. Note that (22) should not be solved via ordinary numerical optimization, which may return a feasible or locally optimal candidate without establishing the global bounds needed in (23). A feasible objective value bounds the minimum from above and the maximum from below; (23) requires the opposite. Instead, for fixed T,mT,m, we perform a binary search on a threshold value, θ. Real quantifier elimination [25, Ch. 14] decides feasibility with ST,m<θS_T,m<θ or ST,m>θS_T,m>θ, and these answers update rational enclosures of the minimum and maximum. Repeating at rational midpoints encloses both extrema to any demanded accuracy. Round jT,m−T2−mj^-_T,m-T2^-m downward and jT,m++T2−mj^+_T,m+T2^-m upward to rational numbers LT,mL_T,m and UT,mU_T,m, respectively, with each rounding error at most 2−m2^-m. This gives (14). Applying Lemma B with g(t,z)=qT(z)g(t,z)=q_T(z) and gm(t,z)=q^T,m(z)g_m(t,z)= q_T,m(z) shows that these bounds satisfy UT,m−LT,m→0U_T,m-L_T,m→ 0 for each fixed T. With c0=0c_0=0, define mT m_T =minm≥1:UT,m−LT,m<2−T, = \m≥ 1:U_T,m-L_T,m<2^-T \, (24) (ℓT,uT) ( _T,u_T) =(LT,mT,UT,mT),cT=maxcT−1,0,ℓT. = (L_T,m_T,U_T,m_T ),\;\;c_T= \c_T-1,0, _T \. For fixed T,mT,m, the binary search above refines the numerical bounds on the two extrema. If LT,m,UT,mL_T,m,U_T,m do not satisfy (24), increase m and solve (22) again. The gap convergence above proves that this loop terminates and yields ℓT≤AT≤uT _T≤ A_T≤ u_T with uT−ℓT<2−Tu_T- _T<2^-T. By construction, [0,T](t)qT(z(t))1_[0,T](t)q_T(z(t)) increases pointwise to h(q(t))h(q(t)). Monotone convergence theorem then gives AT↑CA_T C, also when C=+∞C=+∞. If S≤TS≤ T, then ℓS≤AS≤AT _S≤ A_S≤ A_T. Hence the running maximum in (24) is non-decreasing and satisfies AT−2−T<ℓT≤cT≤AT≤CA_T-2^-T< _T≤ c_T≤ A_T≤ C. Therefore cT↑Cc_T C. By this point, one completes the proof to the following result: Theorem 4. Assume that we can calculate rational approximations of the model data and the safe initial state z0z_0 to any prescribed accuracy, with explicit error bounds. For every integer T≥1T≥ 1, the construction (24) produces a rational number cTc_T in finitely many steps from these data alone, and cT↑C=∫0∞h(q(t))t,c_T C= _0^∞h(q(t))\,dt, (25) with C=+∞C=+∞ allowed. By Theorem 4 and (25), C<∞C<∞ exactly when cT\c_T\ is bounded. If C<∞C<∞, then for every fixed Δ>0 >0, 0≤1Δ∫t+Δh(q(s))s≤1Δ∫t∞h(q(s))s→0.0≤ 1 _t^t+ h(q(s))\,ds≤ 1 _t^∞h(q(s))\,ds→ 0. Pointwise convergence of h(q(t))h(q(t)) to zero, however, is not guaranteed unless h(q(t))h(q(t)) is uniformly continuous, by Barbalat’s lemma. V Numerical Demo This section shows how the controller behaves on a 60-agent configuration with two initial holes. The snapshots show the agent motion and the numerical Čech complex, while the area plot follows the two holes that are present at the initial time. t=0.0t=0.0 s t=5.0t=5.0 s t=10.0t=10.0 s t=15.0t=15.0 s Fig. 4: The two-hole peanut with state-dependent damping. Each column shows the configuration and numerical Čech complex from the same accepted trajectory; the terminal saved frame has β1=0 _1=0. Fig. 5: Areas (coordinate units squared) of the two tracked initial holes over the 15 s run, computed using open disks. Both curves read zero by the end of the run. The initial positions form the two-hole peanut shown on the left of Fig. 4. Table I lists the parameters and integrator settings, where the Chebyshev polynomial degree d=343d=343 is determined by the stated polynomial-tail tolerance. The run integrates the controller in Section I-B. The implementation evaluates the same finite-degree local operators with centralized single-process sparse operations. It emulates the algebra of the distributed recurrence but does not measure message delays or execute a networked message-passing protocol. TABLE I: Parameters for the numerical demo Setting Value Agents; initial number of holes N=60N=60; 2 Sensing-disk radius rs=0.234r_s=0.234 Edge/triangle regularization width 2.5×10−42.5× 10^-4 Local ridge λR _R; gain κ; force cap cHc_H 0.050.05; 88; 0.80.8 Number of columns in Z; random seed K=48K=48; 419419 τ; d 6060; 343343 Polynomial error tolerance 5×10−35× 10^-3 Repulsion distances a=0.18a=0.18; b=0.24b=0.24 Repulsion gain; smoothing 0.50.5; τU=0.25 _U=0.25 Damping gain α=1.2α=1.2 Integrator DOP853 Absolute; relative tolerance 10−910^-9; 10−610^-6 Maximum step 0.050.05 Horizon; saved-frame grid 1515s; 0.10.1s V-A Results For the 15-second demo, the saved frames had β1=2 _1=2 initially, decreasing to 1 by t=8t=8 s. The first observed zero occurred at t=9.9t=9.9 s, and every saved frame from that time through the terminal t=15t=15 s frame had β1=0 _1=0. Fig. 4 should be read from left to right: the agents move inward around the two uncovered regions, and the displayed complex has no holes by the final saved frame. Fig. 5 reports the raw areas of the two tracked initial holes. The two tracked curves first read zero respectively at 7.87.8 and 9.89.8s and remain zero afterwards. VI Conclusion and Future Work We developed a harmonic-covariance controller that turns harmonic edge information into bounded forces on the agents and admits distributed evaluation. Under the stated assumptions, the regularized closed loop has a unique global trajectory from every safe initial state. For this controlled dynamics, the initial state and model determine an offline sequence of rational numbers cTc_T. Each stage is a finite calculation and does not query trajectory data. The sequence is non-decreasing and converges to the cumulative normalized hole area. When the cost is finite, the average hole area on every fixed-length window tends to zero as the window moves forward. Yet this average statement implies neither pointwise convergence to zero nor finite-time hole closure. The numerical demo in Section V illustrates the controller architecture, with both initial holes closing in finite time. Despite the results in this paper, a key difficulty persists: construction of qTq_T and quantifier elimination in Section IV-D4 are still heavy in computation, prohibiting an efficient algorithm that tests whether the holes of a given initial state will close. Apart from addressing this problem, the next steps include weakening assumptions (A2), (A3), and (A4) and testing an actual distributed implementation rather than the current synchronous emulation. Acknowledgments This work was supported by the Air Force Office of Scientific Research (AFOSR) under award number FA9550-23-1-0046. In writing this paper, large language models (LLMs), specifically, the Codex models [26], were used for polishing human-written text, and for enhancing logical continuity and consistency. LLMs were used for verifying the correctness of the mathematical proofs and for increasing clarity of the writing. LLMs were also used in implementation of the numerical demos, including generation of parts of the code for simulation and visualization. All edits suggested by AI were manually reviewed by the authors before implementation. All ideas and technical details originated purely from the authors. References [1] C. W. Reynolds (1987) Flocks, herds and schools: a distributed behavioral model. ACM SIGGRAPH Computer Graphics 21 (4), p. 25–34. Cited by: §I. [2] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet (1995) Novel type of phase transition in a system of self-driven particles. Physical Review Letters 75 (6), p. 1226–1229. Cited by: §I. [3] V. de Silva and R. Ghrist (2007) Coverage in sensor networks via persistent homology. Algebraic & Geometric Topology 7, p. 339–358. Cited by: §I. [4] R. Ghrist (2008) Barcodes: the persistent topology of data. Bulletin of the American Mathematical Society 45 (1), p. 61–75. Cited by: §I. [5] H. Edelsbrunner and J. Harer (2010) Computational topology: an introduction. American Mathematical Society. Cited by: §I, §I-A. [6] A. Tahbaz-Salehi and A. Jadbabaie (2010) Distributed coverage verification in sensor networks without location information. IEEE Transactions on Automatic Control 55 (8), p. 1837–1849. External Links: Document Cited by: §I. [7] A. Muhammad and A. Jadbabaie (2007) Dynamic coverage verification in mobile sensor networks via switched higher order Laplacians. In Proc. Robotics: Science and Systems I (RSS), External Links: Document Cited by: §I. [8] J. Derenick, V. Kumar, and A. Jadbabaie (2010) Towards simplicial coverage repair for mobile robot teams. In Proc. IEEE International Conference on Robotics and Automation (ICRA), p. 5472–5477. External Links: Document Cited by: §I. [9] S. Yadokoro and S. Bhattacharya (2023) Weighted combinatorial Laplacian and its application to coverage repair in sensor networks. Note: arXiv:2312.04825 Cited by: §I. [10] M. Rysz, P. M. Pardalos, and S. S. Mehta (2025) A simplicial homology approach for assessing and rectifying coverage of sensor networks for improved crop management. European Journal of Operational Research 325 (1), p. 204–218. External Links: Document Cited by: §I. [11] X. Feng and C. Zhai (2026) Distributed coverage control on poriferous surface via poly-annulus conformal mapping. Note: arXiv:2601.13688 Cited by: §I. [12] K. Elamvazhuthi, C. Adams, and S. Berman (2016) Coverage and field estimation on bounded domains by diffusive swarms. In Proc. IEEE Conference on Decision and Control (CDC), p. 2867–2874. External Links: Document Cited by: §I. [13] J. Hansen and R. Ghrist (2021) Opinion dynamics on discourse sheaves. SIAM Journal on Applied Mathematics 81 (5), p. 2033–2060. External Links: Document Cited by: §I. [14] T. Hanks, H. Riess, S. Cohen, T. Gross, M. Hale, and J. Fairbanks (2025) Distributed multi-agent coordination over cellular sheaves. In Proc. IEEE Conference on Decision and Control (CDC), p. 2461–2468. External Links: Document Cited by: §I. [15] M. Yang, E. Isufi, M. T. Schaub, and G. Leus (2022) Simplicial convolutional filters. IEEE Transactions on Signal Processing 70, p. 4633–4648. External Links: Document Cited by: §I, §I-B. [16] M. T. Schaub, A. R. Benson, P. Horn, G. Lippner, and A. Jadbabaie (2020) Random walks on simplicial complexes and the normalized Hodge 1-Laplacian. SIAM Review 62 (2), p. 353–391. Cited by: §I, §I-A. [17] J. B. Lasserre, D. Henrion, C. Prieur, and E. Trélat (2008) Nonlinear optimal control via occupation measures and LMI-relaxations. SIAM Journal on Control and Optimization 47 (4), p. 1643–1666. Cited by: §I. [18] D. Henrion, M. Korda, M. Kružík, and R. Rios-Zertuche (2024) Occupation measure relaxations in variational problems: the role of convexity. SIAM Journal on Optimization 34 (2), p. 1708–1731. External Links: Document Cited by: §I. [19] A. Jadbabaie, J. Lin, and A. S. Morse (2003) Coordination of groups of mobile autonomous agents using nearest neighbor rules. IEEE Transactions on Automatic Control 48 (6), p. 988–1001. Cited by: §I-B. [20] R. Olfati-Saber (2006) Flocking for multi-agent dynamic systems: algorithms and theory. IEEE Transactions on Automatic Control 51 (3), p. 401–420. Cited by: §I-B. [21] F. Cucker and S. Smale (2007) Emergent behavior in flocks. IEEE Transactions on Automatic Control 52 (5), p. 852–862. Cited by: §I-B. [22] D. I. Shuman, P. Vandergheynst, and P. Frossard (2011) Chebyshev polynomial approximation for distributed signal processing. In 2011 International Conference on Distributed Computing in Sensor Systems and Workshops, p. 1–8. External Links: Document, 1105.1891 Cited by: §I-B. [23] L. Zhang, B. M. Sadler, R. S. Blum, and S. Bhattacharya (2021) Inter-cluster transmission control using graph modal barriers. IEEE Transactions on Signal and Information Processing over Networks 7, p. 275–293. External Links: Document Cited by: §I-B. [24] NIST Digital Library of Mathematical Functions Jacobi–anger expansions, equation 10.35.2. Note: https://dlmf.nist.gov/10.35.E2 Cited by: §I-B. [25] S. Basu, R. Pollack, and M. Roy (2006) Algorithms in real algebraic geometry. 2nd edition, Algorithms and Computation in Mathematics, Vol. 10, Springer, Berlin, Heidelberg. External Links: Document Cited by: §IV-D4. [26] OpenAI Codex. Note: AI coding assistant External Links: Link Cited by: Acknowledgments. Appendix A Constructions A-A Construction of FkF_k Fix T and k, and let dz=4Nd_z=4N. Choose a rational 0<γT<rT−a0< _T<r_T-a. On a rational grid, retain every cube that meets KTK_T and refine until its diameter is below γT/(22) _T/(2 2). These cubes lie in the safe set, and a consistent subdivision gives a compact rational polyhedral neighborhood KT⊂PT⊂ΩK_T⊂ P_T⊂ with a finite conforming rational triangulation kT_k. At each distinct vertex ζ of kT_k, evaluate the known formula F(ζ)F(ζ) directly and store one shared rational value cζ,k∈ℚdzc_ζ,k ^d_z satisfying ∥cζ,k−F(ζ)∥2≤ηk c_ζ,k-F(ζ) _2≤ _k. For a simplex σ=convζ0,…,ζdzσ=conv\ _0,…, _d_z\, let λσ,j _σ,j be its rational barycentric coordinates and define Fk(z)=∑j=0dzλσ,j(z)cζj,k,z∈σ.F_k(z)= _j=0^d_z _σ,j(z)c_ _j,k,\;\;z∈σ. (26) The same model data and positive safety margin give a rational Lipschitz upper bound LT≥1L_T≥ 1 for F on the retained simplices. With hkh_k denoting the maximum simplex diameter, refine until hk≤2−k−1/LT,ηk=2−k−1.h_k≤ 2^-k-1/L_T,\;\; _k=2^-k-1. (27) Finite refinement achieves these bounds. Shared nodal values make FkF_k continuous across simplex faces. For each simplex σ and z∈σ∩KTz∈σ∩ K_T, barycentric nonnegativity and (27) give ∥Fk(z)−F(z)∥2≤ηk+LThk≤2−k F_k(z)-F(z) _2≤ _k+L_Th_k≤ 2^-k. Taking the supremum over KTK_T proves Proposition 3. A-B Construction of qTq_T We construct qTq_T satisfying (12). Start with a fixed graph ℰE on [N][N] and draw each edge as the straight segment between its endpoints. Let ℰ(q)W_E(q) be the set of points about which some polygonal cycle in ℰE has nonzero winding number. The graph segments themselves have zero area. The normalized area enclosed by the graph but left uncovered by the sensing disks is hℰ(q)=|ℰ(q)∖A(q)|/(8Nrs)2.h_E(q)=|W_E(q) A(q)|/(8Nr_s)^2. (28) Take two configurations q,q~q, q of the same labeled nodes, ordered so that ∥q~∥∞≤∥q∥∞ q _∞≤ q _∞. Set Q=∥q∥∞Q= q _∞ and δ=maxi∥qi−q~i∥δ= _i q_i- q_i , so both configurations lie in [−Q,Q]2N[-Q,Q]^2N. Take the translation direction of one disk as horizontal. At each height, its left and right boundary points each move by at most δ. These give an area bound (2rs)(2δ)=4rsδ(2r_s)(2δ)=4r_sδ. Moving all N disks therefore changes coverage over at most 4rsNδ4r_sNδ of area. Next set q(s)=(1−s)q+sq~q(s)=(1-s)q+s q for 0≤s≤10≤ s≤ 1. The graph ℰE is fixed throughout the position box, including configurations where its edges exceed the sensing range 2rs2r_s. Each edge is bounded by the diagonal of [−Q,Q]2[-Q,Q]^2, namely 22Q≤3Q2 2Q≤ 3Q, and every point on it has speed at most δ. Thus each edge sweeps area at most ∫013Qδs=3Qδ _0^13Qδ\,ds=3Qδ. Cycle winding numbers are unchanged outside the swept regions, so the graph-enclosed region changes over at most 3Q|ℰ|δ3Q|E|δ of area. Adding the two area bounds and dividing by (8Nrs)2(8Nr_s)^2 gives |hℰ(q)−hℰ(q~)|≤Lℰ,Qδ,Lℰ,Q=4rsN+3Q|ℰ|(8Nrs)2.|h_E(q)-h_E( q)|≤ L_E,Qδ,\;\;L_E,Q= 4r_sN+3Q|E|(8Nr_s)^2. (29) It remains to approach the sensing graph without making a discontinuous edge decision. For every complete-graph edge e=i,je=\i,j\, define be(q)=4rs2−∥qi−qj∥2b_e(q)=4r_s^2- q_i-q_j ^2 and set pe,T(q)=min(1,max(0,Tbe(q))).p_e,T(q)= (1, (0,Tb_e(q)) ). (30) The number pe,Tp_e,T acts as the inclusion weight of edge e; its exclusion weight is 1−pe,T1-p_e,T. Averaging over the finitely many graphs gives qT(q)=∑ℰ[∏e∈ℰpe,T(q)][∏e∉ℰ(1−pe,T(q))]hℰ(q).q_T(q)= _E [ _e p_e,T(q) ] [ _e (1-p_e,T(q)) ]h_E(q). (31) For the sensing graph ℰ(q)E(q), a point outside A(q)A(q) is encircled by a cycle exactly when it lies in a hole, so hℰ(q)(q)=h(q)h_E(q)(q)=h(q). Since hℰh_E is nondecreasing under edge addition, the convex combination in (31), supported on sensing subgraphs, gives 0≤qT≤h0≤ q_T≤ h. Continuity follows from (29) and (30). For fixed q, since pe,T(q)p_e,T(q) is non-decreasing with T, qTq_T is also non-decreasing. Every nonsensing edge has zero weight. Once Tbe(q)≥1Tb_e(q)≥ 1 for every sensing edge, (31) reduces to hℰ(q)(q)=h(q)h_E(q)(q)=h(q). Hence qT(q)↑h(q)q_T(q) h(q), proving (12). A-C Construction of q^T,m q_T,m Fix T and m. Choose a rational Q+>0Q_+>0 so that P=[−Q+,Q+]2NP=[-Q_+,Q_+]^2N contains all q-components of KTK_T, and take a rational Lipschitz upper bound ΛT≥1 _T≥ 1 for qTq_T on P; its construction is given below. Take a rational triangulation such that any two configurations q,q~q, q in the same simplex satisfy maxi‖qi−q~i‖2≤ρm≔2−m−1/ΛT. _i\|q_i- q_i\|_2≤ _m 2^-m-1/ _T. At each shared vertex qνq^ν, compute yν∈ℚy_ν satisfying |yν−qT(qν)|≤2−m−1|y_ν-q_T(q^ν)|≤ 2^-m-1. For each graph, divide a box containing the region in (28) into small squares. Discard squares verified to lie entirely outside the region. The total area of squares verified to lie entirely inside is a lower bound; adding the area of every square that intersects the boundary gives an upper bound. Normalize these bounds and evaluate (31) using rational input bounds. Refine the grid until the resulting gap between upper and lower bounds falls below 2−m2^-m, and take its midpoint as yνy_ν. For q∈Pq∈ P, let λν(q)≥0 _ν(q)≥ 0 be its barycentric coordinate with respect to the vertices qνq^ν of the containing simplex and define the piecewise affine approximation q^T,m(q)=∑νλν(q)yν. q_T,m(q)= _ν _ν(q)y_ν. Then it follows that |q^T,m(q)−qT(q)|≤∑νλν(q)(2−m−1+ΛTρm)=2−m.| q_T,m(q)-q_T(q)|≤ _ν _ν(q) (2^-m-1+ _T _m )=2^-m. Thus this construction satisfies (19). To choose ΛT _T, note that every graph-enclosed region lies in [−Q+,Q+]2[-Q_+,Q_+]^2, so 0≤hℰ≤H0≤ h_E≤ H for a rational H≥4Q+2/(8Nrs)2H≥ 4Q_+^2/(8Nr_s)^2. For q,q~∈Pq, q∈ P, set δ=maxi∥qi−q~i∥δ= _i q_i- q_i . The edge lengths are at most 3Q+3Q_+, and |pe,T(q) |p_e,T(q) −pe,T(q~)|≤T|∥qi−qj∥2−∥q~i−q~j∥2| -p_e,T( q)|≤ T | q_i-q_j ^2- q_i- q_j ^2 | ≤T(3Q++3Q+)(δ+δ)=12TQ+δ. ≤ T(3Q_++3Q_+)(δ+δ)=12TQ_+δ. (31) then gives |qT(q)−qT(q~)| |q_T(q)-q_T( q)| ≤maxℰ|hℰ(q)−hℰ(q~)|+H∑e|pe,T(q)−pe,T(q~)| ≤ _E|h_E(q)-h_E( q)|+H _e|p_e,T(q)-p_e,T( q)| ≤(maxℰLℰ,Q++12TQ+H(N2))δ. ≤ ( _EL_E,Q_++12TQ_+H N2 )δ. Any rational ΛT≥1 _T≥ 1 not smaller than this coefficient suffices. Appendix B Proofs Summing the pair forces gives ∑i(qi−q¯)⊤FiH _i(q_i- q) F_i^H =λ∑e((qi−q¯)⊤Fie+(qj−q¯)⊤Fje) =λ _e ((q_i- q) F_i^e+(q_j- q) F_j^e ) =−λκ∑eaese∥qj−qi∥2≤0. =-λκ _ea_e\,s_e\,\|q_j-q_i\|^2≤ 0. The last inequality uses λ,κ,ae,se≥0λ,κ,a_e,s_e≥ 0. ∎ Lemma A. Fix T≥1T≥ 1 and use KTK_T and XTX_T defined in Section IV-D. Among the finite nonnegative measures on XTX_T, the equations ∫XT(∂tψ+F⊤∇zψ)μ=−ψ(0,z0) _X_T ( _tψ+F _zψ )dμ=-ψ(0,z_0) for every C1C^1 function ψ with ψ(T,⋅)=0ψ(T,·)=0 have the single solution μT=dt⊗δz(t) _T=dt _z(t) with μT(XT)=T _T(X_T)=T. By definition, for g∈C(XT)g∈ C(X_T), one has that ∫XTgdμT=∫0Tg(t,z)t _X_Tgd _T= _0^Tg(t,z)dt. Taking g=∂tψ+F⊤∇zψg= _tψ+F _zψ and applying (16) gives the stated equations. Conversely, for any solution μ, taking ψ(t,z)=∫tTζ(s)sψ(t,z)= _t^Tζ(s)\,ds gives ∫XTζ(t)μ=∫0Tζ(t)t _X_Tζ(t)\,dμ= _0^Tζ(t)\,dt, with ζ∈C([0,T])ζ∈ C([0,T]). Hence πtμ=dt _tμ=dt and μ(XT)=Tμ(X_T)=T. Meanwhile, local Lipschitzness of F gives a constant L such that ∥F(x)−F(y)∥≤L∥x−y∥ F(x)-F(y) ≤ L x-y for x,y∈KTx,y∈ K_T. Let z(t)z(t) be the trajectory under F and take Ψ(t,x)=(T−t)e−2Lt‖x−z(t)‖2 (t,x)=(T-t)e^-2Lt x-z(t) ^2. Then Ψ∈C1 ∈ C^1 and Ψ(T,⋅)=0 (T,·)=0. Since z˙=F(z) z=F(z), Cauchy–Schwarz and the Lipschitz bound give ∂tΨ+F(x)⊤∇xΨ≤−e−2Lt‖x−z(t)‖2. _t +F(x) _x ≤-e^-2Lt x-z(t) ^2. The integral identity for Ψ has zero right-hand side. Its nonpositive integrand must therefore vanish μ-almost everywhere. Hence ∥x−z(t)∥=0 x-z(t) =0 almost everywhere. This finally gives μ=dt⊗δz(t)μ=dt _z(t). ∎ Lemma B. Fix a positive integer T and a bounded continuous function g:XT→ℝg:X_T , with uniform rational approximations on XTX_T with any prescribed error bound and corresponding trajectory integral in Lemma A. For every rational ε>0 >0, without querying the trajectory, a terminating finite calculation involving the model, the initial state z0z_0, the compact set KTK_T, and these approximation data returns rational numbers L≤UL≤ U with U−L<ε,L≤∫XTgdμT≤U.U-L< ,\;\;L≤ _X_Tg\,d _T≤ U. Use the terminal-zero functions ψℓ _ defined in Section IV-D. They are dense in the required C1C^1 class. For each k, Proposition 3 supplies the field FkF_k used in aℓka_ k, while Γℓ _ and βℓk _ k are chosen as in Section IV-D. Equation (20) supplies the tolerance in the finite constraints. Let T,mP_T,m be the set of nonnegative measures μ on XTX_T satisfying μ(XT)=Tμ(X_T)=T and |∫XTaℓkμ+βℓk|≤(TΓℓ+1)2−k, 1≤ℓ,k≤m. | _X_Ta_ k\,dμ+ _ k |≤ (T _ +1 )2^-k,\;\;1≤ ,k≤ m. The true measure belongs to every T,mP_T,m by (20). These sets are nested. They are also weakly compact because XTX_T is compact, every measure is nonnegative with the same finite mass T, and the defining inequalities are weakly closed. For each fixed ℓ , letting k→∞k→∞ recovers its exact integral identity; density of the ψℓ _ family and Lemma A therefore give ⋂m≥1T,m=μT _m≥ 1P_T,m=\ _T\. Thus supμ∈T,m|∫XTgμ−∫XTgdμT|→0. _μ _T,m | _X_Tg\,dμ- _X_Tg\,d _T |→ 0. (32) We now make a fixed stage finite. Choose a continuous rational piecewise-polynomial gmg_m on XTX_T with ∥gm−g∥∞,XT≤2−m g_m-g _∞,X_T≤ 2^-m, and define Hm(x)=((aℓk(x))1≤ℓ,k≤m,gm(x))∈ℝm2+1H_m(x)=( (a_ k(x) )_1≤ ,k≤ m,g_m(x)) ^m^2+1. For every μ∈T,mμ _T,m, the normalized integral of HmH_m lies in convHm(XT)convH_m(X_T). Carathéodory’s theorem in convex geometry therefore gives points xr∈XTx_r∈ X_T, coefficients λr≥0 _r≥ 0 with ∑rλr=1 _r _r=1, and at most m2+2m^2+2 terms such that 1T∫XTHmμ=∑r=1MλrHm(xr). 1T _X_TH_m\,dμ= _r=1^M _rH_m(x_r). Set wr=Tλrw_r=T _r, padding to M=m2+2M=m^2+2 with zero weights. The extrema jm±j_m^± of ∑rwrgm(xr) _rw_rg_m(x_r) under (22) equal those of ∫gmμ g_m\,dμ over T,mP_T,m. Since μ(XT)=Tμ(X_T)=T, |∫(gm−g)μ|≤T2−m| (g_m-g)\,dμ|≤ T2^-m. Together with (32), this gives jm±→∫gdμTj_m^±→ g\,d _T and jm−T2−m≤∫XTgdμT≤jm++T2−m.j_m^--T2^-m≤ _X_Tg\,d _T≤ j_m^++T2^-m. (33) The procedure in Section IV-D4 converts (33) into rational bounds with vanishing gap, so refinement terminates. ∎