We present GPU-SLS, a GPU-parallelized framework for safe, robust nonlinear model predictive control (MPC) that scales to high-dimensional uncertain robotic systems and long planning horizons. Our method jointly optimizes an inequality-constrained, dynamically-feasible nominal trajectory, a tracking controller, and a closed-loop reachable set under disturbance, all in real-time. To efficiently compute nominal trajectories, we develop a sequential quadratic programming procedure with a novel GPU-accelerated quadratic program (QP) solver that uses parallel associative scans and adaptive caching within an alternating direction method of multipliers (ADMM) framework. The same GPU QP backend is used to optimize robust tracking controllers and closed-loop reachable sets via system level synthesis (SLS), enabling reachability-constrained control in both fixed- and receding-horizon settings. We achieve substantial performance gains, reducing nominal trajectory solve times by 97.7% relative to state-of-the-art CPU solvers and 71.8% compared to GPU solvers, while accelerating SLS-based control and reachability by 237x. Despite large problem scales, our method achieves 100% empirical safety, unlike high-dimensional learning-based reachability baselines. We validate our approach on complex nonlinear systems, including whole-body quadrupeds (61D) and humanoids (75D), synthesizing robust control policies online on the GPU in 20 milliseconds on average and scaling to problems with 2 x 10^5 decision variables and 8 x 10^4 constraints. The implementation of our method is available at https://github.com/Jeff300fang/gpu_sls.
Figures & tables
Fig. 1: (a) : GPU-accelerated nonlinear constrained whole-body control (61 states, 12 controls) executed on hardware with a Unitree Go2 EDU quadruped, navigating through an obstacle field in real-time at 50 Hz. (b) : Overhead view of simulated positions reconstructed from hardware-experiment data showing successful navigation around the obstacles. (c) : Graph illustrating minimum distance to any obstacle, demonstrating zero obstacle collisions.
Fig. 2: Schematic of our method, GPU-SLS. At each SQP iteration (s) , we compute a nominal trajectory update (δx,δu) using a GPU-parallelized ADMM-based QP solver that exploits caching and parallel associative scans (PASs) for acceleration. The resulting dual variables are post-processed into SLS-compatible duals τ (App. B ), which inform the objective in the SLS controller optimization problem ( 18a ). This optimization is likewise GPU-parallelized via PASs, yielding closed-loop response matrices Φ that implicitly define the controller. These matrices define constraint tightenings ( 18e ) and ( 18f ), i.e., closed-loop reachable sets, which are used to tighten the constraints of the subsequent nominal trajectory solve, ensuring robust feasibility under closed-loop control. The next SQP iteration then proceeds after fully parallelizable Jacobian and Hessian updates, along with state, input, and function evaluations. Green denotes nominal trajectory computations; blue denotes feedback controller and reachable set computations. GEMM is a O(logN) matrix multiplication routine.
Fig. 3: (a) : Comparison of average solve-time scaling with horizon length for various NMPC solvers on a torque-constrained 10-link pendulum stabilization task under random external disturbances. Our method consistently outperforms all baselines for long-horizon problems. (b) : Solve-time scaling of GPU-SLS with a family of n -link torque-constrained pendulums.
Fig. 4: DeepReach rollouts (a) and GPU-SLS rollouts (b) for an obstacle avoidance Dubins car system subject to both adversarial and random disturbance. Deepreach is shown to fail to consistently certify safety, while our method certifies safety 100% of all rollouts.
Fig. 5: Runtime comparison of FastSLS vs GPU-SLS for solving the system in 4 . GPU-SLS consistently outperforms FastSLS due to its logarithmic scaling vs FastSLS’s cubic scaling.
Fig. 6: Dubins car experiment using GPU-SLS to calculate a robust trajectory for 20 meters through a dense constraint field of 30 obstacles. This demonstrates our solver’s ability to handle large-scale trajectory optimization and robust constraint satisfaction at scale.
Fig. 7: Rollouts for uncertain nonlinear dynamics on a 6D planar quadrotor using (a) : Classical Closed-Loop MPC and (b) : robust MPC via GPUSLS. Classical MPC fails to reliably avoid the obstacle, whereas our method consistently achieves safe navigation.
Fig. 8: (a) : Whole-body humanoid (75D, 19C) navigating around an obstacle using ( 14 ) calculated from GPU-SLS. (b) : A visualization of the robust forward tubes (blue squares) generated by ( 18 ) and the humanoid’s rolled out trajectory remaining within the tubes.
Fig. 9: Simulated rollouts of 8 under adversarial disturbances. Using the controller ( 14 ) the humanoid is able to safely remain within the tubes, demonstrating our method’s ability to develop robust control policies for complex, high-dimensional nonlinear systems.
Fig. 10: Runtime comparisons between our GPU solver and the CPU-based solver OSQP for varying horizons on a whole rigid-body quadruped (61D) navigating around an obstacle. Across all horizons, our method outperforms OSQP, achieving up to a 95% reduction in solve time at the longest horizon.
Fig. 11: (a) : GPU-accelerated robust whole-body control (61 states, 12 controls) using RTI GPU-SLS executed on hardware with a Unitree Go2 EDU quadruped. The robot successfully navigates robustly around an obstacle in real-time at 50 Hz. (b) : Graph illustrating the robot’s distance to the obstacle and the maximum tube size for each time step along the trajectory. The tubes are used to tighten the constraints in the nominal trajectory generator, ensuring that the nominal policy remains safe despite disturbances.
Appendix figures & tables5 assets
Supplementary material from the paper’s appendix.
Appendix
Fig. 12: Upsweep phase of the parallel associative scan. Leaf nodes correspond to the ordered elements in S and parent nodes store partial reductgions computed bottom-up using ⊗ .
Fig. 13: Downsweep phase of the parallel associative scan. Using the subtree aggregations during the upsweep, prefixes are propagated top-down along the same binary tree following the rules ( 52 ).
Fig. 14: Whole-body humanoid (75D, 19C) robustly navigating through two obstacles using GPU-SLS in a Real-Time Iteration scheme. Our formulation is able to solve for robust control policies on average in 22 ms on a NVIDIA RTX 4090 GPU, demonstrating the real-time capabilities of our formulation.
Fig. 15: Visualization of (px,py) -slices of the robust reachable tubes at certain time steps for the whole-body humanoid rollout in Fig. 14 . At each time step, the calculated robust policy guarantees that our system stays within the calculated tubes, ensuring that the constraints are never violated even under external disturbances. Each obstacle is inflated to account for the size of the humanoid’s torso.
Fig. 16: Simulated rollouts of Fig. 8 , showing disturbed states and robust tubes under adversarial disturbances.
School of Aerospace Engineering, Tsinghua University, 100084 Beijing, China. · Department of Aerospace Science and Technology, Politecnico di Milano, 20156 Milan, Italy.