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?