Numerical Trajectory Optimization Method

• Numerical trajectory optimization is a technique that transforms continuous optimal control problems into finite mathematical programs using occupation measures and moment relaxations.
• It employs inverse problem formulations and atomic approximations to reconstruct optimal trajectories via linear programming, ensuring global convergence as grid resolution improves.
• The method enhances classical approaches by enabling coordinate-wise reconstruction and quantifying residual errors, ultimately supporting certified optimal solutions for complex systems.

A numerical trajectory optimization method refers to a suite of algorithmic techniques that solve optimal control problems by transforming continuous or infinite-dimensional trajectory planning tasks—subject to system dynamics, path constraints, and other operational restrictions—into finite or countably infinite mathematical programming problems, which are then solved using numerical optimization. Modern approaches combine direct transcription with advanced solvers, moment hierarchies, and inverse problems, enabling global optimality guarantees, systematic trade-offs between accuracy and computational expense, and principled reconstruction of state–control trajectories from moment relaxations.

1. Occupation Measures, Moments, and Atomic Approximation

The occupation measure framework elevates the classical trajectory optimization setting to one in which the (state, control, time) trajectory is encoded as a measure μ over the product space $\mathcal{Z} \subset \mathbb{R}^{1+n+m}$, where $n$ is the state dimension and $m$ the control dimension. The moments of μ are defined by

$y_\alpha = \langle z^\alpha, \mu \rangle,$

for each multi-index α. The hierarchy of moment relaxations is constructed by solving successive semidefinite programs (SDPs), where moments up to some degree $2r$ are computed. Under mild conditions (such as the monotonicity and compact support of the state–control process), the sequence of all moments uniquely defines μ.

Since only a finite number of moments can be computed (due to computational constraints and truncation), exact recovery of the underlying continuous trajectory is infeasible. To circumvent this, the occupation measure is approximated by an atomic measure: $$\tilde{\mu} = \sum_{\beta=1}^q w_\beta \, \delta_{z_\beta},$$ where $\{z_\beta\}_{\beta=1}^q$ is a user-specified grid in $\mathcal{Z}$ with mesh size $\varepsilon$, and $w \in \mathbb{R}_+^q$ are nonnegative weights to be determined. This reduces the original infinite-dimensional problem into a finite, tractable optimization.

2. Inverse Problem: Trajectory Reconstruction via Linear Programming

Given the sequence of truncated moments $b = (y_\alpha)_\alpha$ from SDP relaxation, the numerical trajectory reconstruction problem is posed as a minimum-norm inverse problem: $$\lambda_\varepsilon^* = \min_{w \in \mathbb{R}_+^q} \| b - A w \|_\infty,$$ where $A_{\alpha\beta} = z_\beta^\alpha$ and the minimization is performed over the nonnegative orthant. This LP is always feasible: for any finite moment data and grid, there exists at least one choice of $w$ for which the atomic moments approximate the given vector $b$ as closely as possible.

Importantly, the parameter $\lambda_\varepsilon^*$ quantifies the moment-reconstruction residual. As the grid is refined (ε → 0) and higher-order moments are included (r → ∞),

$$\lim_{\varepsilon \rightarrow 0,\, r \rightarrow \infty} \lambda_\varepsilon^* = 0.$$

This guarantees that in the limit, the reconstructed atomic measure converges weakly* to the true occupation measure.

The reconstructed trajectory is then formed by extracting the grid points $z_\beta$ for which $w_\beta^*$ is significant—these points define sampled locations on the optimal (state, control, time) path.

3. Computational Scalability and Coordinate-Wise Reconstruction

Directly matching the full moment sequence in high-dimensional $\mathcal{Z}$ is computationally prohibitive due to exponential growth in both the number of moments and grid points. The method addresses this by supporting coordinate-wise or marginal reconstruction: only moments involving specific (time, state, or control) coordinates are used, greatly reducing $A$'s size and the complexity of the resulting LP.

For instance, by reconstructing time versus a single state or control coordinate (i.e., working with moments $y_{i, 0, ..., 0, j, 0, ..., 0}$), one may efficiently recover 1D or 2D projections of the trajectory, avoiding the full curse of dimensionality while enabling accurate initialization and verification.

This flexibility is critical in practical applications and enables trade-offs between mesh resolution, computational complexity, and the fidelity of trajectory reconstruction.

4. Comparison with Classical Trajectory Optimization Methods

Traditional numerical trajectory optimization approaches—such as direct transcription, multiple shooting, and Hamilton–Jacobi–Bellman grid methods—require discretizing the trajectory or state–control space and formulating an NLP or DP recursion. These methods typically:

• Impose feasibility and continuity via explicit constraints and often require hand-tuned tolerances or ad hoc penalty weights.
• Rely on good initialization to ensure local convergence and often only guarantee local optimality or feasibility.
• Suffer from mesh-dependent artifacts and the curse of dimensionality.

In contrast, the moment-based approach:

• Leverages convex SDP relaxations, offering global convergence guarantees and lower bounds on cost.
• The subsequent LP for atomic measure reconstruction is always feasible, with no need to set or adjust feasibility tolerances.
• Provides explicit control of grid size and LP problem dimensions, with recovery residuals that directly reflect (and can be systematically reduced by adjusting the mesh).
• Produces trajectory reconstructions that can "hot start" traditional local solvers: if the local trajectory matches the relaxed cost within a prescribed tolerance, global optimality is certified.

5. Practical Implementation and Algorithmic Steps

The overall method consists of the following structured algorithm:

1. SDP Moment Relaxation: Solve a sequence of SDPs to obtain truncated moment vectors $b = (y_\alpha)_\alpha$ and lower bounds on the optimal cost.
2. Grid Definition: Enumerate a finite set $\mathcal{Z}_\varepsilon$ of points in (time × state × control) with mesh-size $\varepsilon$.
3. Linear System Construction: Build matrix $A$ with entries $A_{\alpha\beta} = z_\beta^\alpha$ for the chosen multi-indices and grid points.
4. LP Solution: Solve

$$\min_{w \in \mathbb{R}_+^q} \left\| b - A w \right\|_\infty$$

for the optimal atomic weights $w^*$.

1. Support Extraction: Identify grid points $z_\beta$ with significant $w_\beta^*$ as the approximate trajectory points.
2. Local Solver Refinement and Verification: Use the reconstructed trajectory to initialize a local (direct) optimal control solver. If the refined trajectory has cost within tolerance of the SDP's optimal value, global optimality is certified.

6. Case Studies and Extensions

The method's broad applicability is demonstrated on several complex optimal control archetypes:

• Double Integrator (Minimum-Time): The atomic measure approach yields time-series for both state and control that closely match the analytical optimum, even at moderate grid and moment resolutions.
• Non-Convex Integrator with Obstacles: Even with non-convex state constraints and challenging path topology, the method generates feasible atomic trajectories that can efficiently initialize a local method, leading to rapid convergence to the global optimum even where standard methods stall or require meticulous initialization.
• Nonlinear System Limit Sets (Van der Pol Oscillator): By formulating a measure feasibility problem for invariant sets, the atomic measure reconstruction extracts information about global attractors and equilibria—illustrating that the framework extends beyond classic trajectory optimization to invariant set extraction and nonlinear system analysis.

7. Mathematical Summary and Guarantees

Key equations central to the method are:

Name	Formula	Description
Linear Program (LP)	$$\lambda_\varepsilon^* = \min_{w \ge 0} \|b - A w\|_\infty$$	Recovers weights of atomic measure
Atom-to-Moment Map	$A_{\alpha\beta} = z_\beta^\alpha$	Monomial basis at mesh points
Atomic Representation	$$\tilde{\mu} = \sum_{\beta=1}^q w_\beta \delta_{z_\beta}$$	Dirac sum approximates occupation measure
Convergence	$$\lim_{\varepsilon \to 0,\ r \to \infty} \lambda_\varepsilon^* = 0$$	Weak* convergence to true solution

This systematic approach guarantees that, as computational resources allow, the reconstructed trajectory converges to the true globally optimal trajectory and control.

8. Impact and Limitations

The moment-LP inversion method fundamentally advances numerical trajectory optimization by:

• Providing a pathway to globally optimal trajectory computation for complex, nonlinear, or non-convex control problems.
• Avoiding manual tolerance tuning and initialization heuristics found in classical direct/discrete methods.
• Guaranteeing LP feasibility and offering direct quantification and control of numerical error via $\lambda_\varepsilon^*$.
• Enabling certified optimality by matching local refinement solutions to relaxed global bounds.

Limitations include the computational cost of high-order SDPs and the exponential scaling of the moment cone and support grid size, which restricts practical application to moderate system dimensions unless specialized exploitation of structure or coordinate-wise projections are employed. Nonetheless, for low- to moderate-dimensional systems and as an initializer or global validator for local solvers, the approach is both powerful and flexible.

---

This trajectory optimization methodology represents a rigorous, globally convergent alternative to traditional numerical methods, particularly impactful in problems where global guarantees and initialization of complex local solvers are required, as well as in dynamical systems analysis beyond classical optimal control (Claeys, 2014).