Continuous-Time Transcription and Collocation

The preceding chapters formulated and optimized trajectories through discrete transition maps. Physical models are often given instead by ordinary differential equations, so their states and controls are functions of time. A finite-dimensional optimizer cannot choose an entire function directly. Continuous-time trajectory optimization therefore begins by replacing those functions and their differential equations with finitely many variables and constraints on those variables.

How can that replacement preserve enough of the differential equation to make the resulting trajectory meaningful between the stored time points?

This replacement is called transcription, and the resulting finite optimization problem is a nonlinear program (NLP).

The direct-collocation formulation used here stores state and control values at selected times and interpolates between them with low-degree polynomials. The differential equation is enforced at selected points on each time interval. This chapter develops that construction from one simple example through Euler, trapezoidal, and Hermite--Simpson transcriptions, then applies it to the motion of an overhead crane.

A One-Interval Example¶

How can a constraint on four endpoint values ensure that the chosen control produces the required displacement?

Consider a point that must move from $x(0)=0$ to $x(1)=1$. Its velocity is the control, so

and the objective penalizes squared control effort,

Both $x$ and $u$ are unknown functions. As a first finite approximation, retain only their endpoint values $x_0,x_1,u_0,u_1$ and let the control vary linearly between $u_0$ and $u_1$. Integrating $\dot x=u$ means that the change in state equals the area under this control. The shaded region below is a unit-width trapezoid. Its area is the average of its two endpoint heights, $u_0$ and $u_1$, multiplied by the width.

from pathlib import Path
import sys

from IPython.display import HTML, display

code_dir = Path.cwd() / "code"
if str(code_dir) not in sys.path:
    sys.path.insert(0, str(code_dir))

from collocation_widgets import render_linear_control_area

control_area = render_linear_control_area()
display(HTML(control_area))

The optimizer must choose endpoint states and controls that agree with this area calculation. Otherwise, it could choose $x_0=0$ and $x_1=1$ while setting both controls to zero: the endpoint conditions would hold, but zero velocity could not produce the displacement. To exclude such choices, the stored state change $x_1-x_0$ must equal the change $(u_0+u_1)/2$ predicted by integrating the control. Moving both quantities to the left gives the constraint

The left-hand side measures the mismatch between these two changes and is called the defect. Requiring that mismatch to vanish gives a defect constraint, which enforces the dynamics within the chosen approximation. The equation is algebraic because it relates the four numbers $x_0,x_1,u_0,u_1$ through arithmetic operations; there is no unknown function to differentiate or integrate when evaluating it. An optimizer can therefore check the constraint directly for each candidate set of endpoint values.

Approximating the running-cost integral by the average of its endpoint values gives the finite nonlinear program

The defect integrates the linear control exactly. The objective uses another approximation: $u(t)^2$ is generally quadratic, so averaging its endpoint values need not give its exact integral. Both calculations become exact for a constant control.

The boundary conditions reduce the defect to $u_0+u_1=2$. Substituting $u_1=2-u_0$ into the objective gives

The squared term is minimized at $u_0=1$, which also gives $u_1=1$. The resulting interpolation is the exact solution $u(t)=1$ and $x(t)=t$.

This example already contains the main ingredients of direct collocation. The function values became optimization variables, integration became a weighted sum, and the differential equation became an equality constraint. The remaining sections construct these operations systematically for nonlinear vector dynamics and higher-degree polynomials.

From a Continuous Problem to a Finite NLP¶

Which quantities and constraints must a transcription retain from the continuous-time control problem?

We retain the earlier notation $\mathbf x\in\mathbb R^n$ for the state, $\mathbf u\in\mathbb R^m$ for the control, and $c$ for cost. Here $t$ is continuous physical time, while $k$ will index mesh intervals. The running cost $c(\mathbf x,\mathbf u,t)$ is a cost per unit time; integrating it over an interval produces the counterpart of the discrete stage cost $c_k$. The terminal cost is written $c_f$ because the final physical time is $t_f$. Scalar examples omit boldface.

With inequalities $\mathbf g\leq\mathbf0$ and equalities $\mathbf h=\mathbf0$, the continuous-time Bolza problem is

The equality $\mathbf h=\mathbf0$ imposes endpoint conditions, while $\mathbf g\leq\mathbf0$ imposes constraints along the path, such as actuator limits or a bound on the state. A fixed final time is supplied as data and omitted from the decision variables. Setting $c=0$ gives the Mayer special case, while setting $c_f=0$ gives the Lagrange special case. All three forms use the same transcription.

Two transcription strategies differ in which values become decision variables and in how they impose the differential equation:

Shooting and collocation are therefore distinct transcription strategies. A collocation method may reproduce a familiar integration formula, but it exposes the state values to the NLP rather than hiding every state update inside a simulation. The term direct indicates that the continuous problem is converted directly into an optimization problem, without first deriving necessary conditions such as the Pontryagin equations.

To apply a local approximation repeatedly over the horizon, divide time into intervals with endpoints

This sequence is a mesh, with $k=0,\ldots,N-1$ indexing its intervals. We write $\mathbf x_k$ and $\mathbf u_k$ for values at mesh time $t_k$, as in the earlier discrete-time chapters; the subscript is an index, not a physical time. On each interval, the normalized coordinate is $\tau=(t-t_k)/h_k$, or equivalently

Every physical interval is thereby mapped to the same reference interval $[0,1]$. Differentiation and integration formulas can be constructed once on this reference interval; the length $h_k$ then supplies the physical-time scaling. If the final time is also optimized, the interval lengths become variables or fixed fractions of the variable horizon; the reference-interval operators remain unchanged.

Polynomial pieces on a time mesh¶

A single polynomial can represent a trajectory over the entire horizon; the opening example’s solution $x(t)=t$ already does so. Why introduce several pieces? Consider the same scalar dynamics $\dot x=u$, starting from $x(0)=0$, but now moving right at unit speed until $t=1/2$ and then left at unit speed:

The state is continuous, but its slope changes abruptly at the switching time. Two line segments represent it exactly. A single polynomial has a continuous derivative, so it cannot reproduce this corner exactly, although it can approximate the trajectory. The ODE here holds on either side of the switch; the control’s value at that one instant does not affect its integral.

The comparison below uses the quadratic $q(t)=2t(1-t)$, which passes through the same states at $t=0,1/2,1$. Matching these three values does not reproduce the motion between them: the quadratic has a continuously varying slope, while the prescribed velocity switches from +1 to -1.

Even when the trajectory is smooth, some portions may change much faster than others. Separate polynomial pieces let us shorten the mesh intervals near a rapid change while retaining longer intervals elsewhere. Increasing the degree of one global polynomial adds flexibility across the whole horizon. The local representation also gives the optimizer useful structure: each interval’s dynamics constraints involve only values on that interval and neighboring endpoint states. A global polynomial couples values across the horizon through its differentiation matrix.

For a trajectory that is smooth throughout the horizon, a single polynomial of sufficiently high degree can be an efficient choice. Piecewise polynomials give us local control over resolution, allow changes in slope at joins, and keep the constraint derivatives sparse. They require us to connect adjacent pieces explicitly so that the state remains continuous.

We therefore represent the trajectory by a separate polynomial on each mesh interval $[t_k,t_{k+1}]$. The full approximation $\mathbf{x}_h$ is piecewise polynomial. To construct the piece on interval $k$, use its local coordinate $\tau=(t-t_k)/h_k$, where $h_k=t_{k+1}-t_k$, and write

Thus $\mathbf{p}_k(0)$ is the state at the start of this interval and $\mathbf{p}_k(1)$ is the state at its end. On the next interval, the local coordinate starts again at zero, and a different polynomial $\mathbf{p}_{k+1}$ describes the trajectory. The reference interval and its integration weights can be reused even when the physical intervals have different lengths.

One Polynomial, Two Coordinate Systems¶

Should a polynomial trajectory be stored by its coefficients or by its values at selected points, and how are the two descriptions related?

The opening example stored a linear control by its two endpoint values. Each piece of a collocation trajectory needs a similar reconstruction: given the values the optimizer stores, what curve do they define between nodes? Start with a scalar polynomial $p$ on the reference interval $[0,1]$. The construction applies to each state and control component.

Suppose $p(0)=y_0$ and $p(1)=y_1$, and restrict $p$ to degree at most one. A line has the form $p(\tau)=a_0+a_1\tau$. The first endpoint condition gives $a_0=y_0$, and the second gives $a_0+a_1=y_1$, so $a_1=y_1-y_0$. Substitution and regrouping give

The weights $1-\tau$ and $\tau$ enforce the two endpoint conditions: at 0, the formula returns $y_0$, and at 1, it returns $y_1$. Between the endpoints, the weights vary continuously; at $\tau=1/4$, for example, the value is $3y_0/4+y_1/4$. This is the unique line through the two prescribed values. Without the degree restriction, other curves could pass through them. Even for this same line, $y_0+(y_1-y_0)\tau$ is an equivalent representation. The weighted form is convenient because each stored value has its own function that determines its contribution throughout the interval. The functions $1-\tau$ and $\tau$ form a basis for linear polynomials, just as 1 and $\tau$ do. In the first basis, the coefficients are the endpoint values $y_0,y_1$; in the monomial basis, they are $a_0=y_0$ and $a_1=y_1-y_0$.

Higher-degree pieces use the same idea with more stored values. To specify which curves those values may define, restrict the degree to at most $r$. The resulting polynomial space is

This space, denoted by $\mathcal P_r$, has dimension $r+1$. Choosing a basis $\{\phi_0,\ldots,\phi_r\}$ gives coefficient coordinates

The monomial choice $\phi_j(\tau)=\tau^j$ gives the familiar expression $a_0+a_1\tau+\cdots+a_r\tau^r$. A different basis changes the coefficients used to describe the same function.

The same polynomial can instead be identified by its values. Choose $r+1$ distinct points

and record

These points are called support nodes, or interpolation nodes: once the degree is restricted to at most $r$, the $r+1$ stored values determine the whole polynomial. Storing values is useful for trajectory optimization because a bound on the state at a node then becomes a bound on a decision variable. To evaluate the dynamics between nodes, however, we need a formula that reconstructs the polynomial from those values.

With more nodes, the same construction needs one weight function $\ell_j$ per stored value $y_j$. To recover $y_j$ at its own node without altering the values at the other nodes, $\ell_j$ must equal one at $\sigma_j$ and zero at every other support node. A polynomial with those zeros contains the factors $(\tau-\sigma_m)$ for all $m\ne j$. Dividing their product by its value at $\sigma_j$ makes the value there equal to one. This constructs the Lagrange basis function, also called a cardinal function,

At a support node $\sigma_i$, one factor in the numerator is zero unless $i=j$. When $i=j$, every numerator equals its corresponding denominator. Consequently,

where the Kronecker delta $\delta_{ij}$ is shorthand for one when $i=j$ and zero otherwise. This identity specifies the values at the nodes only. At any other $\tau$, the product formula gives a smoothly varying polynomial weight, which can be negative or exceed one for higher degrees. Thus, unlike a line between two endpoint values, a higher-degree interpolant can take values outside the range of its stored values. The cardinal functions reconstruct the polynomial throughout the interval by

At a node $\sigma_i$, all terms except $y_i\ell_i(\sigma_i)=y_i$ vanish. Between nodes, the weighted sum supplies the intervening values. This expression is the Lagrange form of $p$. When the values $y_j$ are prescribed and $p$ is constructed to pass through them, $p$ is called the Lagrange interpolating polynomial. Each $\ell_j$ is one basis function; $p$ is their weighted sum. Changing one stored value $y_j$ by $\Delta y_j$ changes the curve by $\Delta y_j\ell_j(\tau)$, so the same function describes how that variable affects the entire polynomial.

The uniqueness of this reconstruction follows from a basic root-counting argument. If two polynomials in $\mathcal P_r$ have the same $r+1$ nodal values, their difference has $r+1$ distinct roots. A nonzero polynomial of degree at most $r$ cannot have that many roots, so the two polynomials must be identical.

Consequently, every polynomial in $\mathcal P_r$ has a unique expansion in $\ell_0,\ldots,\ell_r$. These functions form the Lagrange basis associated with the chosen support nodes. Its coefficients are exactly the nodal values $y_j=p(\sigma_j)$. Storing a polynomial by its values is therefore also a basis representation, with the basis chosen so that its coefficients have a direct interpretation as function values.

For a scalar polynomial, collect the coefficients and values into column vectors $\mathbf a=(a_0,\ldots,a_r)^\top$ and $\mathbf y=(y_0,\ldots,y_r)^\top$. The coefficient and nodal descriptions are related by the evaluation matrix

Distinct support nodes make $V$ invertible by the same root-counting argument. Thus $\mathbf{a}$ and $\mathbf{y}$ are two coordinate vectors for one polynomial, rather than two different approximations.

For example, take $p(\tau)=1+2\tau-\tau^2$ and the nodes $0,\tfrac12,1$. In monomial coordinates,

whereas the nodal coordinates are

In the Lagrange basis for these three nodes, the same quadratic is

Here $\ell_0,\ell_1,\ell_2$ select the nodes $0,\tfrac12,1$, respectively. Both coefficient vectors describe the same quadratic. The transcription used here chooses coordinates like $\mathbf y$ for the state and control: values that can be inserted directly into the dynamics, cost, and bounds. Monomial coefficients remain useful for deriving formulas that act on those values.

Polynomial spaces, bases, and nodes¶

The construction separates three decisions that are easy to conflate:

• The polynomial space $\mathcal P_r$ specifies which functions are available.

• The basis specifies coordinates for a member of that space. Monomial and Lagrange bases span the same $\mathcal P_r$.

• The nodes specify where values, residuals, or integrals are evaluated.

Changing the basis does not change the exact polynomial space, although it can change numerical conditioning. Changing the nodes changes the interpolation and the operators built from it. The node families introduced later can still use Lagrange nodal coordinates in the NLP.

Polynomial Interpolation and Least-Squares Regression¶

When should a polynomial pass through every supplied value, and when should it approximate those values with a smaller number of coefficients?

Interpolation and polynomial regression impose different requirements on the same data. Six values at six distinct nodes determine one polynomial of degree at most five that passes through every point. If the values are noisy observations and the aim is to estimate a quadratic trend, a quadratic will generally be unable to pass through all six. Least-squares regression then chooses its three coefficients to minimize the sum of squared discrepancies. The figure below applies both choices to the same six points. Write $V$ for the basis evaluation matrix, as in $\mathbf y=V\mathbf a$. It has six rows in both cases, but six columns for the degree-five interpolant and three for the quadratic fit.

If $V$ is square and invertible, minimizing the squared residual also returns the exact interpolant. The two procedures coincide in that case because the chosen polynomial space can satisfy every value condition.

import numpy as np
import matplotlib.pyplot as plt

observation_nodes = np.linspace(0.0, 1.0, 6)
observation_values = (
    0.25
    + 0.95 * observation_nodes
    - 0.25 * observation_nodes**2
    + np.array([0.00, 0.10, -0.07, 0.08, -0.09, 0.02])
)
plot_nodes = np.linspace(0.0, 1.0, 401)

def evaluate_lagrange(nodes, values, points):
    result = np.zeros_like(points)
    for j, node in enumerate(nodes):
        cardinal = np.ones_like(points)
        for m, other_node in enumerate(nodes):
            if m != j:
                cardinal *= (points - other_node) / (node - other_node)
        result += values[j] * cardinal
    return result

interpolated_values = evaluate_lagrange(
    observation_nodes, observation_values, plot_nodes
)
regression_matrix = np.vander(observation_nodes, 3, increasing=True)
regression_coefficients, *_ = np.linalg.lstsq(
    regression_matrix, observation_values, rcond=None
)
regression_curve = np.vander(plot_nodes, 3, increasing=True) @ regression_coefficients
regression_at_nodes = regression_matrix @ regression_coefficients

with plt.rc_context({
    "font.family": "serif",
    "font.serif": ["STIX Two Text", "Times New Roman", "DejaVu Serif"],
    "mathtext.fontset": "stix",
    "font.size": 9,
    "axes.labelsize": 9,
    "xtick.labelsize": 8,
    "ytick.labelsize": 8,
}):
    figure, axis = plt.subplots(figsize=(5.5, 2.8), constrained_layout=True)
    blue = "#0072B2"
    orange = "#D55E00"
    gray = "#4D4D4D"

    axis.plot(plot_nodes, interpolated_values, color=blue, linewidth=1.8,
              label="Degree-five interpolation")
    axis.plot(plot_nodes, regression_curve, color=orange, linewidth=1.8,
              linestyle="--", label="Quadratic least squares")
    axis.vlines(observation_nodes, regression_at_nodes, observation_values,
                color="#777777", linewidth=1.2, zorder=2)
    axis.scatter(observation_nodes, observation_values, color=gray,
                 edgecolor="white", linewidth=0.6, s=30, zorder=3,
                 label="Same six points")
    axis.set_xlabel(r"$\tau$")
    axis.set_ylabel("value")
    axis.spines[["top", "right"]].set_visible(False)
    axis.grid(alpha=0.18, linewidth=0.6)
    axis.set_xlim(-0.025, 1.025)
    axis.set_ylim(0.15, 1.25)
    axis.legend(loc="upper left", frameon=False, fontsize=8)

display(figure)
plt.close(figure)

In a nodal state representation of direct collocation, the stored states are unknown decision variables rather than observations. Each candidate vector of support values defines one interpolating polynomial exactly. The optimizer selects a candidate whose ODE residual vanishes at the collocation nodes. A numerical solver may temporarily work with a scalar measure of constraint violation, but the collocation conditions remain equality constraints rather than a statistical regression loss.

Fixed Operators from Nodal Values¶

How can the stored values give us both the slopes and the integral of every candidate polynomial?

The nodal representation also turns differentiation and integration into matrix and vector products. Once the nodes are chosen, the required arrays can be computed once and reused for every candidate trajectory. The optimizer can then evaluate slopes and integrals as it changes the nodal values, without recomputing those arrays. Each cardinal function is a known polynomial, so its derivative and antiderivative can be computed in closed form by changing its coefficients. Constructing the arrays therefore does not require finite-difference approximations or numerical quadrature of the basis functions. The resulting operators are exact for the represented polynomial, apart from floating-point error in their construction and application.

This is possible because the cardinal functions stay fixed while only their coefficients $y_j$ change. In the representation

each $y_j$ is constant with respect to $\tau$, so differentiation acts only on $\ell_j$. To obtain the slopes needed by the ODE, differentiate and then evaluate at a node $\sigma_i$:

while integration over the reference interval gives

The resulting differentiation matrix $D$ and integration weights $w$ depend only on the chosen nodes. Their entries record the derivatives and areas of the fixed cardinal functions. Multiplication by the current nodal values combines these contributions into the derivative or integral of the current polynomial. This is how the transcription supplies the calculus needed by the dynamics and cost through arithmetic on the decision variables.

For a concrete construction, suppose a quadratic is stored by its values $y_0,y_m,y_1$ at $0,\tfrac12,1$, where $m$ denotes the midpoint. The desired outputs are its slope at each of these nodes and its integral over $[0,1]$. To build the arrays that return those outputs for any choice of the three values, first construct the three cardinal functions multiplying them:

These are the three contributions in $p(\tau)=y_0\ell_0(\tau)+y_m\ell_m(\tau)+y_1\ell_1(\tau)$. Their derivatives are $4\tau-3$, $4-8\tau$, and $4\tau-1$. At the left endpoint they give the weights $-3,4,-1$, so $p'(0)=-3y_0+4y_m-y_1$. Evaluation at the midpoint and right endpoint supplies the other two rows of the differentiation matrix:

Thus a matrix-vector product turns the three stored values into the three nodal derivatives. Integrating the same cardinal functions gives

The row $w^\mathsf T$ turns the same three values into the exact integral of their quadratic interpolant. The matrix $D$ differentiates every polynomial in $\mathcal P_2$ exactly, and $w$ integrates every polynomial in that space exactly. “Fixed operator” means that these arrays remain constant while the NLP changes the nodal values.

For the earlier polynomial $p(\tau)=1+2\tau-\tau^2$, the stored values are $\mathbf y=(1,\tfrac74,2)^\top$. Applying the arrays gives

The first result agrees with $p'(\tau)=2-2\tau$ at the three nodes; the second agrees with the integral of $1+2\tau-\tau^2$ over $[0,1]$. The figure collects the coordinate change and these two operations for that same polynomial.

A transcription may use nodes for three distinct purposes:

A method often reuses one set of points for two or three roles. Support nodes belong to a particular polynomial: the state, control, and slope polynomials need not use the same interpolation nodes. For example, a state polynomial can be supported at one set of nodes and differentiated at different collocation nodes. In the integration construction below, the collocation nodes also serve as support nodes for the slope polynomial. Both integration and differentiation forms have collocation nodes. Their algebraic residuals differ: one matches state values after integrating slopes, while the other matches slopes after interpolating state values.

From Nodal Slopes to Collocation Constraints¶

How do those fixed maps convert differential equations into algebraic constraints at the chosen collocation nodes?

The preceding operators apply to any polynomial. A trajectory must also satisfy the dynamics. On each interval, the construction will interpolate the slopes supplied by the ODE and integrate them from the stored left state. Constraints will then require the stored states to agree with the resulting polynomial.

Interval $k$ extends from physical time $t_k$ to $t_{k+1}$ and contains one polynomial piece. Its local time coordinate runs from $\tau=0$ at the left endpoint to $\tau=1$ at the right endpoint. Choose $s$ distinct collocation nodes $\tau_1,\ldots,\tau_s$ on $[0,1]$. A node $\tau_j$ specifies one time within this interval, $t_k+h_k\tau_j$; for example, $\tau_j=1/2$ specifies its midpoint. The state and control values at that time are $\mathbf{x}_{k,j}$ and $\mathbf{u}_{k,j}$.

This evaluation point, together with its state and control values, is called a stage. A stage is associated with a single time, so it has no separate beginning or end. The index $k$ identifies the interval and $j$ identifies a stage in that interval. Depending on the method, stages may be interior points or may coincide with an endpoint: $\tau_j=0$ places a stage at the start and $\tau_j=1$ places one at the end. At each stage, the differential equation prescribes the physical-time slope

for $j=1,\ldots,s$. The figure below shows one such evaluation for the scalar dynamics $\dot x=u$. The state and control are read at the same physical time; the ODE then gives the slope that the state polynomial must match there.

As before, $\mathbf p_k(\tau)=\mathbf x_h(t_k+h_k\tau)$ denotes the state piece in normalized time. Here $\ell_j$ is the Lagrange cardinal function for the slope nodes $\tau_1,\ldots,\tau_s$; the $s$ slope values determine a polynomial of degree at most $s-1$. Interpolating them gives the physical-time derivative on this one interval:

This polynomial agrees with each sampled slope $\mathbf{f}_{k,j}$ at its collocation node. Each sample was evaluated at a stored stage state $\mathbf x_{k,j}$; the state curve recovered by integration need not yet pass through those stored states. The stage equations below will enforce that agreement. To recover state values from these slopes, we need to reverse differentiation. For example, both $\tau+2\tau^2$ and $1+\tau+2\tau^2$ have derivative $1+4\tau$. Each is an antiderivative of $1+4\tau$: differentiating it returns that function. Adding any constant gives another antiderivative, because the derivative of a constant is zero. The slopes therefore determine changes in state, but a starting value is needed to determine the state itself. This starting value fixes the integration constant.

The fundamental theorem of calculus connects antiderivatives to definite integrals. For a scalar function $p$ with a continuous derivative, it states that integrating the derivative between two points gives the difference between the function values:

Here $\eta$ is a dummy variable that runs over the integration interval; $\tau$ specifies where that interval ends. Over a short interval of width $\Delta\eta$, the change in $p$ is approximately $p'(\eta)\Delta\eta$. The definite integral is the limit of sums of these signed changes as the interval widths approach zero. Adding $p(0)$ to this accumulated change recovers the value at $\tau$. Polynomial derivatives are continuous, so the theorem applies to the pieces used here. For a vector state, the same identity applies to each component.

The sampled slopes $\mathbf f_{k,j}$ are derivatives with respect to physical time, whereas $\mathbf p_k$ is expressed in local time. Since $t=t_k+h_k\tau$, the chain rule gives

The factor $h_k$ accounts for the duration of the physical interval: a change of $\Delta\tau$ in local time corresponds to $h_k\Delta\tau$ units of physical time. For example, a constant physical slope $\mathbf f$ acting from 0 to $\tau_i$ changes the state by $h_k\tau_i\mathbf f$.

Applying the fundamental theorem with the left-endpoint condition $\mathbf p_k(0)=\mathbf x_k$ now gives the state piece:

The final expression follows by substituting the slope polynomial and integrating its terms separately; the stage slopes $\mathbf f_{k,j}$ are constant with respect to $\eta$. At $\tau=0$, every integral vanishes, leaving $\mathbf p_k(0)=\mathbf x_k$. Differentiating the expression returns the prescribed slope polynomial. Thus the stored left state selects the unique antiderivative with that starting value. Integrating a polynomial of degree at most $s-1$ gives a state piece of degree at most $s$.

For the scalar example above, $h_k=2$ and the physical-time slope is $1/2+2\tau$, so $p_k'(\tau)=1+4\tau$. Starting from $x_k=1$ gives

The brackets mean evaluating the antiderivative at the upper limit and subtracting its value at the lower limit. At $\tau_i=1/2$, the accumulated change is 1, so $p_k(1/2)=2$. At $\tau=1$, the accumulated change is 3, so $p_k(1)=4$. These are two evaluations of the same polynomial: both integrals begin at the left endpoint, but one stops at a stage and the other at the right endpoint.

This integration is exact for the interpolated slope polynomial. It does not generally give the exact solution of the original ODE: between collocation nodes, the interpolated slope can differ from $\mathbf f(\mathbf p_k(\tau),\mathbf u_h(t_k+h_k\tau),t_k+h_k\tau)$.

The integration construction already determines the derivative at each collocation node. Because $\ell_j(\tau_i)=\delta_{ij}$, any candidate stage values give

This identity holds by construction, before the NLP constraints have been satisfied. It does not yet guarantee that the reconstructed curve satisfies the ODE: the dynamics on the right are evaluated at the stored state $\mathbf x_{k,i}$, which may differ from the curve’s value $\mathbf p_k(\tau_i)$.

The remaining condition is agreement between these two state values. The point $(\tau_i,\mathbf{x}_{k,i})$ must lie on the graph of $\mathbf{p}_k$. Imposing $\mathbf{x}_{k,i}=\mathbf{p}_k(\tau_i)$ at every stage of every interval gives the stage equations

Their coefficients are

Each stage equation therefore sets a state-matching residual to zero:

Once this residual vanishes, the stored state used to evaluate $\mathbf f_{k,i}$ is the state on the reconstructed curve. Substituting $\mathbf x_{k,i}=\mathbf p_k(\tau_i)$ into the derivative identity above then gives the ODE at that node. Thus the integration form enforces the ODE at collocation nodes even though its explicit stage residuals compare state values rather than slopes.

For example, consider $\dot x=x$ on an interval of length $h_k=1$, starting from $x_k=1$, with one collocation node at $\tau_1=1/2$. A trial stage value $x_{k,1}=3$ supplies the slope $f_{k,1}=3$ and hence the state piece $p_k(\tau)=1+3\tau$. Its midpoint value is $5/2$, not the stored value 3. The polynomial’s slope is 3, but the ODE evaluated on the curve at the midpoint asks for $5/2$. The state equation $x_{k,1}=1+\tfrac12 x_{k,1}$ instead gives $x_{k,1}=2$. The reconstructed polynomial is then $p_k(\tau)=1+2\tau$, whose midpoint state and slope both equal 2. The state-matching constraint has made the ODE hold at the collocation node.

Row $i$ of $A$ integrates each slope basis function from 0 to $\tau_i$. Here $i$ selects the destination stage and $j$ runs over all slope stages that contribute to the polynomial. All these integrals start at 0, including when $i>1$; they do not start at the preceding stage. The stage equations are simultaneous constraints, since the slopes on the right depend on the state and control variables being chosen.

For example, take three slope stages at $(\tau_1,\tau_2,\tau_3)=(0,\tfrac12,1)$ and fix the destination $i=2$. Every integral in row 2 of $A$ now ends at $\tau_2=\tfrac12$. Varying $j$ selects the basis function being integrated: $\ell_1$, $\ell_2$, or $\ell_3$. Figure 9 shows these three signed areas. Although the third slope is evaluated at the right endpoint $\tau_3=1$, its basis function $\ell_3$ is nonzero between 0 and $\tfrac12$. That slope therefore contributes to the midpoint state as well. The sum runs over all $j=1,\ldots,s$, not only stages with $j\leq i$.

The right endpoint has local coordinate 1. Evaluating the same polynomial there uses the full integration span $[0,1]$, in place of the span $[0,\tau_i]$ used for stage $i$. Requiring this value $\mathbf{p}_k(1)$ to equal the stored endpoint state $\mathbf{x}_{k+1}$ gives

The coefficients $b_j$ integrate the cardinal functions over the full reference interval. In code, both $A_{ij}$ and $b_j$ are obtained by forming an antiderivative polynomial and evaluating it at the integration limits. No repeated integration of the nonlinear dynamics is involved in computing these coefficients. The worked calculation in A reusable coefficient generator and residual evaluator shows the coefficient arithmetic. If a stage lies at $\tau_i=1$, then $A_{ij}=b_j$ for every $j$, so its stage equation gives the same polynomial value as the endpoint equation. Its state must therefore equal $\mathbf{x}_{k+1}$; an implementation can share that variable.

The endpoint equation is a defect constraint: it requires the stored endpoint to match the endpoint obtained by integrating the slopes of piece $k$. The next piece starts from that same stored value, so the two pieces meet:

This shared variable enforces state continuity across the join. If an implementation stores separate endpoint variables for the two pieces, it must impose an equality between them. Choosing an endpoint as a collocation node does not by itself join the pieces. Continuity of the state also does not require continuity of its derivative; neighboring pieces can meet with different slopes.

The running-cost integral can be treated by exactly the same construction. As in the earlier Bolza-to-Mayer reduction, introduce a scalar state $y(t)$ that records the cost accumulated since $t_0$:

The original objective is now $c_f(\mathbf{x}(t_f),t_f)+y(t_f)$. The augmented state $(\mathbf{x}^{\top},y)^{\top}$ has two rates of change: $\mathbf{f}$ supplies the physical state rate and $c$ supplies the cost rate. Applying the same collocation method to this augmented ODE evaluates both rates at each stage’s state, control, and time. In particular, the cost rate at stage $j$ of interval $k$ is

Let $y_k$ and $y_{k,i}$ denote the accumulated-cost values at the interval’s left endpoint and at stage $i$. Integrating the cost-rate interpolant to a stage and to the right endpoint gives

These are the same stage and endpoint equations as for $\mathbf{x}$, with $c_{k,j}$ in place of $\mathbf{f}_{k,j}$. The increase in accumulated cost across interval $k$ therefore supplies its contribution to the discrete objective:

Here $\mathbf{u}_h$ is the piecewise control approximation. The single-index $c_k$ is an interval cost, whereas $c_{k,j}$ is a sampled cost rate; the factor $h_k$ gives them different units. The weighted sum integrates the polynomial interpolating the sampled cost rates exactly. It approximates the running-cost integral because the composed function $c(\mathbf{x}_h(t),\mathbf{u}_h(t),t)$ need not itself be that polynomial. Such a weighted-sum approximation is called a quadrature rule.

Since $y_0=0$, summing the interval increments gives $y_N=\sum_{k=0}^{N-1}c_k$. The accumulated-cost variables can therefore be eliminated and the sum used directly in the objective. Introducing $y$ explains why applying one collocation scheme to both rates uses the same nodes and weights; it does not require adding cost variables to the implementation. A different quadrature rule is also possible, with the state and control polynomials evaluated at that rule’s nodes.

Choosing a polynomial for the control¶

The state and control play different roles. The ODE constrains the state’s derivative, so integrating an interpolant of $s$ slopes gives a state polynomial of degree at most $s$. There is no corresponding ODE for the control in this problem. Its representation is a separate choice: a constant command on each interval, a line between endpoint commands, or a higher-degree polynomial.

For example, with two endpoint slopes the state polynomial is quadratic. The control can still be the line

These degrees need not match. Also, substituting polynomial state and control curves into a nonlinear $\mathbf f$ need not produce a polynomial. Collocation matches the state derivative to $\mathbf f$ at the chosen nodes, rather than requiring the ODE to hold identically between them.

To make the control choice explicit, select distinct control support nodes $\rho_0,\ldots,\rho_{d_u}\in[0,1]$ and construct their cardinal functions $\psi_0,\ldots,\psi_{d_u}$. The decision variables $\widehat{\mathbf u}_{k,r}$ are control values at these support nodes:

Thus the stage controls $\mathbf u_{k,j}$ are computed from the control variables; they are independent variables only when the representation allows independent values at all those stages. For a constant control, $d_u=0$ and every entry of $B$ is one. For a linear control supported at $0,1$, row $j$ of $B$ is $(1-\tau_j,\tau_j)$.

Using more control support values than collocation stages leaves some control variations invisible to those stage samples. The examples below use $d_u+1\leq s$; a richer control representation needs additional sampling or other constraints to account for those variations.

The ODE requires continuous states, but it can admit controls that jump between intervals. Independent control polynomials therefore need no continuity constraint unless the model or chosen parameterization requires one. Sharing endpoint controls between neighboring linear pieces enforces continuity. With separate local variables, the equivalent equality is

At a jump, endpoint stages use the control belonging to their own interval. Finally, bounds on support values guarantee bounds throughout a constant or linear control segment, because its values are convex combinations of its endpoints. A higher-degree polynomial can overshoot those values, so support bounds alone no longer give that guarantee.

Assembling the nonlinear program¶

The state and control representations now specify which quantities the optimizer chooses and which it evaluates from those choices. The mesh states $\mathbf x_k,\mathbf x_{k+1}$ are shared decision variables at the ends of interval $k$. The stage states $\mathbf x_{k,j}$ are decision variables at the collocation nodes; a stage at an endpoint can use the mesh state itself. For the control, the optimizer chooses the support values $\widehat{\mathbf u}_{k,r}$ at $\rho_r$. Evaluating the control polynomial at $\tau_j$ gives $\mathbf u_{k,j}$. Together with $\mathbf x_{k,j}$ and the time $t_k+h_k\tau_j$, this evaluated control supplies the inputs to $\mathbf f$ and $c$.

Collect the independent state and control values across all intervals into the decision vector

Here $\operatorname{col}$ stacks the listed vectors into one column, omitting duplicate variables when a stage shares a mesh endpoint. Stage controls are evaluated using $B$, and any chosen control-continuity equalities are added to the constraints. If $t_f$ is optimized, it is included in $\mathbf z$ as well. The resulting finite optimization problem has the form

In the preceding chapter’s NLP notation, $F(\mathbf z)$ is this scalar objective, $H(\mathbf z)=\mathbf0$ stacks the stage, endpoint, and boundary equalities, and $G(\mathbf z)\leq\mathbf0$ stacks the sampled path and bound inequalities. The slope $\mathbf f_{k,j}$ is a vector, distinct from $F$.

The original path and bound constraints apply at every continuous time. This NLP imposes them at selected points, usually its support or collocation nodes. Feasibility at those nodes does not rule out a violation between them. Checking the polynomial between nodes and simulating the chosen control provide complementary diagnostics, developed below.

Every interval constraint touches only local stage variables and neighboring endpoint states. Consequently, most derivatives of the constraints with respect to the decision variables are zero. With variables ordered by time, the nonzero blocks lie near the diagonal of the Jacobian, the matrix of constraint derivatives. NLP solvers can exploit this sparse, block-banded pattern. The arrays $A$ and $b$ are numerical constants computed before optimization.

Equivalent differentiation form¶

The differentiation form uses the same collocation nodes $\tau_i$ as the integration form. It changes which consistency condition is built into the representation and which must be enforced by a residual. The integration form constructs a curve whose nodal derivatives match the evaluated stage slopes, then constrains its values to match the stored stage states. The differentiation form constructs a curve from state-support values, evaluates its stage states directly on that curve, and constrains its derivatives to match the ODE.

To represent the same state polynomials, retain the same degree. The integral of the degree-$(s-1)$ slope polynomial has degree at most $s$, so take $d=s$ and choose $d+1$ support nodes $\sigma_0,\ldots,\sigma_d$ for the state polynomial. Write $\ell_r^{\mathrm{state}}$ for their cardinal functions to distinguish this support-node basis from the slope-node basis $\ell_j$. Denote the support values by $\widehat{\mathbf x}_{k,r}$; these need not be the stage values $\mathbf x_{k,j}$ at $\tau_j$:

At a collocation node $\tau_i$, the state is a weighted sum of its support values, and its derivative with respect to $\tau$ is another weighted sum. The fixed arrays containing these weights are

The $i$th row of $E$ computes $\mathbf p_k(\tau_i)$, the point on the state polynomial that we supply to the dynamics function. The $i$th row of $D$ computes $\mathbf p_k'(\tau_i)$, the tangent slope measured per unit of local time $\tau$. Both calculations use the same stored support states $\widehat{\mathbf x}_{k,r}$; they do not introduce new decision variables. In particular, the stage state is defined as $\mathbf x_{k,i}=\mathbf p_k(\tau_i)=\sum_r E_{ir}\widehat{\mathbf x}_{k,r}$. State-value consistency therefore holds for every candidate support array; there is no separate state-matching residual at a stage.

To compare this tangent slope with the rate returned by $\mathbf f$, both must be measured in physical time. A local-time increment $\Delta\tau$ corresponds to a physical-time increment $\Delta t=h_k\Delta\tau$. We therefore divide the local-time slope by $h_k$ to obtain the physical-time slope: this is the chain rule $d/dt=(1/h_k)d/d\tau$ used earlier. The remaining constraint sets the physical-time slope residual to zero: