Paper deep dive
Riemannian Density-Driven Optimal Control: Tangent-Space LQR for Second-Order Multi-Agent Systems on Curved Manifolds
Kooktae Lee, Ruchika Singh
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 9/22/2026, 4:18:47 AM
Summary
This paper introduces Riemannian Density-Driven Optimal Control (R-D2OC), an extension of Euclidean D2OC for second-order multi-agent systems on curved Riemannian manifolds. The method maps local distribution objectives to tangent spaces via logarithmic maps, solves a finite-horizon Linear Quadratic Regulator (LQR) problem, and executes control via parallel transport and intrinsic dynamics. The authors provide theoretical bounds on curvature-induced approximation errors and prove a conditional discrete-descent property for the closed-loop system, validated by simulations on a 3D ellipsoidal manifold.
Entities (10)
Relation Signals (7)
Ruchika Singh → affiliatedwith → Texas Tech University
confidence 99% · Kooktae Lee and Ruchika Singh are with the Department of Mechanical and Aerospace Engineering, Texas Tech University
Kooktae Lee → affiliatedwith → Texas Tech University
confidence 99% · Kooktae Lee and Ruchika Singh are with the Department of Mechanical and Aerospace Engineering, Texas Tech University
R-D2OC → uses → LQR
confidence 97% · ...uses its weighted center as the reference for a finite-horizon LQR.
R-D2OC → uses → Logarithmic Map
confidence 96% · The proposed Riemannian D2OC (R-D2OC) constructs a local distribution objective in the tangent space of each agent through logarithmic maps...
R-D2OC → extends → D2OC
confidence 95% · This paper extends D2OC to second-order multi-agent systems evolving on Riemannian manifolds. The proposed Riemannian D2OC (R-D2OC)...
R-D2OC → executescontrolvia → Parallel Transport
confidence 94% · The resulting control is executed on the manifold through intrinsic second-order dynamics and parallel transport...
R-D2OC → validatedon → 3D Ellipsoidal Manifold
confidence 92% · Numerical simulations on a 3D ellipsoidal manifold demonstrate distribution-level control...
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Density-Driven Optimal Control (D2OC) provides an effective framework for steering multi-agent systems toward prescribed spatial distributions. However, existing D2OC formulations are primarily developed for Euclidean domains and do not directly account for intrinsic manifold geometry. This paper extends D2OC to second-order multi-agent systems evolving on Riemannian manifolds. The proposed Riemannian D2OC (R-D2OC) constructs a local distribution objective in the tangent space of each agent through logarithmic maps and uses its weighted center as the reference for a finite-horizon LQR. The resulting control is executed on the manifold through intrinsic second-order dynamics and parallel transport within a receding-horizon scheme. We establish local curvature-dependent bounds that quantify the approximation introduced by the tangent-space reduction and characterize the resulting target bias. Furthermore, we derive a conditional discrete-descent result showing that the closed-loop objective decreases when a local velocity-alignment condition is satisfied. Numerical simulations on a 3D ellipsoidal manifold demonstrate distribution-level control and empirically support the proposed approximation and descent results.
Tags
Links
- Source: https://arxiv.org/abs/2609.22678v1
- Canonical: https://arxiv.org/abs/2609.22678v1
Trouble viewing inline? Open PDF directly →
Full Text
42,715 characters extracted from source content.
Expand or collapse full text
Riemannian Density-Driven Optimal Control: Tangent-Space LQR for Second-Order Multi-Agent Systems on Curved Manifolds Kooktae Lee Affiliation: Kooktae Lee and Ruchika Singh are with the Department of Mechanical and Aerospace Engineering, Texas Tech University, Lubbock, TX 79409, USA, email: kooktae.lee@ttu.edu, ruchisin@ttu.edu Ruchika Singh †thanks: *This work was supported by NSF CAREER Grant CMMI-DCSD-2638508. Affiliation: Kooktae Lee and Ruchika Singh are with the Department of Mechanical and Aerospace Engineering, Texas Tech University, Lubbock, TX 79409, USA, email: kooktae.lee@ttu.edu, ruchisin@ttu.edu Abstract Density-Driven Optimal Control (D2OC) provides an effective framework for steering multi-agent systems toward prescribed spatial distributions. However, existing D2OC formulations are primarily developed for Euclidean domains and do not directly account for intrinsic manifold geometry. This paper extends D2OC to second-order multi-agent systems evolving on Riemannian manifolds. The proposed Riemannian D2OC (R-D2OC) constructs a local distribution objective in the tangent space of each agent through logarithmic maps and uses its weighted center as the reference for a finite-horizon LQR. The resulting control is executed on the manifold through intrinsic second-order dynamics and parallel transport within a receding-horizon scheme. We establish local curvature-dependent bounds that quantify the approximation introduced by the tangent-space reduction and characterize the resulting target bias. Furthermore, we derive a conditional discrete-descent result showing that the closed-loop objective decreases when a local velocity-alignment condition is satisfied. Numerical simulations on a 3D ellipsoidal manifold demonstrate distribution-level control and empirically support the proposed approximation and descent results. I Introduction Steering multi-agent systems toward prescribed spatial distributions while minimizing control effort is a fundamental challenge in robotics [1]. Density-Driven Optimal Control (D2OC) addresses this via a Lagrangian framework penalizing spatial distribution mismatch along control effort [1, 2, 3, 4, 5]. Existing D2OC methods, however, rely on Euclidean spaces where quadratic density objectives yield standard linear quadratic tracking representations. When agent configurations evolve on curved manifolds, these Euclidean formulations break down because linear vector operations violate intrinsic geometric constraints. To address these geometric limitations in multi-agent control, prior efforts have broadly explored three avenues. First, geometric control on Riemannian manifolds and Lie groups intrinsically resolves rotational singularities for attitude and rigid-body tracking [6, 7, 8]. Yet, these methods address internal system configuration constraints rather than agent navigation over environmentally constrained curved surfaces or collective spatial distribution control. Second, geodesic path-planning methods leverage intrinsic distances to generate valid or shortest paths on curved spaces [9, 10], while primarily focusing on geometric path generation for individual agents without accounting for system dynamics or collective distribution control. Third, swarm and coverage-control methods shape agent configurations via Voronoi or density-based objectives [11, 12, 13], but are predominantly Euclidean and do not preserve intrinsic manifold geometry. Thus, an optimal distribution-control framework accounting for intrinsic geometry, system dynamics, and computational efficiency remains undeveloped. To resolve these limitations, we develop a local Riemannian realization of D2OC by constructing the distribution target in the tangent space at the current agent state. The logarithmic map converts local reference samples of target distribution into tangent vectors, whose weighted average defines a virtual target for a finite-horizon LQR. The resulting control is then executed through the intrinsic second-order dynamics using parallel transport, while the receding-horizon update continuously shifts the tangent-space reference with the agent state. This construction retains the fixed-dimensional structure of LQR while providing explicit control over the geometric approximation introduced by the tangent-space reduction. We establish a geometric alignment condition under which a local discrete-time descent property holds, providing a formal criterion for manifold-based tracking control. The main contributions of this paper are thus threefold: 1) a local Riemannian formulation of D2OC (R-D2OC) for second-order agents on curved spaces; 2) a rigorous analysis of approximation errors introduced by the tangent-space reduction under local curvature; and 3) a conditional discrete-descent guarantee for the resulting receding-horizon controller under a local velocity-alignment condition. I Preliminaries and Problem Formulation I-A Density-Driven Optimal Control Framework The Euclidean density-driven optimal control framework steers a multi-agent system toward a prescribed spatial distribution represented by weighted reference samples. As described in [1, 2, 3, 5], the process consists of three iterative stages. For a given coverage mission, a reference distribution is represented by sample points qjj=1N\q_j\_j=1^N with uniform weights 1/N1/N. In Stage A, each agent identifies a local subset of reference samples and computes a control input driving its motion toward the corresponding weighted spatial target. Stage B decreases the sample weights based on coverage progress, while Stage C exchanges updated weights among neighboring agents within communication range. Repeated execution enables progressive coverage of the reference distribution. Due to space limits, further details are omitted and can be found in the cited literature. This paper focuses on Stage A. At each time step, an agent selects a local subset of reference samples and constructs a local target represented by their weighted barycenter. The agent then solves a finite-horizon optimal control problem to track this target while minimizing control effort. In Euclidean D2OCD^2OC, this problem reduces to standard linear-quadratic control since agent dynamics and targets are represented in Euclidean coordinates. However, this construction relies on linear vector operations that are not intrinsically defined on a general Riemannian manifold. Consequently, local target construction and its LQR realization cannot be directly applied to agents moving on curved spaces. In what follows, we develop an intrinsic counterpart of Stage A using logarithmic maps and tangent-space dynamics. I-B Density Representation on Manifold To lift the aforementioned Euclidean D2OCD^2OC framework onto generic curved domains, we now formulate the density representation and optimal transport metrics on a smooth Riemannian manifold. Consider agent i with configuration xi(t)∈ℳx_i(t) , where (ℳ,g)(M,g) is a smooth d-dimensional Riemannian manifold. Let vi(t)∈Txi(t)ℳv_i(t)∈ T_x_i(t)M denote its velocity. At planning step k, the agent considers its local subset of reference points ikS_i^k. The local desired density for agent i is represented by the finite positive measure νik=∑j∈ikπjδqj,qj∈ℳ,πj≥0, _i^k= _j _i^k _j _q_j,\;q_j , _j≥ 0, with total local mass defined as αi=∑j∈ikπj=1Mi>0, _i= _j _i^k _j= 1M_i>0, where πj _j represents the mass transport weight associated with reference point qjq_j, δ(⋅) _(·) denotes a Dirac measure, and MiM_i stands for the fixed total agent capacity assigned to agent i. The symbol MiM_i is determined by the ratio of total operation time to the sampling interval (Ttotal/ΔtT_total/ t), reflecting a uniform per-step energy consumption profile. Extending the Wasserstein distance concept to the Riemannian setting, the optimal transport cost between the agent’s current state (modeled as a Dirac measure at xix_i with mass αi _i) and the local reference measure νik _i^k is given by the squared Riemannian Wasserstein distance: 22(αiδxi,νik)=∑j∈ikπjdistg(xi,qj)2, _2^2( _i _x_i, _i^k)= _j _i^k _jdist_g(x_i,q_j)^2, (1) where distg(xi,qj)dist_g(x_i,q_j) denotes the Riemannian distance between xix_i and qjq_j induced by the metric tensor g on the configuration manifold M. Thus, the density mismatch naturally translates into a weighted sum of squared Riemannian distances over the agent’s local sample set. To characterize the geometric center toward which the agent is steered, the corresponding local Riemannian barycenter for agent i is defined as the minimizer of this transport cost: q¯R,i∈argminz∈12∑j∈ikπjdistg(z,qj)2, q_R,i∈ *arg\,min_z 12 _j _i^k _jdist_g(z,q_j)^2, (2) where U is a strongly convex neighborhood containing the relevant local sample points. I-C Intrinsic Second-Order Dynamics For notational simplicity, the agent index i is omitted hereafter. Each agent obeys x˙=v,∇x˙v=u,u(t)∈Tx(t)ℳ, x=v, _ xv=u, u(t)∈ T_x(t)M, (3) where ∇ is the Levi–Civita connection, thus u is the covariant acceleration. We consider the Bolza problem J J =∫t0tf[V(x(t))+12‖u(t)‖x(t)2]dt+Φ(x(tf)), = _t_0^t_f [V(x(t))+ 12 u(t) _x(t)^2 ]dt+ (x(t_f)), (4) V(x) V(x) =12∑j∈kπjdistg(x,qj)2, = 12 _j ^k _jdist_g(x,q_j)^2, Φ(x) (x) =wf2∑j∈kπjdistg(x,qj)2,wf≥0, = w_f2 _j ^k _jdist_g(x,q_j)^2, 18.49988ptw_f≥ 0, subject to (3) and a fixed initial state. The constant wf≥0w_f≥ 0 penalizes the terminal distribution error relative to the running effort. I Tangent-Space Reduction I-A Normal Coordinates To map global geometric entities from the Riemannian manifold ℳM into a local Euclidean vector space suitable for linear control design, we utilize exponential and logarithmic maps. The exponential map Expxk:Txkℳ→ℳExp_x_k T_x_kM maps a tangent vector from the tangent space at xkx_k to a point along the geodesic on ℳM. Its inverse, the logarithmic map Logxk:ℳ→TxkℳLog_x_k → T_x_kM, projects points from the manifold back onto the tangent space TxkℳT_x_kM, where the norm of the tangent vector corresponds to the Riemannian distance along the minimizing geodesic. At replanning step k, we define the local coordinates in the tangent space via the logarithmic map as ξ=Logxk(x),ηj=Logxk(qj).ξ=Log_x_k(x), _j=Log_x_k(q_j). (5) By choosing an orthonormal basis of TxkℳT_x_kM, we explicitly identify this tangent space with ℝdR^d. For each x in the local normal neighborhood, let x→xkP_x→ x_k denote parallel transport along the unique radial geodesic s↦Expxk(sξ),s∈[0,1].s _x_k(sξ), s∈[0,1]. Define the transported velocity ζ=x→xkv∈Txkℳ.ζ=P_x→ x_kv∈ T_x_kM. Assume the sectional curvature satisfies |Sec|≤κ|Sec|≤κ on a compact normal neighborhood containing the local trajectory and the relevant target samples. To characterize the linearization error on the manifold, consider the radial geodesic γ(s)=Expxk(sξ)γ(s)=Exp_x_k(sξ) for s∈[0,1]s∈[0,1]. For any direction v∈TxkMv∈ T_x_kM, the action of the differential DExpxk(ξ)vDExp_x_k(ξ)v is represented by a Jacobi field J(s)J(s) along γ(s)γ(s) with J(0)=0J(0)=0 and ∇γ˙J(0)=v _ γJ(0)=v. Integrating the Jacobi equation ∇γ˙2J+R(J,γ˙)γ˙=0 _ γ^2J+R(J, γ) γ=0 yields the integral form J(1)=v−∫01(1−s)R(J(s),γ˙)γ˙s.J(1)=v- _0^1(1-s)R(J(s), γ) γ\,ds. (6) Noting that J(1)=DExpxk(ξ)vJ(1)=DExp_x_k(ξ)v, applying parallel transport x→xkP_x→ x_k along γ and subtracting v=Iv=Iv gives (x→xkDExpxk(ξ)−I)v= (P_x→ x_kDExp_x_k(ξ)-I )v= (7) −∫01(1−s)γ(s)→xkR(J(s),γ˙)γ˙ds. - _0^1(1-s)P_γ(s)→ x_kR(J(s), γ) γ\,ds. Taking the norm on both sides and using the sectional curvature bound |Sec|≤κ|Sec|≤κ, the geodesic speed ‖γ˙‖=‖ξ‖\| γ\|=\|ξ\|, and the Jacobi field estimate ‖J(s)‖≤C0‖v‖\|J(s)\|≤ C_0\|v\| (where C0>0C_0>0 is a local constant bounded near unity), the integral term is bounded by C1κ‖ξ‖2‖v‖C_1κ\|ξ\|^2\|v\| for a constant C1>0C_1>0 depending on C0C_0. Taking the supremum over ‖v‖=1\|v\|=1 leads directly to the operator norm estimate ‖x→xkDExpxk(ξ)−I‖≤C1κ‖ξ‖2. _x→ x_kDExp_x_k(ξ)-I ≤ C_1κ ξ ^2. (8) To understand the dynamic behavior in the tangent space TxkMT_x_kM, recall that ξ=Logxk(x)ξ=Log_x_k(x) and its velocity variable in the tangent space is defined via parallel transport as ζ=x→xkx˙ζ=P_x→ x_k x. Differentiating ξ with respect to time yields ξ˙=DLogxk(x)x˙=DLogxk(x)xk→xζ ξ=DLog_x_k(x) x=DLog_x_k(x)P_x_k→ xζ. Since the logarithmic map is the inverse of the exponential map, its differential satisfies DLogxk(x)=(DExpxk(ξ))−1DLog_x_k(x)=(DExp_x_k(ξ))^-1. Applying the operator norm estimate from (8) to this inverse differential and utilizing Neumann series expansion, the leading-order linear mapping reduces to the identity operator I, while the curvature-induced perturbation is collected in the residual term rξr_ξ. Consequently, after possibly shrinking the neighborhood, ξ˙=ζ+rξ, ξ=ζ+r_ξ, (9) with ‖rξ‖≤C2κ‖ξ‖2‖ζ‖. r_ξ ≤ C_2κ ξ ^2 ζ . Similarly, consider the acceleration dynamics under the intrinsic covariant derivative ∇x˙x˙=u _ x x=u. Differentiating the parallel-transported velocity ζ=x→xkx˙ζ=P_x→ x_k x along the state trajectory x(t)x(t) requires accounting for the curvature-induced non-commutativity of parallel transport along distinct paths (the path along the trajectory versus the radial geodesic connecting x(t)x(t) to xkx_k). Evaluating this geometric drift via standard Jacobi-field and curvature tensor estimates along the geodesic triangle formed by xkx_k, x(t)x(t), and x(t+Δt)x(t+ t) yields ζ˙=u0+rζ, ζ=u_0+r_ζ, (10) where u0=x→xku_0=P_x→ x_ku represents the parallel-transported control input, and the residual acceleration error rζr_ζ is bounded by ‖rζ‖≤C3(κ‖ξ‖‖ζ‖2+κ‖ξ‖2‖u0‖). r_ζ ≤ C_3 (κ ξ ζ ^2+κ ξ ^2 u_0 ). The constants C1,C2,C3C_1,C_2,C_3 depend only on the selected compact neighborhood and uniform bounds on the relevant geometric quantities, confirming that the Euclidean double integrator (ξ˙=ζ,ζ˙=u0)( ξ=ζ, ζ=u_0) serves as the exact leading-order local model in the fixed tangent space TxkMT_x_kM. I-B Local Expansion of the Density Objective For ξ,ηjξ, _j in a sufficiently small common normal neighborhood, the squared Riemannian distance admits the local expansion 12distg(Expxk(ξ),Expxk(ηj))2 12dist_g (Exp_x_k(ξ),Exp_x_k( _j) )^2 =12‖ξ−ηj‖xk2+Rj(ξ,ηj). = 12 ξ- _j _x_k^2+R_j(ξ, _j). (11) In normal coordinates centered at xkx_k, the squared Riemannian distance has a Euclidean quadratic leading term and a fourth-order curvature correction. Thus, for ξ and ηj _j in a sufficiently small common normal neighborhood, there exists C4>0C_4>0, depending only on the neighborhood, such that |Rj(ξ,ηj)|≤C4κ(‖ξ‖+‖ηj‖)4,|R_j(ξ, _j)|≤ C_4κ ( ξ + _j )^4, (12) with C4C_4 depending only on the chosen neighborhood. The remainder vanishes identically when the neighborhood is flat. Therefore, V(x)=12∑j∈kπj‖ξ−ηj‖2+RV(ξ), V(x)= 12 _j ^k _j ξ- _j ^2+R_V(ξ), (13) where |RV(ξ)|≤C4κ∑j∈kπj(‖ξ‖+‖ηj‖)4. |R_V(ξ)|≤ C_4κ _j ^k _j ( ξ + _j )^4. (14) Defining the weighted tangent-space center ηlin=1α∑jπjηj _lin= 1α _j _j _j, we have 12∑jπj‖ξ−ηj‖2=α2‖ξ−ηlin‖2+ck, 12 _j _j ξ- _j ^2= α2 ξ- _lin ^2+c_k, where ck=12∑jπj‖ηj−ηlin‖2c_k= 12 _j _j _j- _lin ^2 is independent of ξ. The resulting tangent-space virtual target is ηlin _lin, which is mapped back to the manifold as q~R=Expxk(ηlin). q_R=Exp_x_k( _lin). (15) Remark 1. Note that ηlin _lin serves as a tangent-space virtual target, which coincides with the true Riemannian barycenter in flat space but differs under nonzero curvature. IV Tangent-Space Target Aggregation and LQR Control IV-A Local LQR Controller Neglecting the local geometric remainders in (9)–(14) gives the tangent-space approximation minu0(⋅)Jloc _u_0(·) J_loc =12∫tktf[α‖ξ−ηlin‖2+‖u0‖2]t = 12 _t_k^t_f [α ξ- _lin ^2+ u_0 ^2 ]dt (16) +wfα2‖ξ(tf)−ηlin‖2, + w_fα2 ξ(t_f)- _lin ^2, s.t.ξ˙ .t. ξ =ζ,ζ˙=u0, =ζ, ζ=u_0, where the terminal weight follows directly from (4). Define the tracking error e=ξ−ηline=ξ- _lin and state vector y=[e⊤,ζ⊤]⊤y=[e ,ζ ] . Since ηlin _lin is fixed during one replanning horizon, the resulting nominal linear error dynamics satisfy y˙=Ay+Bu0, y=Ay+Bu_0, where A=[0I00],B=[0I]. A= bmatrix0&I\\ 0&0 bmatrix,\,B= bmatrix0\\ I bmatrix. An exact zero-order-hold discretization with step Δt t gives ym+1=Adym+Bdumy_m+1=A_dy_m+B_du_m with Ad=[IΔtI0I],Bd=[12Δt2IΔtI] A_d= bmatrixI& tI\\ 0&I bmatrix,\,B_d= bmatrix 12 t^2I\\ tI bmatrix at discrete time index m. Approximating the running integral gives Jd=12yT⊤QfyT+∑m=0T−112(ym⊤Q¯ym+um⊤R¯um),J_d= 12y_T Q_fy_T+ _m=0^T-1 12 (y_m Qy_m+u_m Ru_m ), (17) with Q¯=Δtdiag(αI,0) Q= t\,diag(α I,0), R¯=ΔtI R= tI, and Qf=diag(wfαI,0)Q_f=diag(w_fα I,0), thus PT=QfP_T=Q_f. Backward dynamic programming gives um∗=−Kmymu_m^*=-K_my_m with Sm S_m =R¯+Bd⊤Pm+1Bd,Km=Sm−1Bd⊤Pm+1Ad, = R+B_d P_m+1B_d, K_m=S_m^-1B_d P_m+1A_d, (18) Pm P_m =Q¯+Ad⊤Pm+1Ad−Ad⊤Pm+1BdSm−1Bd⊤Pm+1Ad. = Q+A_d P_m+1A_d-A_d P_m+1B_dS_m^-1B_d P_m+1A_d. Because R¯≻0 R 0 and Pm+1⪰0P_m+1 0, Sm≻0S_m 0, the backward recursion is well-defined. Remark 2 (Fixed-Dimensional Computational Complexity). For a horizon T, solving the discrete Riccati recursion requires (Td3)O(Td^3) operations, as the state dimension of y depends solely on the manifold dimension d. Meanwhile, computing the tangent-space target ηlin _lin via weighted summation requires (|k|d)O(|S^k|d) operations. Consequently, the target sample size |k||S^k| affects only the linear vector aggregation, leaving the matrix Riccati state dimension completely invariant to sample size. IV-B Receding-Horizon Implementation for R-D2OC At replanning step k, the local control command is synthesized in the tangent space TxkℳT_x_kM and updated recursively. Given the current manifold state xk∈ℳx_k and a set of discrete target points qj\q_j\, we first compute the logarithmic map coordinates ηjk=Logxk(qj)∈Txkℳ _j^k=Log_x_k(q_j)∈ T_x_kM. The weighted spatial centroid in TxkℳT_x_kM is evaluated as ηlink=1α∑jπjηjk _lin^k= 1α _j _j _j^k with total weight α=∑jπjα= _j _j. Since xkx_k serves as the origin of the chosen orthonormal basis for TxkℳT_x_kM, the initial state vector for the local LQR policy is formulated as yk=[(−ηlink)⊤,vk⊤]⊤y_k=[(- _lin^k) ,v_k ] , where vk∈Txkℳv_k∈ T_x_kM is the current velocity vector. The feedback acceleration command is then obtained as uk=−K0yk∈Txkℳu_k=-K_0y_k∈ T_x_kM. To execute this control over the interval t∈[tk,tk+Δt]t∈[t_k,t_k+ t], uku_k is continuously pushed forward along the trajectory x(t)x(t) via parallel transport, u(t)=xk→x(t)uku(t)=P_x_k→ x(t)u_k, driving the intrinsic system dynamics x˙=v x=v and ∇x˙v=u(t) _ xv=u(t). A local second-order numerical update over the time step Δt t is then used to advance the state on ℳM: xk+1=Expxk(Δtvk+12Δt2uk), x_k+1=Exp_x_k ( tv_k+ 12 t^2u_k ), (19) vk+1=xk→xk+1(vk+Δtuk). v_k+1=P_x_k→ x_k+1 (v_k+ tu_k ). This update provides a local approximation of the intrinsic dynamics over each sampling interval. The tangent-space origin is then shifted from xkx_k to xk+1x_k+1 at each step to complete the receding-horizon control loop. V Theoretical Analysis: Approximation Bias and Local Conditional Descent V-A Objective Approximation Error Define the local objective approximation error at an arbitrary trajectory point E(x,xk)=V(x)−12∑jπj‖Logxk(x)−ηjk‖xk2, E(x;x_k)=V(x)- 12 _j _j _x_k(x)- _j^k _x_k^2, (20) where ηjk=Logxk(qj) _j^k=Log_x_k(q_j). This is exactly the gap between the true running cost and the quadratic surrogate the local LQR controller actually minimizes. Proposition 1 (Local objective error). Suppose x and all target samples belong to a common sufficiently small normal neighborhood of xkx_k. Then |E(x,xk)|≤C4κα(‖ξ‖xk+σktar)4,|E(x;x_k)|≤ C_4κα ( ξ _x_k+ _k^tar )^4, (21) where ξ=Logxk(x)ξ=Log_x_k(x) for the local coordinate of x and σktar=maxj‖ηjk‖xk _k^tar= _j _j^k _x_k. Proof. From (14), |E(x,xk)|≤C4κ∑jπj(‖ξ‖+‖ηjk‖)4|E(x;x_k)|≤ C_4κ _j _j ( ξ + _j^k )^4 ≤C4κα(‖ξ‖+σktar)4.≤ C_4κα ( ξ + _k^tar )^4. ∎ This is a pointwise approximation, not a cumulative convergence result. A cumulative bound requires additional assumptions controlling target spread and prediction error along a trajectory. V-B Bias of the Tangent-Space Center To characterize the Riemannian barycenter q¯R q_R in (2), define the weighted Fréchet energy [14] as ℱ(x)=12∑jπjdistg(x,qj)2F(x)= 12 _j _jdist_g(x,q_j)^2, whose critical point satisfies gradℱ(q¯R)=0gradF( q_R)=0. Although ℱF is identical to the state cost V in (4), we use ℱF here to emphasize its role in the barycenter analysis. To quantitatively compare q~R q_R in (15) with q¯R q_R, we require ℱF to be locally strongly convex near the barycenter. Otherwise, their intrinsic distance cannot be stably controlled. Assumption 1 (Nondegenerate local barycenter). There exists a strongly convex neighborhood BU_B containing the reference point x0x_0, all target samples qjq_j, the tangent-space virtual target q~R q_R, and the Riemannian barycenter q¯R q_R, such that the minimizing geodesics used below remain in BU_B and Hessℱ(x)⪰λminI,x∈B,HessF(x) _ I, x _B, (22) for some λmin>0 _ >0. Remark 3. Assumption 1 is standard in Riemannian optimization [14]. As ℱF is a weighted sum of squared distance functions, Hessℱ(x)=αI+(καρmax2)HessF(x)=α I+O(κα _ ^2) with α=∑jπjα= _j _j. Thus, strong convexity naturally holds within a sufficiently small normal neighborhood. Lemma 1 (Local change-of-base-point estimate). Let x0x_0 and q belong to a sufficiently small strongly convex normal neighborhood, η=Logx0(q)η=Log_x_0(q), δ∈Tx0ℳδ∈ T_x_0M, x~=Expx0(δ) x=Exp_x_0(δ). Assume |Sec|≤κ|Sec|≤κ on the neighborhood. Then, there exists C5>0C_5>0, independent of ‖η‖+‖δ‖ η + δ , such that x~→x0Logx~(q) _ x→ x_0Log_ x(q) =η−δ+r(η,δ), =η-δ+r(η,δ), (23) ‖r(η,δ)‖x0 r(η,δ) _x_0 ≤C5κ(‖η‖+‖δ‖)3. ≤ C_5κ ( η + δ )^3. Proof. This is the local third-order change-of-base-point expansion for the Riemannian logarithmic map, following from the Jacobi-field expansion of the exponential map and the variation of radial parallel transport. The Euclidean relation η−δη-δ is exact when curvature vanishes. The first nonzero curvature correction is cubic in the local displacement, with uniform remainder on a compact strongly convex neighborhood. ∎ This local expansion serves as a key technical tool in our framework, enabling us to compare gradℱgradF evaluated at two distinct base points within the same vector space, up to a controlled cubic error. Theorem 1 (Local curvature bias). Let ρmax=maxj‖Logx0(qj)‖x0 _ = _j _x_0(q_j) _x_0, ηlin=1α∑jπjLogx0(qj) _lin= 1α _j _jLog_x_0(q_j), q~R=Expx0(ηlin) q_R=Exp_x_0( _lin). Suppose the conditions of Assumption 1 hold with some λmin>0 _min>0, |Sec|≤κ|Sec|≤κ on BU_B, and ρmax≤ρ∗ _ ≤ _* for a sufficiently small ρ∗>0 _*>0. Then, there exists C6>0C_6>0, independent of ρmax _ , such that distg(q~R,q¯R)≤C6καλminρmax3.dist_g( q_R, q_R)≤ C_6 κα _ _ ^3. (24) Proof. Let ηj=Logx0(qj) _j=Log_x_0(q_j). By definition, ‖ηj‖≤ρmax _j ≤ _ and ‖ηlin‖≤ρmax _lin ≤ _ . Applying Lemma 1 with δ=ηlinδ= _lin, x~=q~R x= q_R, for each j, we have q~R→x0Logq~R(qj)=ηj−ηlin+rj,‖rj‖x0≤8C5κρmax3.P_ q_R→ x_0Log_ q_R(q_j)= _j- _lin+r_j,\; r_j _x_0≤ 8C_5κ _ ^3. (25) By the same first-variation formula [14] applied at q~R q_R, gradℱ(q~R)=−∑jπjLogq~R(qj)gradF( q_R)=- _j _jLog_ q_R(q_j). Parallel transporting to Tx0ℳT_x_0M and using (25), q~R→x0gradℱ(q~R)=−∑jπj(ηj−ηlin+rj)=−∑jπjrj,P_ q_R→ x_0gradF( q_R)=- _j _j ( _j- _lin+r_j )=- _j _jr_j, due to the fact ∑jπj(ηj−ηlin)=0 _j _j( _j- _lin)=0. Since the parallel transport preserves the Riemannian norm, we have ‖gradℱ(q~R)‖q~R≤8C5καρmax3.\|gradF( q_R)\|_ q_R≤ 8C_5κα _ ^3. (26) Let γ:[0,1]→ℳγ:[0,1] be the minimizing geodesic from q¯R q_R to q~R q_R with constant speed, remaining in BU_B by Assumption 1. Since gradℱ(q¯R)=0gradF( q_R)=0, parallel transport along γ and the fundamental theorem of calculus give q~R→q¯Rgradℱ(q~R)=∫01γ(s)→q¯R[Hessℱ(γ(s))γ˙(s)]s.P_ q_R→ q_RgradF( q_R)= _0^1P_γ(s)→ q_R [HessF(γ(s)) γ(s) ]ds. (27) Because γ˙(s) γ(s) is parallel along the geodesic, Assumption 1 gives ⟨q~R→q¯Rgradℱ(q~R),γ˙(0)⟩q¯R _ q_R→ q_RgradF( q_R), γ(0) _ q_R ≥λmin∫01‖γ˙(s)‖2s ≥ _ _0^1 γ(s) ^2ds (28) =λmindistg(q¯R,q~R)2. = _ dist_g( q_R, q_R)^2. By Cauchy–Schwarz, λmindistg(q¯R,q~R)≤‖gradℱ(q~R)‖ _ dist_g( q_R, q_R)≤ ( q_R) , which combined with (26) gives (24) with C6=8C5C_6=8C_5. ∎ V-C Local Descent under Tangent-Space LQR In the receding-horizon implementation, the target q¯ q is instantiated at each replanning step as the tangent-space target q~Rk q_R^k. The following result directly connects the proposed LQR controller to local descent. In particular, the required alignment condition is imposed on the velocity generated by the tangent-space LQR and explicitly accounts for the local target spread. For a fixed target q¯∈ℳ q , define Ψ(x)=12distg(x,q¯)2,r(x)=Logx(q¯). (x)= 12dist_g(x, q)^2, r(x)=Log_x( q). (29) The first variation formula [14] gives dΨx[v]=−⟨r(x),v⟩x.d _x[v]=- r(x),v _x. (30) Assumption 2 (LQR alignment and local boundedness). For the tangent-space LQR controller in Section IV, there exist constants c1,c2,Γ,Lv,Lu>0c_1,c_2, ,L_v,L_u>0 and a sufficiently small neighborhood such that, at each replanning step, ⟨rk,vk⟩xk≥c1‖rk‖xk2−c2κρk3, r_k,v_k _x_k≥ c_1 r_k _x_k^2-c_2κ _k^3, (31) where rk=Logxk(q¯),ρk=max‖rk‖xk,σktar,σktar=maxj‖ηjk‖xk,r_k=Log_x_k( q),\; _k= \ r_k _x_k, _k^tar \,\; _k^tar= _j _j^k _x_k, and σktar≤Γ‖rk‖xk. _k^tar≤ r_k _x_k. (32) Moreover, along the corresponding LQR closed-loop segment t∈[tk,tk+Δt]t∈[t_k,t_k+ t], ‖v(t)‖x(t)≤Lvρk,‖u(t)‖x(t)≤Luρk. v(t) _x(t)≤ L_v _k, u(t) _x(t)≤ L_u _k. (33) The constants are uniform over the selected local neighborhood. Theorem 2 (Local descent under tangent-space LQR). Suppose the trajectory generated by the tangent-space LQR remains in a compact strongly convex neighborhood on which Ψ has uniformly bounded Hessian. Under Assumption 2, there exist ρ∗>0ρ^*>0, Δt∗>0 t^*>0, and C7>0C_7>0 such that, whenever ρk≤ρ∗ _k≤ρ^* and 0<Δt≤Δt∗0< t≤ t^*, Ψ(xk+1)−Ψ(xk)≤−C7Δt‖rk‖xk2. (x_k+1)- (x_k)≤-C_7 t r_k _x_k^2. (34) Thus, the proposed tangent-space LQR yields local descent whenever its closed-loop velocity satisfies the gain-dependent alignment condition (31). Proof. Let ρk=max‖rk‖xk,σktar _k= \ r_k _x_k, _k^tar\. By (32), ρk≤Γρ‖rk‖xk,Γρ:=max1,Γ. _k≤ _ρ r_k _x_k, _ρ:= \1, \. (35) Along the LQR closed-loop trajectory, d2dt2Ψ(x(t))=HessΨ(x(t))[v(t),v(t)]+⟨gradΨ(x(t)),u(t)⟩x(t). d^2dt^2 (x(t))=Hess (x(t))[v(t),v(t)]+ (x(t)),u(t) _x(t). (36) Since the trajectory remains in a compact strongly convex neighborhood, there exists CH>0C_H>0 such that ‖HessΨ‖≤CH ≤ C_H. Moreover, ‖gradΨ(x(t))‖=distg(x(t),q¯) (x(t)) =dist_g(x(t), q) is locally bounded by CrρkC_r _k for some uniform Cr>0C_r>0. Using (33) therefore gives |d2dt2Ψ(x(t))|≤(CHLv2+CrLu)ρk2=:C8ρk2. | d^2dt^2 (x(t)) |≤ (C_HL_v^2+C_rL_u ) _k^2=:C_8 _k^2. (37) Taylor’s formula consequently yields Ψ(xk+1)−Ψ(xk) (x_k+1)- (x_k) ≤−Δt⟨rk,vk⟩xk+12C8Δt2ρk2. ≤- t r_k,v_k _x_k+ 12C_8 t^2 _k^2. (38) Using (31) and (35), Ψ(xk+1)−Ψ(xk) (x_k+1)- (x_k) ≤−c1Δt‖rk‖xk2+c2κΔtρk3+12C8Δt2ρk2 ≤-c_1 t r_k _x_k^2+c_2κ t _k^3+ 12C_8 t^2 _k^2 (39) ≤[−c1+c2κΓρ3ρk+12C8Γρ2Δt]Δt‖rk‖xk2. ≤ [-c_1+c_2κ _ρ^3 _k+ 12C_8 _ρ^2 t ] t r_k _x_k^2. Choose ρ∗>0ρ^*>0 and Δt∗>0 t^*>0 such that c2κΓρ3ρ∗≤c14,12C8Γρ2Δt∗≤c14.c_2κ _ρ^3ρ^*≤ c_14, 12C_8 _ρ^2 t^*≤ c_14. (40) Then, for ρk≤ρ∗ _k≤ρ^* and 0<Δt≤Δt∗0< t≤ t^*, Ψ(xk+1)−Ψ(xk)≤−c12Δt‖rk‖xk2. (x_k+1)- (x_k)≤- c_12 t r_k _x_k^2. (41) Hence (34) holds with C7=c1/2C_7=c_1/2. ∎ Remark 4. The alignment condition (31) is imposed on the velocity generated by the proposed tangent-space LQR, while (33) captures the corresponding local closed-loop bounds. The additional inequality (32) ensures that the target-spread term introduced by the tangent-space approximation is of the same local scale as the tracking error. Therefore, the result provides a direct sufficient condition under which the proposed LQR inherits local descent on the nonlinear manifold. VI Simulation Results To validate the theoretical guarantees, we implement the proposed R-D2OC framework on a 3D ellipsoidal manifold, where multiple second-order agents coordinate to match a target distribution over the curved geometry. Fig. 1 depicts the overall system behavior and spatial trajectories. The agents successfully track the reference distribution along the Riemannian manifold while maintaining smooth coordination paths. Originating from tightly clustered initial positions (black crosses), trajectories evolve smoothly across the surface to match their time-averaged distribution with the targets (blue dots). This confirms that the second-order tangent-space LQR formulation respects geometric constraints while driving the team toward the desired distribution. To evaluate the geometric approximation properties and empirically assess Theorem 1, we analyze the log-log scaling relationship in Fig. 2(a). The empirical approximation bias exhibits a fitted slope of 2.732.73, which is consistent with the theoretical cubic-order (ρmax3)O( _ ^3) upper bound. This observation supports the predicted local scaling of the tangent-space center bias as the target spread varies. We next evaluate the local descent property in Theorem 2 using Fig. 2(b). The closed-loop trajectories satisfy the alignment condition of Assumption 2 with c1=1.40c_1=1.40, yielding C7=c1/2=0.70C_7=c_1/2=0.70 in (34). Within the empirical local neighborhood ρk≤ρ∗≈4.1 _k≤ρ^*≈ 4.1 (shaded region), the observed objective decrease aligns well with the theoretical bound −C7Δt‖rk‖2-C_7 t\|r_k\|^2. Beyond ρ∗ρ^*, this local bound no longer holds, which is consistent with the conditional and local nature of Theorem 2. With an average per-agent Riccati solve time of 0.49ms0.49~ms relative to the 0.1s0.1~s replanning interval Δt t, the computation remains well within the budget for onboard execution. Fig. 1: R-D2OC result: Three-dimensional multi-agent trajectories and reference distribution on the Riemannian ellipsoid manifold. (a) (b) Fig. 2: Empirical verification of main theoretical results. (a) Approximation bias distg(q~R,q¯R)dist_g( q_R, q_R) versus target spread, validating the third-order scaling in Theorem 1. (b) Empirical descent rate −ΔΨk/Δt- _k/ t versus ρk2 _k^2, with the shaded region marking the neighborhood ρk≤ρ∗≈4.1 _k≤ρ^*≈ 4.1 where Theorem 2 holds. VII Conclusion This paper presented an intrinsic R-D2OC framework for second-order multi-agent systems on Riemannian manifolds. By mapping the intrinsic transport problem to local tangent spaces, the framework reduces online optimization to linear-quadratic subproblems with computational complexity depending strictly on the manifold dimension. We provided rigorous theoretical guarantees, including third-order barycentric bias bounds and a conditional discrete-time descent property. Numerical simulations on the sphere validated these theoretical results, confirming efficient distribution control on curved spaces. References [1] S. Seo and K. Lee (2025) Density-driven optimal control for efficient and collaborative multiagent nonuniform coverage. IEEE Transactions on Systems, Man, and Cybernetics: Systems 55 (12), p. 9340–9354. External Links: Document Cited by: §I, §I-A. [2] S. Seo and K. Lee (2026) Density-driven multi-drone coordination for efficient farm coverage and management in smart agriculture. IEEE Transactions on Control Systems Technology 34 (2), p. 711–724. External Links: Document Cited by: §I, §I-A. [3] K. Lee and E. Brook (2025) Connectivity-preserving multi-agent area coverage via density-driven optimal control (d2^2oc). IEEE Control Systems Letters 9, p. 2723–2728. External Links: Document Cited by: §I, §I-A. [4] K. Lee (2026) Optimal transport-based decentralized multi-agent distribution matching. IEEE Transactions on Automatic Control 71 (8), p. 5478–5485. External Links: Document Cited by: §I. [5] K. Lee (2026) Density-driven optimal control: convergence guarantees for stochastic lti multi-agent systems. Automatica 193, p. 113231. External Links: ISSN 0005-1098, Document Cited by: §I, §I-A. [6] F. Bullo and A. D. Lewis (2005) Geometric control of mechanical systems. Cited by: §I. [7] T. Lee, M. Leok, and N. H. McClamroch (2008) Computational geometric optimal control of rigid bodies. Communications in Information and Systems 8 (4), p. 445–472. Cited by: §I. [8] T. Lee (2015) Global exponential attitude tracking controls on so(3). IEEE Transactions on Automatic Control 60 (10), p. 2837–2842. Cited by: §I. [9] K. Wu, T. Ho, S. A. Huang, K. Lin, Y. Lin, and J. Liu (2016) Path planning and replanning for mobile robot navigation on 3d terrain: an approach based on geodesic. Mathematical Problems in Engineering 2016 (1), p. 2539761. Cited by: §I. [10] Y. Chen, L. Li, and X. Ji (2014) Smooth and accurate trajectory planning for industrial robots. Advances in Mechanical Engineering 6, p. 342137. Cited by: §I. [11] J. Cortés, S. Martínez, T. Karatas, and F. Bullo (2004) Coverage control for mobile sensing networks. IEEE Transactions on Robotics and Automation 20 (2), p. 243–255. External Links: Document Cited by: §I. [12] F. Bullo, J. Cortés, and S. Martinez (2009) Distributed control of robotic networks. Applied Mathematics Series. Princeton University Press 3 (7), p. 4–4. Cited by: §I. [13] V. Krishnan and S. Martínez (2022) A multiscale analysis of multi-agent coverage control algorithms. Automatica 146, p. 110516. External Links: Document Cited by: §I. [14] P. Absil, R. Mahony, and R. Sepulchre (2008) Optimization algorithms on matrix manifolds. Princeton University Press. Cited by: §V-B, §V-B, §V-C, Remark 3.