Error Checking and Mesh Refinement

Why check a collocation solution?

Pockit represents states and controls with piecewise interpolation polynomials. A nonlinear-programming solution can satisfy equations at collocation nodes while still having excessive interpolation error or a between-node constraint violation. Mesh checks estimate whether the current polynomial representation has enough resolution; they do not replace physical validation.

Two error mechanisms require different treatment:

• A genuine jump should coincide with a Radau mesh boundary. Moving or adding a boundary is more useful than raising polynomial degree across the jump.
• A continuous but under-resolved function needs more intervals, more points per interval, or both.

Check a discontinuous control

Lobatto interpolation is continuous across shared mesh nodes, so Pockit's discontinuity workflow applies only to Radau phases. Mark the relevant bounded expression as bang-bang in Phase.set_phase_constraint() before solving. After a successful single-phase solve with no static parameters, check it with:

discontinuous_ok = system.check_discontinuous(
    solution,
    tolerance_discontinuous=1.0e-3,
    tolerance_mesh=1.0e-4,
)

Within each interval, the marked expression is scaled to ([0,1]). The check passes when it stays within tolerance_discontinuous of one bound throughout that interval. Intervals shorter than tolerance_mesh are skipped.

If the check fails, update the mesh and use the interpolated result as the next initial guess:

from pockit.optimizer import ipopt

solution = system.refine_discontinuous(
    solution,
    tolerance_discontinuous=1.0e-3,
    num_point_min=6,
    num_point_max=12,
    mesh_length_min=1.0e-3,
    mesh_length_max=1.0,
)
solution, info = ipopt.solve(system, solution)

For a system with multiple phases and static parameters, pass one phase result per phase followed by the static values. The returned list has the same layout:

values = [solution_0, solution_1, static_values]
values = system.refine_discontinuous(values)

The names in this short layout example refer to results already returned by a multi-phase solve.

Singular-arc limitation

This check assumes that every marked expression approaches its lower or upper bound. A true singular control can remain strictly inside both bounds, so check_discontinuous() and refine_discontinuous() cannot identify, solve, or validate a singular arc. See Singular Arcs: Identification and Current Limitations.

Check a continuous approximation

For continuous regions, Pockit follows the ph mesh-refinement method of Patterson, Hager, and Rao:

Patterson, M. A., W. W. Hager, and A. V. Rao (2015). "A mesh refinement method for optimal control." Optimal Control Applications and Methods, 36(4), 398-421.

Check a single-phase result with explicit absolute and relative tolerances:

continuous_ok = system.check_continuous(
    solution,
    absolute_tolerance_continuous=1.0e-8,
    relative_tolerance_continuous=1.0e-8,
    tolerance_mesh=1.0e-4,
)

If more resolution is needed, refine and solve again:

solution = system.refine_continuous(
    solution,
    absolute_tolerance_continuous=1.0e-8,
    relative_tolerance_continuous=1.0e-8,
    num_point_min=6,
    num_point_max=12,
    mesh_length_min=1.0e-3,
    mesh_length_max=1.0,
)
solution, info = ipopt.solve(system, solution)

num_point_min and num_point_max limit the points assigned to an interval. mesh_length_min and mesh_length_max limit normalized interval length. The refinement method does not reduce the number of mesh points.

Handle both error types

System.check() and System.refine() combine the two workflows. At most one refinement is performed per call: unresolved discontinuities take priority; otherwise, the continuous-error estimate controls the update.

if not system.check(
    solution,
    absolute_tolerance_continuous=1.0e-8,
    relative_tolerance_continuous=1.0e-8,
    tolerance_discontinuous=1.0e-3,
):
    solution = system.refine(
        solution,
        absolute_tolerance_continuous=1.0e-8,
        relative_tolerance_continuous=1.0e-8,
        tolerance_discontinuous=1.0e-3,
        num_point_min=6,
        num_point_max=12,
        mesh_length_min=1.0e-3,
        mesh_length_max=1.0,
    )
    solution, info = ipopt.solve(system, solution)

Repeat only after checking solver status at each iteration, and impose a clear maximum number of updates. The hyper-sensitive example shows a complete adaptive loop for narrow endpoint boundary layers.

Validate after refinement

A passing mesh check is evidence about the interpolation, not a certificate of optimality or physical feasibility. Reconstruct states, controls, and path constraints on a dense physical-time grid; inspect left and right limits at jumps; and compare objectives and engineering diagnostics across at least two discretizations. Independently integrate the original dynamics when practical.