GPU-Parallel Linearization Error Bounds for Real-Time Robust Optimal Control of Nonlinear and Neural Network Dynamics Keyi Shen⋆
Abstract— This paper studies real-time robust optimal control for uncertain nonlinear systems, where linear time-varying (LTV) approximations make planning tractable but require sound linearization error bounds (LEBs) to guarantee robust constraint satisfaction. We develop tight, differentiable, GPUparallel LEBs for LTV approximations of nonlinear and neural network (NN) dynamics. For analytic dynamics, we introduce path-based Hessian bounds that are tighter than standard interval methods. For NN dynamics, we derive certified LEBs using NN verifier-generated affine relaxations and local Jacobian corrections. We adapt a GPU-parallel system-level synthesis LTV-based robust control solver to be compatible with these LEBs by extending it to handle right-invertible disturbance matrices and non-zero-centered disturbance sets for tight zonotopic uncertainty propagation. Our method, GPUSLS-LEO, enables online optimization of robust feedback policies that account for linearization error, producing tight, formally verified reachable tubes. On complex nonlinear and NN dynamics up to 168 state dimensions, our method can compute robust control policies on the GPU at rates up to 67 Hz, reducing solve times and conservativeness relative to baselines while preserving formal guarantees and real-time performance. Project Website
Anutam Srinivasan
Glen Chou
(a)
(b)
(c)
arXiv:2607.01203v1 [eess.SY] 1 Jul 2026
Jeffrey Fang⋆
Fig. 1. (a): Robust tubes from GPUSLS-LEO for a neural T pusher system and disturbed rollouts for a rotation and a push trajectory against random disturbances. (b): MPC rollout using our method on real dynamics, successfully moving the T to the goal. (c): Robust tubes of x and y position showing tight tubes and all simulated rollouts staying within the tubes.
§ Code (GitHub)
I. I NTRODUCTION Safe real-time nonlinear control under uncertainty is essential for the resilient operation of robots, drones, and spacecraft. This motivates methods for solving the robust nonlinear optimal control problem (RNOCP), which seeks a control policy that robustly guarantees constraint satisfaction under worst-case disturbances. In practice, RNOCP solvers typically optimize (1) a nominal trajectory and (2) a stabilizing feedback controller to track it under disturbance. Robust constraint satisfaction is enforced by tightening the nominal constraints using an overapproximation of the closed-loop tracking-error reachable tube, ensuring that the closed-loop dynamics remain in a safe subset of the state and input space. However, exact reachability via error propagation through nonlinear dynamics [1]–[3] is generally intractable in real time. To address this, [4]–[7] approximate the nonlinear dynamics along a nominal trajectory with a linear time-varying (LTV) system and compute reachable sets using the LTV dynamics. System-level synthesis (SLS) [4], [8] is one such LTV robust control framework that can be solved in real time over long horizons for high-dimensional systems [9], [10]. Despite strong performance, LTV-based methods only guarantee robust constraint satisfaction given a bound on the linearization error, which captures the deviation between the ⋆ Equal contribution. All authors are with the Georgia Institute of Technology, Atlanta, GA, USA, 30308. {jfang301, kshen84,
asrinivasan350, chou}@gatech.edu
true nonlinear dynamics and their LTV approximation. Existing methods often use global bounds [4], [11], which yield conservative (loose) reachable tubes. Tighter local bounds around the nominal trajectory are often possible, but must be computed efficiently for real-time planning and informative enough to guide the optimizer toward safe trajectories with minimal linearization error. Thus, the challenge is to develop bounds that are simultaneously (1) tight, (2) real-time computable, and (3) differentiable, so an RNOCP solver can use them to reshape the nominal trajectory and tighten the reachable tubes. As existing methods do not achieve all three, many LTV-based RNOCP solvers must omit linearization error, forgoing formal guarantees for performance [12], [13]. To close this gap, we propose a family of tight, differentiable linearization error bounds (LEBs) for RNOCP solvers with nonlinear dynamics described by analytic functions and neural networks (NNs), with a GPU-parallel implementation in JAX. To make these LEBs compatible with the realtime SLS solver of [9], [10], we extend the solver 1) to accommodate right-invertible disturbance matrices, enabling direct zonotopic propagation of the LEBs while avoiding the conservative overapproximations in [9]; and 2) to handle disturbance sets that are not zero-centered. Finally, we show empirically that these bounds enable real-time, GPU-based computation of formally verified RNOCP solutions, with limited conservativeness for systems with up to 168 states and solve times as low as 15 ms. Our contributions are: • A set of tight, provably-sound linearization error bounding
methods, including path-based interval bounds for analytic systems and NN verifier-based bounds for NN dynamics. • A novel nonlinear SLS formulation that includes linearization error through tight zonotopic propagation, including treatment of right-invertible disturbance matrices and nonzero-centered disturbance sets, guaranteeing robust constraint satisfaction for the original nonlinear system. • A GPU-accelerated nonlinear SLS solver, GPUSLS-LEO, that uses a formulation of these linearization bounds in JAX to enable differentiable, GPU-parallel evaluation and gradient-based robust nonlinear optimal control. • Evaluation on analytic and NN dynamics models, maintaining real-time control rates (up to 67 Hz) for systems up to 168D and problems with ≈ 2 × 105 decision variables. II. R ELATED W ORK While linearization-free reachability methods exist [1]–[3], direct error propagation through nonlinear dynamics is generally intractable in real time, especially beyond ≈10 states [1], [3]. Thus, RNOCP solvers often rely on linearized error propagation for efficiency [13], [14], but introduce linearization error that is difficult to tightly bound, and is often ignored [10], [13]–[15], underapproximated via sampling [4], [11], or handled by assuming access to a linearization bound oracle [16]–[20]. RNOCP solvers that explicitly compute LEBs typically rely on global Lipschitz estimates [11] or interval arithmetic [5]–[7], [21]. Interval arithmetic has been used both to compute reachable sets [22], [23] and to support reachability by constructing linear differential inclusion bounds (LDI) for nonlinear dynamics [24]. In particular, [25] uses Jacobian overbounding to obtain tighter LDIs for reachability, but is not used for control design or directly compute LEBs. Computationally, existing linearization bounding methods are offline [4], slow, or non-differentiable [7], limiting their use for planning. To close these gaps, our method uses interval arithmetic and JAX to compute tight differentiable LEBs for direct use in optimal control. Moreover, we also adapt neural-network (NN) verification tools [26] to derive affine bounds on linearization error for NN dynamics. As LEBs can be slow to compute [4], they are often generated once offline [4], [11] but are loose if we wish to bound error around one trajectory, as in RMPC. Thus, practitioners often omit linearization error, yielding heuristic reachable tubes that trade formal guarantees for performance [12]. To address this, we show that tight, local, GPU-parallel LEBs can be computed online while providing both real-time performance on high-dimensional systems and hard guarantees. III. P RELIMINARIES AND P ROBLEM S TATEMENT Notation: We denote In ∈ Rn×n as the identity matrix, m×n the as ∥A∥F := p Frobenius norm of a matrix A ∈ R Trace(A⊤ A). Let [N ] := {0, ..., N − 1}, and [M, N ] := {M, ..., N } for M, N ∈ N. For a matrix A ∈ Rn×p we define the row-wise ℓ1 -norm ∥A∥1,r ∈ Rn by ∥A∥1,r := ⊤ [∥A1,: ∥1 , . . . , ∥An,: ∥1 ] . For a ∈ Rn , we define diag(a) ∈ n×n R as its diagonal matrix. Given A ∈ Rm×n with linearly independent rows, we define A‡ ∈ Rn×m such
that AA‡ = Im , where A‡ := A⊤ (AA⊤ )−1 . We denote the Minkowski sum of two sets A and B as A ⊕ B := {a + b | a ∈ A, b ∈ B}. For a function p(x) with x ∈ Rn (or ∂p |x=α and explicitly p ((x1 , . . . , xn ))), we define ∂j p(α) := ∂x j 2
p ∂jk p(α) := ∂x∂j ∂x |x=α as the first and second partial derivak
tives with respect to the j-th and k-th components of the vector input. We define In as the set of all n-dimensional inˆ ∈ In , then ∆· ∆ ˆ := {δ · δ̂ | δ ∈ ∆, δ̂ ∈ ∆}. ˆ tervals. Let, ∆, ∆ Definitions: We consider uncertain nonlinear dynamics xk+1 = f (xk , uk ) + E(xk )wk ,
(1)
nx
where xk ∈ X ⊆ R is the system state at time step k, uk ∈ U ⊆ Rnu denotes the control input, f : Rnx × Rnu → Rnx the dynamics function, E : Rnx → Rnx ×nx the disturbance scaling function, and wk ∈ Enx := {w ∈ Rnx , ∥w∥∞ ≤ 1} the disturbance, normalized to be contained in a unit ℓ∞ -ball. We aim to design a controller π(·) that solves the following robust nonlinear optimal control problem (RNOCP): min
J (x̄, π(·))
(2a)
π(·)
s.t.
∀k ∈ [N ],
xk+1 = f (xk , uk ) + E(xk )wk , x0 = x̄0 ,
(2c) ∀k ∈ [N ],
uk = πk (x0 , . . . , xk ), gk (xk , uk ) ≤ 0, f
g (xN ) ≤ 0,
(2b)
∀wk ∈ Enx ,
∀wk ∈ Enx ,
∀k ∈ [N ],
∀k ∈ [N ].
(2d) (2e) (2f)
−1 where π := {πi }N i=0 is a sequence of causal control policies, x̄0 ∈ Rnx is the initial state, functions gk : Rnx × Rnu → f nc
denote stagewise state-input constraints, and g : R Rnx → Rnf is the terminal constraint. We define the closedloop reachable set Ωk := (Ωxk , Ωuk ) of states and controls as
xj+1 = f (xj , uj ) + E(xj )wj , ∀j ∈ [k], ∀wj ∈ Enx , ∀j ∈ [k], . Ωk := (xk , uk ) uj = πj (x0 , . . . , xj ), ∀j ∈ [k + 1]. (3)
Because (2) is an intractable infinite-dimensional problem, RNOCP solvers often simplify it by computing (A) a nominal state-input trajectory and (B) a tracking controller around that trajectory. In this paper, we use one such method, GPUSLS [10], which adapts [9] for real-time solution on the GPU. We then extend [10] to account for linearization error in Sec. V. Our approach uses SLS, which we now review. A. System Level Synthesis (SLS) Consider the following uncertain LTV dynamics xk+1 = Ak xk + Bk uk + Ek wk .
(4)
SLS is a control design framework that optimizes over causal disturbance feedback controllers uk = vk +
k−1 X
Φuk,j Ej wj ,
(5)
j=0
with nominal control vk ∈ Rnu , i.e., we assign a disturbance feedback matrix Φuk,j ∈ Rnu ×nx for each disturbance Ej wj
and control uk for k > j. For LTV dynamics (4), it can be shown using algebraic manipulations [8] that the resulting closed-loop state sequence can be expressed as xk = zk +
k−1 X
Φxk,j Ej wj , z0 = x̄0 ,
(6)
j=0
where zk ∈ Rnx is the nominal state and the closed-loop response Φxk,j ∈ Rnx ×nx captures the effect of disturbance wj on state xk . Starting with Φxj+1,j = Inx , SLS propagates the disturbance via Φxk+1,j = Ak Φxk,j + Bk Φuk,j , for all j ∈ [N ] and k ∈ [j + 1, N − 1]. Since (5)–(6) are the true closedloop trajectory under a given disturbance sequence, Lk−1 the exact closed-loop reachable sets are Ωkx = zk + j=0 Φxk,j Ej Enx Lk−1 and Ωuk = vk + j=0 Φuk,j Ej Enx . To apply the methods of Sec. III-A to the nonlinear dynamics (1), we approximate (1) with the nominal dynamics zk+1 = f (zk , vk ) and LTV error dynamics along a nominal state-control trajectory −1 nx ×(N +1) nu ×N z := {zk }N , v := {vk }N . k=0 ∈ R k=0 ∈ R nu nx For a nominal (zk , vk ) ∈ R × R , we construct an LTV approximation of the form (4) by defining Ak :=
∂f ∂f , Bk := , Ek := E(zk ). (7) ∂x (zk ,vk ) ∂u (zk ,vk )
Denote δzk := xk − zk and δvk := uk − vk . Then, the residual error dk (xk , uk , wk ) between (1) and its first-order Taylor approximation about (zk , vk ), i.e. f (zk , vk )+Ak δzk + Bk δvk , can be written as, dk (xk ,uk , wk ) := Ek wk + rlin + rdist , |k {z k } rkdist := E(xk ) − Ek wk , :=r(xk ,uk ,wk ) rklin := f (xk , uk ) − f (zk , vk ) − Ak δzk − Bk δvk .
(8a) (8b) (8c)
Specifically, (8c) captures the Taylor remainder of the nonlinear dynamics f , while (8b) accounts for the error induced by the state dependence of the disturbance matrix E. Given this setup, we study the following problems. B. Problem Statement Problem 1 (Linearization Error Bound). For all k ∈ [N ], find a set-valued overapproximation Ēk ⊆ Rnx and associated scaling matrix R̄k ∈ Rnx ×nx and offset c̄k ∈ Rnx where rk (xk , uk , wk ) ∈ Ēk ⊆ R̄k Enx ⊕ c̄k , (9) (3) ∀(xk , uk ) ∈ Ωk , ∀wk ∈ Enx . Problem 2 (Real-time RNOCP). Using Ēk to account for the discrepancy between (1) and its LTV approximation, efficiently solve the RNOCP (2) using the LTV surrogate model (4) to robustly satisfy the constraints (2e) and (2f). IV. GPU-PARALLEL L INEARIZATION E RROR B OUNDS To streamline exposition we consider a vector-valued function p(α) : Rn → Rm , over the augmented state α := [x⊤ , u⊤ ]⊤ . Let A ⊂ Rn be a compact set and α̃ ∈ A be a nominal point. The linearization error, i.e., the Lagrange remainder of the first-order Taylor expansion of p at α̃ is, r(α; α̃) := p(α) − p(α̃) − J(α̃)(α − α̃),
α ∈ A. (10)
For the specific case when p(α) = f (x, u), (10) coincides with (8c). Our goal is to compute a sound componentwise interval enclosure R ∈ Im of r(A; α̃), i.e., r(A; α̃) ⊆ R. A. Linearization Error Bounds for Analytic Dynamics By Taylor’s theorem in Lagrange form, for each component ℓ ∈ {1, . . . , m}, there exists ξ ∈ A such that rℓ (α; α̃) = 12 δ ⊤ Hℓ (ξ)δ,
δ := α − α̃.
(11)
where Hℓ (ξ) is the ℓ-th item of the Hessian H at ξ. To evaluate the interval enclosure of rℓ (A; α̃), previous work [5]–[7], [21], separately evaluate δ ⊤ , Hℓ (ξ) and δ in interval arithmetic, and then perform interval products as in (12), yielding the enclosure Rℓ,classic : rℓ (A; α̃) ⊆ Rℓ,classic := 21 ∆⊤ [Hℓ ](A)∆,
(12)
where ∆ := A ⊖ α̃ ⊂ Rn denotes the displacement set and [Hℓ ](A) is an interval enclosure of the Hessian over A with entries [Hℓ,jk ](A) (where Hℓ,jk (α) := ∂jk pℓ (α)). In practice, α̃ is usually chosen in the center of A [7], leading to a symmetric ∆ for smaller enclosure of rℓ (A; α̃). Retaining the tight quadratic form, (11) can also be first expanded componentwise, followed by performing interval arithmetic on the expansion. It yields (13), where interval ∆j denotes the j-th slice of ∆, after the symmetry of the interval Hessian is applied. When ∆ is symmetric, Rℓ,expand ⊆ Rℓ,classic is guaranteed because (13) preserves the nonnegativity of δk2 , for all δk ∈ ∆k and for all k ∈ {1, . . . , n}: P rℓ (A; α̃) ⊆ Rℓ,expand = 1≤j<k≤n [Hℓ,kj ](A) ∆j · ∆k Pn + 21 k=1 [Hℓ,kk ](A) ∆2k . (13) However, the classical bounds in (12) and (13) evaluate all interval enclosures of the Hessian entries over the full set A, causing conservativeness due to dependency and wrapping effects. To mitigate this, we use a coordinate-wise path from α̃ to α that enables tighter bounds on second-order terms. For any α ∈ A, define the sequence of intermediate points (k)
αα̃ := (α1 , . . . , αk , α̃k+1 , . . . , α̃n ), (0)
(n)
k ∈ [n + 1], (14)
with αα̃ = α̃ and αα̃ = α. This defines a coordinatewise path [24], [25] that incrementally moves from α̃ to α along each coordinate. Correspondingly, define the prefix sets A(k) (α̃) := {α ∈ A | αℓ = α̃ℓ , ∀ℓ ∈ [k + 1, n]}. Theorem 1 (Path-Based Hessian Bound): Let p : Rn → m R be twice continuously differentiable on an interval A ⊂ Rn , and let α̃ ∈ A. Define the displacement set ∆ := A ⊖ α̃. For each α ∈ A, define the coordinate-wise path α(k) and the corresponding prefix sets A(k) (α̃). Then, for each component ℓ ∈ {1, . . . , m}, the linearization error satisfies P rℓ (A; α̃) ⊆ Rℓ,path:= 1≤j<k≤n [Hℓ,kj ](A(j) (α̃)) ∆j · ∆k Pn + 12 k=1 [Hℓ,kk ](A(k) (α̃)) ∆2k . (15) Proof: Fix α ∈ A and define δ := α − α̃. Using telescoping along the coordinate-wise path, we have
Pn (k) (k−1) pℓ (α) − pℓ (α̃) = k=1 pℓ (αα̃ ) − pℓ (αα̃ ) . For each k, applying Taylor’s theorem in 1D yields (k)
(k−1)
(k−1)
)δk + 12 ∂kk pℓ (ξk ) δk2 , (16) (k−1) (k) for some ξk ∈ [αα̃ , αα̃ ] ⊆ A(k) (α̃). Summing (16) over k and subtracting ∇pℓ (α̃)⊤ δ yields pℓ (α) − pℓ (α̃) − ∇p (α̃)⊤ δ on the LHS which is rℓ (α; α̃) and the RHS is in Pn Pnℓ (k−1) ) − ∂k pℓ (α̃) δk + 12 k=1 ∂kk pℓ (ξk ) δk2 . k=1 ∂k pℓ (αα̃ Expanding the gradient difference along the same path for k, Pk−1 (k−1) (j) we have ∂k pℓ (αα̃ ) − ∂k pℓ (α̃) = j=1 ∂k pℓ (αα̃ ) − (j−1) ∂k pℓ (αα̃ ) . By the Mean Value Theorem, for each (j−1) (j) j < k, there exists ηj,k ∈ [αα̃ , αα̃ ] ⊆ A(j) such that (j) (j−1) ∂k pℓ (αα̃ ) − ∂k pℓ (αα̃ ) P = ∂kj pℓ (ηj,k ) δj . Substituting yields rℓ (α; α̃) = 1≤j<k≤n ∂kj pℓ (ηj,k ) δj δk + Pn 1 2 (j) (k) ∂ p (ξ ) δ . Since η (α̃), j,k ∈ A (α̃), ξk ∈ A k k=1 kk ℓ k 2 and δ ∈ ∆, taking interval enclosures yields (15). Thm. 1 implies that our path-based Hessian bound is never looser than the expanded interval Hessian bound, as each Hessian entry is evaluated over a subset of the full domain. Corollary 1 (Path-Based Hessian is Never Looser): When ∆ is symmetric, the enclosure in (15) is guaranteed to be no looser than (13) and (12). That is, pℓ (αα̃ ) − pℓ (αα̃
) = ∂k pℓ (αα̃
Rℓ,path ⊆ Rℓ,expand ⊆ Rℓ,classic . (17) Proof: For each j ∈ {1, . . . , n}, the prefix set satisfies A(j) ⊆ A by construction. Then, for every Hessian entry, [Hℓ,kj ](A(j) ) ⊆ [Hℓ,kj ](A), [Hℓ,kk ](A(k) ) ⊆ [Hℓ,kk ](A). Multiplying by ∆j ∆k and ∆2k preserves inclusion and summing the resulting terms yields Rℓ,path ⊆ Rℓ,expand . Unlike Rℓ,classic , Rℓ,expand and Rℓ,path preserve componentwise expansion to produce tight, off-centered intervals (via the quadratic form rather than symmetric bounds). We exploit the added tightness in the RNOCP setting (Sec. V).
r(α) := (W − J(α̃))α + J(α̃)α̃ − p(α̃) + b, r(α) := (W − J(α̃))α + J(α̃)α̃ − p(α̃) + b.
(20)
Then, an interval enclosure of r(A; α̃) ⊆ [r, r] is given by r i = min r i (α), α∈A
r i = max r i (α), α∈A
∀i ∈ [1, m], (21)
Proof: Substituting the CROWN lower bound (18) into the expression for r(α; α̃) (10), we yield r(α; α̃) ≥ Wα + b − p(α̃) − J(α̃)(α − α̃), which simplifies to r(α; α̃) ≥ (W − J(α̃))α + J(α̃)α̃ − p(α̃) + b = r(α). Similarly, we apply the CROWN upper bound. Thus, ∀α ∈ A: r(α) ≤ r(α; α̃) ≤ r(α). For bounds valid over A, we take the componentwise minimum of r(α) and maximum of r(α) over A, as in (21). Thm. 3 reconciles the CROWN slopes and the local linearization slope J(α̃), thereby producing a sound remainder bound around the nominal point. Since r(α) and r(α) are affine in α, problem (21) is inexpensive when A is a box set. C. Zonotope Bounds of rk (xk , uk , wk ) To obtain the zonotope bound for rk (xk , uk , wk ) (9), we lin compute the enclosure [rlin k , r k ] of the linearization error dist rklin as in Sec. IV-A or IV-B and the enclosure [rdist k , rk ] of the disturbance error rkdist via interval arithmetic in (22), where [E](Ωxk ) is the interval enclosure of E over Ωxk : rkdist = E(xk ) − Ek wk ∈ [E](Ωxk ) ⊖ Ek Enx . (22) Finally, we obtain rk (xk , uk , wk ) ∈ [rk , rk ], where rk = dist dist rlin and rk = rlin k + r k . Its zonotope bound is k + rk given by R̄k Enx ⊕ c̄k , where c̄k = 12 (rk + rk ) and R̄k = 1 2 diag(r k − r k ). We denote this mapping compactly as R̄k = Γr (Ωk , zk , vk ; f, E), c̄k = Γc (Ωk , zk , vk ; f, E). (23) D. Differentiable, Parallel Bound Implementation
B. Linearization Error Bounds for Neural Dynamics For NN dynamics f , computing a tight interval of the Hessian can be difficult. Instead, we use CROWN [26], [27] to obtain affine bounds of an NN over a box input set, and convert them into an interval enclosure of r(α; α̃). Theorem 2 (CROWN bounds (Thm. 3.2 [26])): Given an NN, p : Rn → Rm , and a bounded input α ∈ A where A is an interval, there exist linear lower/upper bounds on p(α): Wα + b ≤ p(α) ≤ Wα + b,
(18)
where W, W ∈ Rm×n and b, b ∈ Rm . The bound in (18) does not directly bound r(α; α̃), since the slopes W, W generally do not match the Jacobian J(α̃). We address this via further relaxation in Thm. 3. Theorem 3 (Linearization error bound from CROWN): Let p : Rn → Rm be an NN, A ⊂ Rn be a box set, α̃ ∈ A, J(α̃) be the Jacobian of p at α̃, and let (W, b, W, b) satisfy the CROWN bounds (18) on A. Then, ∀α ∈ A: r(α) ≤ r(α; α̃) ≤ r(α),
(19)
We implement our LEBs in JAX, with the PBH using immrax for interval arithmetic and jax verify providing CROWN-based bounds for NN dynamics. JAX offers differentiability and GPU parallelization, enabling gradientbased refinement of the nominal point α̃ and scalable input domain partitioning to reduce conservativeness for SLS. V. U SING L INEARIZATION E RROR B OUNDS IN GPUSLS We discuss the use of our LEBs in SLS (Sec. V-A) and describe an efficient GPU-based implementation (Sec. V-B). A. Nonlinear SLS Formulation Following [9, App. A], SLS can be extended to nonlinear systems by planning a nominal trajectory for the disturbancefree dynamics zk+1 = f (zk , vk ) and modeling tracking error dynamics via the LTV system of (4). To solve (2) approximately, we modify the nonlinear SLS formulation of [9, App. A] to be compatible with 1) non-zero-centered
disturbance sets and 2) the LEBs given in Sec. IV-A-IV-B: min J(z, v) + H̃0 (Φ, D)
(24a)
z, v,Φ
s.t.
zk+1 = f (zk , vk ), ∀k ∈ [N ], z0 = x̄0 , (24b) Φxk+1,j = Ak Φxk,j + Bk Φuk,j , Φxj+1,j = Inx , (24c) ∀j ∈ [N ], ∀k ∈ [j + 1, N − 1], Ēk := Zr (zk , vk , µk , σk ), ∀k ∈ [N ], (24d) c̄k := Zc (zk , vk , µk , σk ), σk = hσk (Φ, D),
∀k ∈ [N ],
(24e)
µk = hµk (Φ, D),
gk (zk , vk ) + hk (Φ, D) ≤ 0, g f (zN ) + hf (Φ, D) ≤ 0,
∀k ∈ [N ], (24f) ∀k ∈ [N ], (24g) (24h)
where Ak , Bk are the linearized dynamics (7) at time k and functions Zr (·), Zc (·) (24d)–(24e) are the Minkowski sums of the exogenous disturbance and linearization error as a function of the nominal (zk ,vk ), reachable tube center offset µk ∈ Rnx +nu , and widths σk ∈ Rnx +nu respectively [4], −1 [12], and D := {(c̄j , Ēj )}N j=0 . Formally, for all k ∈ [N ] in nx (24d)–(24e), c̄k ∈ R and Ēk ∈ Rnx ×2nx are given as Ak = [zk⊤ , vk⊤ ]⊤ + µk ⊕ diag(σk )Enx +nu , (25a) Ēk := [E(zk ) R̄k ] = E(zk ) Γr (Ak , zk , vk ; f, E) , (25b) {z } | c̄k := Γc (Ak , zk , vk ; f, E), Zr (zk ,vk ,µk ,σk ) (25c) {z } |
⊤ where Φk,j := Φxk,j ⊤ Φuk,j ⊤ , Φ collects Φk,j for all u x . The constraint , R̃ ∈ Sn++ k, j ∈ [N ], and Q̃, Q̃N ∈ Sn++ tightenings and tube widths in (28)-(30) are derived for the worst case disturbances using the dual norm property [28]. The cost H̃0 (Φ, D) (31) penalizes uncertainty (the tube widths (28)). Following [9], [10], we solve (24) via an iterative process that alternates between finding a (A) nominal trajectory and a (B) robust controller. To solve (A), we find (z, v) by solving a constraint-tightened NOCP (32), min J(z, v)
Robust constraint satisfaction is enforced through constraint tightenings hk (Φ, D) and hf (Φ, D). In short, these tightenings capture the propagation of unit-normalized disturbances in (26), which are then scaled by the disturbance bounds Ēj and shifted by offset c̄j . For simplicity, we derive the tightenings linear constraints, i.e., assuming ⊤ gk (zk , vk ) = Gk zk , vk and g f (zN ) = Gf zN , for Gk ∈ Rnc ×(nx +nu ) , GN ∈ Rnc ×nx . Valid tightenings can be obtained for nonlinear gk and g f by including additional error terms [12, Eq. 22]. Formally, we define the tube center (27) and width (28), stacked (terminal) constraint tightenings (29)-(30), and surrogate system-level response cost (31) as, Pk−1 hµk (Φ, D) = j=0 Φk,j c̄j (27) Pk−1 σ hk (Φ, D) = j=0 Φk,j Ēj 1,r (28) Pk−1 hk (Φ, D) = j=0 (Gk Φk,j c̄j + Gk Φk,j Ēj 1,r ) (29) PN −1 f h (Φ, D) = j=0 (GN ΦxN,j c̄j + GN ΦxN,j Ēj 1,r ) (30) PN −1 PN −1 H̃0 (Φ, D) = j=0 ( k=j (∥Q̃1/2 Φxk,j Ēj ∥2F (31) 1/2
+ ∥R̃1/2 Φuk,j Ēj ∥2F ) + ∥Q̃N ΦxN,j Ēj ∥2F ) ,
zk+1 = f (zk , vk ),
s.t.
∀k ∈ [N ],
gk (zk , vk ) + hk (Φ, D) ≤ 0, f
z0 = x̄0 , (32b) ∀k ∈ [N ],
f
g (zN ) + h (Φ, D) ≤ 0.
(32c) (32d)
Then, (32) mirrors the structure of (2) and can be solved using the GPU-parallel method of [10, Sec. IV]. After solving (A) (32), we linearize around the resulting (z, v) and solve for a robust controller (B) by finding a Φ that optimizes (31): min x u
Φ ,Φ
N −1 N −1 X X
∥Qk,j Φk,j Ēj ∥2F + ∥QN,j ΦxN,j Ēj ∥2F (33a)
j=0 k=j
s.t. Φxk+1,j = Ak Φxk,j + Bk Φuk,j , Φxj+1,j = Inx ,
Zc (zk ,vk ,µk ,σk )
where c̄k is the offset and Ēk is the generator matrix of a zonotope for the Minkowski sum of the exogenous disturbance and linearization error: E(zk )Enx ⊕ (c̄k ⊕ R̄k Enx ) ≡ c̄k ⊕ Ēk E2nx . Using (25b) and (25c), the resulting reachable tube overapproximation is, Pk−1 Lk−1 Ω̄xk := zk + j=0 Φxk,j c̄j + j=0 Φxk,j Ēj E2nx , (26) Pk−1 Lk−1 Ω̄uk := vk + j=0 Φuk,j c̄j + j=0 Φuk,j Ēj E2nx .
(32a)
z, v
(33b)
∀j ∈ [N ], ∀k ∈ [j + 1, N − 1], (33c)
where Qk,j , QN,j are defined according to [9, Eq. 20, 23] (see App. V) for consistency between (32) and (33). Notably, the right-multiplication by Ēj prevents the cost in (33a) from being expressed as a standard LQR objective in Φx and Φu . Thus, it is not directly solvable using the Riccati-based methods of [9], [10]. Instead, we solve (34), which has an LQR cost and can be solved via the efficient solvers [9], [10]: min x u
Φ ,Φ
s.t.
N −1N −1 X X
∥Qk,j Φk,j ∥2F + ∥QN,j ΦxN,j ∥2F
(34a)
j=0 k=j Φxk+1,j = Ak Φxk,j + Bk Φuk,j , Φxj+1,j = Inx , ∀j ∈ [N ], ∀k ∈ [j + 1, N − 1].
(34b) (34c)
Note that (34) coincides with (33) if Ēj = I for all j ∈ [N ]. In Prop. 1, we prove that the set of minimizers of (33) and (34) are equivalent if Ēj is right-invertible for all j ∈ [N ] (see App. III for the proof). Proposition 1: Let (Φ̃x , Φ̃u ) be an optimizer of (33), and let (Φ̂x , Φ̂u ) be an optimizer of (34). If Ēj is right-invertible for all j ∈ [N ], i.e., there exists Ēj‡ such that Ēj Ēj‡ = Inx , then the two optimizers coincide: Φ̃x = Φ̂x , Φ̃u = Φ̂u . When implementing the method of Sec. IV, we add a small positive diagonal padding to R̄k so that it is always invertible, and thus Ēj (25) is always right-invertible for all j ∈ [N ]. After solving (34), we compute the constraint tightening terms in (27)–(30). This updates the reachable tube centers µk and widths σk : σk ← hσk (Φ, D) and µk ← hµk (Φ, D), which are used as constraint tightenings in the next iteration of nominal trajectory optimization (32). This iterates until convergence (as in [9, Alg. 1]); we terminate after a fixed
B. Efficient SLS Implementation on the GPU We propose GPUSLS-Linearization Error Optimization (GPUSLS-LEO), which builds on GPUSLS [10] to explicitly account for linearization error propagation. We first replace the ellipsoidal disturbance representation of [10] with a zonotope, modifying the constraint tightenings and cost terms via (29)-(31). At each controller update, we modify GPUSLS by including a zonotopic bound on the linearization error via (25). We then modify the initialization of Φx via (34c). All other steps closely follow the GPUSLS procedure. In particular, GPUSLS-LEO solves (24) iteratively via sequential quadratic programming (SQP). At each iteration, GPUSLS-LEO quadraticizes the cost and linearizes xk+1 = f (xk , uk ) around the current (z, v) to obtain Ak , Bk (4). It then solves (A) a local QP approximation of (32) to update (z, v), (B) (34) to get Φ, and updates (C) the tightenings via (24f). Since our LEBs are differentiable, we can penalize error accumulation during step (A) by modifying the nominal ¯ v) := J(z, v) + λJerr (z, v), where objective (32a) to J(z, Jerr (z, v) =
N −1 N −1 X X
∥Φk,j Zr (zj , vj , µj , σj )∥2F ,
A. Linearization Error We compare the LEB from our path-based Hessian (PBH) method on the Satellite (nx = 7, nu = 3) and Quadrotor (nx = 12, nu = 4) systems (App. I) with several baselines. Global sampling [4] estimates worst-case curvature offline (taking over 7 minutes [4]) over the global set X × U, yielding an empirical over-approximation. Interval Hessian (IH) (classic) [6], [7] evaluates (12) on a local set T ⊂ (X × U). IH (classic)-CORA denotes the standard IH implementation in CORA [6], [29]. Random sampling uses 10K samples from the local set and gives an under-approximation. We sample nominal points (z, v) from X × U and define local sets T := (z, v) ⊕ ϵEnx +nu with ϵ = 0.01, 0.05, 0.1, 0.15, 0.2. Fig. 2 shows the mean interval width of (23). For Satellite, where nonlinearity is mild, all over-approximation methods remain close to the randomsampling lower bound as ϵ grows, though PBH is still the tightest. For Quadrotor, PBH consistently gives the tightest bounds. By contrast, global sampling is loose at small ϵ, suggesting global curvature misses local behavior, while the interval Hessian becomes overly conservative at ϵ = 0.20. Table I compares online runtime. PBH runs in under 1 ms, enabling real-time use. Although slower than the global bound, which is efficiently calculated online after slow offline computation, PBH is 60% and 27% (satellite and quadrotor respectively) faster than dense random sampling and negligible running time compared to CORA IH.
number of iterations or if the change in the iterates is below a threshold (see App. II for an algorithm block). Formulation Novelty: Here, we summarize the changes made to prior nonlinear SLS formulations [9, App. A] [4]: 1) Zonotopes enable exact tightening computation (27)– (30) by concatenating the generators of the exogenous disturbance and linearization error, unlike prior ellipsoidal methods [9, App. A], but this yields right-invertible Ēj incompatible with existing solvers [9], [10]. We address this by setting Φxj+1,j = Inx and including the disturbance scaling directly into the tightenings and cost (31), with a compatible reformulation (34) justified by Prop. 1. 2) Finally, we extend [4], [9] to handle non-zero-centered disturbances, enabling direct use of interval bounds from (23) without extra conservativeness due to zero-centering.
!!
(a)
(b)
Fig. 2. (a): LEBs on the satellite system. (b): LEBs on the quadrotor system. In both systems, our method achieves the tightest over-approximation.
(35)
j=0 k=j
penalizes the propagated residual error (8a). Here, Φ, µj , σj are from the previous iteration and λ is a weighting term. Notably, (35) biases (z, v) toward regions of low propagated exogenous and linearization error disturbances. Finally, we show that by solving (24), GPUSLS-LEO ensures containment in the reachable tubes (26), i.e., xk ∈ Ω̄xk and uk ∈ Ω̄uk , for all k ∈ [N ]. Formally, we have (proof in App. IV): Theorem 4: Let z, v, and Φ be a feasible solution for (24). Then, (26) is guaranteed to overapproximate the true closed-loop reachable set (3), i.e., Ωk ⊆ Ω̄k for all k ∈ [N ]. VI. R ESULTS We evaluate LEB tightness (Sec. VI-A) and the robust controllers given by GPUSLS-LEO (Sec. VI-B). Runtime tests use an NVIDIA RTX 4090, except the multi-quadrotor and long-horizon experiments, which require an NVIDIA H200. System definitions appear in App. I.
TABLE I RUNTIME ( IN MS ) COMPARISON OF LINEARIZATION BOUND METHODS . System
Random
Global
IH-classic
IH-classic-CORA
Path-based
Satellite Quadrotor
0.894 1.278
0.143 0.156
0.218 0.453
11570 227200
0.361 0.936
TABLE II L INEARIZATION ERROR UNDER NOMINAL OPTIMIZATION VIA PGD. System
Initial
Sampling
PGD
PGD better (%)
Satellite Quadrotor
0.2680 0.02253
0.2365 0.01639
0.2241 0.01590
100 100
To quantify the benefit of differentiability in our pathbased Hessian implementation, we search for a nominal point within a fixed perturbation radius ϵ that minimizes the LEB size. Starting from random nominal points, we run projected
(a)
(b)
)
)
)
(
)
(a) )
Fig. 3. Satellite (7D). (a): Tube sizes for our method, GPUSLS-LEO, compared to baselines and the variant without linearization error. Our method is the least conservative, while minimally increasing tube size (b): Tubes for the first four state dimensions; our tubes capture 100% of sampled rollouts.
gradient descent (PGD) within a local neighborhood of radius ϵopt to reduce the error interval width. We compare PGD with random sampling of nominal points under the same time budget, where both methods use the PBH to evaluate linearization error. Table II shows that across 100% of test cases, PGD yields smaller error interval widths than random sampling across both systems, showing the value of gradientbased optimization for finding nominal points with tighter LEBs, which can benefit downstream SLS. B. Using Linearization Bounds for RNOCP In this section, we evaluate GPUSLS-LEO on a suite of analytic and neural dynamical systems against various RNOCP baselines. We further assess the impact of key ablations, including the effect of optimizing over the linearization error gradients and the use of nonzero-centered zonotopes. All benchmarked methods are evaluated under the same disturbance set for the given dynamical system and, where applicable, are run for a maximum of 100 SQP iterations. 1) Satellite: We evaluate our method on a 7D satellite attitude control task, where the objective is to stabilize the system to a desired orientation under bounded disturbances. We adopt the experimental setup of [4] and benchmark against their approach. In Fig. 3 (a) we compare the tube widths σk for: GPUSLS-LEO, GPUSLS-E (ellipsoidal propagation as in [9, App. A]), GPUSLS (no linearization error propagation) [10], and NL-SLS ( [4]). Our method achieves the tightest tubes among all baselines. We have on average 20% smaller tubes than NL-SLS, which relies on a global maximum Hessian bound and is therefore locally overlyconservative along a trajectory. Compared to GPUSLS-E, our method achieves 45% smaller tubes, as it involves overapproximating the Minkowski sum of ellipsoids, whereas our zonotope-based combination is exact. We observe that incorporating our LEB pipeline does not introduce significant conservativeness, as the tube widths increase by only 4% compared to GPUSLS. In Fig. 3 (b), we show the rollouts of 20 random and 64 adversarial disturbances for the quaternion states, showing 100% containment under our method. 2) Planar quadrotor: We evaluate our method on a 6D planar quadrotor against control contraction metrics (CCM) [30]. Both methods must navigate through a dense
Fig. 4. Planar quadrotor (6D). (a): Tube area comparison against CCM and immrax, showing our method is less conservative. (b): Tubes and rollouts for GPUSLS-LEO and CCM. Both methods achieve 100% safety, but the conservativeness of the CCM forces the system to take a suboptimal path.
)
(b)
(a)
(b) Fig. 5. Planar quadrotor (6D). (a): System navigating through an obstacle field over a 20 m trajectory with a planning horizon length of 7500. (b): Tube widths along the horizon. Despite the long horizon, the tubes remain bounded due to the GPUSLS-LEO controller’s optimization of tube sizes.
obstacle field over a horizon of 40 steps. Fig. 4(a) shows our method achieving 91% smaller tube area than the CCM for the (x, y) position states. This arises because the CCM certifies safety over the entire state space, introducing inherent conservativeness. GPUSLS-LEO, on the other hand, certifies safety along a specific trajectory, enabling tighter, trajectorydependent bounds. We also evaluate the closed-loop system on an interval-based reachability analysis tool immrax [24], showing its calculated tube area diverging after a few steps. In Fig. 4(b) we compare the calculated trajectories and tubes for our method and CCM. Our approach is less conservative, allowing a more direct, task-optimal path, while the CCM is forced to take a longer route. In Fig. 5(a), we show the trajectory and tubes over a long horizon of 7500 steps, which involves optimizing a state-control trajectory with 1.95×105 decision variables. Due to the tightness of our LEBs and the optimization of the closed-loop controller in SLS, our tubes remain well-bounded and finite over the horizon (Fig. 5(b)), demonstrating its scalability to long horizons. 3) Coupled 3D Quadrotor: To demonstrate scalability to high dimensional systems, we evaluate our method on a coupled system of 14 3D quadrotors (168D) over a horizon of 25. The system is tasked with navigating through an obstacle field while the quadrotors are coupled via a spring force toward the group centroid. We show the planned trajectories and the tubes in Fig. 6(a), and observe all rollouts remain in the tight tubes despite the high dimensional state space,
Runtime (ms)
(b)
(a)
10 2 GPUSLS-LEO GPUSLS PBH
10 1 10 0
2
3
4
11
12
13
14
Fig. 6. System of 14 coupled 3D quadrotors (168D). (a): Computed trajectories and tubes in an obstacle field. Even at high dimensions, our method is able to compute tight tubes. (b): Tubes and 100 rollouts for states py and vy for 4 quadrotors. 100% of rollouts stay safely in the tubes. (b)
6 7 8 9 10 Number of Quadrotors
Fig. 8. Coupled Quadrotor. Per-iteration runtime evaluation across increasing numbers of coupled quadrotors. GPUSLS scales favorably with state dimension, growing approximately logarithmically with state and control space. Path-Based Hessians (PBH) achieves low runtimes at small dimensions but scales more poorly as the state dimension increases.
5
dynamics, where the pusher successfully completes the task. We note that we do not formally consider learning error; however, it can be incorporated as in [31], [32]. TABLE III P ER ITERATION RUNTIME BREAKDOWN ( MS ) ACROSS SYSTEMS . System
(c)
Satellite Planar Quadrotor Coupled Quadrotor (96D) Neural T-Dynamics
GPUSLS Lin. Error Grad. Remainder Total 11 18 40 630
3 1 – 100
<1 <1 48 32
15 20 88 762
DT 1000 150 100 1000
(a)
Fig. 7. Neural T-pusher (5D). (a): Tube size comparison across ablations of our method, showing nonzero-centered tubes and linearization error gradients reduce tube widths. (b): Rollouts of our method; all trajectories remain in the computed robust tubes. (c): Rollouts of GPUSLS. Without formally considering linearization error, rollouts leave the tube.
demonstrating the method’s scalability to state dimensionality. In Fig. 6(b), we plot the tubes and the corresponding rollouts for py and vy for 100 adversarial disturbances. 4) Neural T-pusher: We demonstrate that our formulation can account for linearization error propagation through learned neural dynamics. We train an MLP-based model for the T-pusher system and task it with rotating and pushing the object to the goal configuration (Fig. 1). We evaluate our method on several ablations (Fig. 7(a)), including GPUSLS, GSL (shorthand for GPUSLS-LEO), GSL-C (zero-centered zonotopes) and GSL (λ = 1), which incorporates a penalty on the linearization error as defined in V-B. Compared to GPUSLS, GPUSLS-LEO produces tubes that are on average 99% larger. However, Fig. 7(b) shows that our method achieves 100% containment of all rollouts, while Fig. 7(c) illustrates how trajectories from GPUSLS fall out of the tube, showing the need to formally consider linearization error. Fig. 7(a) shows that adding linearization error penalties can reduce conservativeness, with λ = 1 decreasing tube widths by 8%. We also show the importance of nonzerocentered zonotopes, as using zero-centered zonotopes leads to 25% larger tubes. We present rollouts and disturbance deviations in Fig. 1(a). In Fig. 1(b), we demonstrate the method’s ability for closed-loop MPC on simulated contact
5) Runtime: We evaluate the runtime of our method across all systems. In Tab. III, we show the per iteration breakdown in milliseconds of the different components of our methodology. In all systems aside from the coupled quadrotor, the main computation load is in the GPUSLS, showing that our proposed modifications do not drastically increase run time. We also show that all total runtimes fall below the discretization time step DT, demonstrating the viability of our method under a real-time iteration MPC scheme [10]. Finally, we evaluate the method’s scaling with growing state and control dimensions using an increasing number of coupled quadrotors. Fig. 8 shows that GPUSLS incurs higher initial runtime but scales favorably with state and control dimension due to its logarithmic scaling [10]. PBH on the other hand, exhibits low initial runtime, but can slow down in higher dimensions due to the cost of Hessian computations. VII. C ONCLUSION This paper presents a framework for computing tight, differentiable, and sound linearization error bounds (LEBs) for nonlinear analytic and NN dynamics and for using them in GPUSLS-LEO, a real-time GPU-parallel RNOCP solver based on SLS. The approach reduces conservativeness, outperforming state-of-the-art reachability analysis tools, while maintaining formal guarantees and achieving performance on systems with up to 168 states, real-time control rates (up to 67 Hz), and problems with ≈ 2 × 105 decision variables. Future work will investigate how our LEBs can be extended to bound linearization errors for contact-rich planning [28] and nonlinear output-feedback control, including through extensions of the SLS-based methods described in [33]–[35].
R EFERENCES [1] A. Majumdar and R. Tedrake, “Funnel libraries for real-time robust feedback motion planning,” IJRR, vol. 36, no. 8, 2017. [2] M. Althoff, “Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets,” in HSCC, 2013. [3] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-jacobi reachability: A brief overview and recent advances,” in CDC, 2017. [4] A. Leeman, J. Köhler, A. Zanelli, S. Bennani, and M. Zeilinger, “Robust nonlinear optimal control via system level synthesis,” TAC, 2025. [5] L. Schäfer and M. Althoff, “Robust optimal control using set-based reachability analysis,” in Eur. Control Conf. (ECC). IEEE, 2025. [6] D. D. Leister and J. P. Koeln, “Robust model predictive control for nonlinear discrete-time systems using iterative time-varying constraint tightening,” in American Control Conference (ACC). IEEE, 2025. [7] M. Althoff, O. Stursberg, and M. Buss, “Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization,” in CDC. IEEE, 2008, pp. 4042–4048. [8] J. Anderson, J. C. Doyle, S. H. Low, and N. Matni, “System level synthesis,” Annu. Rev. Control., vol. 47, pp. 364–393, 2019. [9] A. P. Leeman, J. Kohler, F. Messerer, A. Lahr, M. Diehl, and M. N. Zeilinger, “Fast system level synthesis: Robust model predictive control using riccati recursions,” IFAC, vol. 58, no. 18, 2024. [10] J. Fang and G. Chou, “Safe large-scale robust nonlinear mpc in milliseconds via reachability-constrained system level synthesis on the gpu,” in Robotics: Science and Systems (RSS), 2026. [11] T. Kim, P. Elango, and B. Açıkmeşe, “Joint synthesis of trajectory and controlled invariant funnel for discrete-time systems with locally lipschitz nonlinearities,” IJRNC, vol. 34, no. 6, 2024. [12] S. Zhan, C.-Y. Chiu, A. Leeman, and G. Chou, “Robustly constrained dynamic games for uncertain nonlinear dynamics,” in ICRA, 2026. [13] F. Messerer and M. Diehl, “An efficient algorithm for tube-based robust nonlinear optimal control with optimal linear feedback,” 2021. [14] Z. Manchester and S. Kuindersma, “Dirtrel: Robust nonlinear direct transcription with ellipsoidal disturbances and lqr feedback,” 2017. [15] M. Diehl, H. G. Bock, and E. Kostina, “An approximation technique for robust nonlinear optimization,” Math. Prog., vol. 107, no. 1, 2006. [16] B. Houska, F. Logist, J. Van Impe, and M. Diehl, “Robust optimization of nonlinear dynamic systems with application to a jacketed tubular reactor,” Journal of Process Control, vol. 22, no. 6, 2012. [17] M. Villanueva, R. Quirynen, M. Diehl, B. Chachuat, and B. Houska, “Robust mpc via min–max differential inequalities,” Automatica, 2017. [18] M. Cannon, D. Ng, and B. Kouvaritakis, “Successive linearization nmpc for a class of stochastic nonlinear systems,” in NMPC, 2009. [19] M. Cannon, J. Buerger, B. Kouvaritakis, and S. Rakovic, “Robust tubes in nonlinear model predictive control,” TAC, vol. 56, no. 8, 2011. [20] A. Richards, “Robust model predictive control for time-varying systems,” in CDC. IEEE, 2005, pp. 3747–3752. [21] M. Rungger and M. Zamani, “Accurate reachability analysis of uncertain nonlinear systems,” in HSCC, 2018, pp. 61–70. [22] D. Limon, J. Bravo, T. Alamo, and E. Camacho, “Robust mpc of constrained nonlinear systems based on interval arithmetic,” IEE Proceedings-Control Theory and Applications, vol. 152, no. 3, 2005. [23] K. Shen and G. Chou, “Parallel differentiable reachability for learning and planning with certified neural dynamics and controllers,” in Robotics: Science and Systems (RSS), 2026. [24] A. Harapanahalli, S. Jafarpour, and S. Coogan, “immrax: A parallelizable and differentiable toolbox for interval analysis and mixed monotone reachability in jax,” ADHS, vol. 58, no. 11, 2024. [25] A. Harapanahalli and S. Coogan, “A linear differential inclusion for contraction analysis to known trajectories,” IEEE TAC, 2025. [26] H. Zhang, T.-W. Weng, P.-Y. Chen, C.-J. Hsieh, and L. Daniel, “Efficient neural network robustness certification with general activation functions,” NeurIPS, vol. 31, 2018. [27] K. Xu, Z. Shi, H. Zhang, Y. Wang, K.-W. Chang, M. Huang, B. Kailkhura, X. Lin, and C.-J. Hsieh, “Automatic perturbation analysis for scalable certified robustness and beyond,” NeurIPS, 2020. [28] W.-C. Li and G. Chou, “Certified gradient-based contact-rich manipulation via smoothing-error reachable tubes,” in Robotics: Science and Systems (RSS), 2026. [29] M. Althoff, “An introduction to CORA 2015,” in Workshop on Applied Verification for Continuous and Hybrid Systems, 2015. [30] A. Sasfi, M. N. Zeilinger, and J. Köhler, “Robust adaptive mpc using control contraction metrics,” Automatica, vol. 155, p. 111169, 2023.
[31] A. Srinivasan, A. Leeman, and G. Chou, “Safety beyond the training data: Robust out-of-distribution mpc via conformalized system level synthesis,” Learning for Dynamics and Control (L4DC), 2026. [32] D. Nath, A. Srinivasan, H. Yin, R. Jiang, J. Fang, and G. Chou, “Pixels to proofs: Probabilistically-safe latent world model control via parallel conformal robust mpc,” arXiv preprint arXiv:2606.15594, 2026. [33] A. P. Leeman, S. Zhan, M. N. Zeilinger, and G. Chou, “Vision-sls: Safe perception-based control from learned visual representations via system level synthesis,” Robotics: Science and Systems (RSS), 2026. [34] H. Zhou and V. Tzoumas, “Safe control of partially-observed linear time-varying systems with minimal worst-case dynamic regret,” in 2023 62nd IEEE Conference on Decision and Control (CDC). IEEE, 2023, pp. 8781–8787. [35] W.-C. Li, J. Fang, S. Polisetti, Y. Song, and G. Chou, “Robustness without wrinkles: Parallel simulation and robust mpc for certified deformable manipulation,” arXiv preprint arXiv:2606.14188, 2026.
A PPENDIX
C. Coupled 3D Quadrotor
In the appendix, we first discuss details on the system dynamics used in the experimental results (App. I). We then provide an algorithm block describing how sequential quadratic programming is used to solve GPUSLS-LEO (App. II). Next, we discuss and prove Proposition 1 and a supporting lemma in App. III. We then provide proofs for the remaining theoretical results that were omitted from the main manuscript in App. IV. Finally, we discuss in more detail the Lagrange multiplier variable τ used in our framework in App. V. A PPENDIX I DYNAMICAL S YSTEM D EFINITIONS In this section, we give an overview of the discrete-time dynamical system definitions for each system. The discretetime dynamics are obtained via forward Euler discretization with timestep ∆t, xk+1 = xk + ∆t(ẋk ). A. Satellite We use the same satellite dynamics and experimental setup as [4]. Specifically, we define the dynamics as,
ζ(ω)q IS−1 (v − w × (IS ω)) 0 −ω1 −ω2 −ω3 1 ω1 0 ω3 −ω2 ζ(ω) := 0 ω1 2 ω2 −ω3 ω3 ω2 −ω1 0 ẋ =
(36a)
(36b)
with states x := (q, ω) ∈ R7 , where q ∈ R4 is the attitude quaternion, ω ∈ R3 the angular rotation rate, v ∈ R3 the input control torque, and IS = diag(5, 2, 1) the symmetric inertia matrix. We define a constant disturbance scaling of E = 5 · 10−3 · diag(0, 0, 0, 0, 1, 1, 1).
We consider a system of N quadrotors. For a given quadrotor i, we define its dynamics as vx,i ṗx,i vy,i ṗy,i v z,i ṗz,i pi + (qi sin ϕi + ri cos ϕi ) tan θi ϕ̇i q cos ϕ − r sin ϕ i i i i θ̇i qi sin ϕi +ri cos ϕi cos θi ˙ξi = ψ̇i = Ti . v̇ m (cos ψi sin θi cos ϕi + sin ψi sin ϕi ) x,i Ti v̇ m (sin ψi sin θi cos ϕi − cos ψi sin ϕi ) y,i Ti v̇ m (cos θi cos ϕi ) − g z,i τ J −J y z ϕ,i ṗ Jx qi ri + Jx i τ Jz −Jx θ,i q̇i Jy pi ri + Jy τ Jx −Jy ṙi pi qi + ψ,i Jz
Jz
(38a) The state of quadrotor i is defined as ξi := (px,i , py,i , pz,i , ϕi , θi , ψi , vx,i , vy,i , vz,i , pi , qi , ri ) ∈ R12 where (px,i , py,i , pz,i ) denotes position, (ϕi , θi , ψi ) the roll, pitch, and yaw angles, (vx,i , vy,i , vz,i ) the translational velocities, and (pi , qi , ri ) denote body angular velocities. The input is defined as υi := (Ti , τϕ,i , τθ,i , τψ,i ) ∈ R4 , where Ti is the collective thrust and (τϕ,i , τθ,i , τψ,i ) are the body torques. For a system of N coupled quadrotors, we stack the states and inputs, giving the overall state and control vectors, x := (ξ1 , . . . , ξN ) ∈ R12N ,
u := (υ1 , . . . , υN ) ∈ R4N . (39) We additionally incorporate a centroid spring-damper coupling in the translational dynamics. We define the centroid position and velocity as N
1 X p̄k = pi,k , N i=1
N
1 X v̄k = vi,k . N i=1
(40)
The coupling acceleration applied to quadrotor i is B. Planar Quadrotor We consider a planar quadrotor with the following dynamics vx vy ϕ̇ ẋ = (37a) − 1 (u1 + u2 ) sin(ϕ) 1 m (u1 + u2 ) cos(ϕ) − g m L J (u2 − u1 ) The state is defined as x := (px , py , ϕ, vx , vy , ϕ̇) ∈ R6 and input u := (u1 , u2 ) ∈ R2 . The position is denoted as (px , py ), ϕ the pitch angle, (vx , vy ) the translational velocities, and ϕ̇ the angular velocity. The inputs (u1 , u2 ) are the individual rotor thrusts. We define the mass m = 2.0576, gravitational acceleration g = 9.81, arm length L = 0.25 and moment of inertia J = 0.01. We define the disturbance scaling as E = 5 · 10−2 · diag(0, 0, 0, 1, 1, 0).
acoup i,k = −αN (pi,k − p̄k ) − γ(vi,k − v̄k )
(41)
where α and γ are the coupling gains, with α specifying the stiffness of the spring term and γ controlling the damping coefficient. In our experiments, we set α = 0.75 and γ = 0.25. Therefore, the discrete-time dynamics become ṗx,i,k ṗy,i,k ṗz,i,k ϕ̇i,k θ̇ i,k . ξi,k+1 = ξi,k + ∆t (42) ψ̇i,k coup v̇i,k + a i,k ṗi,k q̇ i,k
ṙi,k We use a disturbance scaling matrix of E = 2·10−2 ·I12N
A PPENDIX II A LGORITHM B LOCK
Proof: Observe that (33) can be rewritten into (43) with the variable change Ψk,j = Φk,j Ēj , since by using (33b) and (33c) we have that,
Algorithm 1 GPUSLS-LEO Require: x̄0 , J, f (x, u), E(x), gk (x, u), g f (x) 1: while SQP NOT CONVERGED do 2: A, B, G ← SQP: L INEARIZE S YSTEM 3: Q, M, R, q, r ← SQP: Q UADRATICIZE C OST 4: while SLS NOT CONVERGED do 5: z, v ← O PTIMIZE N OMINAL T RAJECTORY 6: σ, µ ← C OMPUTE INTERVAL 7: Ē, c̄ ← C OMPUTE RESIDUAL ERROR ZONOTOPE 8: τ , Q ← U PDATE DUALS AND COST 9: Φx , Φu ← O PTIMIZE C ONTROLLER 10: h(Φ, D), hf (Φ, D) ← U PDATE T IGHTENINGS 11: end while 12: end while 13: return z, v, Φx , Φu Here, Q, M, R, q, r correspond to the quadratic approximation of the nominal cost. A PPENDIX III O N S OLVING FOR Φ u
Solving (33) yields optimal values for Φ and Φ which penalize the size of the tightenings (31). However, this leads to terms in (33a) that scale the response matrices Φx and Φu by Ēj in the cost (33a). This remains a quadratic cost in Φ but no longer maintains the structure of an LQR cost function. To address this, we instead solve (34) which has an LQR structure and can be solved via Riccati recursions. For self-containment, we restate the original problem (33) and the reformulated problem (34), and the equivalent optimizer Proposition (Prop. 1) below: Original Problem:
Φ ,Φ
N −1 N −1 X X
∥Qk,j Φk,j Ēj ∥2F + ∥QN,j ΦxN,j Ēj ∥2F (33a)
j=0 k=j
s.t. Φxk+1,j = Ak Φxk,j + Bk Φuk,j ,
Φxj+1,j Ēj = Inx Ēj
min x u
Ψ ,Ψ
N −1N −1 X X
⇐⇒ Eq. (43c).
∥Qk,j Ψk,j ∥2F + ∥QN,j ΨxN,j ∥2F
(43a)
j=0 k=j s.t. Ψxk+1,j = Ak Ψxk,j + Bk Ψuk,j , (43b) x Ψj+1,j = Ēj , ∀j ∈ [N ], ∀k ∈ [j + 1, N − 1] (43c)
From [10], (43) can be solved via the parallel Riccati recursions, GN,j = QxN,j ,
Kk,j = −Gk,j Bk,j , ⊤ ⊤ Pk,j = Qxk,j + Ak Pk+1,j Ak + Kk,j Bk,j , ⊤ −1 Gk,j = Quk,j + Bk Pk+1,j Bk , ⊤ Bk,j = Qux Pk+1,j Ak , k,j + Bk
(44)
Ψxj+1,j = Ēj ,
Ψuk,j = Kk,j Ψxk,j ,
Ψxk+1,j = (Ak + Bk Kk,j )Ψxk,j .
(45)
Then, since Ēj is right-invertible and (43) represents the system level constraint (SLC) assumed in Lemma 1, the controller gains Ψ obey, Ψxj+1,j = Ēj ,
Ψuk,j Ēj‡ = Kk,j Ψxk,j Ēj‡
Ψxk+1,j Ēj‡ = (Ak + Bk Kk,j ) Ψxk,j Ēj‡ .
(46)
Then, by unraveling the recursion (46), we have for k > j, k−1 Y I Ψk,j Ēj‡ = nx (Am + Bm Km,j ) Ψxj+1,j Ēj‡ , Kk,j m=j+1
Reformulated Problem: PN −1PN −1 2 x 2 min j=0 k=j ∥Qk,j Φk,j ∥F + ∥QN,j ΦN,j ∥F (34a) x u Φ ,Φ
Φxj+1,j = Inx ,
and
(33b)
Φxj+1,j = Inx , ∀j ∈ [N ], ∀k ∈ [j + 1, N − 1], (33c)
s.t. Φxk+1,j = Ak Φxk,j + Bk Φuk,j ,
⇐⇒ Eq. (43b),
and parallel forward propagations x
min x u
Φxk+1,j Ēj = Ak Φxk,j Ēj + Bk Φuk,j Ēj
(34b)
∀j ∈ [N ], ∀k ∈ [j + 1, N − 1] (34c)
In Prop. 1, we prove that the optimizer sets of (33) and (34) are equivalent. Proposition 1 (Equivalent Optimizers): Let (Φ̃x , Φ̃u ) be an optimizer of (33), and let (Φ̂x , Φ̂u ) be an optimizer of the reformulated problem (34). Suppose that Ēj is rightinvertible for all j ∈ [N ], i.e., there exists Ēj‡ such that Ēj Ēj‡ = Inx . Then the two optimizers coincide: Φ̃x = Φ̂x , Φ̃u = Φ̂u .
(47) while the remaining gains (Ψk,j where k ≤ Q j) are zero, n where we define the ordered matrix product i=1 Ri := Rn Rn−1 · · · R1 . To recover the optimal gains for (33), we use the rightinvertibility property of Ēj and substitute Ψk,j Ēj‡ = Φ̃k,j into (47): k−1 Y I Φ̃k,j = nx (Am + Bm Km,j ) Φ̃xj+1,j Kk,j m=j+1 k−1 Y I = nx (Am + Bm Km,j ) Inx . (48) Kk,j m=j+1
The second equality stems from Φ̃xj+1,j satisfying (33c). Now, observe (34) has the same Riccati recursions as (44),
thus resulting in an optimal gain, similarly structured to (47) where Ej‡ is replaced with Inx , Φ̂: k−1 I nx Y Φ̂k,j = (Am + Bm Km,j ) Φ̂xj+1,j , (49) Kk,j
Φ̂u from (58) we arrive at the first equality in (52). The second equality arises from algebraic manipulation using the definition Ē ‡ = Ē ⊤ (Ē Ē ⊤ )−1 . Lastly, manipulating (52) to see K(Φx Ē ‡ ) = Φu Ē ‡ , and expanding the system level constraint in (51) we have that
m=j+1
where Φ̂xj+1,j = Inx , as constrained by (34c). Hence ∀ k, j ∈ [N ] Φ̃k,j = Φ̂k,j , so Φ̃ = Φ̂. Lemma 1: Let, Ē := blk diag(Inx , Ē0 , . . . , Ēn−1 ), where (i) (i) for all i ∈ [N ], Ēi ∈ Rnx ×nw is right-invertible, and nw ≥ nx . Enforcing the uncertain LTV dynamics xk+1 = Ak xk + Bk uk + Ēk wk , wk ∈ En(k) , w
(50)
is equivalent to enforcing the system level constraint (51), i.e., (50) ⇐⇒ (51), Φx (IN nx − ZA) −ZB = Ē. (51) Φu Furthermore, the controller gain mapping for the SLC in (51) is, K = (Φu E ‡ )(Φx Ē ‡ )−1 = (Φu Ē ⊤ )(Φx Ē ⊤ )−1 , ‡
(52)
‡ ). Lastly, given K, the = blk diag(I, Ē0‡ , . . . , Ēn−1
where Ē forward propagation of disturbance gains obeys: Φxj+1,j = Ēj ,
Φuk,j Ēj‡ = Kk,j Φxk,j Ēj‡
Φxk+1,j Ēj‡ = (Ak + Bk Kk,j ) Φxk,j Ēj‡ . Proof: First, let δk = Ek wk be such that xk+1 = Ak xk + Bk uk + δk .
(53)
(54)
⊤ ⊤ Now, let x = [x⊤ u = 0 , . . . , xn−1 ] , ⊤ ⊤ ⊤ ⊤ ⊤ ]⊤ . Then, [u0 , . . . , un−1 ] , and δ = [δ0 , . . . , δn−1 defining the gain mappings Φ̂x and Φ̂u such that, x x Φ̂ = δ, (55) u Φ̂u
we have from [8]: (IN nx − ZA)
−ZB
Φ̂x = IN nx Φ̂u
(56)
⊤ Defining w = [w0⊤ , . . . , wn−1 ]⊤ , we observe:
δ = Ēw
(57)
Then, the gain mapping from w to x and u via Φx and Φu follows from (55) and (57): x x x x Φ Φ̂ Ē Φ̂ = w= w= Ēw. (58) u u Φu Φ̂ Ē Φ̂u This establishes that Φx = Φ̂x Ē and Φu = Φ̂u Ē. Hence, by multiplying (56) by Ē, we arrive at (51). Thus we have shown (51) =⇒ (50). To show the equivalence between (50) and (51), we can show (50) =⇒ (51), by directly recovering the matrices Ak , Bk , Ēk ∀k ∈ [N ] from (51) to produce (50). To prove (52), we observe that, K = Φ̂u (Φ̂x )−1 , [8]. Then using the substitutions Φx Ē ‡ = Φ̂x and Φu Ē ‡ =
Φxj+1,j = Ēj ,
Φuk,j Ēj‡ = Kk,j Φxk,j Ēj‡
Φxk+1,j Ēj‡ = (Ak + Bk Kk,j ) Φxk,j Ēj‡ .
(59)
A PPENDIX IV P ROOFS Theorem 4: Let z, v, and Φ be a feasible solution for (24). Then, (26) is guaranteed to overapproximate the true closed-loop reachable set (3), i.e., Ωk ⊆ Ω̄k for all k ∈ [N ]. Proof: We will first prove (26) is correct assuming valid disturbances. Let Wk specify the disturbance ∀k ∈ [N ]. Then, from SLS [4], we have that the overapproximation of the reachable set is, Lk−1 Ω̄xk := zk + j=0 Φxk,j Wj , (60) Lk−1 Ω̄uk := vk + j=0 Φuk,j Wj . Now, suppose Wk can be expressed as a zonotope, c̄k ⊕ Ēk E2nx . Making this substitution in (60) we arrive at, Lk−1 Ω̄xk := zk + j=0 Φxk,j (c̄k ⊕ Ēk E2nx ), (61) Lk−1 Ω̄uk := vk + j=0 Φuk,j (c̄k ⊕ Ēk E2nx ). Using the associativity and commutativity of the Minkowski sum (61) can be rewritten into, Lk−1
Lk−1 x x j=0 Φk,j c̄j ⊕ j=0 Φk,j Ēj E2nx , Lk−1 u Lk−1 u u Ω̄k := vk + j=0 Φk,j c̄j ⊕ j=0 Φk,j Ēj E2nx . Ω̄xk := zk +
(62)
Lk−1 x Then recalling that j=0 Φk,j c̄j is a Minkowski sum of Pk−1 x singleton sets, we can replace it with j=0 Φk,j c̄j (and Lk−1 u similarly do this for j=0 Φk,j c̄j ), and arrive at at the disturbance reachable set (26). Now, we will show the disturbance bounds, Zr (zj , vj , µj , σj ), Zc (zj , vj , µj , σj ), in (24d)-(24e) are valid overapproximations Ēk (9) and certify that the closed-loop reachable sets Ω̄xk and Ω̄uk overapproximate the true reachable set Ωxk and Ωuk (3). Base Case: The disturbance at timestep k = 0 is precisely E(z0 ) = E(x̄0 ), which is captured by (25), and thus (24d)(24e) provides overapproximate disturbance bounds. Inductive Step: Suppose that for all timesteps j ∈ [k], dj (xj , uj , wj ) ∈ Ej Enx ⊕ R̄j Enx ⊕ c̄j . Then, we have a overapproximate reachable sets, Ωxk ⊆ Ω̄xk and Ωuk ⊆ Ω̄uk , at timestep k. Define Ω̄k := {[x⊤ , u⊤ ]⊤ | x ∈ Ω̄xk , u ∈ Ω̄uk }. Then observe that Ω̄k ⊆ Ak := [zk⊤ , vk⊤ ]⊤ + µk ⊕ diag(σk )Enx +nu , as defined in (25a). Then, by construction, ⊤ ⊤ for all [x⊤ ∈ Ω̄k and for all wk ∈ Enx , we have k , uk ] that r(xk , uk , wk ) ∈ c̄k ⊕ R̄k Enx , where (c̄k , R̄k ) is computed using Ak in (25). Similarly, the nominal disturbance Ek wk ∈ E(zk )Enx . Thus, the total disturbance (8a), rklin + rkdist + Ek wk = d(xk , uk , wk ) ∈ c̄k ⊕ R̄k Enx ⊕ E(zk )Enx .
Concatenating the two generator matrices produces (25b). Thus, the disturbance bound at timestep k in (24d)-(24e) is an overapproximation and contains the true residual error, i.e., dk (xk , uk , wk ) (8a). Having established valid disturbance bounds for all steps k ∈ [N ], SLS guarantees that the overapproximation of the closed-loop reachable set (26), Ω̄xk and Ω̄uk , are valid, i.e., Ωxk ⊂ Ω̄xk and Ωuk ⊂ Ω̄uk for all k ∈ [N ]. A PPENDIX V GPUSLS D UALS τ As discussed in Sec. V-A, the RNOCP problem is solved by iteratively alternating between solving a nominal trajectory optimization and a controller update. Since these two optimizations are solved separately, we use the dual variable τ to enforce consistency between the two optimizations as discussed in [10, App. B]. Due to our reformulation, we modify the procedure as follows. As with [10, App. B], we define the auxiliary terms βk,j ∈ Rnc and βN,j ∈ Rnf , where ◦2 βk,j = Gk Φk,j c̄j + ∥Gk Φk,j Ēj ∥1,r
(63a)
∀j ∈ [N ], ∀k ∈ [j, N ] ◦2 x βN,j = GN ΦN,j c̄j + ∥GN ΦxN,j Ēj ∥1,r
(63b)
∀j ∈ [N ]
(63c) (63d)
where [·]◦2 denotes element-wise square. We can then calculate our dual variables τ by λ̄k τk,j = p βk,j + ϵ
∀j ∈ [N ],
∀k ∈ [j, N ]
(64)
where we denote λ̄k as the dual variable of the nominal trajectory optimization (corresponding to the dual variable λk in [10]), introduced here to avoid notation conflict. Finally, we can define our cost terms as 1/2 √ Q̃ 0 , Qk,j = diag τk,j Gk , 0 R̃1/2 (65) √ 1/2 QN,j = diag τN,j GN , Q̃N . Following [9] and [10], this formulation preserves dual consistency and yields an equivalent reformulation of the RNOCP under the proposed decomposition.