Multiple-Interval LGR Pseudospectral Methods

Single-interval pseudospectral methods approximate the state and control over the entire time horizon using a global polynomial. These methods are highly accurate when the solution is smooth. However, many realistic optimal control problems contain:

• bang–bang controls;

• state or control discontinuities;

• boundary layers;

• sharp changes in curvature;

• active-set transitions;

• nonsmooth dynamics;

• path-constraint junctions.

A single global polynomial is often inadequate for such problems. The practical remedy is to divide the horizon into multiple mesh intervals and perform local collocation in each interval.

This section develops the multiple-interval Legendre–Gauss–Radau (LGR) pseudospectral method in detail. The main goals are to:

1. review LG, LGR, and LGL collocation structures;

2. explain why global polynomial approximations may fail;

3. derive the multiple-interval LGR discretization;

4. obtain the quadrature approximation of the objective;

5. derive the local collocation equations;

6. impose state continuity between intervals;

7. describe the block-sparse nonlinear programming structure;

8. explain mesh refinement and practical implementation.

Review of Gaussian Collocation Points¶

Three common sets of orthogonal collocation points are:

1. Legendre–Gauss (LG);

2. Legendre–Gauss–Radau (LGR);

3. Legendre–Gauss–Lobatto (LGL).

All are commonly defined on the reference domain

Legendre–Gauss Points¶

LG points lie strictly inside the interval:

Neither endpoint is included.

For an LG transcription, the state approximation typically includes the initial point, the interior LG nodes, and a separately reconstructed terminal state. The controls are defined only at the LG collocation points.

Legendre–Gauss–Radau Points¶

Forward LGR points include the left endpoint but exclude the right endpoint:

A terminal state value at

is added to the state approximation, but it is not a collocation point.

Thus, an LGR discretization uses:

• $N$ collocation points;

• $N+1$ state approximation points;

• $N$ control approximation points.

Legendre–Gauss–Lobatto Points¶

LGL points include both endpoints:

The state and control are typically represented at all LGL nodes.

Why Multiple Intervals Are Necessary¶

Global Polynomial Approximation¶

Suppose the entire time interval

is mapped to

A global polynomial approximation is then used:

where $\ell_j$ are Lagrange interpolation polynomials.

This approach is spectrally accurate for analytic solutions. The error can decrease extremely rapidly with increasing polynomial degree.

Failure Near Discontinuities¶

If the exact control has a jump, such as

then a global polynomial remains smooth and cannot represent the discontinuity exactly.

Increasing the polynomial degree produces oscillations near the jump. This is the Gibbs phenomenon.

This behavior is analogous to approximating a square wave using a truncated Fourier series.

Local Polynomial Approximation¶

The remedy is to divide the domain into intervals:

where

A separate polynomial is used in each interval. Then:

• discontinuities can be aligned with interval boundaries;

• smooth portions can use high-order polynomials;

• difficult regions can use smaller intervals;

• sparsity is preserved.

Comparison of LG, LGR, and LGL on Multiple Intervals¶

Consider two adjacent intervals:

The state must satisfy continuity:

The control need not be continuous unless the physics or problem formulation requires it.

LG Structure¶

LG points exclude both interval endpoints.

Consequences:

• neither the end of interval 1 nor the beginning of interval 2 is a control collocation point;

• the interface state is reconstructed;

• no control is directly available at the interface.

LGR Structure¶

Forward LGR includes the left endpoint of each interval but excludes the right endpoint.

At the interface:

• the terminal state of interval $k$ equals the initial state of interval $k+1$;

• the control is defined at the start of interval $k+1$;

• only one control value is associated with the interface;

• the control is not artificially double-valued.

This is one reason LGR is attractive for multiple-interval transcription.

LGL Structure¶

LGL includes both endpoints of every interval.

At a shared interface, there are two control values:

These values may differ if the optimal control is discontinuous. If the user incorrectly enforces equality, a physically valid control jump may be suppressed.

Original Optimal Control Problem¶

Consider the Bolza problem

The dynamics are

Endpoint constraints are

Path constraints are

Global Time Transformation¶

Map physical time $t\in[t_0,t_f]$ to

using

The inverse mapping is

The derivative relationships are

Thus,

The cost becomes

Mesh Definition¶

Divide the transformed domain into $K$ mesh intervals:

The $k$th interval is

The quantities

are called mesh points.

For a fixed mesh, they are constants during the nonlinear programming solve.

The mesh width is

A uniform mesh satisfies

A nonuniform mesh allows smaller intervals in difficult regions.

Local Coordinate Mapping¶

Standard LGR points are defined on

For interval $\mathcal{S}_k=[s_{k-1},s_k]$, map $\tau$ to $s$ using

Equivalently,

The derivative scaling is

If the standard LGR nodes are

then the nodes in interval $k$ are

The noncollocated terminal state point is

Local State Approximation¶

In interval $k$, approximate the state using

where

The Lagrange basis polynomial is

It satisfies the cardinal property

Control Approximation¶

The control is represented only at the $N_k$ LGR collocation points:

A polynomial control interpolant may be formed if desired:

However, in the NLP, only the nodal control values are usually required.

Local Differentiation Matrix¶

Differentiate the state approximation:

Evaluate at the collocation points:

where

The local differentiation matrix has dimensions

Scaling from the Standard LGR Matrix¶

Let

be the LGR differentiation matrix on $\tau\in[-1,1]$.

Because

the interval differentiation matrix is

Collocation Equations¶

The transformed dynamics are

At node $i$ of interval $k$:

Thus,

These equations hold for

Physical Time at Each Node¶

The global transformed coordinate $s_i^{(k)}$ maps to physical time through

If $t_0$ or $t_f$ is free, these nodal times depend on the NLP decision variables.

Quadrature over Multiple Intervals¶

For an arbitrary integrand $G(s)$:

In interval $k$, LGR quadrature gives

The interval weights scale from the standard LGR weights:

Therefore,

Discrete Cost Functional¶

The objective becomes

Hence,

State Continuity between Mesh Intervals¶

For a continuous state trajectory:

There are two implementation choices.

Explicit Continuity Constraints¶

Keep both interface state variables and impose

Shared State Variables¶

Use the same NLP variable for both sides:

This eliminates explicit continuity constraints and reduces the number of decision variables.

Control Behavior at Interfaces¶

Forward LGR has no control at the terminal point of interval $k$, but it does have a control at the initial point of interval $k+1$.

Thus, there is one interface control value:

This is especially useful when the control is discontinuous. The method does not force left and right control limits to be equal.

Endpoint Constraints¶

The discrete endpoint constraints are

The initial state comes from the first node of the first interval. The final state comes from the additional terminal state point of the last interval.

Path Constraints¶

At every collocation point:

These constraints are imposed for

Because constraints are enforced only at collocation nodes, the continuous-time solution should be checked between nodes after solving.

Discrete Nonlinear Programming Problem¶

The NLP decision vector contains:

• state values in every mesh interval;

• control values at every LGR node;

• free initial and final times;

• static parameters.

A generic decision vector is

The NLP is

subject to

Global Matrix Structure¶

Let

The global differentiation matrix is block diagonal:

Unlike a single global pseudospectral differentiation matrix, this matrix is sparse.

The global collocation equations can be written schematically as

Sparsity Pattern¶

Each local collocation defect depends only on:

• state variables in the same interval;

• controls in the same interval;

• global time variables;

• static parameters.

Therefore, the Jacobian has a block-banded structure.

This sparsity is one of the principal computational advantages of multiple-interval methods.

Number of Decision Variables¶

Suppose:

• $n_x$ states;

• $n_u$ controls;

• $n_p$ static parameters;

• $N_k$ LGR points in interval $k$.

Without shared interface states, the state decision variables are

The controls contribute

With shared interface states, the number of unique state nodes is

Thus,

Defect Constraints¶

Define the defect at node $i$ of interval $k$:

The collocation conditions are

For all intervals, the number of defect equations is

Integral Consistency Check¶

LGR quadrature also gives an integral form of the state update:

This relation can be used as a consistency check.

Mesh Refinement¶

The mesh should adapt to the solution.

Two basic refinement operations are:

$h$-Refinement¶

Divide an interval into smaller intervals.

This is effective near:

• discontinuities;

• sharp control switches;

• path-constraint junctions;

• localized high-curvature regions.

$p$-Refinement¶

Increase polynomial degree $N_k$ within an interval.

This is effective when the solution is smooth and the approximation error is distributed across the interval.

$hp$-Refinement¶

Combine interval subdivision and degree elevation.

Modern pseudospectral optimal control software often uses adaptive $hp$ refinement.

Error Estimation¶

A common error estimator compares:

1. the polynomial state approximation;

2. a state reconstructed by integrating the dynamics using the approximate control.

Let

be the collocation polynomial and

be an integrated state reconstruction.

Define

If

the interval is refined.

Detecting Nonsmoothness¶

A smooth solution is typically characterized by rapidly decaying polynomial coefficients.

Suppose the local state approximation is written in an orthogonal basis:

If

decreases rapidly with $j$, $p$-refinement is appropriate.

If the coefficients decay slowly or stagnate, the interval may contain a nonsmooth feature and should be subdivided.

Switch Detection¶

For bang–bang problems, switch locations can be inferred from:

• rapid changes in the computed control;

• large defect errors;

• active-bound transitions;

• switching-function sign changes;

• concentration of polynomial oscillation.

Placing a mesh boundary near the switch dramatically improves accuracy.

Algorithmic Workflow¶

Algorithm: Multiple-Interval LGR Pseudospectral Method

Choose an initial mesh \{s_k\}_{k=0}^{K}.
Choose polynomial degrees N_k.
Compute standard LGR nodes and weights.
Map nodes and weights into each mesh interval.
Form local differentiation matrices.
Construct the NLP decision vector.
Impose collocation defects.
Impose state continuity.
Impose endpoint and path constraints.
Solve the NLP.
Estimate interval-wise discretization error.
if all interval errors are below tolerance:
    Accept the solution.
else:
    Apply h-, p-, or hp-refinement.
    Interpolate the previous solution onto the new mesh.
    Resolve the NLP.
end if

Practical Implementation Details¶

Initial Guess¶

The initial guess may be generated using:

• linear state interpolation;

• constant control;

• analytical reduced-order solutions;

• simulation-based trajectories;

• a previously solved nearby case;

• continuation in parameters or constraints.

Scaling¶

States, controls, constraints, and time should be scaled so that typical values are near unity.

For example:

Poor scaling can cause:

• inaccurate derivatives;

• slow convergence;

• misleading feasibility tolerances;

• failure of the NLP solver.

Derivative Computation¶

Derivatives may be obtained using:

• analytic formulas;

• automatic differentiation;

• complex-step differentiation;

• finite differences.

Automatic differentiation is usually preferred for large sparse problems.

MATLAB-Like Assembly Skeleton¶

function data = buildMultipleIntervalLGR(mesh, degree)

K = length(mesh)-1;

Dblocks = cell(K,1);
nodes   = cell(K,1);
weights = cell(K,1);

for k = 1:K

    Nk = degree(k);

    % Standard forward LGR data on [-1,1]
    [tau, wHat, DHat] = lgrData(Nk);

    sL = mesh(k);
    sR = mesh(k+1);
    ds = sR-sL;

    % Map collocation nodes to current mesh interval
    s = 0.5*ds*tau + 0.5*(sR+sL);

    % Scale quadrature weights and differentiation matrix
    w = 0.5*ds*wHat;
    D = (2/ds)*DHat;

    nodes{k}   = s;
    weights{k} = w;
    Dblocks{k} = D;
end

data.nodes   = nodes;
data.weights = weights;
data.D       = blkdiag(Dblocks{:});

end

Defect Assembly Example¶

for k = 1:K

    Nk = degree(k);

    Xk = X{k};        % (Nk+1)-by-nx
    Uk = U{k};        % Nk-by-nu
    tk = timeMap(nodes{k}, t0, tf);

    Fk = dynamics(Xk(1:Nk,:), Uk, tk, p);

    defect{k} = D{k}*Xk ...
              - 0.5*(tf-t0)*Fk;
end

Advantages of Multiple-Interval LGR¶

1. Handles nonsmooth solutions more effectively than a global polynomial.

2. Preserves high-order accuracy within smooth intervals.

3. Produces sparse Jacobian and Hessian structures.

4. Avoids double-valued interface controls.

5. Supports adaptive mesh refinement.

6. Provides a natural framework for large-scale trajectory optimization.

7. Works well with sparse nonlinear programming solvers.

Limitations¶

1. Requires mesh-selection and error-estimation logic.

2. Can generate large NLPs.

3. Path constraints are enforced only at nodes unless additional checks are used.

4. Poorly selected meshes may miss short-duration events.

5. Very low-order intervals may reduce accuracy.

6. Very high-order intervals may become ill-conditioned.

7. Discontinuous states require multiphase formulations rather than simple continuity.

Multiple Phases versus Multiple Mesh Intervals¶

Multiple mesh intervals are used when the dynamics and state remain continuous.

Multiple phases are used when there are:

• discontinuous states;

• changes in dynamics;

• staging or separation;

• impact events;

• different control sets;

• phase-specific constraints.

For multiple phases, event constraints connect the terminal state of one phase to the initial state of another:

Verification and Validation¶

A converged NLP solution is not automatically a trustworthy continuous-time solution.

Verification should include:

1. mesh refinement until the objective stabilizes;

2. independent integration of the dynamics;

3. endpoint residual checks;

4. path-constraint checking between nodes;

5. comparison across mesh strategies;

6. sensitivity to tolerances;

7. examination of control oscillations;

8. Hamiltonian or costate checks when available.

Common Mistakes¶

• forgetting the factor $(t_f-t_0)/2$;

• forgetting the interval scaling $2/\Delta s_k$ in $\boldsymbol{D}^{(k)}$;

• using standard LGR weights without interval scaling;

• imposing control continuity when it is not required;

• failing to impose state continuity;

• confusing global coordinate $s$ with local coordinate $\tau$;

• treating mesh points as NLP variables in a fixed-mesh transcription;

• applying one global differentiation matrix across discontinuities;

• failing to check path constraints between nodes;

• using excessive polynomial order in a nonsmooth interval.

Summary¶

The multiple-interval LGR method extends single-interval pseudospectral collocation to realistic optimal control problems.

The principal ideas are:

1. transform the physical time domain to a fixed global domain;

2. divide the transformed domain into mesh intervals;

3. map standard LGR points into each interval;

4. approximate the state locally using Lagrange polynomials;

5. enforce the dynamics at local LGR points;

6. approximate the cost using composite LGR quadrature;

7. enforce state continuity between intervals;

8. exploit the resulting block-sparse NLP structure;

9. refine the mesh until the discretization error is acceptable.

This method combines spectral accuracy in smooth regions with the flexibility required for nonsmooth controls and localized trajectory features.

Connection. Local interval equations become a useful computational method only after their states, controls, differentiation blocks, and continuity conditions are assembled globally.