Paper deep dive
GPU-Parallel Linearization Error Bounds for Real-Time Robust Optimal Control of Nonlinear and Neural Network Dynamics
Jeffrey Fang, Keyi Shen, Anutam Srinivasan, Glen Chou
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 7/5/2026, 10:03:28 AM
Summary
The paper introduces GPUSLS-LEO, a GPU-parallelized, differentiable framework for real-time robust optimal control (RNOCP) of nonlinear and neural network (NN) dynamics. It addresses the challenge of computing tight, real-time linearization error bounds (LEBs) to guarantee robust constraint satisfaction. For analytic dynamics, the authors propose a path-based Hessian bound that is tighter than standard interval methods. For NN dynamics, they utilize NN verifier-generated affine relaxations. The method is implemented in JAX and optimized for GPU, enabling real-time control rates (up to 67 Hz) for high-dimensional systems (up to 168 states) by handling right-invertible disturbance matrices and non-zero-centered disturbance sets via zonotopic uncertainty propagation.
Entities (6)
Relation Signals (4)
Jeffrey Fang → affiliatedwith → Georgia Institute of Technology
confidence 100% · All authors are with the Georgia Institute of Technology
GPUSLS-LEO → implements → Linearization Error Bounds
confidence 95% · Our method, GPUSLS-LEO, enables online optimization of robust feedback policies that account for linearization error
GPUSLS-LEO → addresses → Robust Nonlinear Optimal Control Problem
confidence 90% · Our method, GPUSLS-LEO, enables online optimization of robust feedback policies
GPUSLS-LEO → builton → JAX
confidence 90% · a GPU-parallel implementation in JAX
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
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, GPU-parallel 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.
Tags
Links
- Source: https://arxiv.org/abs/2607.01203v1
- Canonical: https://arxiv.org/abs/2607.01203v1
Trouble viewing inline? Open PDF directly →
Full Text
94,265 characters extracted from source content.
Expand or collapse full text
GPU-Parallel Linearization Error Bounds for Real-Time Robust Optimal Control of Nonlinear and Neural Network Dynamics Jeffrey Fang⋆ Keyi Shen⋆ Anutam Srinivasan Glen Chou ⋆ Equal contribution. All authors are with the Georgia Institute of Technology, Atlanta, GA, USA, 30308. jfang301, kshen84, asrinivasan350, chou@gatech.edu 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, GPU-parallel LEBs for LTV approximations of nonlinear and neural network (N) dynamics. For analytic dynamics, we introduce path-based Hessian bounds that are tighter than standard interval methods. For N dynamics, we derive certified LEBs using N 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 N 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 Code (GitHub) I Introduction 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, 2, 3] is generally intractable in real time. To address this, [4, 5, 6, 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) [8, 4] 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 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]. Figure 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. 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 real-time SLS solver of [10, 9], 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 N verifier-based bounds for N dynamics. • A novel nonlinear SLS formulation that includes linearization error through tight zonotopic propagation, including treatment of right-invertible disturbance matrices and non-zero-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 N dynamics models, maintaining real-time control rates (up to 67 Hz) for systems up to 168D and problems with ≈2×105≈ 2× 10^5 decision variables. I Related Work While linearization-free reachability methods exist [1, 2, 3], direct error propagation through nonlinear dynamics is generally intractable in real time, especially beyond ≈ 10 states [3, 1]. Thus, RNOCP solvers often rely on linearized error propagation for efficiency [14, 13], but introduce linearization error that is difficult to tightly bound, and is often ignored [14, 15, 13, 10], underapproximated via sampling [11, 4], or handled by assuming access to a linearization bound oracle [16, 17, 18, 19, 20]. RNOCP solvers that explicitly compute LEBs typically rely on global Lipschitz estimates [11] or interval arithmetic [7, 5, 6, 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 (N) verification tools [26] to derive affine bounds on linearization error for N 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. I Preliminaries and Problem Statement Notation: We denote In∈ℝn×nI_n ^n× n as the identity matrix, the Frobenius norm of a matrix A∈ℝm×nA ^m× n as ‖A‖ℱ:=Trace(A⊤A)\|A\|_F:= Trace(A A). Let [N]:=0,…,N−1[N]:=\0,...,N-1\, and [M,N]:=M,…,N[M,N]:=\M,...,N\ for M,N∈ℕM,N . For a matrix A∈ℝn×pA ^n× p we define the row-wise ℓ1 _1-norm ‖A‖1,r∈ℝn\|A\|_1,r ^n by ‖A‖1,r≔[‖A1,:‖1,…,‖An,:‖1]⊤\|A\|_1,r [\|A_1,:\|_1,…,\|A_n,:\|_1 ] . For a∈ℝna ^n, we define diag(a)∈ℝn×ndiag(a) ^n× n as its diagonal matrix. Given A∈ℝm×nA ^m× n with linearly independent rows, we define A‡∈ℝn×mA ^n× m such that AA‡=ImAA =I_m, where A‡:=A⊤(AA⊤)−1A :=A (A )^-1. We denote the Minkowski sum of two sets A and B as A⊕B:=a+b∣a∈A,b∈BA B:=\a+b a∈ A,b∈ B\. For a function p(x)p(x) with x∈ℝnx ^n (or explicitly p((x1,…,xn))p ((x_1,…,x_n) )), we define ∂jp(α):=∂p∂xj|x=α _jp(α):= ∂ p∂ x_j|_x=α and ∂jkp(α):=∂2p∂xj∂xk|x=α _jkp(α):= ∂^2p∂ x_j∂ x_k|_x=α as the first and second partial derivatives with respect to the j-th and k-th components of the vector input. We define nI_n as the set of all n-dimensional intervals. Let, Δ,Δ^∈n , _n, then Δ⋅Δ^:=δ⋅δ^∣δ∈Δ,δ^∈Δ · :=\δ· δ δ∈ , δ∈ \. Definitions: We consider uncertain nonlinear dynamics xk+1=f(xk,uk)+E(xk)wk,x_k+1=f(x_k,u_k)+E(x_k)w_k, (1) where xk∈⊆ℝnxx_k ^n_x is the system state at time step k, uk∈⊆ℝnuu_k ^n_u denotes the control input, f:ℝnx×ℝnu→ℝnxf:R^n_x×R^n_u ^n_x the dynamics function, E:ℝnx→ℝnx×nxE:R^n_x ^n_x× n_x the disturbance scaling function, and wk∈ℰnx≔w∈ℝnx,‖w‖∞≤1w_k _n_x \w ^n_x,\|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¯,π(⋅)) J ( x,π(·) ) (2a) s.t. xk+1=f(xk,uk)+E(xk)wk,∀k∈[N], x_k+1=f(x_k,u_k)+E(x_k)w_k, ∀ k∈[N], (2b) x0=x¯0, x_0= x_0, (2c) uk=πk(x0,…,xk),∀k∈[N], u_k= _k(x_0,…,x_k), ∀ k∈[N], (2d) gk(xk,uk)≤0,∀wk∈ℰnx,∀k∈[N], g_k(x_k,u_k)≤ 0, ∀ w_k _n_x, ∀ k∈[N], (2e) gf(xN)≤0,∀wk∈ℰnx,∀k∈[N]. g^f(x_N)≤ 0, ∀ w_k _n_x, ∀ k∈[N]. (2f) where π:=πii=0N−1π:=\ _i\_i=0^N-1 is a sequence of causal control policies, x¯0∈ℝnx x_0 ^n_x is the initial state, functions gk:ℝnx×ℝnu→ℝncg_k:R^n_x×R^n_u ^n_c denote stagewise state-input constraints, and gf:ℝnx→ℝnfg^f:R^n_x ^n_f is the terminal constraint. We define the closed-loop reachable set Ωk:=(Ωkx,Ωku) _k:=( _k^x, _k^u) of states and controls as Ωk:=(xk,uk)|xj+1=f(xj,uj)+E(xj)wj,∀j∈[k],∀wj∈ℰnx,∀j∈[k],uj=πj(x0,…,xj),∀j∈[k+1].. -7.0pt _k:= \(x_k,u_k)\; |\; aligned &x_j+1=f(x_j,u_j)+E(x_j)w_j,\ \ ∀ j∈[k],\\ & 18.49988pt∀ w_j _n_x, ∀ j∈[k],\\ &u_j= _j(x_0,…,x_j), ∀ j∈[k+1]. aligned \. (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. I-A System Level Synthesis (SLS) Consider the following uncertain LTV dynamics xk+1=Akxk+Bkuk+Ekwk.x_k+1=A_kx_k+B_ku_k+E_kw_k. (4) SLS is a control design framework that optimizes over causal disturbance feedback controllers uk=vk+∑j=0k−1k,juEjwj,u_k=v_k+ _j=0^k-1 ^u_k,jE_jw_j, (5) with nominal control vk∈ℝnuv_k ^n_u, i.e., we assign a disturbance feedback matrix k,ju∈ℝnu×nx ^u_k,j ^n_u× n_x for each disturbance EjwjE_jw_j and control uku_k for k>jk>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+∑j=0k−1k,jxEjwj,z0=x¯0,x_k=z_k+ _j=0^k-1 ^x_k,jE_jw_j,~z_0= x_0, (6) where zk∈ℝnxz_k ^n_x is the nominal state and the closed-loop response k,jx∈ℝnx×nx ^x_k,j ^n_x× n_x captures the effect of disturbance wjw_j on state xkx_k. Starting with j+1,jx=Inx ^x_j+1,j=I_n_x, SLS propagates the disturbance via k+1,jx=Akk,jx+Bkk,ju ^x_k+1,j=A_k ^x_k,j+B_k ^u_k,j, for all j∈[N]j∈[N] and k∈[j+1,N−1]k∈[j+1,N-1]. Since (5)–(6) are the true closed-loop trajectory under a given disturbance sequence, the exact closed-loop reachable sets are Ωxk=zk+⨁j=0k−1k,jxEjℰnx _x^k=z_k+ _j=0^k-1 ^x_k,jE_jE_n_x and Ωku=vk+⨁j=0k−1k,juEjℰnx _k^u=v_k+ _j=0^k-1 ^u_k,jE_jE_n_x. To apply the methods of Sec. I-A to the nonlinear dynamics (1), we approximate (1) with the nominal dynamics zk+1=f(zk,vk)z_k+1=f(z_k,v_k) and LTV error dynamics along a nominal state-control trajectory :=zkk=0N∈ℝnx×(N+1)z:=\z_k\_k=0^N ^n_x×(N+1), :=vkk=0N−1∈ℝnu×Nv:=\v_k\_k=0^N-1 ^n_u× N. For a nominal (zk,vk)∈ℝnx×ℝnu(z_k,v_k) ^n_x×R^n_u, we construct an LTV approximation of the form (4) by defining Ak:=∂f∂x|(zk,vk),Bk:=∂f∂u|(zk,vk),Ek:=E(zk).A_k:= . ∂ f∂ x |_(z_k,v_k),B_k:= . ∂ f∂ u |_(z_k,v_k),E_k:=E(z_k). (7) Denote δzk:=xk−zkδ z_k:=x_k-z_k and δvk:=uk−vkδ v_k:=u_k-v_k. Then, the residual error dk(xk,uk,wk)d_k(x_k,u_k,w_k) between (1) and its first-order Taylor approximation about (zk,vk)(z_k,v_k), i.e. f(zk,vk)+Akδzk+Bkδvkf(z_k,v_k)+A_kδ z_k+B_kδ v_k, can be written as, dk(xk, d_k(x_k, uk,wk):=Ekwk+rklin+rkdist⏟:=r(xk,uk,wk), u_k,w_k):=E_kw_k+ r_k^lin+r_k^dist_ :=r(x_k,u_k,w_k), (8a) rkdist r_k^dist :=(E(xk)−Ek)wk, := (E(x_k)-E_k )w_k, (8b) rklin r_k^lin :=f(xk,uk)−f(zk,vk)−Akδzk−Bkδvk. :=f(x_k,u_k)-f(z_k,v_k)-A_kδ z_k-B_kδ v_k. (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. I-B Problem Statement Problem 1 (Linearization Error Bound). For all k∈[N]k∈[N], find a set-valued overapproximation ℰ¯k⊆ℝnx E_k ^n_x and associated scaling matrix R¯k∈ℝnx×nx R_k ^n_x× n_x and offset c¯k∈ℝnx c_k ^n_x where rk(xk,uk,wk) r_k(x_k,u_k,w_k) ∈ℰ¯k⊆(R¯kℰnx⊕c¯k), ∈ E_k ( R_kE_n_x c_k ), (9) ∀(xk,uk)∈(3)Ωk,∀wk∈ℰnx. ∀(x_k,u_k) eq:reach∈ _k,\ ∀ w_k _n_x. Problem 2 (Real-time RNOCP). Using ℰ¯k E_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 Linearization Error Bounds To streamline exposition we consider a vector-valued function p():ℝn→ℝmp( α):R^n ^m, over the augmented state :=[x⊤,u⊤]⊤ α:=[x ,u ] . Let ⊂ℝnA ^n be a compact set and ~∈ α be a nominal point. The linearization error, i.e., the Lagrange remainder of the first-order Taylor expansion of p at ~ α is, (;~):=p()−p(~)−J(~)(−~),∈. r( α; α):=p( α)-p( α)-J( α)( α- α), α . (10) For the specific case when p()=f(x,u)p( α)=f(x,u), (10) coincides with (8c). Our goal is to compute a sound componentwise interval enclosure ℛ∈mR _m of (;~) r(A; α), i.e., (;~)⊆ℛ r(A; α) . IV-A Linearization Error Bounds for Analytic Dynamics By Taylor’s theorem in Lagrange form, for each component ℓ∈1,…,m ∈\1,…,m\, there exists ξ∈ξ such that ℓ(;~)=12δ⊤Hℓ(ξ)δ,δ:=−~. r_ ( α; α)= 12δ H_ (ξ)δ, δ:= α- α. (11) where Hℓ(ξ)H_ (ξ) is the ℓ -th item of the Hessian H at ξ. To evaluate the interval enclosure of ℓ(;~) r_ (A; α), previous work [7, 5, 6, 21], separately evaluate δ⊤,Hℓ(ξ)δ ,H_ (ξ) and δ in interval arithmetic, and then perform interval products as in (12), yielding the enclosure ℛℓ,classicR_ ,classic: ℓ(;~)⊆ℛℓ,classic:=12Δ⊤[Hℓ]()Δ, r_ (A; α) _ ,classic:= 12 [H_ ](A) , (12) where Δ:=⊖~⊂ℝn :=A α ^n denotes the displacement set and [Hℓ]()[H_ ](A) is an interval enclosure of the Hessian over A with entries [Hℓ,jk]()[H_ ,jk](A) (where Hℓ,jk():=∂jkpℓ()H_ ,jk( α):= _jkp_ ( α)). 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 _j denotes the j-th slice of Δ , after the symmetry of the interval Hessian is applied. When Δ is symmetric, ℛℓ,expand⊆ℛℓ,classicR_ ,expand _ ,classic is guaranteed because (13) preserves the non-negativity of δk2 _k^2, for all δk∈Δk _k∈ _k and for all k∈1,…,nk∈\1,…,n\: ℓ(;~)⊆ℛℓ,expand= r_ (A; α) \;R_ ,expand= ∑1≤j<k≤n[Hℓ,kj]()Δj⋅Δk _ subarrayc1≤ j<k≤ n subarray[H_ ,kj](A)\, _j· _k + + 12∑k=1n[Hℓ,kk]()Δk2. 12 _k=1^n[H_ ,k](A)\, _k^2. (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 ∈ α , define the sequence of intermediate points ~(k):=(1,…,k,~k+1,…,~n),k∈[n+1], α_ α^(k):=( α_1,…, α_k, α_k+1,…, α_n), k∈[n+1], (14) with ~(0)=~ α_ α^(0)= α and ~(n)= α_ α^(n)= α. This defines a coordinate-wise path [24, 25] that incrementally moves from ~ α to α along each coordinate. Correspondingly, define the prefix sets (k)(~):=∈∣ℓ=~ℓ,∀ℓ∈[k+1,n]A^(k)( α):=\ α α_ = α_ ,∀ ∈[k+1,n]\. Theorem 1 (Path-Based Hessian Bound) Let p:ℝn→ℝmp:R^n ^m be twice continuously differentiable on an interval ⊂ℝnA ^n, and let ~∈ α . Define the displacement set Δ:=⊖~ :=A α. For each ∈ α , define the coordinate-wise path (k) α^(k) and the corresponding prefix sets (k)(~)A^(k)( α). Then, for each component ℓ∈1,…,m ∈\1,…,m\, the linearization error satisfies ℓ(;~)⊆ℛℓ,path -6.0pt r_ (A; α) _ ,path :=∑1≤j<k≤n[Hℓ,kj]((j)(~))Δj⋅Δk -3.0pt:= -3.0pt _1≤ j<k≤ n[H_ ,kj](A^(j)( α))\, _j· _k +12∑k=1n[Hℓ,kk]((k)(~))Δk2. + 12 _k=1^n[H_ ,k](A^(k)( α))\, _k^2. (15) Proof: Fix ∈ α and define δ:=−~δ:= α- α. Using telescoping along the coordinate-wise path, we have pℓ()−pℓ(~)=∑k=1n(pℓ(~(k))−pℓ(~(k−1)))p_ ( α)-p_ ( α)= _k=1^n (p_ ( α_ α^(k))-p_ ( α_ α^(k-1)) ). For each k, applying Taylor’s theorem in 1D yields pℓ(~(k))−pℓ(~(k−1))=∂kpℓ(~(k−1))δk+12∂kkpℓ(ξk)δk2,p_ ( α_ α^(k))-p_ ( α_ α^(k-1))= _kp_ ( α_ α^(k-1)) _k+ 12 _kp_ ( _k)\, _k^2, (16) for some ξk∈[~(k−1),~(k)]⊆(k)(~) _k∈[ α_ α^(k-1), α_ α^(k)] ^(k)( α). Summing (16) over k and subtracting ∇pℓ(~)⊤δ∇ p_ ( α) δ yields pℓ()−pℓ(~)−∇pℓ(~)⊤δp_ ( α)-p_ ( α)-∇ p_ ( α) δ on the LHS which is ℓ(;~) r_ ( α; α) and the RHS is in ∑k=1n(∂kpℓ(~(k−1))−∂kpℓ(~))δk+12∑k=1n∂kkpℓ(ξk)δk2 _k=1^n ( _kp_ ( α_ α^(k-1))- _kp_ ( α) ) _k+ 12 _k=1^n _kp_ ( _k)\, _k^2. Expanding the gradient difference along the same path for k, we have ∂kpℓ(~(k−1))−∂kpℓ(~)=∑j=1k−1(∂kpℓ(~(j))−∂kpℓ(~(j−1))) _kp_ ( α_ α^(k-1))- _kp_ ( α)= _j=1^k-1 ( _kp_ ( α_ α^(j))- _kp_ ( α_ α^(j-1)) ). By the Mean Value Theorem, for each j<kj<k, there exists ηj,k∈[~(j−1),~(j)]⊆(j) _j,k∈[ α_ α^(j-1), α_ α^(j)] ^(j) such that ∂kpℓ(~(j))−∂kpℓ(~(j−1))=∂kjpℓ(ηj,k)δj _kp_ ( α_ α^(j))- _kp_ ( α_ α^(j-1))= _kjp_ ( _j,k)\, _j. Substituting yields ℓ(;~)=∑1≤j<k≤n∂kjpℓ(ηj,k)δjδk+12∑k=1n∂kkpℓ(ξk)δk2 r_ ( α; α)= _1≤ j<k≤ n _kjp_ ( _j,k)\, _j _k+ 12 _k=1^n _kp_ ( _k)\, _k^2. Since ηj,k∈(j)(~) _j,k ^(j)( α), ξk∈(k)(~) _k ^(k)( α), 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, ℛℓ,path⊆ℛℓ,expand⊆ℛℓ,classic.R_ ,path _ ,expand _ ,classic. (17) Proof: For each j∈1,…,nj∈\1,…,n\, the prefix set satisfies (j)⊆A^(j) by construction. Then, for every Hessian entry, [Hℓ,kj]((j))⊆[Hℓ,kj](),[Hℓ,kk]((k))⊆[Hℓ,kk]().[H_ ,kj](A^(j)) [H_ ,kj](A),[H_ ,k](A^(k)) [H_ ,k](A). Multiplying by ΔjΔk _j _k and Δk2 _k^2 preserves inclusion and summing the resulting terms yields ℛℓ,path⊆ℛℓ,expandR_ ,path _ ,expand. ∎ Unlike ℛℓ,classicR_ ,classic, ℛℓ,expandR_ ,expand and ℛℓ,pathR_ ,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). IV-B Linearization Error Bounds for Neural Dynamics For N dynamics f, computing a tight interval of the Hessian can be difficult. Instead, we use CROWN [26, 27] to obtain affine bounds of an N over a box input set, and convert them into an interval enclosure of (;~) r( α; α). Theorem 2 (CROWN bounds (Thm. 3.2 [26])) Given an N, p:ℝn→ℝmp:R^n ^m, and a bounded input ∈ α where A is an interval, there exist linear lower/upper bounds on p()p( α): ¯+¯≤p()≤¯+¯, W α+ b≤ p( α)≤ W α+ b, (18) where ¯,¯∈ℝm×n W, W ^m× n and ¯,¯∈ℝm b, b ^m. The bound in (18) does not directly bound (;~) r( α; α), since the slopes ¯,¯ W, W generally do not match the Jacobian J(~)J( α). We address this via further relaxation in Thm. 3. Theorem 3 (Linearization error bound from CROWN) Let p:ℝn→ℝmp:R^n ^m be an N, ⊂ℝnA ^n be a box set, ~∈ α , J(~)J( α) be the Jacobian of p at ~ α, and let (¯,¯,¯,¯)( W, b, W, b) satisfy the CROWN bounds (18) on A. Then, ∀∈∀ α : ¯()≤(;~)≤¯(), r( α)≤ r( α; α)≤ r( α), (19) ¯() r( α) :=(¯−J(~))+J(~)~−p(~)+¯, =( W-J( α)) α+J( α) α-p( α)+ b, (20) ¯() r( α) :=(¯−J(~))+J(~)~−p(~)+¯. =( W-J( α)) α+J( α) α-p( α)+ b. Then, an interval enclosure of (;~)⊆[¯,¯] r(A; α) [ r, r] is given by ¯i=min∈¯i(),¯i=max∈¯i(),∀i∈[1,m], r_i= _ α \; r_i( α), r_i= _ α \; r_i( α), ∀ i∈[1,m], (21) Proof: Substituting the CROWN lower bound (18) into the expression for (;~) r( α; α) (10), we yield (;~)≥¯+¯−p(~)−J(~)(−~) r( α; α)≥ W α+ b-p( α)-J( α)( α- α), which simplifies to (;~)≥(¯−J(~))+J(~)~−p(~)+¯=¯(). r( α; α)≥( W-J( α)) α+J( α) α-p( α)+ b= r( α). Similarly, we apply the CROWN upper bound. Thus, ∀∈∀ α : ¯()≤(;~)≤¯() 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(~)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. IV-C Zonotope Bounds of rk(xk,uk,wk)r_k(x_k,u_k,w_k) To obtain the zonotope bound for rk(xk,uk,wk)r_k(x_k,u_k,w_k) (9), we compute the enclosure [r¯klin,r¯klin][ r_k^lin, r_k^lin] of the linearization error rklinr_k^lin as in Sec. IV-A or IV-B and the enclosure [r¯kdist,r¯kdist][ r_k^dist, r_k^dist] of the disturbance error rkdistr_k^dist via interval arithmetic in (22), where [E](Ωkx)[E]( _k^x) is the interval enclosure of E over Ωkx _k^x: rkdist=(E(xk)−Ek)wk∈([E](Ωkx)⊖Ek)ℰnx.r_k^dist= (E(x_k)-E_k )w_k∈ ([E]( _k^x) E_k )E_n_x. (22) Finally, we obtain rk(xk,uk,wk)∈[r¯k,r¯k]r_k(x_k,u_k,w_k)∈[ r_k, r_k], where r¯k=r¯klin+r¯kdist r_k= r_k^lin+ r_k^dist and r¯k=r¯klin+r¯kdist r_k= r_k^lin+ r_k^dist. Its zonotope bound is given by (R¯kℰnx⊕c¯k) ( R_kE_n_x c_k ), where c¯k=12(r¯k+r¯k) c_k= 12( r_k+ r_k) and R¯k=12diag(r¯k−r¯k) R_k= 12diag( 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). -1.0pt R_k -1.0pt= -1.0pt _r( _k,z_k,v_k;f,E), 6.0pt c_k -1.0pt= -1.0pt _c( _k,z_k,v_k;f,E). (23) IV-D Differentiable, Parallel Bound Implementation We implement our LEBs in JAX, with the PBH using immrax for interval arithmetic and jax_verify providing CROWN-based bounds for N dynamics. JAX offers differentiability and GPU parallelization, enabling gradient-based refinement of the nominal point ~ α and scalable input domain partitioning to reduce conservativeness for SLS. V Using Linearization Error Bounds in GPUSLS We discuss the use of our LEBs in SLS (Sec. V-A) and describe an efficient GPU-based implementation (Sec. V-B). V-A Nonlinear SLS Formulation Following [9, App. A], SLS can be extended to nonlinear systems by planning a nominal trajectory for the disturbance-free dynamics zk+1=f(zk,vk)z_k+1=f(z_k,v_k) 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,, -13.0pt _ subarraycz,\,v, subarray\ \ J(,)+H~0(,) J(z,v)+ H_0( , D) (24a) s.t. zk+1=f(zk,vk),∀k∈[N],z0=x¯0, z_k+1=f(z_k,v_k), ∀ k∈[N], z_0= x_0, (24b) k+1,jx=Akk,jx+Bkk,ju,j+1,jx=Inx, ^x_k+1,j=A_k ^x_k,j+B_k ^u_k,j,\ \ ^x_j+1,j=I_n_x, (24c) ∀j∈[N],∀k∈[j+1,N−1], ∀ j∈[N], 5.0pt∀ k∈[j+1,N-1], E¯k:=r(zk,vk,k,k),∀k∈[N], E_k:=Z_r(z_k,v_k, μ_k, σ_k), ∀ k∈[N], (24d) c¯k:=c(zk,vk,k,k),∀k∈[N], c_k:=Z_c(z_k,v_k, μ_k, σ_k), ∀ k∈[N], (24e) k=hkσ(,),k=hkμ(,),∀k∈[N], σ_k=h^σ_k( , D),~ μ_k=h^μ_k( , D),\,\,∀ k∈[N], (24f) gk(zk,vk)+hk(,)≤0,∀k∈[N], g_k(z_k,v_k)+h_k( , D)≤ 0, ∀ k∈[N], -5.0pt (24g) gf(zN)+hf(,)≤0, g^f(z_N)+h^f( , D)≤ 0, (24h) where Ak,BkA_k,B_k are the linearized dynamics (7) at time k and functions r(⋅),c(⋅)Z_r(·),Z_c(·) (24d)–(24e) are the Minkowski sums of the exogenous disturbance and linearization error as a function of the nominal (zkz_k,vkv_k), reachable tube center offset k∈ℝnx+nu μ_k ^n_x+n_u, and widths k∈ℝnx+nu σ_k ^n_x+n_u respectively [12, 4], and :=(c¯j,E¯j)j=0N−1 D:=\( c_j, E_j)\_j=0^N-1. Formally, for all k∈[N]k∈[N] in (24d)–(24e), c¯k∈ℝnx c_k ^n_x and E¯k∈ℝnx×2nx E_k ^n_x× 2n_x are given as k -4.0ptA_k =([zk⊤,vk⊤]⊤+k)⊕diag(k)ℰnx+nu, = ([z_k ,v_k ] + μ_k ) diag( σ_k)E_n_x+n_u, (25a) E¯k -4.0pt E_k :=[E(zk)R¯k]=[E(zk)Γr(k,zk,vk;f,E)]⏟r(zk,vk,k,k), :=[E(z_k)\ \ R_k]= bmatrixE(z_k)\ \ _r\! (A_k,z_k,v_k;f,E ) bmatrix_Z_r(z_k,v_k, μ_k, σ_k), (25b) c¯k -4.0pt c_k :=Γc(k,zk,vk;f,E)⏟c(zk,vk,k,k), := _c\! (A_k,z_k,v_k;f,E )_Z_c(z_k,v_k, μ_k, σ_k), (25c) where c¯k c_k is the offset and E¯k E_k is the generator matrix of a zonotope for the Minkowski sum of the exogenous disturbance and linearization error: E(zk)ℰnx⊕(c¯k⊕R¯kℰnx)≡c¯k⊕E¯kℰ2nxE(z_k)E_n_x ( c_k R_kE_n_x)≡ c_k E_kE_2n_x. Using (25b) and (25c), the resulting reachable tube overapproximation is, Ω¯kx:=zk+∑j=0k−1k,jxc¯j+⨁j=0k−1k,jxE¯jℰ2nx, _k^x=z_k+ _j=0^k-1 ^x_k,j c_j+ _j=0^k-1 ^x_k,j E_jE_2n_x, (26) Ω¯ku:=vk+∑j=0k−1k,juc¯j+⨁j=0k−1k,juE¯jℰ2nx. _k^u=v_k+ _j=0^k-1 ^u_k,j c_j+ _j=0^k-1 ^u_k,j E_jE_2n_x. Robust constraint satisfaction is enforced through constraint tightenings hk(,)h_k( , D) and hf(,)h^f( , D). In short, these tightenings capture the propagation of unit-normalized disturbances in (26), which are then scaled by the disturbance bounds E¯j E_j and shifted by offset c¯j c_j. For simplicity, we derive the tightenings assuming linear constraints, i.e., gk(zk,vk)=Gk[zk,vk]⊤g_k(z_k,v_k)=G_k bmatrixz_k,\ v_k bmatrix and gf(zN)=GfzNg^f(z_N)=G^fz_N, for Gk∈ℝnc×(nx+nu),GN∈ℝnc×nxG_k ^n_c×(n_x+n_u),G_N ^n_c× n_x. Valid tightenings can be obtained for nonlinear gkg_k and gfg^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, hkμ(,) h_k^μ( , D) =∑j=0k−1k,jc¯j = _j=0^k-1 _k,j c_j (27) hkσ(,) h_k^σ( , D) =∑j=0k−1‖k,jE¯j‖1,r = _j=0^k-1 \| _k,j E_j \|_1,r (28) hk(,) h_k( , D) =∑j=0k−1(Gkk,jc¯j+‖Gkk,jE¯j‖1,r) = _j=0^k-1(G_k _k,j c_j+ \|G_k _k,j E_j \|_1,r) (29) hf(,) h^f( , D) =∑j=0N−1(GNN,jxc¯j+‖GNN,jxE¯j‖1,r) = _j=0^N-1(G_N ^x_N,j c_j+ \|G_N ^x_N,j E_j \|_1,r) -2.0pt (30) H~0(,) H_0( , D) =∑j=0N−1(∑k=jN−1(∥Q~1/2k,jxE¯j∥ℱ2 = _j=0^N-1\!( _k=j^N-1\!(\| Q^1/2 ^x_k,j E_j\|_F^2 (31) +∥R~1/2k,juE¯j∥ℱ2)+∥Q~N1/2N,jxE¯j∥ℱ2), 11.38109pt+\| R^1/2 ^u_k,j E_j\|_F^2)+\| Q_N^1/2 ^x_N,j E_j\|_F^2)\,, where k,j:=[k,jx⊤k,ju⊤]⊤ _k,j:= bmatrix ^x_k,j & -5.69054pt ^u_k,j bmatrix , collects k,j _k,j for all k,j∈[N]k,j∈[N], and Q~,Q~N∈++nx,R~∈++nu Q, Q_N _++^n_x, R _++^n_u. The constraint tightenings and tube widths in (28)-(30) are derived for the worst case disturbances using the dual norm property [28]. The cost H~0(,) 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, -5.0pt _ subarraycz,\,v subarray J(,) J(z,v) (32a) s.t. zk+1=f(zk,vk),∀k∈[N],z0=x¯0, z_k+1=f(z_k,v_k), ∀ k∈[N], z_0= x_0, (32b) gk(zk,vk)+hk(,)≤0,∀k∈[N], g_k(z_k,v_k)+h_k( , D)≤ 0, ∀ k∈[N], (32c) gf(zN)+hf(,)≤0. g^f(z_N)+h^f( , D)≤ 0. (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): minx,u -3.0pt _ ^x, ^u\ ∑j=0N−1∑k=jN−1‖k,jk,jE¯j‖ℱ2+‖N,jN,jxE¯j‖ℱ2 _j=0^N-1 _k=j^N-1\|Q_k,j _k,j E_j\|^2_F+\|Q_N,j ^x_N,j E_j\|^2_F (33a) s.t. k+1,jx=Akk,jx+Bkk,ju, ^x_k+1,j=A_k ^x_k,j+B_k ^u_k,j, (33b) j+1,jx=Inx,∀j∈[N],∀k∈[j+1,N−1], ^x_j+1,j=I_n_x,~∀ j∈[N],~∀ k∈[j+1,N-1], -4.0pt (33c) where k,j,N,jQ_k,j,Q_N,j are defined according to [9, Eq. 20, 23] (see App. E) for consistency between (32) and (33). Notably, the right-multiplication by E¯j E_j prevents the cost in (33a) from being expressed as a standard LQR objective in x ^x and u ^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]: minx,u -8.0pt _ ^x, ^u\ ∑j=0N−1∑k=jN−1‖k,jk,j‖ℱ2+‖N,jN,jx‖ℱ2 _j=0^N-1 -3.0pt _k=j^N-1\|Q_k,j _k,j\|^2_F+\|Q_N,j ^x_N,j\|^2_F (34a) s.t. k+1,jx=Akk,jx+Bkk,ju, ^x_k+1,j=A_k ^x_k,j+B_k ^u_k,j, (34b) j+1,jx=Inx,∀j∈[N],∀k∈[j+1,N−1]. ^x_j+1,j=I_n_x,~∀ j∈[N],~∀ k∈[j+1,N-1]. (34c) Note that (34) coincides with (33) if E¯j=I E_j=I for all j∈[N]j∈[N]. In Prop. 1, we prove that the set of minimizers of (33) and (34) are equivalent if E¯j E_j is right-invertible for all j∈[N]j∈[N] (see App. C for the proof). Proposition 1 Let (~x,~u)( ^x, ^u) be an optimizer of (33), and let (^x,^u)( ^x, ^u) be an optimizer of (34). If E¯j E_j is right-invertible for all j∈[N]j∈[N], i.e., there exists E¯j‡ E_j such that E¯jE¯j‡=Inx E_j E_j =I_n_x, then the two optimizers coincide: ~x=^x,~u=^u ^x= ^x, ^u= ^u. When implementing the method of Sec. IV, we add a small positive diagonal padding to R¯k R_k so that it is always invertible, and thus E¯j E_j (25) is always right-invertible for all j∈[N]j∈[N]. After solving (34), we compute the constraint tightening terms in (27)–(30). This updates the reachable tube centers k μ_k and widths k σ_k: k←hkσ(,)andk←hkμ(,), σ_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 number of iterations or if the change in the iterates is below a threshold (see App. B 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 E¯j E_j incompatible with existing solvers [9, 10]. We address this by setting j+1,jx=Inx ^x_j+1,j=I_n_x 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. V-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 ^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)x_k+1=f(x_k,u_k) around the current (z,v)(z,v) to obtain Ak,BkA_k,B_k (4). It then solves (A) a local QP approximation of (32) to update (z,v)(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 objective (32a) to J¯(z,v):=J(z,v)+λJerr(z,v) J(z,v):=J(z,v)+λ J_ err(z,v), where Jerr(z,v)=∑j=0N−1∑k=jN−1‖k,jr(zj,vj,j,j)‖ℱ2, -2.0ptJ_ err(z,v)= _j=0^N-1 _k=j^N-1\| _k,jZ_r(z_j,v_j, μ_j, σ_j)\|_F^2, (35) penalizes the propagated residual error (8a). Here, ,j,j , μ_j, σ_j are from the previous iteration and λ is a weighting term. Notably, (35) biases (z,v)(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∈Ω¯kxx_k∈ _k^x and uk∈Ω¯kuu_k∈ _k^u, for all k∈[N]k∈[N]. Formally, we have (proof in App. D): 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 _k _k for all k∈[N]k∈[N]. VI Results 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. A. VI-A Linearization Error We compare the LEB from our path-based Hessian (PBH) method on the Satellite (nx=7n_x=7, nu=3n_u=3) and Quadrotor (nx=12n_x=12, nu=4n_u=4) systems (App. A) 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)(z,v) from ×X×U and define local sets :=(z,v)⊕ϵℰnx+nuT:=(z,v) _n_x+n_u with ϵ=0.01,0.05,0.1,0.15,0.2ε=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 random-sampling 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ε=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. Figure 2: (a): LEBs on the satellite system. (b): LEBs on the quadrotor system. In both systems, our method achieves the tightest over-approximation. TABLE I: Runtime (in ms) comparison of linearization bound methods. System Random Global IH-classic IH-classic-CORA Path-based Satellite 0.894 0.143 0.218 11570 0.361 Quadrotor 1.278 0.156 0.453 227200 0.936 TABLE I: Linearization error under nominal optimization via PGD. System Initial Sampling PGD PGD better (%) Satellite 0.2680 0.2365 0.2241 100 Quadrotor 0.02253 0.01639 0.01590 100 To quantify the benefit of differentiability in our path-based 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 gradient descent (PGD) within a local neighborhood of radius ϵopt _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 I shows that across 100% of test cases, PGD yields smaller error interval widths than random sampling across both systems, showing the value of gradient-based optimization for finding nominal points with tighter LEBs, which can benefit downstream SLS. VI-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. Figure 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. VI-B1 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 σ_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 overly-conservative 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. Figure 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. Figure 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. VI-B2 Planar quadrotor We evaluate our method on a 6D planar quadrotor against control contraction metrics (CCM) [30]. Both methods must navigate through a dense 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)(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, trajectory-dependent 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×1051.95× 10^5 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. Figure 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 pyp_y and vyv_y for 4 quadrotors. 100% of rollouts stay safely in the tubes. VI-B3 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, demonstrating the method’s scalability to state dimensionality. In Fig. 6(b), we plot the tubes and the corresponding rollouts for pyp_y and vyv_y for 100 adversarial disturbances. Figure 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. VI-B4 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λ=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λ=1 decreasing tube widths by 8%. We also show the importance of nonzero-centered 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 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 I: Per iteration runtime breakdown (ms) across systems. System GPUSLS Lin. Error Grad. Remainder Total DT Satellite 11 3 <1<1 15 1000 Planar Quadrotor 18 1 <1<1 20 150 Coupled Quadrotor (96D) 40 – 48 88 100 Neural T-Dynamics 630 100 32 762 1000 Figure 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. VI-B5 Runtime We evaluate the runtime of our method across all systems. In Tab. I, 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 Conclusion This paper presents a framework for computing tight, differentiable, and sound linearization error bounds (LEBs) for nonlinear analytic and N 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≈ 2× 10^5 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, 34, 35]. References [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, p. 4042–4048. [8] J. Anderson, J. C. Doyle, S. H. Low, and N. Matni, “System level synthesis,” Annu. Rev. Control., vol. 47, p. 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, p. 3747–3752. [21] M. Rungger and M. Zamani, “Accurate reachability analysis of uncertain nonlinear systems,” in HSCC, 2018, p. 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, p. 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. In the appendix, we first discuss details on the system dynamics used in the experimental results (App. A). We then provide an algorithm block describing how sequential quadratic programming is used to solve GPUSLS-LEO (App. B). Next, we discuss and prove Proposition 1 and a supporting lemma in App. C. We then provide proofs for the remaining theoretical results that were omitted from the main manuscript in App. D. Finally, we discuss in more detail the Lagrange multiplier variable τ used in our framework in App. E. Appendix A Dynamical System Definitions In this section, we give an overview of the discrete-time dynamical system definitions for each system. The discrete-time dynamics are obtained via forward Euler discretization with timestep Δt t, xk+1=xk+Δt(x˙k)x_k+1=x_k+ t( x_k). A-A Satellite We use the same satellite dynamics and experimental setup as [4]. Specifically, we define the dynamics as, x˙ x =[ζ(ω)qIS−1(v−w×(ISω))] = bmatrixζ(ω)q\\ I_S^-1(v-w×(I_Sω)) bmatrix (36a) ζ(ω) ζ(ω) :=12[0−ω1−ω2−ω3ω10ω3−ω2ω2−ω30ω1ω3ω2−ω10] := 12 bmatrix0&- _1&- _2&- _3\\ _1&0& _3&- _2\\ _2&- _3&0& _1\\ _3& _2&- _1&0 bmatrix (36b) with states x:=(q,ω)∈ℝ7x:=(q,ω) ^7, where q∈ℝ4q ^4 is the attitude quaternion, ω∈ℝ3ω ^3 the angular rotation rate, v∈ℝ3v ^3 the input control torque, and IS=diag(5,2,1)I_S=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)E=5· 10^-3·diag(0,0,0,0,1,1,1). A-B Planar Quadrotor We consider a planar quadrotor with the following dynamics x˙ x =[vxvyϕ˙−1m(u1+u2)sin(ϕ)1m(u1+u2)cos(ϕ)−gLJ(u2−u1)] = bmatrixv_x\\ v_y\\ φ\\ - 1m(u_1+u_2) (φ)\\ 1m(u_1+u_2) (φ)-g\\ LJ(u_2-u_1) bmatrix (37a) The state is defined as x:=(px,py,ϕ,vx,vy,ϕ˙)∈ℝ6x:=(p_x,p_y,φ,v_x,v_y, φ) ^6 and input u:=(u1,u2)∈ℝ2u:=(u_1,u_2) ^2. The position is denoted as (px,py)(p_x,p_y), ϕφ the pitch angle, (vx,vy)(v_x,v_y) the translational velocities, and ϕ˙ φ the angular velocity. The inputs (u1,u2)(u_1,u_2) are the individual rotor thrusts. We define the mass m=2.0576m=2.0576, gravitational acceleration g=9.81g=9.81, arm length L=0.25L=0.25 and moment of inertia J=0.01J=0.01. We define the disturbance scaling as E=5⋅10−2⋅diag(0,0,0,1,1,0)E=5· 10^-2·diag(0,0,0,1,1,0). A-C Coupled 3D Quadrotor We consider a system of N quadrotors. For a given quadrotor i, we define its dynamics as ξ˙i ξ_i =[p˙x,ip˙y,ip˙z,iϕ˙iθ˙iψ˙iv˙x,iv˙y,iv˙z,ip˙iq˙ir˙i]=[vx,ivy,ivz,ipi+(qisinϕi+ricosϕi)tanθiqicosϕi−risinϕiqisinϕi+ricosϕicosθiTim(cosψisinθicosϕi+sinψisinϕi)Tim(sinψisinθicosϕi−cosψisinϕi)Tim(cosθicosϕi)−gJy−JzJxqiri+τϕ,iJxJz−JxJypiri+τθ,iJyJx−JyJzpiqi+τψ,iJz]. = bmatrix p_x,i\\ p_y,i\\ p_z,i\\ φ_i\\ θ_i\\ ψ_i\\ v_x,i\\ v_y,i\\ v_z,i\\ p_i\\ q_i\\ r_i bmatrix= bmatrixv_x,i\\ v_y,i\\ v_z,i\\ p_i+(q_i _i+r_i _i) _i\\ q_i _i-r_i _i\\ q_i _i+r_i _i _i\\ T_im ( _i _i _i+ _i _i )\\ T_im ( _i _i _i- _i _i )\\ T_im ( _i _i )-g\\ J_y-J_zJ_xq_ir_i+ _φ,iJ_x\\ J_z-J_xJ_yp_ir_i+ _θ,iJ_y\\ J_x-J_yJ_zp_iq_i+ _ψ,iJ_z bmatrix. (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)∈ℝ12 _i:=(p_x,i,p_y,i,p_z,i, _i, _i, _i,v_x,i,v_y,i,v_z,i,p_i,q_i,r_i) ^12 where (px,i,py,i,pz,i)(p_x,i,p_y,i,p_z,i) denotes position, (ϕi,θi,ψi)( _i, _i, _i) the roll, pitch, and yaw angles, (vx,i,vy,i,vz,i)(v_x,i,v_y,i,v_z,i) the translational velocities, and (pi,qi,ri)(p_i,q_i,r_i) denote body angular velocities. The input is defined as υi:=(Ti,τϕ,i,τθ,i,τψ,i)∈ℝ4 _i:=(T_i, _φ,i, _θ,i, _ψ,i) ^4, where TiT_i is the collective thrust and (τϕ,i,τθ,i,τψ,i)( _φ,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)∈ℝ12N,u:=(υ1,…,υN)∈ℝ4N.x:=( _1,…, _N) ^12N, u:=( _1,…, _N) ^4N. (39) We additionally incorporate a centroid spring-damper coupling in the translational dynamics. We define the centroid position and velocity as p¯k=1N∑i=1Npi,k,v¯k=1N∑i=1Nvi,k. p_k= 1N _i=1^Np_i,k, v_k= 1N _i=1^Nv_i,k. (40) The coupling acceleration applied to quadrotor i is ai,kcoup=−αN(pi,k−p¯k)−γ(vi,k−v¯k)a_i,k^coup=-α N(p_i,k- p_k)-γ(v_i,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α=0.75 and γ=0.25γ=0.25. Therefore, the discrete-time dynamics become ξi,k+1=ξi,k+Δt[p˙x,i,kp˙y,i,kp˙z,i,kϕ˙i,kθ˙i,kψ˙i,kv˙i,k+ai,kcoupp˙i,kq˙i,kr˙i,k]. _i,k+1= _i,k+ t bmatrix p_x,i,k\\ p_y,i,k\\ p_z,i,k\\ φ_i,k\\ θ_i,k\\ ψ_i,k\\ v_i,k+a^coup_i,k\\ p_i,k\\ q_i,k\\ r_i,k bmatrix. (42) We use a disturbance scaling matrix of E=2⋅10−2⋅I12NE=2· 10^-2· I_12N Appendix B Algorithm Block Algorithm 1 GPUSLS-LEO 1:x¯0 x_0, J, f(x,u)f(x,u), E(x)E(x), gk(x,u)g_k(x,u), gf(x)g^f(x) 2:while SQP not converged do 3: A, B, G ←SQP: Linearize System← SQP: Linearize System 4: ,,,,←SQP: Quadraticize CostQ,M,R,q,r← SQP: Quadraticize Cost 5: while SLS not converged do 6: z, v ←Optimize Nominal Trajectory← Optimize Nominal Trajectory 7: ,←Compute interval σ, μ← Compute interval 8: ¯ E, ¯ c ←Compute residual error zonotope← Compute residual error zonotope 9: ,←Update duals and cost τ, Q← Update duals and cost 10: x,u←Optimize Controller x, u← Optimize Controller 11: h(,),hf(,)←Update Tighteningsh( , D),h^f( , D)← Update Tightenings 12: end while 13:end while 14:return ,,x,uz,v, ^x, ^u Here, ,,,,Q,M,R,q,r correspond to the quadratic approximation of the nominal cost. Appendix C On Solving for Solving (33) yields optimal values for x ^x and u ^u which penalize the size of the tightenings (31). However, this leads to terms in (33a) that scale the response matrices x ^x and u ^u by E¯j E_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: minx,u -3.0pt _ ^x, ^u\ ∑j=0N−1∑k=jN−1‖k,jk,jE¯j‖ℱ2+‖N,jN,jxE¯j‖ℱ2 _j=0^N-1 _k=j^N-1\|Q_k,j _k,j E_j\|^2_F+\|Q_N,j ^x_N,j E_j\|^2_F (33a) s.t. k+1,jx=Akk,jx+Bkk,ju, ^x_k+1,j=A_k ^x_k,j+B_k ^u_k,j, (33b) j+1,jx=Inx,∀j∈[N],∀k∈[j+1,N−1], ^x_j+1,j=I_n_x,~∀ j∈[N],~∀ k∈[j+1,N-1], -4.0pt (33c) Reformulated Problem: minx,u -5.0pt _ ^x, ^u\ ∑j=0N−1∑k=jN−1‖k,jk,j‖ℱ2+‖N,jN,jx‖ℱ2 _j=0^N-1 -3.0pt _k=j^N-1\|Q_k,j _k,j\|^2_F+\|Q_N,j ^x_N,j\|^2_F (34a) s.t. k+1,jx=Akk,jx+Bkk,ju, ^x_k+1,j=A_k ^x_k,j+B_k ^u_k,j, (34b) j+1,jx=Inx,∀j∈[N],∀k∈[j+1,N−1] ^x_j+1,j=I_n_x,~∀ 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)( ^x, ^u) be an optimizer of (33), and let (^x,^u)( ^x, ^u) be an optimizer of the reformulated problem (34). Suppose that E¯j E_j is right-invertible for all j∈[N]j∈[N], i.e., there exists E¯j‡ E_j such that E¯jE¯j‡=Inx E_j E_j =I_n_x. Then the two optimizers coincide: ~x=^x,~u=^u ^x= ^x, ^u= ^u. Proof: Observe that (33) can be rewritten into (43) with the variable change k,j=k,jE¯j _k,j= _k,j E_j, since by using (33b) and (33c) we have that, k+1,jxE¯j=Akk,jxE¯j+Bkk,juE¯j ^x_k+1,j E_j=A_k ^x_k,j E_j+B_k ^u_k,j E_j ⇔Eq. (43b), . eq:fastsls_qp_x_prop_var, and j+1,jxE¯j=InxE¯j ^x_j+1,j E_j=I_n_x E_j ⇔Eq. (43c). . eq:fastsls_qp_x_init_var. minx,u -3.0pt _ ^x, ^u\ ∑j=0N−1∑k=jN−1‖k,jk,j‖ℱ2+‖N,jN,jx‖ℱ2 _j=0^N-1 -3.0pt _k=j^N-1\|Q_k,j _k,j\|^2_F+\|Q_N,j ^x_N,j\|^2_F (43a) s.t. k+1,jx=Akk,jx+Bkk,ju, ^x_k+1,j=A_k ^x_k,j+B_k ^u_k,j, (43b) j+1,jx=E¯j,∀j∈[N],∀k∈[j+1,N−1] ^x_j+1,j= E_j,~∀ j∈[N],~∀ k∈[j+1,N-1] (43c) From [10], (43) can be solved via the parallel Riccati recursions, N,j _N,j =N,jx,k,j=−k,jℬk,j, =Q^x_N,j, _k,j=-G_k,jB_k,j, (44) k,j _k,j =k,jx+(Ak)⊤k+1,jAk+(k,j)⊤ℬk,j, =Q^x_k,j+ (A_k ) P_k+1,jA_k+ (K_k,j ) B_k,j, k,j _k,j =(k,ju+(Bk)⊤k+1,jBk)−1, = (Q^u_k,j+ (B_k ) P_k+1,jB_k )^-1, ℬk,j _k,j =k,jux+(Bk)⊤k+1,jAk, =Q^ux_k,j+ (B_k ) P_k+1,jA_k, and parallel forward propagations j+1,jx=E¯j,k,ju=k,jk,jx, ^x_j+1,j= E_j, ^u_k,j=K_k,j ^x_k,j, (45) k+1,jx=(Ak+Bkk,j)k,jx. ^x_k+1,j=(A_k+B_kK_k,j) ^x_k,j. Then, since E¯j E_j is right-invertible and (43) represents the system level constraint (SLC) assumed in Lemma 1, the controller gains obey, j+1,jx=E¯j,k,juE¯j‡=k,jk,jxE¯j‡k+1,jxE¯j‡=(Ak+Bkk,j)k,jxE¯j‡. -6.0pt aligned & ^x_j+1,j= E_j,~~~~~~~~~ ^u_k,j E_j =K_k,j ^x_k,j E_j \\ & ^x_k+1,j E_j = (A_k+B_kK_k,j ) ^x_k,j E_j .\\ aligned (46) Then, by unraveling the recursion (46), we have for k>jk>j, k,jE¯j‡=[Inxk,j](∏m=j+1k−1(Am+Bmm,j))j+1,jxE¯j‡, _k,j E_j = bmatrixI_n_x\\ K_k,j bmatrix ( _m=j+1^k-1 (A_m+B_mK_m,j ) ) ^x_j+1,j E_j , (47) while the remaining gains (Ψk,j _k,j where k≤jk≤ j) are zero, where we define the ordered matrix product ∏i=1nRi:=RnRn−1⋯R1 _i=1^nR_i:=R_nR_n-1·s R_1. To recover the optimal gains for (33), we use the right-invertibility property of E¯j E_j and substitute k,jE¯j‡=~k,j _k,j E _j= _k,j into (47): ~k,j _k,j =[Inxk,j](∏m=j+1k−1(Am+Bmm,j))~j+1,jx = bmatrixI_n_x\\ K_k,j bmatrix ( _m=j+1^k-1 (A_m+B_mK_m,j ) ) ^x_j+1,j =[Inxk,j](∏m=j+1k−1(Am+Bmm,j))Inx. = bmatrixI_n_x\\ K_k,j bmatrix ( _m=j+1^k-1 (A_m+B_mK_m,j ) )I_n_x. (48) The second equality stems from ~j+1,jx ^x_j+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‡E_j is replaced with InxI_n_x, : ^k,j=[Inxk,j](∏m=j+1k−1(Am+Bmm,j))^j+1,jx, _k,j= bmatrixI_n_x\\ K_k,j bmatrix ( _m=j+1^k-1 (A_m+B_mK_m,j ) ) ^x_j+1,j, (49) where ^j+1,jx=Inx ^x_j+1,j=I_n_x, as constrained by (34c). Hence ∀k,j∈[N]~k,j=^k,j∀\,k,j∈[N]~ _k,j= _k,j, so ~= = . ∎ Lemma 1 Let, ¯:=blk_diag(Inx,E¯0,…,E¯n−1) E:= blk\_diag(I_n_x, E_0,…, E_n-1), where for all i∈[N],E¯i∈ℝnx×nw(i)i∈[N], E_i ^n_x× n_w^(i) is right-invertible, and nw(i)≥nxn_w^(i)≥ n_x. Enforcing the uncertain LTV dynamics xk+1=Akxk+Bkuk+E¯kwk,wk∈ℰnw(k),x_k+1=A_kx_k+B_ku_k+ E_kw_k,~w_k _n_w^(k), (50) is equivalent to enforcing the system level constraint (51), i.e., (50)⇔ (51), [(INnx−)−][xu]=¯. bmatrix(I_Nn_x- Z A)&- Z B bmatrix bmatrix ^x\\ ^u bmatrix= E. (51) Furthermore, the controller gain mapping for the SLC in (51) is, =(u‡)(x¯‡)−1=(u¯⊤)(x¯⊤)−1,K=( ^u E )( ^x E )^-1=( ^u E )( ^x E )^-1, (52) where ¯‡=blk_diag(I,E¯0‡,…,E¯n−1‡) E = blk\_diag(I, E_0 ,…, E_n-1 ). Lastly, given K, the forward propagation of disturbance gains obeys: j+1,jx=E¯j,k,juE¯j‡=k,jk,jxE¯j‡k+1,jxE¯j‡=(Ak+Bkk,j)k,jxE¯j‡. -6.0pt aligned & ^x_j+1,j= E_j,~~~~~~~~~ ^u_k,j E_j =K_k,j ^x_k,j E_j \\ & ^x_k+1,j E_j = (A_k+B_kK_k,j ) ^x_k,j E_j .\\ aligned (53) Proof: First, let δk=Ekwk _k=E_kw_k be such that xk+1=Akxk+Bkuk+δk.x_k+1=A_kx_k+B_ku_k+ _k. (54) Now, let =[x0⊤,…,xn−1⊤]⊤,=[u0⊤,…,un−1⊤]⊤,and=[δ0⊤,…,δn−1⊤]⊤ x=[x_0 ,…,x_n-1 ] ,~ u=[u_0 ,…,u_n-1 ] ,~and~ δ=[ _0 ,…, _n-1 ] . Then, defining the gain mappings ^x ^x and ^u ^u such that, []=[^x^u], bmatrix x\\ u bmatrix= bmatrix ^x\\ ^u bmatrix δ, (55) we have from [8]: [(INnx−)−][^x^u]=INnx bmatrix(I_Nn_x- Z A)&- Z B bmatrix bmatrix ^x\\ ^u bmatrix=I_Nn_x (56) Defining =[w0⊤,…,wn−1⊤]⊤ w=[w_0 ,…,w_n-1 ] , we observe: =¯ δ= E w (57) Then, the gain mapping from w to x and u via x ^x and u ^u follows from (55) and (57): []=[xu]=[^x¯^u¯]=[^x^u]¯. bmatrix x\\ u bmatrix= bmatrix ^x\\ ^u bmatrix w= bmatrix ^x E\\ ^u E bmatrix w= bmatrix ^x\\ ^u bmatrix E w. (58) This establishes that x=^x¯ ^x= ^x E and u=^u¯ ^u= ^u E. Hence, by multiplying (56) by ¯ E, we arrive at (51). Thus we have shown (51)⟹(50) eq:slc eq:ltv_proof. To show the equivalence between (50) and (51), we can show (50)⟹(51) eq:ltv_proof eq:slc, by directly recovering the matrices Ak,Bk,E¯kA_k,B_k, E_k ∀k∈[N]∀ k∈[N] from (51) to produce (50). To prove (52), we observe that, K=^u(^x)−1K= ^u( ^x)^-1, [8]. Then using the substitutions x¯‡=^x ^x E = ^x and u¯‡=^u ^u E = ^u from (58) we arrive at the first equality in (52). The second equality arises from algebraic manipulation using the definition ¯‡=¯⊤(¯¯⊤)−1 E = E ( E E )^-1. Lastly, manipulating (52) to see (x¯‡)=u¯‡K( ^x E )= ^u E , and expanding the system level constraint in (51) we have that j+1,jx=E¯j,k,juE¯j‡=k,jk,jxE¯j‡k+1,jxE¯j‡=(Ak+Bkk,j)k,jxE¯j‡. -6.0pt aligned & ^x_j+1,j= E_j,~~~~~~~~~ ^u_k,j E_j =K_k,j ^x_k,j E_j \\ & ^x_k+1,j E_j = (A_k+B_kK_k,j ) ^x_k,j E_j .\\ aligned (59) ∎ Appendix D Proofs 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 _k _k for all k∈[N]k∈[N]. Proof: We will first prove (26) is correct assuming valid disturbances. Let kW_k specify the disturbance ∀k∈[N]∀ k∈[N]. Then, from SLS [4], we have that the overapproximation of the reachable set is, Ω¯kx:=zk+⨁j=0k−1k,jxj, _k^x=z_k+ _j=0^k-1 ^x_k,jW_j, (60) Ω¯ku:=vk+⨁j=0k−1k,juj. _k^u=v_k+ _j=0^k-1 ^u_k,jW_j. Now, suppose kW_k can be expressed as a zonotope, c¯k⊕E¯kℰ2nx c_k E_kE_2n_x. Making this substitution in (60) we arrive at, Ω¯kx:=zk+⨁j=0k−1k,jx(c¯k⊕E¯kℰ2nx), _k^x=z_k+ _j=0^k-1 ^x_k,j( c_k E_kE_2n_x), (61) Ω¯ku:=vk+⨁j=0k−1k,ju(c¯k⊕E¯kℰ2nx). _k^u=v_k+ _j=0^k-1 ^u_k,j( c_k E_kE_2n_x). Using the associativity and commutativity of the Minkowski sum (61) can be rewritten into, Ω¯kx:=zk+⨁j=0k−1k,jxc¯j⊕⨁j=0k−1k,jxE¯jℰ2nx, _k^x=z_k+ _j=0^k-1 ^x_k,j c_j _j=0^k-1 ^x_k,j E_jE_2n_x, (62) Ω¯ku:=vk+⨁j=0k−1k,juc¯j⊕⨁j=0k−1k,juE¯jℰ2nx. _k^u=v_k+ _j=0^k-1 ^u_k,j c_j _j=0^k-1 ^u_k,j E_jE_2n_x. Then recalling that ⨁j=0k−1k,jxc¯j _j=0^k-1 ^x_k,j c_j is a Minkowski sum of singleton sets, we can replace it with ∑j=0k−1k,jxc¯j _j=0^k-1 ^x_k,j c_j (and similarly do this for ⨁j=0k−1k,juc¯j _j=0^k-1 ^u_k,j c_j), and arrive at at the disturbance reachable set (26). Now, we will show the disturbance bounds, r(zj,vj,j,j)Z_r(z_j,v_j, μ_j, σ_j), c(zj,vj,j,j)Z_c(z_j,v_j, μ_j, σ_j), in (24d)-(24e) are valid overapproximations ℰ¯k E_k (9) and certify that the closed-loop reachable sets Ω¯kx ^x_k and Ω¯ku ^u_k overapproximate the true reachable set Ωkx ^x_k and Ωku ^u_k (3). Base Case: The disturbance at timestep k=0k=0 is precisely E(z0)=E(x¯0)E(z_0)=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]j∈[k], dj(xj,uj,wj)∈Ejℰnx⊕R¯jℰnx⊕c¯jd_j(x_j,u_j,w_j)∈ E_jE_n_x R_jE_n_x c_j. Then, we have a overapproximate reachable sets, Ωkx⊆Ω¯kx _k^x ^x_k and Ωku⊆Ω¯ku _k^u ^u_k, at timestep k. Define Ω¯k:=[x⊤,u⊤]⊤∣x∈Ω¯kx,u∈Ω¯ku _k:=\[x ,u ] x∈ ^x_k,u∈ ^u_k\. Then observe that Ω¯k⊆k:=([zk⊤,vk⊤]⊤+k)⊕diag(k)ℰnx+nu _k _k:= ([z_k ,v_k ] + μ_k ) diag( σ_k)E_n_x+n_u, as defined in (25a). Then, by construction, for all [xk⊤,uk⊤]⊤∈Ω¯k[x_k ,u_k ] ∈ _k and for all wk∈ℰnxw_k _n_x, we have that r(xk,uk,wk)∈c¯k⊕R¯kℰnxr(x_k,u_k,w_k)∈ c_k R_kE_n_x, where (c¯k,R¯k)( c_k, R_k) is computed using kA_k in (25). Similarly, the nominal disturbance Ekwk∈E(zk)ℰnxE_kw_k∈ E(z_k)E_n_x. Thus, the total disturbance (8a), rklin+rkdist+Ekwk=d(xk,uk,wk)∈c¯k⊕R¯kℰnx⊕E(zk)ℰnxr_k lin+r_k^dist+E_kw_k=d(x_k,u_k,w_k)∈ c_k R_kE_n_x E(z_k)E_n_x. 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)d_k(x_k,u_k,w_k) (8a). Having established valid disturbance bounds for all steps k∈[N]k∈[N], SLS guarantees that the overapproximation of the closed-loop reachable set (26), Ω¯kx ^x_k and Ω¯ku ^u_k, are valid, i.e., Ωkx⊂Ω¯kx ^x_k⊂ ^x_k and Ωku⊂Ω¯ku ^u_k⊂ ^u_k for all k∈[N]k∈[N]. ∎ Appendix E GPUSLS Duals τ 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∈ℝnc _k,j ^n_c and βN,j∈ℝnf _N,j ^n_f, where βk,j _k,j =[Gkk,jc¯j+‖Gkk,jE¯j‖1,r]∘2 = [G_k _k,j c_j+\|G_k _k,j E_j\|_1,r ] 2 (63a) ∀j∈[N],∀k∈[j,N] ∀ j∈[N], ∀ k∈[j,N] (63b) βN,j _N,j =[GNN,jxc¯j+‖GNN,jxE¯j‖1,r]∘2 = [G_N ^x_N,j c_j+\|G_N ^x_N,j E_j\|_1,r ] 2 (63c) ∀j∈[N] ∀ j∈[N] (63d) where [⋅]∘2[·] 2 denotes element-wise square. We can then calculate our dual variables τ by τk,j=λ¯kβk,j+ϵ∀j∈[N],∀k∈[j,N] _k,j= λ_k _k,j+ε ∀ j∈[N], ∀ k∈[j,N] (64) where we denote λ¯k λ_k as the dual variable of the nominal trajectory optimization (corresponding to the dual variable λk _k in [10]), introduced here to avoid notation conflict. Finally, we can define our cost terms as k,j=(diag(τk,j)Gk,[Q~1/200R~1/2]), _k,j= (diag ( _k,j )G_k, bmatrix Q^1/2&0\\ 0& R^1/2 bmatrix ), (65) N,j=(diag(τN,j)GN,Q~N1/2). _N,j= (diag ( _N,j )G_N, Q_N^1/2 ). Following [9] and [10], this formulation preserves dual consistency and yields an equivalent reformulation of the RNOCP under the proposed decomposition.