Paper deep dive
Nonlinear Density-Driven Optimal Control (D2OC) for Multi-Agent Spatial Coverage via Sequential Convex Programming
Julian Martinez, Kooktae Lee
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 10/6/2026, 9:47:25 AM
Summary
This paper presents a nonlinear extension of Density-Driven Optimal Control (D2OC) for multi-agent spatial coverage. It utilizes Sequential Convex Programming (SCP) to handle discrete-time control-affine nonlinear systems by locally linearizing dynamics over a finite horizon. This approach transforms the problem into a strictly convex Quadratic Program (QP) that preserves the Wasserstein-based distribution-matching objective while incorporating input constraints. The method provides an explicit finite-horizon error bound and demonstrates computational efficiency compared to nonlinear Model Predictive Control (NMPC) in simulations with unicycle and quadrotor agents.
Entities (10)
Relation Signals (6)
Julian Martinez → affiliatedwith → Texas Tech University
confidence 99% · Julian Martinez and Kooktae Lee are with the Department of Mechanical and Aerospace Engineering, Texas Tech University
Kooktae Lee → affiliatedwith → Texas Tech University
confidence 99% · Julian Martinez and Kooktae Lee are with the Department of Mechanical and Aerospace Engineering, Texas Tech University
D2OC → uses → Wasserstein Distance
confidence 96% · D2OC drives the collective spatial distribution of agents toward a desired density through a Wasserstein-based objective.
D2OC → extendedby → SCP
confidence 95% · We extend this framework to multi-step finite-horizon control for discrete-time control-affine nonlinear systems using sequential convex programming.
SCP → yields → QP
confidence 95% · At each control update, the nonlinear dynamics are locally linearized over the prediction horizon, yielding a strictly convex quadratic program
D2OC → comparedto → NMPC
confidence 90% · Simulations with unicycle and quadrotor teams show coverage performance comparable to nonlinear model predictive control, while substantially reducing computation time.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper presents a nonlinear extension of Density-Driven Optimal Control (D2OC) for multi-agent spatial coverage with prescribed density distributions. Rather than assigning individual target locations, D2OC drives the collective spatial distribution of agents toward a desired density through a Wasserstein-based objective. We extend this framework to multi-step finite-horizon control for discrete-time control-affine nonlinear systems using sequential convex programming. At each control update, the nonlinear dynamics are locally linearized over the prediction horizon, yielding a strictly convex quadratic program that preserves the Wasserstein barycentric structure while directly incorporating input constraints. We further characterize the effect of constrained control deviations and nonlinear Taylor remainders on the accuracy of the local linear prediction, establishing an explicit finite-horizon error bound and a two-step specialization for receding-horizon implementation. The resulting method retains the decentralized, distribution-driven nature of D2OC while providing a computationally efficient optimization procedure for nonlinear multi-agent systems. Simulations with unicycle and quadrotor teams show coverage performance comparable to nonlinear model predictive control, while substantially reducing computation time. These results demonstrate a tractable and theoretically characterized framework for density-driven spatial coverage under nonlinear dynamics.
Tags
Links
- Source: https://arxiv.org/abs/2610.04545v1
- Canonical: https://arxiv.org/abs/2610.04545v1
Trouble viewing inline? Open PDF directly →
Full Text
49,101 characters extracted from source content.
Expand or collapse full text
Nonlinear Density-Driven Optimal Control (D2OC) for Multi-Agent Spatial Coverage via Sequential Convex Programming Julian Martinez Affiliation: Julian Martinez and Kooktae Lee are with the Department of Mechanical and Aerospace Engineering, Texas Tech University, Lubbock, TX 79409, USA, email: mar85781@ttu.edu, kooktae.lee@ttu.edu. Kooktae Lee †thanks: *This work was supported by NSF CAREER Grant CMMI-DCSD-2638508. Affiliation: Julian Martinez and Kooktae Lee are with the Department of Mechanical and Aerospace Engineering, Texas Tech University, Lubbock, TX 79409, USA, email: mar85781@ttu.edu, kooktae.lee@ttu.edu. Abstract This paper presents a nonlinear extension of Density-Driven Optimal Control (D2OC) for multi-agent spatial coverage with prescribed density distributions. Rather than assigning individual target locations, D2OC drives the collective spatial distribution of agents toward a desired density through a Wasserstein-based objective. We extend this framework to multi-step finite-horizon control for discrete-time control-affine nonlinear systems using sequential convex programming. At each control update, the nonlinear dynamics are locally linearized over the prediction horizon, yielding a strictly convex quadratic program that preserves the Wasserstein barycentric structure while directly incorporating input constraints. We further characterize the effect of constrained control deviations and nonlinear Taylor remainders on the accuracy of the local linear prediction, establishing an explicit finite-horizon error bound and a two-step specialization for receding-horizon implementation. The resulting method retains the decentralized, distribution-driven nature of D2OC while providing a computationally efficient optimization procedure for nonlinear multi-agent systems. Simulations with unicycle and quadrotor teams show coverage performance comparable to nonlinear model predictive control, while substantially reducing computation time. These results demonstrate a tractable and theoretically characterized framework for density-driven spatial coverage under nonlinear dynamics. I Introduction Multi-agent spatial coverage aims to distribute autonomous agents over a spatial domain according to sensing, monitoring, or mission-dependent requirements. Classical coverage control methods, such as the Voronoi-based framework, partition the workspace into agent-specific Voronoi cells and steer agents toward their cell centroids [1]. While effective for geometric coverage, such methods do not directly address non-uniform agent distributions. This motivates formulating non-uniform coverage as a spatial distribution-matching problem. Optimal transport (OT) provides a principled framework for comparing and matching spatial distributions through Wasserstein distances [2, 3]. This perspective has also been incorporated into multi-agent coverage control, where the objective is formulated directly in terms of the discrepancy between the empirical distribution of agents and a prescribed spatial density, rather than through individually assigned target locations [4]. Building on this distribution-based viewpoint, Density-Driven Optimal Control (D2OC) formulates non-uniform multi-agent coverage as a decentralized optimal-control problem that drives the empirical agent distribution toward a reference density through local sample selection, weight updates, and neighbor information exchange [5, 6, 7, 8]. The D2OC framework was initially developed for linear time-invariant (LTI) dynamics and later extended to linear time-varying (LTV) dynamics [6, 7]. Its extension to nonlinear systems, however, has so far been limited to a one-step prediction horizon. A multi-step nonlinear formulation must account for nonlinear agent dynamics while preserving the Wasserstein-based coverage objective, decentralized structure, and input constraints. NMPC provides a framework for nonlinear dynamics and constraints through repeated finite-horizon optimization [9], while distributed receding-horizon methods address computational and communication requirements in multi-agent systems [10]. Sequential convex programming (SCP), in contrast, locally approximates nonlinear dynamics and solves a sequence of convex subproblems, and has been applied to real-time nonlinear optimal control and successive convexification [11, 12, 13]. This makes SCP a natural candidate for extending D2OC to multi-step nonlinear coverage although the resulting local approximation introduces a discrepancy between predicted and actual trajectories that must be characterized. In this paper, we extend D2OC to discrete-time control-affine nonlinear systems using sequentially linearized local models. Each optimization step yields a strictly convex local Quadratic Program (QP) with input constraints while retaining the Wasserstein-based coverage objective. We further characterize constraint-induced control deviations and derive a finite-horizon bound on the resulting linearization error. Thus, the main contributions of this paper are summarized as follows: 1) We extend the Wasserstein-based D2OC formulation to discrete-time control-affine nonlinear systems while preserving its decentralized coverage mechanism and distribution-matching objective; 2) We formulate each nonlinear D2OC update as a strictly convex local QP via sequential convex programming, retaining the Wasserstein barycenter structure and actuator constraints while enabling efficient computation; 3) We quantify input-saturation effects through the deviation between constrained and unconstrained QP solutions and incorporate this deviation into the linearization-error analysis; 4) We establish a finite-horizon error bound accounting for control deviations and second-order Taylor remainders, with a receding-horizon specialization for first-input application. The resulting framework combines the distribution-level objective of D2OC with efficient local optimization and explicit prediction-error characterization for nonlinear multi-agent dynamics. I Problem Setup Notation: The sets of real numbers are denoted by ℝR. The sets ℤ>0Z_>0 and ℤ≥0:=ℤ>0∪0Z_≥ 0:=Z_>0∪\0\ denote positive and non-negative integers. The space ℝnR^n represents n-dimensional column vectors. The Euclidean norm of a vector x is written as ‖\|x\|, and the transpose of a matrix A as A⊤A . The zero matrix m×n∈ℝm×n0_m× n ^m× n, the identity matrix n∈ℝn×nI_n ^n× n, and the all-ones vector n∈ℝn1_n ^n are denoted with their sizes as subscripts. The operator diag([⋅])diag([·]) constructs a block diagonal matrix from its arguments, while col(⋅)col(·) vertically concatenates its vector or matrix arguments. The Kronecker and Hadamard (elementwise) products are denoted by ⊗ and ⊙ , respectively. We consider a network of nan_a agents, where each agent is modeled by a discrete-time control-affine nonlinear system xi,k+1 x_i,k+1 =fi(xi,k)+Gi(xi,k)ui,k, =f_i(x_i,k)+G_i(x_i,k)u_i,k, (1) yi,k y_i,k =Cixi,k, =C_ix_i,k, where xi,k∈ℝnx_i,k ^n, ui,k∈ℝmu_i,k ^m, and yi,k∈ℝdy_i,k ^d denote the state, control input, and output, respectively. The symbols fi:ℝn→ℝnf_i:R^n ^n and Gi:ℝn→ℝn×mG_i:R^n ^n× m describe the nonlinear drift and input dynamics, respectively, and Ci∈ℝd×nC_i ^d× n maps the state to the coverage output. We assume that fif_i and GiG_i are sufficiently smooth over the region of interest. The objective is to control the agents so that their empirical spatial distribution approaches a prescribed reference distribution while respecting input and communication constraints. We adopt the Density-Driven Optimal Control (D2OC) framework [5, 7, 6, 8] for this purpose. The original D2OC formulation is based on a Wasserstein-distance objective and admits a QP formulation for linear time-invariant dynamics. In this work, this structure is used as the starting point for extending D2OC to control-affine nonlinear systems. In particular, we retain the Wasserstein-based coverage objective while replacing the globally linear prediction model with a sequentially linearized local model. For a discrete probability measure, let δy _y denote the Dirac measure at y. We characterize the desired spatial distribution ν by a finite set of sample points qjj=1N⊂ℝd\q_j\_j=1^N ^d with associated weights βj _j, and the agent distribution by the empirical measure ρk _k at time k as follows: ρk:=1k+1∑t=0k(1na∑i=1naδyi,t),ν:=∑j=1Nβjδqj. _k:= 1k+1 _t=0^k ( 1n_a _i=1^n_a _y_i,t ), ν:= _j=1^N _j _q_j. (2) The coverage objective is to reduce the discrepancy between ρk _k and ν, measured by the quadratic Wasserstein distance 22(μ,ν)=minπℓj∑ℓ=1M∑j=1Nπℓj∥yℓ−qj∥2,W_2^2(μ,ν)= _ _ j _ =1^M _j=1^N _ j\|y_ -q_j\|^2, (3) subject to the standard transport constraints πℓj≥0,∑jπℓj=αℓ,∑ℓπℓj=βj. _ j≥ 0, _j _ j= _ , _ _ j= _j. Rather than solving the global Wasserstein problem centrally, the Density-Driven Optimal Control (D2OC) framework operates through the following decentralized cycle [5, 7, 8]: 1. Local sample selection and control: each agent selects nearby sample points with nonzero remaining weights and computes a feasible control input to reduce its local Wasserstein cost. 2. Weight update: after executing the control action, the agent decreases the remaining weights of the sample points it covers, recording its coverage progress. 3. Neighbor information exchange: when agents communicate, they synchronize their sample weights using the information received from neighboring agents, avoiding redundant coverage. Thus, the first two stages are decentralized, while communication is used only to maintain consistency of the coverage information. In this work, we focus on the optimal-control step in Stage 1 and its extension to nonlinear agent dynamics. I Formulation of the D2OC Cost Function via Quadratic Programming We first recall how the Wasserstein objective is formulated in D2OC. Unlike the standard optimal transport problem, where the transport plan πℓj _ j is the optimization variable, D2OC uses the current agent position to analytically determine the transport plan and identify where the agent should move to match the prescribed density. Over a prediction horizon, these weights evolve with the predicted agent position. For agent i, let i,k+hS_i,k+h denote the selected reference samples at prediction step k+hk+h, with corresponding transport weights πj,k+h _j,k+h. The local D2OC cost is ∑h=1Hi2(k+h):=∑h=1H∑j∈i,k+hπj,k+h‖yi,k+h−qj‖2. _h=1^HW_i^2(k+h):= _h=1^H _j i,k+h _j,k+h||y_i,k+h-q_j||^2. (4) Defining the weighted barycenter at each prediction step as q¯i,k+h:=∑j∈i,k+hπj,k+hqj∑j∈i,k+hπj,k+h, q_i,k+h:= _j _i,k+h _j,k+hq_j _j _i,k+h _j,k+h, (5) the Wasserstein cost becomes a weighted least-squares tracking problem with respect to the time-varying barycenter. Stacking the predicted outputs and barycenters gives Yik|H Y_i^k|H :=[yi,k+1⊤⋯yi,k+H⊤]⊤, := bmatrixy_i,k+1 &·s&y_i,k+H bmatrix , (6) Ξik|H _i^k|H :=[q¯i,k+1⊤⋯q¯i,k+H⊤]⊤. := bmatrix q_i,k+1 &·s& q_i,k+H bmatrix . Using ik|H:=diag(∑j∈i,k+hπj,k+h,d)h=1H, _i^k|H:=diag ( _j _i,k+h _j,k+h,I_d )_h=1^H, the local Wasserstein objective can then be written as a weighted least-squares tracking cost as follows. Proposition 1 (Wasserstein-to-barycenter reduction). (Proposition 1 in [6]) Let the transport weights πj,k+h _j,k+h be fixed over the prediction horizon. Then, ∑h=1Hi2(k+h)=‖ik|H(Yik|H−Ξik|H)‖2+const. _h=1^HW_i^2(k+h)= \| _i^k|H (Y_i^k|H- _i^k|H ) \|^2+const. (7) For the LTI dynamics, the stacked output prediction is affine in the stacked control input Yik|H=ΘiUik|H+Φixi,k,Y_i^k|H= _iU_i^k|H+ _ix_i,k, (8) where Uik|H=[ui,k⊤⋯ui,k+H−1⊤]⊤,U_i^k|H= bmatrixu_i,k &·s&u_i,k+H-1 bmatrix , (9) Θi=[CiBi0⋯0CiAiBiCiBi⋯0⋱CiAiH−1BiCiAiH−2Bi⋯CiBi], _i= bmatrixC_iB_i&0&·s&0\\ C_iA_iB_i&C_iB_i&·s&0\\ & & & \\ C_iA_i^H-1B_i&C_iA_i^H-2B_i&·s&C_iB_i bmatrix, (10) Φi=[(CiAi)⊤(CiAi2)⊤⋯(CiAiH)⊤]⊤. _i= bmatrix(C_iA_i) &(C_iA_i^2) &·s&(C_iA_i^H) bmatrix . (11) Consequently, Proposition 1 establishes that the complex distribution matching problem can be effectively recast into a standard weighted barycentric tracking task together with the weight update stage at each time step introduced in Section I. IV Nonlinear Optimality Conditions and SCP-Based Local QP Formulation Consider the discrete-time control-affine nonlinear dynamics in (1). For clarity, we omit the agent index i throughout this section. Let k denote the current decision time and t∈k,…,k+H−1t∈k,…,k+H-1 the prediction-time index. The dynamics at each prediction step are given by xt+1=f(xt)+G(xt)ut,yt=Cxt,t=k,…,k+H−1.x_t+1=f(x_t)+G(x_t)u_t,\;\;y_t=Cx_t, t=k,…,k+H-1. Let the stage cost at time t be ℓt(xt,ut):=(yt−q¯t)⊤Qt(yt−q¯t)+ut⊤Rut, _t(x_t,u_t):=(y_t- q_t) Q_t(y_t- q_t)+u_t R\,u_t, (12) where Qt∈ℝd×dQ_t ^d× d is symmetric and positive semidefinite and R≻0R 0. Over a prediction horizon of length H, let the terminal cost be ℓk+H(xk+H)=(yk+H−q¯k+H)⊤Qk+H(yk+H−q¯k+H), _k+H(x_k+H)=(y_k+H- q_k+H) Q_k+H(y_k+H- q_k+H), if needed. IV-A PMP Characterization and NMPC Baseline The finite-horizon nonlinear optimal control problem can be characterized by discrete-time Pontryagin’s Minimum Principle (PMP). For the control-affine dynamics (1), any locally optimal solution satisfies the first-order necessary conditions: xt+1 x_t+1 =f(xt)+G(xt)ut, =f(x_t)+G(x_t)u_t, (13a) λt _t =2C⊤Qt(Cxt−q¯t)+∂Ψ∂x(xt,ut)⊤λt+1, =2C Q_t(Cx_t- q_t)+ ∂ ∂ x(x_t,u_t) _t+1, (13b) λk+H _k+H =2C⊤Qk+H(Cxk+H−q¯k+H), =2C Q_k+H(Cx_k+H- q_k+H), (13c) where Ψ(x,u):=f(x)+G(x)u. (x,u):=f(x)+G(x)u. (14) Under box control constraints umin≤ut≤umaxu_min≤ u_t≤ u_max, the control satisfies the pointwise Hamiltonian minimization condition ut⋆∈argminumin≤u≤umaxu⊤Ru+λt+1⊤G(xt)u.u_t ∈ _u_min≤ u≤ u_max \u Ru+ _t+1 G(x_t)u \. (15) Equivalently, this condition follows from the KKT conditions of the box-constrained quadratic problem. When R is diagonal, its solution is available in closed form as ut⋆=sat[umin,umax](−12R−1G(xt)⊤λt+1).u_t =sat_[u_min,u_max] (- 12R^-1G(x_t) _t+1 ). (16) Thus, the saturation in (16) is the exact pointwise minimizer rather than a heuristic clipping operation. These conditions formulate a two-point boundary-value problem (TPBVP) with xkx_k given and λk+H _k+H specified at the terminal horizon. A standard shooting approach treats the initial costate λk=s _k=s as a decision variable, propagating the state-costate dynamics to find a root of the terminal residual F(s):=λk+H(s)−2C⊤Qk+H(Cxk+H(s)−q¯k+H)=0.F(s):= _k+H(s)-2C Q_k+H (Cx_k+H(s)- q_k+H )=0. (17) While this PMP-based shooting method serves as our baseline NMPC reference, the saturation operator in (16) introduces active-set switching and hence makes F(s)F(s) generally only piecewise smooth, which can complicate Newton-type root-finding and Jacobian evaluation. To overcome these limitations, our primary computational strategy relies on Sequential Convex Programming (SCP). Rather than solving the resulting piecewise-smooth shooting problem, we locally linearize the dynamics around a nominal trajectory and solve a sequence of convex QPs, with the input bounds imposed directly as QP constraints. This provides a convenient numerical formulation for handling input constraints while preserving the Wasserstein tracking objective. The resulting local QP formulation is presented next. IV-B SCP-Based Local QP Reduction The preceding PMP formulation provides the necessary optimality conditions for the original nonlinear problem and motivates a shooting-based NMPC procedure as a reference method. The proposed approach instead avoids solving the nonlinear TPBVP directly. We employ a receding-horizon SCP scheme in which the nonlinear dynamics are linearized once at each decision time. After the first optimized input is applied, the horizon is shifted, the dynamics are relinearized about the new operating point, and a new convex QP is solved. This provides a convenient numerical formulation for handling input constraints while preserving the Wasserstein tracking objective. At decision time k, let (xk,uknom)(x_k,u_k^nom) denote the nominal state-input pair used for linearization. Given the definition of Ψ in (14), define Ak A_k :=∂Ψ∂x(xk,uknom)=∂f∂x(xk)+∑j=1m∂Gj∂x(xk)uk,jnom, := ∂ ∂ x(x_k,u_k^nom)= ∂ f∂ x(x_k)+ _j=1^m ∂ G_j∂ x(x_k)u_k,j^nom, (18) Bk B_k :=∂Ψ∂u(xk,uknom)=G(xk). := ∂ ∂ u(x_k,u_k^nom)=G(x_k). A first-order Taylor expansion of Ψ about (xk,uknom)(x_k,u_k^nom) gives xt+1≈Ψ(xk,uknom)+Ak(xt−xk)+Bk(ut−uknom).x_t+1≈ (x_k,u_k^nom)+A_k(x_t-x_k)+B_k(u_t-u_k^nom). (19) Defining δxt:=xt−xk,δut:=ut−uknom,δ x_t:=x_t-x_k, δ u_t:=u_t-u_k^nom, (20) and ck:=Ψ(xk,uknom)−xk,c_k:= (x_k,u_k^nom)-x_k, (21) the local affine model becomes δxt+1=ck+Akδxt+Bkδut,t=k,…,k+H−1,δ x_t+1=c_k+A_kδ x_t+B_kδ u_t, t=k,…,k+H-1, (22) with δxk=0δ x_k=0. In general, ck≠0c_k≠ 0 because the linearization point need not be an equilibrium of the nominal dynamics. Note that the linearization is performed once at the current decision point (xk,uknom)(x_k,u_k^nom). The resulting AkA_k, BkB_k, and ckc_k are then fixed throughout the H-step prediction horizon. This fixed-model construction provides a single consistent linearized prediction model for the local QP and enables explicit analysis of the resulting closed-loop approximation error. To separate the input-independent part of the affine prediction from the input-dependent part, let δxt0δ x_t^0 denote the prediction of (22) with δut=0δ u_t=0. Thus, δxt+10=ck+Akδxt0,δxk0=0.δ x_t+1^0=c_k+A_kδ x_t^0, δ x_k^0=0. (23) Define the input-induced state deviation by δx^t:=δxt−δxt0.δ x_t:=δ x_t-δ x_t^0. (24) Subtracting (23) from (22) gives δx^t+1=Akδx^t+Bkδut,δx^k=0.δ x_t+1=A_kδ x_t+B_kδ u_t, δ x_k=0. (25) Since yt=Cxty_t=Cx_t, the corresponding input-induced output deviation is δy^t=Cδx^t.δ y_t=Cδ x_t. Stacking these output and input deviations over the horizon, define ΔY:=col(δy^k+1,…,δy^k+H), Y:=col (δ y_k+1,…,δ y_k+H ), (26) ΔU:=col(δuk,…,δuk+H−1). U:=col (δ u_k,…,δ u_k+H-1 ). Because AkA_k and BkB_k are fixed over the horizon, the standard finite-horizon input-output relation gives ΔY=ΘΔU, Y= U, (27) where Θ is given by (10). Further, defining Y=col(yk+1,…,yk+H),Y0=col(yk+10,…,yk+H0),Y=col (y_k+1,…,y_k+H ),\;Y^0=col (y_k+1^0,…,y_k+H^0 ), (28) with yt0=C(xk+δxt0),y_t^0=C (x_k+δ x_t^0 ), we obtain the affine prediction Y=Y0+ΘΔU.Y=Y^0+ U. (29) Using Ξ for the stacked Wasserstein barycenter in (6), the local tracking objective is constructed by J(ΔU)=‖(Y0+ΘΔU−Ξ)‖2+ΔU⊤ℛΔU,J( U)= \| (Y^0+ U- ) \|^2+ U R U, (30) where ℛ=IH⊗RR=I_H R, R≻0.R 0. The above cost J corresponds to the least square cost in (7) combined with the input penalty term. Expanding (30) with respect to ΔU U gives J(ΔU)=12ΔU⊤ℋΔU+g⊤ΔU+const.,J( U)= 12 U H U+g U+const., (31) where the constant term is independent of ΔU U, and ℋ =2[(Θ)⊤(Θ)+ℛ], =2 [( ) ( )+R ], (32) g g =2(Θ)⊤(Y0−Ξ). =2( ) (Y^0- ). Since ℛ≻0R 0, it follows that ℋ≻0H 0, establishing strict convexity of the local QP. Theorem 1 (Optimal Solution for Nonlinear D2D^2OC via Strictly Convex Local QP). At each SCP iteration, the local control update is obtained from the quadratic program minΔU _ U 12ΔU⊤ℋΔU+g⊤ΔU 12 U H U+g U (33) s.t. .t. umin(H)≤Unom+ΔU≤umax(H), u_ ^(H)≤ U^nom+ U≤ u_ ^(H), where ℋH and g are defined in (32), and Unom U^nom :=H⊗uknom, :=1_H u_k^nom, umin(H) u_ ^(H) :=H⊗umin,umax(H):=H⊗umax. :=1_H u_ , u_ ^(H):=1_H u_ . Thus, the box constraints are imposed at every prediction stage t=k,…,k+H−1t=k,…,k+H-1. The local QP is strictly convex and therefore admits a unique global minimizer whenever its constraints are feasible. If the unconstrained minimizer is feasible, it is given by ΔUun=−ℋ−1g. U^un=-H^-1g. (34) Otherwise, the unique constrained minimizer is characterized by the KKT conditions. Proof. Let M:=ΩΘM:= . From (32), ℋ=2M⊤M+2ℛ,ℛ=IH⊗R.H=2M M+2R, =I_H R. Since M⊤M⪰0M M 0 and R≻0R 0, we have ℛ≻0R 0, and hence ℋ≻0.H 0. Therefore, the QP objective is strictly convex while the box constraints are affine. Consequently, whenever the feasible set is nonempty, the QP admits a unique global minimizer. If the unconstrained minimizer satisfies the box constraints, the first-order condition ℋΔUun+g=0H U^un+g=0 directly gives (34). For the constrained case, define the Lagrangian ℒ(ΔU,λ+,λ−)= ( U,λ^+,λ^-)= 12ΔU⊤ℋΔU+g⊤ΔU 12 U H U+g U +(λ+)⊤(Unom+ΔU−umax(H)) +(λ^+) (U^nom+ U-u_ ^(H) ) +(λ−)⊤(umin(H)−Unom−ΔU). +(λ^-) (u_ ^(H)-U^nom- U ). (35) The KKT conditions are ℋΔU⋆+g+λ+−λ− U +g+λ^+-λ^- =0, =0, (36) λ+,λ− λ^+,λ^- ≥0, ≥ 0, (37) umin(H)≤Unom+ΔU⋆ u_ ^(H)≤ U^nom+ U ≤umax(H), ≤ u_ ^(H), (38) λ+⊙(Unom+ΔU⋆−umax(H)) λ^+ (U^nom+ U -u_ ^(H) ) =0, =0, (39) λ−⊙(umin(H)−Unom−ΔU⋆) λ^- (u_ ^(H)-U^nom- U ) =0. =0. (40) Since the objective is convex and the constraints are affine, these KKT conditions are necessary and sufficient for global optimality. Strict convexity of the objective further guarantees uniqueness of ΔU⋆ U . ∎ Remark 1 (Significance of the SCP-QP Formulation). While the general iterative solution mechanism aligns with SCP, Theorem 1 provides a non-trivial reduction tailored specifically to the D2OCD^2OC framework. Rather than applying generic Sequential Quadratic Programming (SQP), this formulation explicitly maps the Wasserstein barycenter tracking metrics into the Hessian ℋH and gradient g matrices. This structured reduction guarantees strict convexity (ℋ≻0H 0) and uniqueness of ΔU⋆ U at each step, serving as the precise computational backbone for the rigorous error bound and certificate analysis developed in the sequel. V Linearization Error and Finite-Horizon Error Certification The local QP in Section IV is constructed from a first-order linearization of the nonlinear dynamics. Although this reduction yields a convex optimization problem, the resulting prediction can deviate from the nonlinear trajectory due to both input clipping and the higher-order Taylor remainder of the dynamics. This section quantifies these effects and establishes a finite-horizon bound on the resulting linearization error. The following assumptions are provided to analyze this linearization error. Assumption 1 (Regularity and curvature bounds). f,G∈C2f,G∈ C^2 on a convex set ×⊂ℝn×ℝmX×U ^n×R^m containing the relevant state-input trajectories, with ‖fxx(x)‖≤Mf,‖∂xxGj(x)‖≤MGx,‖∂xGj(x)‖≤MGu,∀(x,u)∈×.\|f_x(x)\|≤ M_f,\;\| _xG_j(x)\|≤ M_Gx,\;\| _xG_j(x)\|≤ M_Gu,\;∀(x,u) ×U. Set M¯:=Mf+mMGx M:=M_f+mM_Gx and M¯′:=mMGu M :=mM_Gu, where m denotes the input dimension. This regularity condition ensures that the second-order Taylor remainder of Ψ is uniformly bounded over X. Furthermore, the convexity of X guarantees that line segments connecting candidate trajectories remain inside X, enabling the application of the mean-value theorem. Assumption 2 (Uniform sensitivity on the operative tube). Let Ψ(x,u) (x,u) be defined as in (14). There exists a finite constant L>0L>0 such that ‖∂xΨ(x,u)‖≤L,∀(x,u)∈×.\| _x (x,u)\|≤ L, ∀(x,u) ×U. Moreover, define κ:=supx∈‖G(x)‖<∞.κ:= _x \|G(x)\|<∞. Then, ‖∂uΨ(x,u)‖=‖G(x)‖≤κ,∀(x,u)∈×.\| _u (x,u)\|=\|G(x)\|≤κ, ∀(x,u) ×U. The symbol κ limits the single-step dynamic perturbation from input clipping, while the Lipschitz bound L permits the propagation of this perturbation across the prediction horizon via a discrete Grönwall inequality [14]. Assumption 3 (Small nominal excursion). There exist constants εx>0 _x>0 and εu>0 _u>0 such that ∥δxtun∥≤εx,∥δutun∥≤εu,t=k,…,k+H−1.\|δ x_t^un\|≤ _x, \|δ u_t^un\|≤ _u, t=k,…,k+H-1. The above assumption ensures that the nominal trajectory remains within the local neighborhood in which the subsequent sensitivity bounds are valid. Lemma 1 (Uniform control deviation bound). Let ΔU⋆ U and ΔUun U^un denote the constrained and unconstrained minimizers of the local QP in Theorem 1. Define δumax:=‖ΔU⋆−ΔUun‖. _u := \| U - U^un \|. (41) Then, for every prediction stage t=k,…,k+H−1t=k,…,k+H-1, ‖ut⋆−utun‖≤δumax. \|u_t -u_t^un \|≤ _u . (42) Proof. Recall that δut⋆:=ut⋆−uknomδ u_t :=u_t -u_k^nom, δutun:=utun−uknom.δ u_t^un:=u_t^un-u_k^nom. Since ΔU⋆ U and ΔUun U^un are stacked as ΔU⋆ U =col(δuk⋆,…,δuk+H−1⋆), =col (δ u_k ,…,δ u_k+H-1 ), (43) ΔUun U^un =col(δukun,…,δuk+H−1un), =col (δ u_k^un,…,δ u_k+H-1^un ), we have ut⋆−utun=δut⋆−δutun.u_t -u_t^un=δ u_t -δ u_t^un. Because δut⋆−δutunδ u_t -δ u_t^un is a block of ΔU⋆−ΔUun U - U^un, the Euclidean norm satisfies ‖ut⋆−utun‖=‖δut⋆−δutun‖≤‖ΔU⋆−ΔUun‖=δumax. \|u_t -u_t^un \|= \|δ u_t -δ u_t^un \|≤ \| U - U^un \|= _u . ∎ The preceding development establishes the local SCP subproblem and quantifies the control deviation introduced by the input constraints. We next use this control deviation as the starting point for analyzing the accuracy of the local linearization along the resulting nonlinear trajectory. Theorem 2 (Linearization error under QP clipping). Under Assumptions 1–3, let A,BA,B denote the local linearization matrices of Theorem 1 at (xk,uknom)(x_k,u_k^nom), with ‖A‖≤α\|A\|≤α. Let ΔUun U^un and ΔU⋆ U be the unconstrained and box-constrained QP solutions, respectively, and let xtunx_t^un and xt⋆x_t denote the corresponding nonlinear trajectories initialized at the common state xkx_k. Define δxt⋆:=xt⋆−xkδ x_t :=x_t -x_k and δut⋆:=ut⋆−uknom,δ u_t :=u_t -u_k^nom, and let xtlinx_t^lin be the state predicted by the local linear model under the applied constrained input, δxt+1lin=ck+Aδxtlin+Bδut⋆,δxklin=0,δ x_t+1^lin=c_k+Aδ x_t^lin+Bδ u_t , δ x_k^lin=0, where ckc_k is described in (21). Set et:=xt⋆−xtline_t:=x_t -x_t^lin, ek=0e_k=0, and let δumax _u be as in (41). Define ε¯x:=εx+κδumax∑j=0N−k−1Lj,ε¯u:=εu+δumax. _x:= _x+κ _u _j=0^N-k-1L^j, _u:= _u+ _u . Then, for every k<N≤k+Hk<N≤ k+H, ‖eN‖≤(∑j=0N−k−1αj)[12M¯ε¯x2+M¯′ε¯xε¯u].\|e_N\|≤ ( _j=0^N-k-1α^j ) [ 12 M _x^2+ M _x _u ]. (44) Equivalently, if α≠1α≠ 1 and L≠1L≠ 1, ε¯x=εx+κδumaxLN−k−1L−1, _x= _x+κ _u L^N-k-1L-1, and ‖eN‖≤αN−k−1α−1[12M¯ε¯x2+M¯′ε¯xε¯u].\|e_N\|≤ α^N-k-1α-1 [ 12 M _x^2+ M _x _u ]. If α=1α=1 or L=1L=1, the corresponding geometric sum is replaced by N−kN-k. Proof. Define the clipping-induced trajectory deviation e^t:=xt⋆−xtun,e^k=0. e_t:=x_t -x_t^un, e_k=0. Since both trajectories remain in ×X×U by Assumption 1, the line segment connecting (xtun,utun)(x_t^un,u_t^un) and (xt⋆,ut⋆)(x_t ,u_t ) lies in ×X×U. Parameterize this segment by xs=xtun+se^t,us=utun+s(ut⋆−utun),s∈[0,1].x_s=x_t^un+s e_t, u_s=u_t^un+s(u_t -u_t^un), s∈[0,1]. By the fundamental theorem of calculus, e^t+1 e_t+1 =Ψ(xt⋆,ut⋆)−Ψ(xtun,utun) = (x_t ,u_t )- (x_t^un,u_t^un) =∫01[∂xΨ(xs,us)e^t+∂uΨ(xs,us)(ut⋆−utun)]ds. = _0^1 [ _x (x_s,u_s) e_t+ _u (x_s,u_s)(u_t -u_t^un) ]ds. Therefore, using ‖∂xΨ‖≤L\| _x \|≤ L, ‖∂uΨ‖≤κ\| _u \|≤κ, and ‖ut⋆−utun‖≤δumax\|u_t -u_t^un\|≤ _u , ‖e^t+1‖≤L‖e^t‖+κδumax.\| e_t+1\|≤ L\| e_t\|+κ _u . Since e^k=0 e_k=0, iteration yields ‖e^t‖≤κδumax∑j=0t−k−1Lj.\| e_t\|≤κ _u _j=0^t-k-1L^j. Hence, for all t≤Nt≤ N, ‖e^t‖≤κδumax∑j=0N−k−1Lj.\| e_t\|≤κ _u _j=0^N-k-1L^j. By Assumption 3, ‖δxtun‖≤εx.\|δ x_t^un\|≤ _x. Thus, ‖δxt⋆‖≤‖δxtun‖+‖e^t‖≤ε¯x.\|δ x_t \|≤\|δ x_t^un\|+\| e_t\|≤ _x. Likewise, ‖δut⋆‖≤‖δutun‖+‖ut⋆−utun‖≤ε¯u.\|δ u_t \|≤\|δ u_t^un\|+\|u_t -u_t^un\|≤ _u. We now expand the nonlinear dynamics about the fixed linearization point (xk,uknom)(x_k,u_k^nom). Since Ψ(x,u)=f(x)+G(x)u (x,u)=f(x)+G(x)u is affine in u, its second derivative with respect to u vanishes. Consequently, the second-order Taylor remainder contains only state-state and state-input terms. By the Hessian bounds in Assumption 1, ‖rt‖≤12M¯‖δxt⋆‖2+M¯′‖δxt⋆‖‖δut⋆‖≤12M¯ε¯x2+M¯′ε¯xε¯u.\|r_t\|≤ 12 M\|δ x_t \|^2+ M \|δ x_t \|\|δ u_t \|≤ 12 M _x^2+ M _x _u. Define η:=12M¯ε¯x2+M¯′ε¯xε¯u.η:= 12 M _x^2+ M _x _u. Because the linear prediction uses the same applied input δut⋆δ u_t as the nonlinear trajectory, subtraction gives the exact error recursion et+1=Aet+rte_t+1=Ae_t+r_t, ek=0.e_k=0. Therefore, eN=∑t=kN−1AN−1−trt,e_N= _t=k^N-1A^N-1-tr_t, and hence ‖eN‖ \|e_N\| ≤∑t=kN−1‖A‖N−1−t‖rt‖ ≤ _t=k^N-1\|A\|^N-1-t\|r_t\| ≤η∑t=kN−1αN−1−t=η∑j=0N−k−1αj. ≤η _t=k^N-1α^N-1-t=η _j=0^N-k-1α^j. Substituting the definition of η proves (44). ∎ Corollary 1 (Two-step specialization of the finite-horizon error bound). Under the conditions of Theorem 2, consider N=k+2N=k+2. Then the linearization error after two prediction steps satisfies ‖ek+2‖ \|e_k+2\| ≤12M¯(εx+κδumax,(0))2 ≤ 12 M ( _x+κ _u ,(0) )^2 (45) +M¯′(εx+κδumax,(0))(εu+δumax). + M ( _x+κ _u ,(0) ) ( _u+ _u ). In particular, since Ψ is affine in u and δxk⋆=0δ x_k =0, the first-step Taylor remainder satisfies rk=0r_k=0, and hence ek+1=0e_k+1=0. Thus, the first nonzero linearization error can occur at N=k+2N=k+2. Proof. At t=kt=k, δxk⋆=0δ x_k =0. Since Ψ is affine in u, the second-order Taylor remainder satisfies rk=0r_k=0, and therefore ek+1=Aek+rk=0.e_k+1=Ae_k+r_k=0. Hence, the error at k+2k+2 is generated only by the Taylor remainder at the second step. Applying the local remainder bound at t=k+1t=k+1 and using the one-step state and input deviations gives (45), where δumax,(0) _u ,(0) denotes the clipping magnitude of the first input block, with δumax,(0)≤δumax _u ,(0)≤ _u . ∎ Remark 2 (Corollary 1 sharpens Theorem 2 at N=k+2N=k+2). Corollary 1 targets the same quantity, ek+2e_k+2, that Theorem 2 bounds at N=k+2N=k+2, but exploits the specific initialization ek=0e_k=0 and the fact that rk=0r_k=0. Consequently, the Grönwall recursion is needed only once, from k to k+1k+1, rather than over the entire horizon. This removes the L-dependent factor 1+L1+L from the state envelope ε¯x _x appearing under direct substitution of N=k+2N=k+2 into (44): only the first input block contributes to e^k+1 e_k+1, which is captured by δumax,(0)≤δumax _u ,(0)≤ _u . The input envelope ε¯u _u, by contrast, still bounds δuk+1⋆δ u_k+1 – the second input block – which δumax,(0) _u ,(0) does not control; it therefore retains the horizon-wide quantity δumax _u from Lemma 1. Thus, (45) provides a bound that is sharper than direct substitution of N=k+2N=k+2 into Theorem 2 in its state-dependent term, while matching it in its input-dependent term. Remark 3 (Implications for Closed-Loop Stability). This bound characterizes local prediction accuracy rather than closed-loop stability or asymptotic convergence. Specifically, (45) quantifies the two-step deviation of the actual nonlinear state from the linearized model prediction under the applied input, without implying state boundedness or target tracking under D2OCD^2OC. Thus, Corollary 1 serves as the foundational disturbance bound required for such analysis rather than a standalone stability proof. VI Simulation We evaluate the proposed controller on two control-affine nonlinear multi-agent systems: a unicycle team and a quadrotor team. The proposed SCP procedure relinearizes the dynamics at each decision step and solves the resulting strictly convex local QP, as characterized in Theorem 1. We compare its coverage performance and computational cost with a nonlinear NMPC baseline under the same objective, horizon, and input constraints. The simulation also evaluates the linearization error bound established in Theorem 2. VI-A Unicycle Model Five agents with state x=[px,py,θ]⊤x=[p_x,p_y,θ] and input u=[v,ω]⊤u=[v,ω] cover a 100×100100× 100 m domain. The sampling time and prediction horizon are set to T=0.1T=0.1 s and H=20H=20, respectively. Only the first input of each horizon is applied before relinearization. The reference distribution ν contains 500500 samples drawn from a five-component Gaussian mixture, as shown in Fig. 1. The initial positions of the five unicycle agents are indicated by cross marks. To evaluate the proposed method, ten simulation runs are performed, with the reference sample points independently redrawn for each of the ten seeds. The sample weights are initially assigned uniformly as 1500 1500 and then decreased according to the rule in [5] as the agents cover the domain. The simulation is terminated when the sum of the remaining sample weights falls below one percent of its initial value. Fig. 1: Five-unicycle team trajectories: (a) NMPC baseline and (b) proposed method. Reference samples are shaded by weight; crosses and circles denote initial and final positions. Fig. 1 shows that both the NMPC and SCP controllers effectively cover the domain according to the prescribed reference density. The Wasserstein distance in Fig. 2(a), which quantifies the discrepancy between the time-averaged agent distribution and the reference distribution, further illustrates the density-driven coverage performance. The mean 2(ρk,ν)W_2( _k,ν) over the final ten evaluations is 5.285.28 for the proposed method and 5.235.23 for NMPC, with no significant difference in variance observed. These results indicate that the linearized update based on Theorem 1 preserves coverage performance comparable to the nonlinear baseline in this simulation. Fig. 2: Unicycle team over ten seeds: (a) Wasserstein distance between the time-averaged agent distribution ρk _k and reference density ν and (b) Ratio of ‖ek+2‖||e_k+2|| to the proposed bound. Fig. 2(b) evaluates the error bound in Theorem2 using (45). Since only the first input of each horizon is applied before relinearization, we evaluate the first nonzero Taylor-remainder case at N=k+2N=k+2, corresponding to Corollary 1. Across the ten seeds, the worst-case and median values are plotted over the simulation steps. The ratio |ek+2|/bound|e_k+2|/bound remains below one for all evaluated steps, consistent with the bound in Corollary 1. At each solve, εx _x and εu _u are instantiated as the largest state and input excursions, respectively, of the unconstrained plan over t=k,…,k+H−1t=k,…,k+H-1. Thus, Assumption 3 is verified on the generated trajectory rather than imposed a priori. The constants κ, L, M¯ M, and M¯′ M are instantiated analogously. Rather than fixing them a priori as uniform bounds over a pre-specified ×X×U, we evaluate κ and L using the maximum Jacobian norms of Ψ sampled along the realized pairs (xt⋆,ut⋆)(x_t ,u_t ) and (xtun,utun)(x_t^un,u_t^un) at each decision time k. The quantities M¯ M and M¯′ M are obtained from finite-difference estimates of the sensitivity of the same Jacobians with respect to x over the corresponding samples. Accordingly, Assumptions 1 and 2 are evaluated on the realized operative tube at each k, rather than as fixed uniform bounds over a pre-specified region. To evaluate the computational advantage of the proposed method, the computation time is measured for each run. The median computation time across the ten runs is 1.611.61 ms for the proposed SCP method and 9.209.20 ms for NMPC, corresponding to a factor of 5.75.7 reduction. This result demonstrates the computational benefit of replacing the nonlinear optimization with the strictly convex QP characterized in Theorem 1. VI-B Quadrotor Model We next evaluate the method on four quadrotors with a full twelve-state rigid-body model in a 100×100×100100× 100× 100 m domain. The input consists of collective thrust and body torques. The sampling time and prediction horizon are set to T=0.05T=0.05 s and H=12H=12, respectively. The reference distribution contains 600600 sample points drawn from a four-component Gaussian mixture, as depicted in Fig. 3, with the samples independently redrawn for each of the ten seeds. The initial positions of the four quadrotors are indicated by cross marks. Although the density-driven coverage trajectories of the NMPC and SCP methods in Fig. 3 differ, the quantitative measure in Fig. 4(a) shows that both controllers achieve essentially the same coverage level. The averages over the final ten evaluations are 12.7412.74 for both methods, with a difference of only 0.0030.003. These results indicate that the proposed SCP-based linearization can drive the nonlinear multi-agent system toward the prescribed reference distribution while maintaining comparable coverage performance to the NMPC baseline. Fig. 3: Four-quadrotor team trajectories: (a) NMPC baseline and (b) proposed method. Crosses and circles denote initial and final positions. Fig. 4: Quadrotor team over ten seeds: (a) Wasserstein distance between ρk _k and ν, and (b) Ratio of ‖ek+2‖||e_k+2|| to the proposed bound. Fig. 4(b) further presents the worst-case and median values of the bound ratio across the ten simulation runs, supporting the error-bound analysis in Theorem2. The parameters κ, L, M¯ M, and M¯′ M are again instantiated locally at each decision time using the same procedure as for the unicycle model. Finally, the computational advantage becomes more pronounced for the quadrotor case due to its higher state dimension. The median computation time is 3.643.64 ms for the proposed SCP method versus 162.2162.2 ms for NMPC, corresponding to a factor of 4545, compared with a factor of 5.75.7 for the unicycle case. VII Conclusion This paper extended D2OC to discrete-time control-affine nonlinear multi-agent systems using sequential convex programming. The proposed approach preserves the Wasserstein-based coverage objective and barycentric structure while reducing each local update to a strictly convex QP with input constraints. We further quantified the deviation between the constrained and unconstrained local QP solutions and incorporated this control deviation into a finite-horizon bound on the error between the nonlinear trajectory and its linearized prediction. Simulations with unicycle and quadrotor teams showed that the proposed method achieves coverage performance comparable to NMPC while substantially reducing computation time. These results demonstrate the potential of the proposed formulation for computationally efficient nonlinear coverage control. Future work will focus on closed-loop stability, recursive feasibility, and extensions to more general nonlinear constraints. 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. External Links: Document Cited by: §I. [2] C. Villani et al. (2009) Optimal transport: old and new. Vol. 338, Springer. Cited by: §I. [3] Y. Chen, T. T. Georgiou, and M. Pavon (2021) Optimal transport in systems and control. Annual Review of Control, Robotics, and Autonomous Systems 4, p. 89–113. External Links: Document Cited by: §I. [4] D. Inoue, Y. Ito, and H. Yoshida (2021) Optimal transport-based coverage control for swarm robot systems: generalization of the voronoi tessellation-based method. IEEE Control Systems Letters 5 (4), p. 1483–1488. External Links: Document Cited by: §I. [5] 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, §VI-A. [6] 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, §I, Proposition 1. [7] 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, §I. [8] K. Lee (20262026) Density-driven optimal control: convergence guarantees for stochastic lti multi-agent systems. Automatica 193. External Links: Document Cited by: §I, §I, §I. [9] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert (2000) Constrained model predictive control: stability and optimality. Automatica 36 (6), p. 789–814. External Links: Document Cited by: §I. [10] W. B. Dunbar and R. M. Murray (2006) Distributed receding horizon control for multi-vehicle formation stabilization. Automatica 42 (4), p. 549–558. External Links: Document Cited by: §I. [11] Q. T. Dinh, C. Savorgnan, and M. Diehl (2011) Real-time sequential convex programming for nonlinear model predictive control and application to a hydro-power plant. In Proceedings of the 50th IEEE Conference on Decision and Control and European Control Conference, p. 583–588. External Links: Document Cited by: §I. [12] Y. Mao, M. Szmuk, and B. Açıkmeşe (2016) Successive convexification of non-convex optimal control problems and its convergence properties. In 2016 IEEE 55th Conference on Decision and Control (CDC), p. 3636–3641. Cited by: §I. [13] Y. Mao, M. Szmuk, and B. Açikmese (2018) A tutorial on real-time convex optimization based guidance and control for aerospace applications. In 2018 American Control Conference (ACC), p. 2410–2416. External Links: Document Cited by: §I. [14] B. G. Pachpatte et al. (1998) Inequalities for differential and integral equations. Technical report Academic press. Cited by: §V.