Adaptive Mesh Refinement for Pseudospectral Optimal Control

Adaptive mesh refinement is the process of repeatedly modifying the discretization until the numerical solution satisfies a prescribed accuracy requirement. In a multiple-interval pseudospectral method, the mesh is defined by:

• the locations of the mesh points;

• the number of mesh intervals;

• the polynomial degree in each interval.

The refinement process is iterative:

1. choose an initial mesh;

2. solve the nonlinear programming problem;

3. estimate the discretization error;

4. modify the mesh;

5. interpolate the previous solution onto the new mesh;

6. resolve the NLP;

7. repeat until the accuracy tolerance is satisfied.

This procedure is analogous to adaptive step-size integration for ordinary differential equations. In smooth regions, relatively large intervals or higher-order polynomials may provide excellent accuracy. Near rapid changes, switching points, singular arcs, boundary layers, or discontinuities, smaller intervals are often necessary.

This section develops an interpolation-and-integration error estimator for multiple-interval Legendre–Gauss–Radau (LGR) collocation and explains how the estimated error can guide $h$-, $p$-, and $hp$-refinement.

Initial Mesh and Polynomial Degrees¶

The transformed time domain is

Choose mesh points

The $k$th mesh interval is

Assign $N_k$ LGR collocation points to interval $k$:

The state is approximated using $N_k+1$ points:

where

Thus, the discretization is determined by:

Mesh Refinement as an Iterative Process¶

Let $\mathcal{M}^{(r)}$ denote the mesh at refinement iteration $r$. The basic process is

The iteration terminates when

where $e_k$ is an interval-wise error indicator and $\epsilon$ is the required tolerance.

State Approximation on One Mesh Interval¶

Consider interval

The LGR state polynomial is

where

Because the state polynomial is known, it can be evaluated at any point in the interval.

Constructing an Enriched LGR Grid¶

To assess accuracy, create a new set of LGR points containing one additional collocation point:

Let the enriched LGR nodes be

The associated additional terminal state point is

Because the new LGR set has a different polynomial degree, most of its interior points differ from the original set.

Interpolated State on the Enriched Grid¶

Evaluate the original state polynomial at the enriched nodes:

Equivalently,

Define the state interpolation matrix

Then

The state interpolation is unique because the state polynomial is explicitly known.

Why Control Interpolation Is Different¶

The control variables in a direct collocation method are nodal samples:

Unlike the state, the control is not necessarily defined by a unique polynomial approximation. Infinitely many continuous functions can pass through the same set of control samples.

Therefore, the control reconstruction is not unique.

Piecewise-Linear Control Reconstruction¶

A simple reconstruction is piecewise-linear interpolation.

Let

denote the control evaluated at the enriched node

Then

Piecewise-linear interpolation is not guaranteed to reproduce the exact control, but it is:

• simple;

• stable;

• local;

• consistent with the control samples;

• usually sufficient for error estimation.

Implicit Integration Form of the LGR Method¶

The LGR collocation method can be written in differentiation form or in an equivalent implicit integration form.

For the original grid in interval $k$:

where:

• $\boldsymbol{A}^{(k)}$ is the LGR integration matrix;

• $\boldsymbol{F}_i^{(k)}$ is the dynamics evaluated at node $i$;

• $\boldsymbol{1}$ replicates the initial state across the rows.

For a scalar state, this becomes

Integrated State on the Enriched Grid¶

On the enriched grid, evaluate the dynamics using the interpolated state and reconstructed control:

The enriched physical time is

Use the enriched LGR integration matrix $\boldsymbol{B}^{(k)}$:

The initial value satisfies

Two Independent State Estimates¶

The enriched grid now contains two state estimates.

Interpolation-Based Estimate¶

This estimate comes directly from the original state polynomial.

Dynamics-Based Estimate¶

is obtained by integrating the dynamics on the enriched grid using the interpolated state and reconstructed control.

If the original solution is sufficiently accurate, the two estimates should be close.

Absolute Error Indicator¶

Define the pointwise error:

For multiple states, the absolute value is taken componentwise.

A scalar interval indicator is

Relative Error Indicator¶

A relative measure is often more useful:

Alternatively, each state can be normalized using a problem-specific scale $S_r$:

The interval is accepted if

Interpretation of the Error Distribution¶

The maximum error alone determines whether refinement is needed. The spatial distribution of the pointwise error helps determine which type of refinement should be used.

Let

Two common patterns are:

Localized Error¶

The tolerance is violated only in a small region.

This may indicate:

• a switching point;

• a boundary layer;

• a sharp curvature change;

• a path-constraint junction;

• a discontinuity.

The preferred response is usually $h$-refinement.

Distributed Error¶

The tolerance is violated throughout most of the interval.

This suggests that the polynomial degree is insufficient over the entire interval.

The preferred response is usually $p$-refinement.

$h$-Refinement¶

In $h$-refinement, an interval is subdivided:

A natural splitting location is the enriched point with maximum error:

and

Safeguards should prevent the new point from being too close to either endpoint.

$p$-Refinement¶

In $p$-refinement, the interval remains unchanged and the polynomial degree is increased:

Typical values are

A small increase is preferred because:

• very high degrees may cause poor conditioning;

• the NLP becomes denser within the interval;

• the error response can be reassessed after the next solve.

Maximum Polynomial Degree¶

Impose an upper bound:

Typical choices might be

If $p$-refinement is indicated but

the interval should be subdivided instead.

Thus:

$hp$-Refinement¶

An $hp$ method permits both:

• interval subdivision ($h$);

• polynomial-degree increase ($p$).

A simple decision rule is:

Localized-Error Metric¶

To automate the decision, define the set

Define the violation fraction:

One simple rule is:

For example,

Then:

• if fewer than 25% of enriched points violate the tolerance, use $h$-refinement;

• otherwise, use $p$-refinement unless $N_{\max}$ has been reached.

Alternative Smoothness Indicator¶

Another strategy is to inspect the decay of polynomial coefficients.

Write the state approximation in a Legendre basis:

For a smooth function, the coefficients should decay rapidly:

If the decay is fast, $p$-refinement is appropriate. If the decay is slow or irregular, $h$-refinement is safer.

Adaptive Refinement Algorithm¶

Algorithm: Adaptive $hp$ Mesh Refinement

Choose an initial mesh and interval degrees.
Set tolerance epsilon, degree increment Delta  p, and N_{\max}.
for r=0,1,2,:
    Solve the NLP on the current mesh.
    for each interval k:
        Set M_k=N_k+1.
        Construct the enriched LGR grid.
        Interpolate the state onto the enriched grid.
        Reconstruct the control using piecewise-linear interpolation.
        Evaluate the dynamics on the enriched grid.
        Integrate the dynamics using the enriched LGR integration matrix.
        Compute the pointwise and interval errors.
        if e_k<=epsilon:
            Keep interval k unchanged.
        else if the error is localized:
            Split interval k near the maximum error.
        else if N_k<N_{\max}:
            Set N_k<- N_k+Delta  p.
        else:
            Split interval k near the maximum error.
        end if
    end for
    if all intervals satisfy the tolerance:
        Stop.
    end if
    Interpolate the previous solution onto the new mesh.
end for

Warm Starting the Refined NLP¶

After changing the mesh, the old solution should be mapped to the new nodes.

For new state node $s_{\mathrm{new}}$:

For new control node $s_{\mathrm{new}}$:

Retain the previous values of:

• initial time;

• final time;

• static parameters;

• phase-link variables.

Warm starts usually reduce the number of NLP iterations significantly.

Pure $h$ Methods¶

A pure $h$ method keeps the polynomial degree fixed:

Only interval widths are adapted.

Advantages:

• simple implementation;

• robust near discontinuities;

• predictable local polynomial degree;

• straightforward error control.

Disadvantages:

• may require many intervals;

• may produce unnecessarily large NLPs;

• can over-refine near isolated switching points.

Pure $p$ Methods¶

A pure $p$ method uses one interval or fixed interval boundaries and increases polynomial degree.

Advantages:

• spectral convergence for smooth solutions;

• compact NLPs;

• high accuracy with relatively few nodes.

Disadvantages:

• performs poorly near discontinuities;

• may exhibit Gibbs oscillations;

• may become ill-conditioned at high degree;

• cannot localize resolution.

Why $hp$ Methods Are Attractive¶

The $hp$ strategy combines the strengths of both methods:

• use high-order approximations in smooth regions;

• use additional intervals near nonsmooth features;

• avoid unnecessary global degree growth;

• avoid excessive interval proliferation;

• keep the NLP smaller for a prescribed accuracy.

The main difficulty is classification: determining whether the interval is smooth enough for $p$-refinement or should instead be divided.

NLP Size Considerations¶

Let

For $n_x$ states and $n_u$ controls, the principal number of nodal variables is

Pure $h$-refinement may increase $K$ rapidly and thereby increase $N$.

Pure $p$-refinement increases the density of local differentiation blocks.

A successful $hp$ method seeks a balance:

Bang–Bang Double Integrator Example¶

Consider

with bounded control

For a minimum-time transfer, the optimal control is typically bang–bang:

Then:

• $x_2$ is piecewise linear;

• $x_1$ is piecewise quadratic.

An ideal mesh places an interval boundary at $t_s$. Then each state segment can be represented exactly by a low-degree polynomial.

A pure $h$ method may repeatedly cluster intervals around the switch. An effective $hp$ method should identify the localized error, place a mesh boundary near the switch, and retain low polynomial degrees on either side.

Why Over-Refinement Is Wasteful¶

Suppose the exact solution is piecewise quadratic. A discretization with many low-order intervals may approximate it accurately, but the NLP can be much larger than necessary.

The goal of adaptive refinement is not merely:

A better goal is:

This is especially important when using general-purpose sparse NLP solvers.

MATLAB-Like Error Estimator¶

function [err,delta,Y,Z] = estimateIntervalError( ...
    X,U,sOld,sNew,t0,tf,Bnew,p)

% Interpolate state polynomial to enriched grid.
Y = interpolateStateLagrange(X,sOld,sNew);

% Reconstruct control using piecewise-linear interpolation.
V = interp1(sOld(1:end-1),U,sNew(1:end-1), ...
            'linear','extrap');

% Map enriched nodes to physical time.
tNew = 0.5*(tf-t0)*sNew(1:end-1) ...
     + 0.5*(tf+t0);

% Evaluate dynamics at enriched collocation points.
Fnew = dynamics(Y(1:end-1,:),V,tNew,p);

% Integrate dynamics from the initial state.
Z = zeros(size(Y));
Z(1,:) = Y(1,:);
Z(2:end,:) = Y(1,:) ...
    + 0.5*(tf-t0)*(Bnew*Fnew);

% Relative pointwise discrepancy.
scale = 1 + abs(Y);
delta = abs(Y-Z)./scale;

% Scalar interval error.
err = max(delta,[],'all');

end

MATLAB-Like Refinement Decision¶

function action = chooseRefinement(delta,tol,Nk,Nmax,rhoCrit)

viol = max(delta,[],2) > tol;
rho = nnz(viol)/numel(viol);

if ~any(viol)
    action.type = 'none';
    return
end

[~,jmax] = max(max(delta,[],2));

if rho <= rhoCrit
    action.type = 'h';
    action.splitIndex = jmax;
elseif Nk < Nmax
    action.type = 'p';
    action.newDegree = Nk + 2;
else
    action.type = 'h';
    action.splitIndex = jmax;
end

end

Safeguards for Mesh Splitting¶

A split point should not be too close to an endpoint.

Define

Choose a minimum fractional distance $\eta$, such as

Then enforce

This avoids creating extremely small intervals.

Stopping Criteria¶

The refinement loop may stop when:

1. every interval satisfies the error tolerance;

2. the maximum number of refinement iterations is reached;

3. the maximum number of nodes is reached;

4. the NLP fails repeatedly;

5. the objective and state trajectories no longer change materially.

A strong convergence test includes both error and solution stabilization:

Possible Failure Modes¶

Inaccurate Control Reconstruction¶

Piecewise-linear control interpolation may be poor near rapid changes or discontinuities.

Possible remedies:

• use one-sided interpolation near switches;

• place mesh points at detected discontinuities;

• use shape-preserving interpolation;

• refine before estimating.

Error Cancellation¶

The interpolated and integrated states may agree by coincidence.

Use multiple checks, including:

• independent forward integration;

• defect residuals on a dense grid;

• objective convergence;

• polynomial coefficient decay.

Repeated Refinement without Convergence¶

Possible causes include:

• incorrect dynamics;

• poor variable scaling;

• inconsistent boundary conditions;

• an under-resolved discontinuity;

• inaccurate derivatives;

• a defective NLP solution.

Practical Recommendations¶

1. Start with a modest number of intervals.

2. Use low-to-moderate polynomial degree.

3. Set a reasonable maximum degree.

4. Normalize state errors using physical scales.

5. Prefer $h$-refinement when nonsmoothness is suspected.

6. Prefer $p$-refinement only in clearly smooth intervals.

7. Warm start every refined solve.

8. Verify the final solution independently.

9. Monitor both accuracy and NLP size.

10. Retain the full refinement history.

Refinement History Table¶

A useful refinement log contains:

This table helps identify stagnation, excessive node growth, and objective instability.

Summary¶

Adaptive mesh refinement converts a fixed-mesh pseudospectral method into an accuracy-controlled computational procedure.

The principal steps are:

1. solve the NLP on a current mesh;

2. create an enriched LGR grid;

3. interpolate the state to that grid;

4. reconstruct the control;

5. integrate the dynamics on the enriched grid;

6. compare interpolated and integrated state values;

7. accept accurate intervals;

8. apply $h$-refinement to localized errors;

9. apply $p$-refinement to distributed smooth errors;

10. enforce a maximum polynomial degree;

11. interpolate the old solution onto the new mesh;

12. repeat until convergence.

Pure $h$ methods are robust but may create large NLPs. Pure $p$ methods are efficient for smooth solutions but perform poorly near discontinuities. An $hp$ method seeks to combine robustness with compactness by placing resolution only where it is needed.