Lagrange multipliers $\lambda$ in Nonholonomic systems and Costates $\lambda$ in Optimal control theory In Structure and Interpretation of Classical Mechanics, Section 1.10.3 by Sussman & Wisdom, talking about systems with non-holonomic constraints, it says the following about the Lagrange equations that we get when we augment the Lagrangian with a non-holonomic constraint:

An interesting feature of these equations is that they involve both $λ$ and $Dλ$. Thus the usual state variables $q$ and $Dq$, with the constraint, are not sufficient to determine a full set of initial conditions for the derived Lagrange equations; we need to specify an initial value for $λ$ as well.

It then goes on to say that this system is "is not fully determined".
These lambdas seem very similar to the "costates" that show up when applying indirect methods for optimal control. In that case, we get costates (also represented by $λ$'s) that have their own ODEs associated with them. Finding the optimal solution involves root-solving for the valid initial conditions of these lambdas along with any unbounded states. While it is not easy, it is certainly possible to start with some trivial initial guess for a short trajectory and then over a series of steps "grow" the trajectory to the one we want.
Can a similar approach work for propagating the dynamics of non-holonomic systems? Can we use something like single-shooting or multiple-shooting method with some initial guess to figure out the "correct" initial value for the lambdas and then numerically integrate the dynamic equations to see how the system evolves?
 A: *

*Sussman & Wisdom are precisely mentioning the 2 possible first principles for affine non-holonomic constraints (=semi-holonomic constraints) in the Lagrangian formalism, namely

*

*d'Alembert's principle (DAP), and


*stationary action principle (SAP),
which happen to be mutually incompatible, cf. e.g. this related Phys.SE post.


*In this answer, we will from now on assume the SAP:
$$\tilde{L}(q,\lambda,\dot{q},t)~=~L(q,\dot{q},t)+\lambda^T \underbrace{f(q,\dot{q},t)}_{\text{non-holonomic}},$$
which makes the Lagrange multipliers $\lambda$ dynamical variables because $\dot{\lambda}$ appear in the corresponding Euler-Lagrange (EL) equations. This means that one must impose adequate boundary conditions (BCs) for the $\lambda$ variables in order for the SAP to be a mathematically well-posed variational principle. The choice of BCs is usually viewed as part of the first principle. Sometimes there's a natural choice of physically relevant BCs.


*As OP already observes the Lagrange multipliers $\lambda$ are indeed closely related to the costates $\lambda$ in optimal control theory (OCT), which in turn play the role of momenta in the Hamiltonian formulation. However, since OCT is a broad research topic, we do not have anything definite to say beyond that.
