From Optimal Control to a Sparse NLP

Pockit solves continuous-time optimal-control problems by direct collocation, a form of direct transcription. It approximates the unknown state and control functions with finitely many nodal values, then expresses the dynamics and other constraints as algebraic equations and inequalities. The result is a sparse nonlinear program (NLP) that Ipopt or SciPy can solve.

The essential map is

\[ \text{continuous OCP} \longrightarrow \text{nodal polynomial approximations and collocation defects} \longrightarrow \text{sparse finite-dimensional NLP}. \]

Continuous-time problem

A single-phase Bolza problem can be written as

\[ \begin{aligned} \min_{x(\cdot),u(\cdot),t_0,t_f,s}\quad & \Phi\!\left(x(t_0),x(t_f),t_0,t_f,s\right) + \int_{t_0}^{t_f} L(x,u,t,s)\,\mathrm{d}t, \\ \text{subject to}\quad & \dot{x}(t)=f(x,u,t,s), \\ & g_L \le g(x,u,t,s) \le g_U, \\ & b_L \le b\!\left(x(t_0),x(t_f),t_0,t_f,s\right) \le b_U. \end{aligned} \]

Here \(x(t)\) is the state, \(u(t)\) is the control, and \(s\) denotes static parameters that do not vary with time and may be shared across phases. Because \(x(\cdot)\) and \(u(\cdot)\) range over functions, the problem is infinite-dimensional. Transcription approximates those functions with a finite set of variables.

Normalize time and form polynomial approximations

Pockit maps physical time to the normalized phase coordinate

\[ \tau=\frac{t-t_0}{\Delta t},\qquad \Delta t=t_f-t_0,\qquad 0\le\tau\le1, \]

so the dynamics become

\[ \frac{\mathrm{d}x}{\mathrm{d}\tau} =\Delta t\,f(x,u,t,s). \]

Pockit divides the normalized phase into mesh intervals. Within each interval, it approximates the states and controls with polynomials defined by their values at Radau or Lobatto nodes. The unknown trajectories are thus represented by arrays of nodal values \(X\) and \(U\); any free phase times and static parameters become additional finite-dimensional variables.

This is a direct method because the nodal state and control values are optimized directly. It does not first derive and solve the costate boundary-value problem from Pontryagin's minimum principle.

Turn dynamics into defect equations

Let \(f_j=f(X_j,U_j,t_j,s)\). For each state component, Pockit enforces integral-form collocation defects of the form

\[ d_r=(T X)_r-\Delta t\sum_j I_{rj}f_j=0. \]

The matrix \(T\) combines state samples to represent state differences within the mesh intervals, while \(I\) integrates the polynomial interpolant of the sampled dynamics. Each row therefore requires the state change implied by the polynomial to agree with the integral of the dynamics. The selected mesh and the Radau or Lobatto basis are encoded in these matrices.

This finite set of defect equations replaces the continuous differential equation \(\dot{x}=f(x,u,t,s)\). The defects are nonlinear algebraic equations whenever the dynamics are nonlinear.

Integral and differentiation-matrix forms

Direct collocation is also often presented in differentiation-matrix form as \(D X=\Delta t F\). Pockit uses the integral-form defects above. Both formulations express the same collocation principle, but their row ordering, endpoint treatment, and multiplier formulas depend on the exact transcription and cannot be mixed directly.

Discretize the objective and constraints

With normalized quadrature weights \(w_j\), Pockit approximates the Bolza objective by

\[ F_N=\Phi+\Delta t\sum_j w_j L(X_j,U_j,t_j,s). \]

The remaining model components become finite-dimensional terms and constraints:

Continuous object	Discrete NLP representation
Dynamics \(\dot{x}=f\)	Integral defect equalities \(d_r=0\)
Running cost or integral	Weighted quadrature sum
State and control bounds	Bounds on nodal decision variables
General path constraint	Algebraic bounds evaluated at applicable nodes
Endpoint condition	Algebraic constraint on endpoint values and times
Phase linkage or reset	System constraint between phase endpoints and static parameters

Collecting all unknown nodal values, free times, and static parameters into a vector \(z\) gives the standard finite problem

\[ \begin{aligned} \min_z\quad & F_N(z),\\ \text{subject to}\quad & c_L\le c(z)\le c_U,\\ & z_L\le z\le z_U. \end{aligned} \]

Apart from shared quantities such as phase times and static parameters, most defect and path-constraint rows depend only on variables in one mesh interval. Pockit preserves this local structure when it constructs symbolic first- and second-order derivatives, producing the sparse Jacobians and Hessians used by the NLP solver.

Radau and Lobatto change the discrete representation

Both families perform the same OCP-to-NLP transformation, but they place nodes and join neighboring interval polynomials differently.

Property	Radau	Lobatto
Endpoint pattern	Asymmetric	Includes both interval endpoints
Shared mesh boundary	Left and right values can be represented separately	Adjacent intervals share a value
Reconstructed function	Can represent a jump at a mesh boundary	Continuous across the phase
Typical reason to choose it	Event-aligned or bang-bang discontinuity	Smooth state and control histories

The original continuous model is unchanged; the node family changes the finite-dimensional approximation space and therefore the NLP seen by the solver. See Radau and Lobatto Interpolation for practical mesh and degree choices.

What solver convergence does and does not establish

Solver convergence applies only to the current finite transcription. It does not, by itself, establish that the reconstructed trajectory accurately solves the continuous problem. Three checks answer different questions:

1. NLP convergence: does the solution satisfy the discrete KKT conditions and nodal constraints to the requested solver tolerances?
2. Discretization convergence: do the objective and trajectory remain stable when the mesh or polynomial degree changes?
3. Continuous-time validation: when the trajectory is evaluated densely, including between nodes, do its dynamics residuals, path constraints, endpoint conditions, and physical invariants meet the required tolerances?

Path constraints are generally enforced at nodes, so a high-order polynomial can overshoot between them. A credible result therefore combines solver convergence with mesh comparison and dense reconstruction. The error checking and mesh refinement guide describes that workflow.

A minimal example of the map

For the double integrator

\[ \dot{x}=v,\qquad \dot{v}=u,\qquad -1\le u\le1, \]

transcription introduces nodal variables for \(x\), \(v\), and \(u\). The position defects match changes in position to the integral of velocity; the velocity defects match changes in velocity to the integral of control. The control limit becomes bounds on the control nodes, the endpoint conditions become algebraic equalities, and a free \(t_f\) becomes one more decision variable. Minimizing \(t_f\) then gives a finite-dimensional sparse NLP.

The getting-started tutorial builds and solves this model. The double-integrator reference derives the analytical bang-bang solution. It is an example and verification benchmark, not a general account of transcription theory.