Shock Capturing and Mesh Refinement

Introduction

Pockit approximates functions using piecewise linear polynomials. Radau nodes allow discontinuities at mesh points (subinterval boundaries); Lobatto nodes always produce continuous approximations within the entire time interval. To accurately fit functions, both discontinuities and continuous parts need to be considered:

For all discontinuities, there should be mesh points as close to them as possible.
For continuous parts, the mesh should be dense enough, and the number of interpolation points should be sufficient to keep the approximation error within acceptable limits.

Pockit can estimate the errors in both discontinuous and continuous parts and automatically adjust the number of mesh points and interpolation points based on the current calculation results.

Shock Capturing

Lobatto nodes cannot accurately approximate discontinuous functions; the methods in this section apply only to Radau nodes.

For Radau nodes, pockit checks the accuracy of discontinuity point estimates as follows: For phase constraints set as bang-bang control in Phase.set_phase_constraint(), it checks whether they are all close to the lower or upper bounds within each subinterval. The relative error tolerance_discontinuous is user-defined, with a default value of (10^{-3}). To test whether the current calculation results meet the above conditions, use the System.check_discontinuous() method:

# The system has n_p Phases, with corresponding Variables v_0, ..., v_{n_p - 1}
# The list or np.array of system variable values is v_sp
S.check_discontinuous(
    [v_0, ..., v_{n_p - 1}, v_sp], # If there is only one phase and no system variables, v can be passed directly
    tolerance_discontinuous,
    tolerance_mesh, # Subintervals shorter than tolerance_mesh are skipped
)

If the check fails, pockit can automatically estimate the positions of discontinuities and add or move mesh points to make them as close to the discontinuities as possible. This is done using the System.refine_discontinuous() method:

# The system has n_p Phases, with corresponding Variables v_0, ..., v_{n_p - 1}
# The list or np.array of system variable values is v_sp
[v2_0, ..., v2_{n_p - 1}, v2_sp] = S.refine_discontinuous(
    [v_0, ..., v_{n_p - 1}, v_sp], # If there is only one phase and no system variables, v can be passed directly
    tolerance_discontinuous,
    num_point_min, # Minimum number of interpolation points per subinterval
    num_point_max, # Maximum number of interpolation points per subinterval
    mesh_length_min, # Minimum subinterval length (assuming total length is 1)
    mesh_length_max # Maximum subinterval length (assuming total length is 1)
)
# The returned v2_* are the interpolations of v_* at the new interpolation points, which can be used as input for the next optimization round

Note: Pockit currently cannot estimate discontinuities in optimal control problems with singular arcs.

For continuous parts, pockit uses the ph method to estimate errors and refine the mesh. The method is referenced from:

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

The current results can be checked using the System.check_continuous() method:

# The system has n_p Phases, with corresponding Variables v_0, ..., v_{n_p - 1}
# The list or np.array of system variable values is v_sp
S.check_continuous(
    [v_0, ..., v_{n_p - 1}, v_sp], # If there is only one phase and no system variables, v can be passed directly
    atol,
    rtol,
    tolerance_mesh, # Subintervals shorter than tolerance_mesh are skipped
)

If the check fails, pockit can automatically estimate the approximation error and increase the number of mesh points or interpolation points. This is done using the System.refine_continuous() method:

# The system has n_p Phases, with corresponding Variables v_0, ..., v_{n_p - 1}
# The list or np.array of system variable values is v_sp
[v2_0, ..., v2_{n_p - 1}, v2_sp] = S.refine_continuous(
    [v_0, ..., v_{n_p - 1}, v_sp], # If there is only one phase and no system variables, v can be passed directly
    atol,
    rtol,
    num_point_min, # Minimum number of interpolation points per subinterval
    num_point_max, # Maximum number of interpolation points per subinterval
    mesh_length_min, # Minimum subinterval length (assuming total length is 1)
    mesh_length_max # Maximum subinterval length (assuming total length is 1)
)
# The returned v2_* are the interpolations of v_* at the new interpolation points, which can be used as input for the next optimization round

Handling Discontinuous and Continuous Parts Simultaneously

Pockit also provides the System.check() and System.refine() methods to handle both discontinuous and continuous parts simultaneously.

These methods check the errors in both discontinuities and continuous parts. If the discontinuity error is not within the acceptable range, the System.refine_discontinuous() method is called, and the mesh refinement for continuous parts is skipped. If the discontinuity error is within the acceptable range, but the continuous part error is not, the System.refine_continuous() method is called. Since ph mesh refinement does not reduce the number of mesh points, continuous refinement will not cause mesh points to move away from the function discontinuities.