Fast and Certifiable Trajectory Optimization

1 Introduction

• The target problems have polynomial losses and dynamics with polynomial constraints, making them generally nonconvex and NP-hard; convex relaxations can provide suboptimality certificates.The introduction focuses on generic nonconvex problems that cannot be losslessly convexified or reduced to combinatorial integer optimization.Source: §1 · p. 2
• Second-order moment-SOS relaxations are costly, while local solvers can get stuck in bad minima and depend on high-quality initial guesses that are difficult to obtain.The paper says first-order relaxations are too loose for trajectory optimization; generic interior-point SDP solvers also scale poorly, cannot be warmstarted, and do not exploit the relaxation structure.Source: §1 · pp. 2–3
• STROM’s C++ conversion generates chain-like multiple-cone SDPs with PSD variables, two orders of magnitude faster than YALMIP and SOSTOOLS and on par with TSSOS.The conversion uses the primal moment perspective rather than the dual SOS perspective used by the Matlab tools.Source: §1 · p. 3
• cuADMM implements sGS-ADMM in CUDA, parallelizing PSD-cone projections and sparse linear solves; it is the first GPU-running ADMM-based SDP solver and achieves up to 10× speedup over other ADMM solvers.The cited comparators are CDCS and SDPNAL+.Source: §1 · p. 3
• Across five applications, cuADMM achieves certificates below 1% suboptimality in minutes or seconds, while other solvers take hours, run out of memory, or take minutes.The applications are inverted pendulum, cart-pole, vehicle landing, flying robot, and car back-in with obstacle avoidance.Source: §1 · pp. 3–4
• 0.66 seconds: data-driven k-nearest-neighbor warmstarting yields a certifiably optimal inverted-pendulum swing-up with a below-1% optimality certificate.The paper reports this as its first certifiably optimal swing-up computed in subseconds.Source: §1 · p. 4

2.1 Polynomial Optimization with Chain-like Sparsity

• Trajectory optimization’s dynamics and set constraints involve variables from each adjacent state-control-state clique, yielding a chain-like sparsity pattern.The objective decomposes into terms associated with individual cliques.Source: §2.1 · p. 5, Figure 1 · p. 5

2.2 The Sparse Moment-SOS Hierarchy

• The sparse moment-SOS hierarchy uses a truncated moment sequence as its convex variable, with matrix constraints providing necessary conditions for a representing measure.The sequence can be interpreted as moments of a probability measure supported on the POP’s feasible set.Source: §2.2 · p. 7
• For bounded clique variables, sparse moment relaxations provide convergent lower bounds and can certify tightness and extract optimal solutions when the rank condition holds.When the condition fails, a spectral-initialization and local-optimization heuristic extracts a near-optimal solution and evaluates its suboptimality certificate.Source: §2.2 · pp. 7–8

2.3 Faster Conversion to Standard SDP

• The conversion produces a chain-structured SDP with only PSD variables, and its package enables real-time conversion in Matlab.The resulting structure mirrors the POP’s chain-like structure.Source: §2.3 · pp. 8, 10

3 Algorithmic Backbone: sGS-ADMM

• sGS-ADMM addresses the poor scalability of interior-point solvers for practical trajectory-optimization SDPs.The method is applied to the Lagrangian dual and has established global convergence for the SDP pair.Source: §3 · p. 10
• Each iteration’s main computational tasks are solving linear systems with fixed coefficient matrix AA* and projecting symmetric matrices onto the PSD cone.The algorithm terminates at the iteration limit or when the maximum KKT residual falls below tolerance.Source: §3 · p. 10
• Because first-order methods achieve only moderate KKT accuracy, the method constructs a valid lower bound to refine the suboptimality certificate.The certificate also requires an upper bound on the trace of each lifted feasible-solution matrix.Source: §3 · pp. 10–11

4 GPU Speedup: cuADMM

• The CUDA implementation assumes relaxation order 2, equal-sized cliques, and inequality constraints of degree 1 or 2, conditions met by all applications considered.Source: §4 · p. 12
• cuADMM parallelizes sparse matrix-vector operations and projections onto the Cartesian product of PSD cones.Clique-indexed storage supports locality in sparse operations, and the implementation uses different eigensolvers for moment and localizing matrices.Source: §4 · pp. 12–13
• The hybrid sparse linear-system solver is about 4 times faster than solving two sparse triangular systems on a GPU and 10% faster than two CPU permutations.It uses GPU vector permutations and CPU CHOLMOD LDLT solving.Source: §4 · p. 12

General Setup

• The experiments covered five trajectory problems, rescaled polynomial variables and coefficients to [−1, 1], and compared MOSEK with CDCS and SDPNAL+ as baselines.Tests ran on a workstation with two NVIDIA RTX 6000 Ada GPUs and 1 TB RAM; reported runtimes exclude local-solver time.Source: §5 · pp. 13–14

Summary Results

• Across the experiments, cuADMM’s mean and median suboptimality gaps were well below 10^-2 except near 10^-2 for car back-in, with up to 10× speedup over other ADMM solvers on large SDPs.MOSEK could not scale beyond pendulum and cart-pole; other ADMM solvers did not produce <1% certificates within the reported times. With warmstarts, pendulum global optimization averaged 0.66 seconds.Source: Table 1 · p. 14, §5 · p. 14

Inverted Pendulum

• With kNN warmstarts, cuADMM solved the pendulum problem in less than 0.7 seconds on average and was 10× faster than MOSEK.The study sampled 10×10 initial states; about 10% were hard for all solvers and formed a spiral-shaped region.Source: §5 · p. 15

Model predictive control

• In the inverted-pendulum MPC simulation, almost all time steps had suboptimality gaps below 10^-2, and average solving time was 0.72 seconds.The figure’s top panel shows MPC with kNN warmstarts.Source: §5 · p. 15, Figure 4 · p. 15

Car Back-in

• In the car back-in tests, cuADMM alone finished all 10 experiments within one hour and produced certifiable results with ξ near 10^-2.The experiments used a separating-hyperplane method for obstacle avoidance; the figure also presents three globally optimal trajectories.Source: §5 · p. 15, Figure 5 · p. 16

Flying Robot

• For the flying-robot problem, cuADMM was the only solver capable of computing optimal trajectories.The robot has four boosters and six continuous-time states, and the figure shows two globally optimal trajectories.Source: §5 · p. 16, Figure 6 · p. 16

6 Conclusion

• STROM combines a C++ sparse-moment-relaxation generator, two orders of magnitude faster than existing Matlab packages, with cuADMM, which solves large-scale SDPs beyond existing solvers’ reach.The paper also demonstrates the potential for real-time certifiable pendulum optimization using data-driven warmstarts.Source: §6 · p. 16

A.1 Trajectory optimization: global solvers

• Prior global-solver work largely emphasized problem reformulation and existing solvers, leaving problem-specific SDP solver design underexplored.The cited work includes sparse Moment-SOS approaches and Shor relaxation; the latter reports a suboptimality gap around 10%.Source: §A.1 · p. 21

A.2 Trajectory optimization: other methods

• Local trajectory solvers remain essential for upper bounds used in global-optimality certificates, and STROM uses them for that purpose rather than replacing them.Local methods include gradient-based and zeroth-order approaches.Source: §A.2 · p. 22

A.3 Polynomial optimization and Moment-SOS hierarchy

• The Moment-SOS hierarchy solves polynomial optimization through a sequence of growing SDPs and can converge finitely under certain constraint qualifications.Correlative, term, and ideal sparsity are among approaches used to improve scalability.Source: §A.3 · p. 22

A.4 Semidefinite programming (SDP) solvers

• Interior-point methods are robust and accurate for small-to-medium SDPs, while first-order ADMM-style methods offer lower per-iteration time and memory costs as scale grows.Existing first-order SDP solvers described here run on CPUs, and specialized low-rank methods can underperform on large multiple-cone SDPs.Source: §A.4 · p. 22

B Insufficiency of First-Order Relaxation

• All sparse second-order pendulum solutions had suboptimality gaps below 10^-2, while none of the first-order moment matrices produced a feasible rounded solution.Across ten randomly selected initial states, first-order SDP objective values had mean 13.7 and median 13.8; second-order values had mean 61.3 and median 59.2.Source: B Insufficiency of First-Order… · pp. 22–23
• First-order relaxation cannot handle polynomial optimization beyond QCQP, including the Car Back-in and Van der Pol problems.Source: B Insufficiency of First-Order… · p. 23

D POP-SDP Conversion Package Comparison

• The sparse moment-relaxation package is the only compared package that converts the pendulum problem to an SDP in real time in MATLAB.The comparison includes SOSTOOLS and YALMIP in MATLAB and TSSOS in Julia.Source: D POP-SDP Conversion Package Comparison · p. 23
• When N = 30, SOSTOOLS takes more than 30 minutes and more than 300GB of memory for pendulum POP-SDP conversion.Source: Table 2 · p. 24

• The pendulum swing-up experiment starts at θ = θ0 and ˙θ = ˙θ0 and targets θ = π and ˙θ = 0 after 30 time steps.The tested dynamics and constraints are discretized from a continuous pendulum model.Source: §E.1 · p. 24
• Three of ten feasible-initialized fmincon attempts failed to swing up the pendulum, and seven of ten noise-perturbed optimal-trajectory attempts became infeasible.The perturbation scales ranged from 10^-1 to 10.Source: §E.1 · pp. 24–25
• Around 10% of initial states were hard for all four SDP solvers, with high-gap states roughly forming a spiral line associated with nonsmooth parts of the optimal value function.The authors note multiple optimal trajectories for these states, which likely explains failure to produce a rank-1 solution.Source: §E.1 · p. 25, Figure 9 · p. 26

E.2 Cart-Pole

• The cart-pole task targets θ = π and zero cart and pole velocities from an initial state, and Figure 10 shows three globally optimal trajectories.The setup uses N = 30 time steps.Source: Figure 10 · p. 27, §E.2 · p. 27
• The simulated cart-pole MPC stabilized the system at 10Hz, although 4 of 47 SDP solutions failed to round to a feasible solution.The controller reused the previous step’s optimal SDP solution as the next initial guess and heuristically applied the extracted first control input, possibly with clipping.Source: §E.2 · p. 27, Figure 11 · p. 28

E.3 Car Back-in

• For car back-in, six lifting variables encode two separating lines, one between the car and each obstacle, making non-collision equivalent to their existence.The car is modeled as a rectangle and must reach its target pose while avoiding two obstacles.Source: §E.3 · pp. 27–29
• Although the car back-in moment matrices are not rank-1, the submatrix for (xk, uk, xk+1) remains nearly rank-1 in practice and supports trajectory extraction.Source: §E.3 · p. 29

E.4 Vehicle Landing

• Figure 12 shows three globally optimal vehicle-landing trajectories.The spacecraft is tasked with landing at x = 0, y = 10, θ = 0, with translational and angular velocities equal to zero.Source: Figure 12 · p. 30, §E.4 · p. 30
• The vehicle-landing discretization includes position updates from velocity over each time step.The stated updates are rx,k+1 − rx,k = ∆t · vx,k and ry,k+1 − ry,k = ∆t · vy,k.Source: §E.4 · p. 30

E.5 Flying Robot

• The flying-robot task is to reach x = 0, y = 0, θ = 0 with translational and angular velocities equal to zero.The robot is described as drifting with four fumaroles.Source: §E.5 · p. 31