I'm trying to solve a trajectory optimization problem for a class of problems like the old Altas-Centaur SLV3 Centaur launch vehicles. Those are a stage and a half design where the 2 LR-87 engines are dropped at an optimized time and the rocket continues on the LR-105 sustainer. Since this is an optimized staging time which has two different burn phases on either side of it rather than a burn-coast transition -- with a mass/thrust/isp discontinuities -- none of the typical mathematical tricks to eliminate the integration of the mass costate applies. I believe the way to solve this (?) is that the mass costate should be integrated and the Weierstrauss-Erdmann corner condition applied to the continuity of the Hamiltonian across the staging time. However without accounting for the discontinuity in the mass costate there is a discontinuity in the Hamiltonian and so this constraint cannot be applied. The question is how to calculate the discontinuity so to be able to apply this condition to solve the optimization problem?
I have solved this via wrapping a fixed time optimization problem with a line search for the optimized time and validated that my problem has a reasonable optimum value. I've also validated that other than the discontinuities the calculation of the full Hamiltonian is stepwise-constant given my integration of the mass costate. For varying fixed staging time trajectories around the optimum solution the discontinuity in the Hamiltonian value changes.
The approach I'm taking is similar to that in e.g. Lu, et al 2008 although I'm only solving the vacuum problem and using an ODE integrator instead of using the analytic solutions of the linearized gravity problem. The same numerical conditioning is applied so that $g_0 = \mu / R_0^2$ and the distances are scaled by $R_0$ the velocities by $\sqrt{R_0 g_0}$ and time by $\sqrt{R_0 / g_0}$. So I am minimizing the integrated thrust:
$$J = - \int_{t_0}^{t_f} \frac{T}{c} dt $$
With the Hamiltonian:
$$ \begin{align} H &= P_r^T V - p_V^T \frac{r}{r^3} + p_V^T 1_T A_T - p_m \frac{T}{c} - \frac{T}{c} \\ &= P_r^T V - p_V^T \frac{r}{r^3} + T \left( \frac{p_V^T 1_T}{m g_0} - \frac{p_m}{c} - \frac{1}{c} \right) := H_0 + T S \end{align} $$
Note this is different from equation 10 in the reference above due to not making the linearized gravity approximation (which should not matter). For most typical burn-coast problems we can write $H_0^{-} + T^{-} S^{-} = H_0^+ + T^+ S^+$ and we can use the constancy of $H_0$ across a coast and that one of $T^{-}$ or $T^{+}$ are zero to simplify the constraints. In this case neither side is a coast so $H_0$ is not constant anywhere and T is also not zero on either side of the corner.
The integration of the mass costate and the terminal constraint for the free final time problem are:
$$ \begin{align} p_m^{'} &= \frac{T \left| p_V \right|}{m^2 g_0} \\ p_m(t_f) &= 0 \end{align} $$
The rest of the problem of the integration of the state and other costates I'll omit, but examples are in the above paper.
If the times $t_0, t_1, t_2, t_f$ correspond to launch, jettison of the engines, jettison of the atlas stage and terminal (free) insertion into the orbital conditions. Then I am trying to use the constraint:
$$ H(t_1)^- = H(t_1)^+ + \Delta H $$
I can solve that via substitution, but that does not constrain the problem, it is just a tautology. I need to find that $\Delta H$ via other means. Note there's also a discontinuity at $t_2$ due to mass jettison as well, but that time is not optimized, and is fixed by the choice of $t_1$ and the constraint of running out of fuel for the sustainer stage.