Paper deep dive
Computationally Efficient Density-Driven Optimal Control via Analytical KKT Reduction and Contractive MPC
Julian Martinez, Kooktae Lee
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 96%
Last extracted: 3/22/2026, 5:03:04 AM
Summary
The paper introduces a computationally efficient Density-Driven Optimal Control (D2OC) framework for multi-agent systems. By deriving an analytical structural reduction of the T-horizon Karush-Kuhn-Tucker (KKT) system, the authors transform the problem into a condensed quadratic program (QP) with O(T) linear scalability, significantly improving upon the conventional O(T^3) complexity. Additionally, the paper establishes Input-to-State Stability (ISS) for the closed-loop system using a contractive Lyapunov constraint, ensuring robust performance in dynamic environments.
Entities (5)
Relation Signals (3)
Analytical structural reduction → transforms → KKT system
confidence 98% · we propose an analytical structural reduction that transforms the T-horizon KKT system into a condensed quadratic program (QP)
Contractive Lyapunov constraint → ensures → Input-to-State Stability
confidence 96% · we incorporate a contractive Lyapunov constraint and prove the Input-to-State Stability (ISS) of the closed-loop system
D2OC → utilizes → KKT system
confidence 95% · implementing this as a predictive controller requires solving a large-scale Karush-Kuhn-Tucker (KKT) system
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Efficient coordination for collective spatial distribution is a fundamental challenge in multi-agent systems. Prior research on Density-Driven Optimal Control (D2OC) established a framework to match agent trajectories to a desired spatial distribution. However, implementing this as a predictive controller requires solving a large-scale Karush-Kuhn-Tucker (KKT) system, whose computational complexity grows cubically with the prediction horizon. To resolve this, we propose an analytical structural reduction that transforms the T-horizon KKT system into a condensed quadratic program (QP). This formulation achieves O(T) linear scalability, significantly reducing the online computational burden compared to conventional O(T^3) approaches. Furthermore, to ensure rigorous convergence in dynamic environments, we incorporate a contractive Lyapunov constraint and prove the Input-to-State Stability (ISS) of the closed-loop system against reference propagation drift. Numerical simulations verify that the proposed method facilitates rapid density coverage with substantial computational speed-up, enabling long-horizon predictive control for large-scale multi-agent swarms.
Tags
Links
- Source: https://arxiv.org/abs/2603.18503v1
- Canonical: https://arxiv.org/abs/2603.18503v1
Trouble viewing inline? Open PDF directly →
Full Text
38,218 characters extracted from source content.
Expand or collapse full text
Computationally Efficient Density-Driven Optimal Control via Analytical KKT Reduction and Contractive MPC Julian Martinez1 and Kooktae Lee1 *This work was supported by NSF CAREER Grant CMMI-DCSD-2145810.1Julian Martinez and Kooktae Lee are with the Department of Mechanical Engineering, New Mexico Institute of Mining and Technology, Socorro, NM 87801, USA, email: julian.martinez@student.nmt.edu, kooktae.lee@nmt.edu. (June 2025) Abstract Efficient coordination for collective spatial distribution is a fundamental challenge in multi-agent systems. Prior research on Density-Driven Optimal Control (D2OC) established a framework to match agent trajectories to a desired spatial distribution. However, implementing this as a predictive controller requires solving a large-scale Karush-Kuhn-Tucker (KKT) system, whose computational complexity grows cubically with the prediction horizon. To resolve this, we propose an analytical structural reduction that transforms the T-horizon KKT system into a condensed quadratic program (QP). This formulation achieves O(T)O(T) linear scalability, significantly reducing the online computational burden compared to conventional O(T3)O(T^3) approaches. Furthermore, to ensure rigorous convergence in dynamic environments, we incorporate a contractive Lyapunov constraint and prove the Input-to-State Stability (ISS) of the closed-loop system against reference propagation drift. Numerical simulations verify that the proposed method facilitates rapid density coverage with substantial computational speed-up, enabling long-horizon predictive control for large-scale multi-agent swarms. I Introduction Multi-agent systems (MAS) have emerged as a transformative framework for executing complex tasks through the coordination of autonomous entities such as UAVs and mobile robots. A fundamental challenge in this coordination is the spatial distribution control problem. Traditionally viewed as a static coverage task, modern missions in precision agriculture [12], environmental monitoring [9], and surveillance [2] now require agents to manage their collective presence as a continuous process. This necessitates a density matching approach where the swarm’s spatial histogram tracks a desired distribution over time. Specifically, many scenarios demand that the swarm maintain a specific time-averaged presence over a given horizon to ensure thorough and continuous coverage of the task domain [11]. Translating macroscopic distribution goals into microscopic control laws is addressed through Eulerian and Lagrangian perspectives. The Eulerian approach treats the MAS as a continuous field governed by Partial Differential Equations (PDEs). Notable methods include mean-field feedback control [14] and backstepping-based formation control [3]. A related class of Eulerian techniques, such as Heat Equation Driven Area Coverage (HEDAC) [4], utilizes virtual potential fields to guide agents. While these frameworks provide rigorous analysis, they often focus on static or terminal matching. Even recent distributed online optimization techniques [5] suffer from slow convergence or increased communication overhead as the predictive horizon extends. Conversely, the Lagrangian approach focuses directly on discrete agent trajectories. While traditional approaches focus on static Voronoi-based coverage [1] or stochastic vehicle routing for persistent service [10], and recent optimal transport-based frameworks [7] primarily address terminal matching, significant progress has also been made in time-averaged behavior through ergodicity. A pioneering Lagrangian contribution is the spectral multiscale coverage framework [8], which uses Fourier-based metrics to track target densities. Although robust, these spectral methods are primarily gradient-driven and often struggle to incorporate high-order dynamics or ensure formal optimality regarding control effort, especially over complex, multi-step predictive horizons. Building upon these principles, the Density-Driven Optimal Control (D2OC) framework provides a rigorous solution for optimal time-averaged distribution matching [11, 12, 6]. However, existing D2OC methodologies are hindered by O(T3)O(T^3) computational complexity and lack robust stability guarantees against reference drift. This paper overcomes these bottlenecks via a Contractive Lyapunov-based MPC framework with the following primary contributions: 1) Analytical KKT Condensation: Unlike prior D2OC works limited by O(T3)O(T^3) complexity, we derive an exact analytical reduction to m-dimensional QPs. This reduces the online cost to O(T)O(T), enabling long-horizon implementation; 2) ISS under Reference Drift: We establish the first rigorous stability framework for D2OC using a contractive Lyapunov constraint, proving the Input-to-State Stability (ISS) of the closed-loop system against reference propagation drift; 3) Real-time Decentralized Scalability: The proposed framework enables independent local optimization while maintaining global distribution matching. Simulations verify that the linear computational growth allows for high-fidelity control of large-scale swarms in real-time. I Problem Description Notations: ℝnR^n and ℝn×mR^n× m denote n-dimensional real vectors and n×mn× m real matrices, respectively. nI_n is the n×n× n identity matrix, and T1_T denotes a T-dimensional vector of ones. The matrix transpose is A⊤A , and ⊗ denotes the Kronecker product. Q⪰Q 0 (R≻R 0) indicates that Q is positive semi-definite (R is positive definite). A discrete distribution is represented by a set of pairs (qj,γj)\(q_j, _j)\, where qjq_j and γj _j are the position and weight of the j-th sample point, respectively. I-A Multi-Agent System Modeling Consider a multi-agent system consisting of N agents. Each agent i∈1,…,Ni∈\1,…,N\ is governed by discrete-time Linear Time-Invariant (LTI) dynamics: ik+1=Aiik+Biuik,yik=Ciik, x_i^k+1=A_i x_i^k+B_iu_i^k, y_i^k=C_i x_i^k, (1) where k is a discrete-time index, i∈ℝn x_i ^n is the state, ui∈ℝmu_i ^m is the control input, and yi∈ℝdy_i ^d is the agent’s position. It is important to note that the proposed D2OC framework inherently supports fully distributed control, as the optimization for each agent is decoupled. While we consider heterogeneous dynamics (Ai,Bi,CiA_i,B_i,C_i) to emphasize the generality of our approach, the subscripts are omitted in the following sections for notational simplicity without loss of generality. I-B Decentralized Density-Driven Optimal Control (D2OC) The objective of D2OC is to drive the empirical time-averaged distribution μk=(yik,αi)i=1Mμ^k=\(y_i^k, _i)\_i=1^M to match a prescribed reference distribution ν=(qj,βj)j=1NSPν=\(q_j, _j)\_j=1^N_SP. Here, yiky_i^k and qjq_j denote the spatial positions, while the weights are uniformly defined as αi=1/M _i=1/M and βj=1/NSP _j=1/N_SP, representing the normalized contribution of each sample point to its respective distribution. Although the details are described in [12, 11], the D2OC framework is composed of a three-stage iterative cycle as follows: • Stage A: Local Sample Selection & Control: Agent i identifies a local set of sample points (SPs), kS^k, near the current agent position yiky_i^k, and solves a T-horizon optimization problem to align its trajectory with the local density. • Stage B: Weight Adjustment: After movement, agent i reduces the weight of SPs at its current location by a certain mass amount to account for coverage progress. • Stage C: Information Exchange: Agents share updated SPs’ weights with neighbors within a communication range threshold to maintain a consistent view of the global coverage state. Iterative execution of these stages ensures that the multi-agent system effectively sweeps high-priority regions encoded by the SPs. This paper focuses on Stage A, where the matching between the agent’s future trajectory and the local SPs is formulated via the Wasserstein distance, which serves as a metric to quantify the dissimilarity between two distributions [13]. The optimization problem for agent i at time k is formulated with horizon T by: minui _u_i J=Φ(ik+T)+∑ℓ=k+T−1L(iℓ,uiℓ) J= ( x_i^k+T)+ _ =k^k+T-1L( x_i ,u_i ) (2) s.t. ik+1=Aik+Buik, x_i^k+1=A x_i^k+Bu_i^k, where the stage cost L(⋅)L(·) and terminal cost Φ(⋅) (·) are defined as: L(iℓ,uiℓ) L( x_i ,u_i ) =12(iℓ|k)2+12(iℓ)⊤Qiℓ+12(uiℓ)⊤Ruiℓ, = 12(W_i |k)^2+ 12( x_i ) Q x_i + 12(u_i ) Ru_i , Φ(ik+T) ( x_i^k+T) =12(ik+T|k)2+12(ik+T)⊤Qik+T. = 12(W_i^k+T|k)^2+ 12( x_i^k+T) Q x_i^k+T. Here, Q⪰Q 0 and R≻R 0 are weighting matrices. The term (iℓ|k)2(W_i |k)^2, referred to as the local Wasserstein distance, is defined as: (iℓ|k)2=minγj≥0∑j∈kγj‖yiℓ−qj‖2,(W_i |k)^2= _ _j≥ 0 _j ^k _j\|y_i -q_j\|^2, (3) where qjq_j denotes the SP position and γj _j the mass transport plan from each yiy_i to qjq_j. Note that index ℓ denotes the predicted time step within the horizon T. Theorem 1 (Analytical Solution [11]) For the LTI system (1), the optimal control sequence u¯ik:=[(uik)⊤,⋯,(uik+T−1)⊤]⊤ u_i^k:=[(u_i^k) ,·s,(u_i^k+T-1) ] minimizing (2) is analytically given by: u¯ik=ℰ13⊤F1+ℰ23⊤F2, u_i^k=E_13 F_1+E_23 F_2, (4) where ℰijE_ij denotes the (i,j)(i,j)-th block submatrix of E−1E^-1 and the constituent blocks of the KKT matrix E are defined as: E11=Q¯⊗T,E23=E32⊤=B⊗T,E33=R⊗T, E_11= Q _T,\ E_23=E_32 =B _T,\ E_33=R _T, E13=E31=E22=,Q¯=(∑j∈kγj)C⊤C+Q, E_13=E_31=E_22=0,\ Q= \ ( _j ^k _j )C C+Q \, E12=E21⊤=[−nA⊤⋯−nA⊤⋱⋮⋱⋮⋱−nA⊤⋯−n]. E_12=E_21 = $ bmatrix-I_n&A &0&·s&0\\ 0&-I_n&A & & \\ & & & &0\\ && &-I_n&A \\ 0&·s&·s&0&-I_n bmatrix$. (5) Moreover, F1F_1 and F2F_2 in (4) are defined by: F1=[(∑j∈kγj)C⊤q¯k⋮(∑j∈kγj)C⊤q¯k],F2=[−Aik⋮], F_1= bmatrix( _j ^k _j)C q^k\\ \\ ( _j ^k _j)C q^k bmatrix,\,F_2= bmatrix-A x_i^k\\ 0\\ \\ 0 bmatrix, where q¯k=(∑j∈kγjqj)/(∑j∈kγj) q^k=( _j ^k _jq_j)/( _j ^k _j) denotes the local barycenter of the sample points within kS^k. Remark 1 By utilizing the analytical block-inversion of the KKT matrix E, the D2OC scheme avoids numerical instabilities associated with indefinite saddle-point problems, enabling real-time decentralized control. However, as the horizon T increases, the dimension of the KKT system expands, naturally increasing the computational burden for direct matrix inversion. I Computational Efficiency Improvement via Multi-Step Reference Structural Reduction I-A From Fixed Target to Multi-Stage Prediction In our previous work [11], the optimal control utilized a stationary reference barycenter q¯k q^k fixed at time k. However, as an agent’s movement progressively reduces nearby sample weights (Stage B), the local target barycenter systematically shifts even under a static global distribution. To anticipate this dynamic drift, we propose a multi-stage MPC framework incorporating a sequence of predicted barycenters q¯k+ℓ+1ℓ=0T−1\ q^k+ +1\_ =0^T-1, allowing the full T-horizon trajectory to align with the evolving local objectives inherent in the coverage process. Definition 1 (Iterative Sample Set and Weight Propagation) For each prediction step ℓ∈0,…,T−1 ∈\0,…,T-1\, the local sample set k+ℓS^k+ is collected around the predicted agent position. Consequently, the time-varying state penalty matrix Q¯k+ℓ Q^k+ and target barycenter q¯k+ℓ+1 q^k+ +1 are recursively defined as: Q¯k+ℓ=(∑j∈k+ℓγj)C⊤C+Q,q¯k+ℓ+1=∑j∈k+ℓγjqj∑j∈k+ℓγj. Q^k+ = ( _j ^k+ _j )C C+Q, q^k+ +1= _j ^k+ _jq_j _j ^k+ _j. (6) I-B Analytical Condensation to m-Dimensional Quadratic Form Although incorporating multi-stage goals increases the theoretical complexity, we prove that the high-dimensional T-horizon problem originally formulated over the stacked input sequence u¯ik∈ℝmT u_i^k ^mT in Theorem 1 can be exactly reduced to an efficient m-dimensional QP in terms of the actual control input uik∈ℝmu_i^k ^m. By incorporating T-step preview information into the condensed objective function, the agent can strategically pursue coverage progress while maintaining a computational complexity that scales only linearly with the horizon T. This ensures that the online burden remains comparable to that of short-horizon controllers, enabling high-fidelity predictive control in real-time. Theorem 2 (Analytical Reduction of Time-Varying MPC) The optimal first-step control uiku_i^k for the T-horizon problem (2) is the unique solution to the condensed QP: minuik _u_i^k 12(uik)⊤Tuik+T⊤uik s.t.umin≤uik≤umax 12(u_i^k) H_Tu_i^k+g_T u_i^k s.t. u_ ≤ u_i^k≤ u_ (7) where T∈ℝm×mH_T ^m× m and T∈ℝmg_T ^m are analytically given by: T _T =R+∑ℓ=0T−1B⊤(Aℓ)⊤Q¯k+ℓAℓB, =R+ _ =0^T-1B (A ) Q^k+ A B, (8) T _T =∑ℓ=0T−1B⊤(Aℓ)⊤Q¯k+ℓ(Aℓ+1ik−q¯k+ℓ+1). = _ =0^T-1B (A ) Q^k+ (A +1 x_i^k- q^k+ +1). (9) Proof: The optimal control law in (4) is derived from the linear optimality conditions, which can be expressed as a KKT system of the form Ez=FEz=F. Here, z=[¯i⊤,λ¯i⊤,u¯i⊤]⊤z=[ x_i , λ_i , u_i ] is the vector of decision variables consisting of the sequences of states, co-states (Lagrange multipliers), and control inputs over the prediction horizon T. The KKT matrix E encapsulates the coupled dynamics and optimality constraints, while the right-hand side vector F=[F1⊤,F2⊤,⊤]⊤F=[F_1 ,F_2 ,0 ] incorporates the initial state and the time-indexed sequence of local barycenters q¯k+ℓ+1\ q^k+ +1\. By exploiting the specific block-sparse structure of E, one can obtain ¯i=E12−⊤(F2−E23u¯ik) x_i=E_12^- (F_2-E_23 u_i^k), λ¯i=E12−1(F1−E11E12−⊤(F2−E23u¯ik)) λ_i=E_12^-1(F_1-E_11E_12^- (F_2-E_23 u_i^k)). Finally, substituting λ¯i λ_i into the third row E23⊤λ¯i+E33u¯ik=E_23 λ_i+E_33 u_i^k=0 yields the condensed system ℋu¯ik=H u_i^k=G, where: ℋ =E33+E23⊤E12−1E11E12−⊤E23, =E_33+E_23 E_12^-1E_11E_12^- E_23, (10) =E23⊤E12−1(E11E12−⊤F2−F1). =E_23 E_12^-1(E_11E_12^- F_2-F_1). (11) Note that E12−1E_12^-1 is an upper block-triangular operator representing the adjoint state transition, which allows for the analytical reduction of the KKT system into a dense QP problem. The control influence is captured by the block vector =[0⊤,…,T−1⊤]⊤K=[K_0 ,…,K_T-1 ] , where each block is defined as ℓ=−AℓBK_ =-A B. The condensed Hessian TH_T, which is the (1,1)(1,1)-block of ℋH, is then derived using E11=diag(Q¯k,…,Q¯k+T−1)E_11=diag( Q^k,…, Q^k+T-1): T=R+∑ℓ=0T−1ℓ⊤Q¯k+ℓℓ=R+∑ℓ=0T−1B⊤(Aℓ)⊤Q¯k+ℓAℓB.H_T=R+ _ =0^T-1K_ Q^k+ K_ =R+ _ =0^T-1B (A ) Q^k+ A B. (12) For the linear term Tg_T, we extract the first m-dimensional block of G defined in (11). By observing the block-structure of E23⊤E12−1E_23 E_12^-1, the first block-row is analytically given by [−B⊤,−B⊤A⊤,…,−B⊤(AT−1)⊤][-B ,-B A ,…,-B (A^T-1) ]. Substituting the initial condition F2=[−Aik,,…]⊤F_2=[-A x_i^k,0,…] and the target-dependent vector F1=[Q¯k+1q¯k+1,…,Q¯k+Tq¯k+T]⊤F_1=[ Q^k+1 q^k+1,…, Q^k+T q^k+T] , the first block Tg_T is evaluated as: T=∑ℓ=0T−1(AℓB)⊤Q¯k+ℓ(Aℓ+1ik−q¯k+ℓ+1)g_T= _ =0^T-1(A B) Q^k+ (A +1 x_i^k- q^k+ +1) (13) Note that for ℓ=0 =0, the term B⊤Q¯k(Aik−q¯k+1)B Q^k(A x_i^k- q^k+1) correctly incorporates the initial state transition. This confirms the exact analytical reduction to an m-dimensional QP. ∎ I-C Computational Efficiency Analysis Proposition 1 (Computational Efficiency of Reduced KKT) The structural reduction in Theorem 2 achieves a complexity transformation from cubic to linear growth with respect to the horizon T as follows: 1. Full KKT Approach: Direct factorization of the (2n+m)T(2n+m)T-dimensional system requires (T3(2n+m)3)O(T^3(2n+m)^3) operations, scaling cubically with the horizon. 2. Proposed Reduced KKT: By utilizing the recursive construction of TH_T and Tg_T, the online cost is reduced to (Tn2m)O(Tn^2m) for matrix updates and (m3)O(m^3) for solving the m-dimensional QP. 3. Acceleration Factor: For T≫1T 1, the computational speedup ratio relative to the full KKT approach is: ℛ=CostKKTCostProposed≈(2n+m)3T3Tn2m+m3≈(T2).R= Cost_ KKTCost_ Proposed≈ (2n+m)^3T^3Tn^2m+m^3 (T^2). (14) Proof: The conventional KKT cost KKTC_KKT depends on the dimension of the stacked vector z∈ℝ(2n+m)Tz ^(2n+m)T, leading to KKT≈(((2n+m)T)3)C_KKT (((2n+m)T)^3). In the proposed method, we exploit the iterative relation Aℓ+1=A⋅AℓA +1=A· A to compute the sequences AℓB\A B\ and (Aℓ)⊤Q¯\(A ) Q\. This recursive preparation allows TH_T and Tg_T to be constructed in T steps, each involving matrix-matrix products of order (n2m)O(n^2m). The total cost propC_prop is then: prop _prop =∑ℓ=0T−1(n2m)⏟Recursive Updates+(m3)⏟QP Solver=(T⋅n2m)+(m3). = _ =0^T-1O(n^2m)_Recursive Updates+ O(m^3)_QP Solver=O(T· n^2m)+O(m^3). (15) As T increases, the cubic term (T3)O(T^3) in the KKT approach is dominated by the linear term (T)O(T) in the proposed method, yielding an asymptotic acceleration of (T2)O(T^2). This ensures real-time feasibility for large-scale systems with long-horizon preview. ∎ IV Contractive Lyapunov Constraint for Receding-Horizon D2OC In the previous section, we showed that the multi-goal D2OC can be analytically reduced to an m-dimensional QP, achieving O(T)O(T) computational complexity. However, since the reference barycenters q¯k+ℓ\ q^k+ \ are dynamically updated via iterative sample selection and the spatial decay of density weights, the closed-loop stability cannot be directly inferred from standard infinite-horizon analysis. In this section, we introduce a contractive Lyapunov constraint designed as a Linear Matrix Inequality (LMI). This approach enforces solver-level stability at each MPC update, rigorously distinguishing the proposed multi-step scheme from a purely reactive horizon-1 QP. IV-A LMI-Based Stability Design via Schur Complement To analyze the tracking performance, we define the time-varying tracking error as ek:=ik−q¯ke^k:= x_i^k- q^k. Based on the discrete-time LTI dynamics (1), the error evolution is governed by: ek+1=Aek+Buik+dk,e^k+1=Ae^k+Bu_i^k+d^k, (16) where dk:=Aq¯k−q¯k+1d^k:=A q^k- q^k+1 represents the reference propagation drift. To guarantee the stability of the error dynamics, we employ a quadratic Lyapunov candidate V(ek):=(ek)⊤PekV(e^k):=(e^k) Pe^k with a symmetric positive definite matrix P≻0P 0. The core of our stability design is to enforce a strict contraction condition, V(ek+1)−V(ek)≤−(ek)⊤QcekV(e^k+1)-V(e^k)≤-(e^k) Q_ce^k, where Qc≻0Q_c 0 determines the convergence rate. To incorporate this requirement into a convex optimization framework, we first formulate it as an LMI. Proposition 2 (Schur-Complement Contractive Constraint) The error dynamics (16) satisfies the contraction condition V(ek+1)−V(ek)≤−(ek)⊤QcekV(e^k+1)-V(e^k)≤-(e^k) Q_ce^k if and only if the control input uiku_i^k satisfies the following LMI: [P−Qc(Aek+Buik+dk)⊤P(Aek+Buik+dk)P]⪰0. bmatrixP-Q_c&(Ae^k+Bu_i^k+d^k) P\\ P(Ae^k+Bu_i^k+d^k)&P bmatrix 0. (17) Proof: Consider a symmetric block matrix M=[XY⊤YZ]M= bmatrixX&Y \\ Y&Z bmatrix. According to the Schur complement lemma, M⪰0M 0 if and only if Z≻0Z 0 and X−Y⊤Z−1Y⪰0X-Y Z^-1Y 0. By mapping the terms in (17) such that X=P−QcX=P-Q_c, Y=P(Aek+Buik+dk)Y=P(Ae^k+Bu_i^k+d^k), and Z=PZ=P, the condition X−Y⊤Z−1Y⪰0X-Y Z^-1Y 0 becomes: (P−Qc)−(Aek+Buik+dk)⊤P⊤P−1P(Aek+Buik+dk)⪰0.(P-Q_c)-(Ae^k+Bu_i^k+d^k) P P^-1P(Ae^k+Bu_i^k+d^k) 0. Since P is symmetric (P⊤=P =P), the term simplifies to PP−1P=P^-1P=P. Substituting this and rearranging the inequality, we obtain: (Aek+Buik+dk)⊤P(Aek+Buik+dk)≤(ek)⊤(P−Qc)ek.(Ae^k+Bu_i^k+d^k) P(Ae^k+Bu_i^k+d^k)≤(e^k) (P-Q_c)e^k. Noting that ek+1=Aek+Buik+dke^k+1=Ae^k+Bu_i^k+d^k and V(ek+1)=(ek+1)⊤Pek+1V(e^k+1)=(e^k+1) Pe^k+1, the above is equivalent to: V(ek+1)≤(ek)⊤Pek−(ek)⊤Qcek,V(e^k+1)≤(e^k) Pe^k-(e^k) Q_ce^k, which yields V(ek+1)−V(ek)≤−(ek)⊤QcekV(e^k+1)-V(e^k)≤-(e^k) Q_ce^k. ∎ IV-B Numerical Reformulation for Real-Time Execution While Proposition 2 provides a robust theoretical foundation for stability via the LMI (17), solving a Semidefinite Program (SDP) in real-time is computationally demanding for high-speed multi-agent coordination. To ensure numerical tractability, we leverage the Cholesky factorization of the Lyapunov matrix, P=L⊤LP=L L, to transform the matrix inequality into a Second-Order Cone (SOC) constraint. By utilizing the properties of the Euclidean norm, the quadratic form (ek+1)⊤Pek+1(e^k+1) Pe^k+1 is expressed as ‖Lek+1‖22\|Le^k+1\|_2^2. Similarly, the upper bound (ek)⊤(P−Qc)ek(e^k) (P-Q_c)e^k can be rewritten using the relation P−Qc=P(I−P−1Qc)P-Q_c=P(I-P^-1Q_c). Taking the square root of both sides, the contraction condition is reformulated as the following SOC constraint: ‖L(Aek+Buik+dk)‖2≤‖LI−P−1Qcek‖2.\|L(Ae^k+Bu_i^k+d^k)\|_2≤\|L I-P^-1Q_c\,e^k\|_2. (18) This reformulation allows the stability-guaranteed controller to be implemented as a Quadratic Constrained Quadratic Program (QCQP), which is significantly faster to solve than general SDPs. IV-C QCQP Formulation Integrating the T-horizon condensed parameters TH_T and Tg_T from Theorem 2, the optimal control input uiku_i^k is obtained by solving the following convex QCQP: minuik _u_i^k 12uik⊤Tuik+T⊤uik 12u_i^k H_Tu_i^k+g_T u_i^k (19a) s.t. ‖L(Aek+Buik+dk)‖2≤ηk, \|L(Ae^k+Bu_i^k+d^k)\|_2≤η^k, (19b) umin≤uik≤umax, u_min≤ u_i^k≤ u_max, (19c) where ηk:=‖LI−P−1Qcek‖2η^k:=\|L I-P^-1Q_c\,e^k\|_2 is a constant scalar at each step k. The SOC constraint (19b) ensures efficient numerical resolution while maintaining the recursive stability proven via the LMI framework. IV-D Closed-Loop Stability Analysis Theorem 3 (Input-to-State Stability of Condensed D2OC) Assume the QCQP (19) is feasible at each step k. Then, the tracking error eke^k is Input-to-State Stable (ISS) with respect to the reference drift dkd^k. Proof: To prove ISS, we establish a dissipation inequality using the quadratic Lyapunov candidate V(ek)=(ek)⊤PekV(e^k)=(e^k) Pe^k. From the quadratic form properties, it holds that λmin(P)‖ek‖2≤V(ek)≤λmax(P)‖ek‖2 _ (P)\|e^k\|^2≤ V(e^k)≤ _ (P)\|e^k\|^2. The feasibility of the contractive constraint (19b) at each step k ensures: V(ek+1)≤V(ek)−λmin(Qc)‖ek‖2.V(e^k+1)≤ V(e^k)- _ (Q_c)\|e^k\|^2. (20) By defining the decay ratio λ=λmin(Qc)/λmax(P)∈(0,1)λ= _ (Q_c)/ _ (P)∈(0,1), (20) can be rewritten as: V(ek+1)≤(1−λ)V(ek).V(e^k+1)≤(1-λ)V(e^k). (21) However, the actual error dynamics (16) involve the reference drift dkd^k. In the presence of dkd^k, the constraint (19b) remains feasible only if the control uiku_i^k can counteract the drift. Iterating (21) from j=0j=0 to k−1k-1 yields: V(ek)≤(1−λ)kV(e0)+∑j=0k−1(1−λ)k−j−1δ(dj),V(e^k)≤(1-λ)^kV(e^0)+ _j=0^k-1(1-λ)^k-j-1δ(d^j), (22) where δ(⋅)δ(·) is a class K function representing the perturbation induced by the reference shift. Using the property of geometric series ∑j=0k−1(1−λ)k−j−1<1/λ _j=0^k-1(1-λ)^k-j-1<1/λ, and taking the square root of the Lyapunov bounds, we obtain: ‖ek‖≤λmax(P)λmin(P)(1−λ)k/2‖e0‖+1λmin(P)⋅λsup0≤j<k‖dj‖. split\|e^k\|≤& _ (P) _ (P)(1-λ)^k/2\|e^0\|+ 1 _ (P)·λ _0≤ j<k\|d^j\|. split (23) The first term corresponds to the class ℒKL function β(‖e0‖,k)β(\|e^0\|,k) which decays exponentially, and the second term corresponds to the class K function γ(sup‖dj‖)γ( \|d^j\|), which defines the ultimate tracking bound. Thus, the system is ISS with respect to the reference propagation drift. ∎ Remark 2 (Predictive Alignment and Feasibility) The feasibility of (19) depends on whether the controller can manage the drift ‖dk‖\|d^k\| within its physical limits. By incorporating future barycenters, the linear term Tg_T allows the agent to preemptively align with future targets before the tracking error becomes critically large. This proactive movement keeps the state well within the safety boundary of the SOC constraint (19b), ensuring recursive feasibility even during rapid density shifts. IV-E Numerical Resolution via Soft-Constrained QCQP To ensure recursive feasibility under strict control saturation constraints uik∈i:=u∣umin≤u≤umaxu_i^k _i:=\u u_ ≤ u≤ u_ \, we reformulate the original QCQP (19) by introducing a slack variable ϵ≥0ε≥ 0. This relaxation transforms the hard contractive constraint (19b) into a soft constraint, guaranteeing that a feasible solution exists even when the required energy decay rate is physically unattainable due to the agent’s limited control authority. The relaxed optimization problem is formulated as follows: minuik,ϵ _u_i^k,ε 12(uik)⊤Tuik+T⊤uik+ρϵ2 12(u_i^k) H_Tu_i^k+g_T u_i^k+ρε^2 (24a) s.t. ‖L(Aek+Buik+dk)‖2≤ℛ(ek)+ϵ, \|L(Ae^k+Bu_i^k+d^k)\|_2 (e^k)+ε, (24b) uik∈,ϵ≥0, u_i^k , ε≥ 0, (24c) where ℛ(ek):=‖LI−P−1Qcek‖2R(e^k):=\|L I-P^-1Q_c\,e^k\|_2 denotes the state-dependent stability radius and ρ≫0ρ 0 is a penalty parameter. In this framework, the agent prioritizes the minimization of the tracking error while satisfying the stability-induced contraction to the greatest extent possible under saturation. IV-F Dual-Newton Solver for Real-Time Execution For resource-constrained on-board processing, solving (24) via generic interior-point methods (IPMs) can be computationally expensive. Exploiting the fact that the control input dimension m is typically much smaller than the state dimension, we derive a specialized Dual-Newton solver that exploits the analytical structure of the condensed QP. By defining the stability residual as f(uik):=‖L(Aek+Buik+dk)‖2−ℛ(ek)f(u_i^k):=\|L(Ae^k+Bu_i^k+d^k)\|_2-R(e^k), the optimal control is obtained by iteratively updating the Lagrange multiplier μ≥0μ≥ 0 associated with (24b). Setting the gradient of the Lagrangian to zero (∇uℒ=0 _uL=0) yields the following control law parameterized by μ: u∗(μ)=proj(−(T+μ⊤P)−1⋅(T+μ⊤P(Aek+dk))). splitu^*(μ)=proj_U (&-(H_T+ PB)^-1\\ &·(g_T+ P(Ae^k+d^k)) ). split (25) The optimal μ∗μ^* is found using the The optimal multiplier is found via the Newton-Raphson update: μν+1=max(0,μν−f(u(μν))∇μf(u(μν))),μ^ν+1= (0,μ^ν- f(u(μ^ν)) _μf(u(μ^ν)) ), (26) where ν is the inner iteration index and ∇μf _μf denotes the sensitivity of the constraint residual with respect to the dual variable. This custom solver achieves sub-millisecond convergence, providing the necessary scalability for dense coverage tasks while maintaining the practical ISS property. (a) Full KKT (b) Reduced KKT (c) Reduced KKT w/ Stability Figure 1: Multi-agent coverage trajectories via D2OC: (a) Full KKT baseline, (b) proposed reduced KKT formulation showing identical optimality, and (c) reduced KKT with recursive stability constraints. V Simulation The performance and robustness of the proposed O(T)O(T) D2OC framework were evaluated against a conventional full KKT-based solver. We considered a multi-agent task with N=10N=10 agents in a 100×100m2100× 100\,m^2 domain, where each agent follows linearized 8th-order quadrotor dynamics on a 2D plane. The communication range for each agent is set to 15m15\,m, denoted by the semi-transparent circles around the agents in Fig. 1. The target field is modeled as a non-convex Gaussian Mixture Model (GMM). To ensure consistent comparison, all simulations for horizons T∈10,…,60T∈\10,…,60\ were executed until reaching 99% global coverage. V-A Trajectory Equivalency and Stability Analysis Multi-agent trajectories are presented in Fig. 1, where full KKT solver and the proposed reduced KKT formulation results are shown in Figs. 1(a) and (b), respectively. The paths are numerically identical, which confirms that the analytical reduction via Theorem 2 preserves the exact KKT optimality conditions of the original problem. Furthermore, Fig. 1(c) demonstrates the results when augmented with stability constraints. Although paths deviate slightly to remain within the Lyapunov-based feasible region, agents successfully converge to the target distribution. This highlights the ability to ensure closed-loop stability without sacrificing the coverage mission. V-B Computational Scalability and Efficiency The computation time is depicted in Fig. 2, where the horizon T is varied from 10 to 60 in increments of 10. While the full KKT solver’s execution time escalates to 26.46 ms at T=60T=60, the proposed method maintains a consistent processing time of 2.70 ms. This 9.78-fold speedup validates the O(T)O(T) complexity of our formulation, providing a substantial computational margin for high-frequency control loops where conventional solvers fail at large horizons. Beyond mean execution time, the proposed method ensures timing predictability. As shown by the shaded regions in Fig. 2, the full KKT solver suffers from numerical jitter with a standard deviation (σ) up to 3.77 ms. In contrast, the proposed method maintains extreme consistency with σ=0.25σ=0.25 ms at T=60T=60. This jitter suppression is a direct consequence of the reduced dimensionality and well-conditioned optimization landscape, which is paramount for safety-critical systems where timing-induced delays must be minimized. Figure 2: Computation time comparison b/w the full KKT and the reduced KKT with stability guarantee. VI Conclusion In this paper, we presented an analytical structural reduction of the T-horizon D2OC problem into a computationally efficient m-dimensional QP. By leveraging the specific block-triangular structure of the KKT matrix, we derived a closed-form solution that incorporates multi-stage target barycenters without the cubic growth in computational complexity. Furthermore, our stability analysis proves that the resulting condensed MPC framework ensures ISS of the tracking error, even under time-varying reference drifts. The synergy between the linear term Tg_T and the contractive Lyapunov constraints ensures that the agent preemptively aligns with future density evolution while maintaining recursive feasibility. Our results demonstrate that the proposed framework achieves both the numerical stability of high-dimensional optimality conditions and the real-time scalability required for complex multi-agent coverage tasks. References [1] J. Cortes, S. Martinez, T. Karatas, and F. Bullo (2004) Coverage control for mobile sensing networks. IEEE Transactions on Robotics and Automation 20 (2), p. 243–255. Cited by: §I. [2] S. Du, X. Sun, M. Cao, and W. Wang (2017) Pursuing an evader through cooperative relaying in multi-agent surveillance networks. Automatica 83, p. 155–161. Cited by: §I. [3] G. Freudenthaler and T. Meurer (2020) PDE-based multi-agent formation control using flatness and backstepping: analysis, design and robot experiments. Automatica 115, p. 108897. Cited by: §I. [4] S. Ivić, B. Crnković, and I. Mezić (2016) Ergodicity-based cooperative multiagent area coverage via a potential field. IEEE transactions on cybernetics 47 (8), p. 1983–1993. Cited by: §I. [5] V. Krishnan and S. Martínez (2025) Distributed online optimization for multi-agent optimal transport. Automatica 171, p. 111880. Cited by: §I. [6] K. Lee and E. Brook (2025) Connectivity-preserving multi-agent area coverage via optimal-transport-based density-driven optimal control (d2oc). IEEE Control Systems Letters. Cited by: §I. [7] K. Lee (2026) Optimal transport-based decentralized multi-agent distribution matching. IEEE Transactions on Automatic Control (), p. 1–8. External Links: Document Cited by: §I. [8] G. Mathew and I. Mezić (2011) Metrics for ergodicity and design of ergodic dynamics for multi-agent systems. Physica D: Nonlinear Phenomena 240 (4-5), p. 432–442. Cited by: §I. [9] G. Notomista, C. Pacchierotti, and P. R. Giordano (2022) Multi-robot persistent environmental monitoring based on constraint-driven execution of learned robot tasks. In 2022 international conference on robotics and automation (icra), p. 6853–6859. Cited by: §I. [10] M. Pavone, E. Frazzoli, and F. Bullo (2010) Adaptive and distributed algorithms for vehicle routing in a stochastic and dynamic environment. IEEE Transactions on automatic control 56 (6), p. 1259–1274. Cited by: §I. [11] 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, §I-B, §I-A, Theorem 1. [12] S. Seo and K. Lee (2026) Density-driven multidrone 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, §I-B. [13] C. Villani et al. (2009) Optimal transport: old and new. Vol. 338, Springer. Cited by: §I-B. [14] T. Zheng, Q. Han, and H. Lin (2021) Transporting robotic swarms via mean-field feedback control. IEEE Transactions on Automatic Control 67 (8), p. 4170–4177. Cited by: §I.