I'm using the constraint equation H0(t_cf) = 0 at the endpoint of the coast as the condition for terminating the coast. This only works when the coast is before the terminal stage and the target conditions are free-final time. The free-final time condition means that H(tf) = 0 and as the burn starts and ends at an optimum time and the mass is constant across the coast-burn transition that condition applies across the coast as well. As H = H0 + TS and S = 0 at the switching point implies that H0 = 0 at the termination of the coast.

This is invalidated if the burntime is fixed and H != 0. This is also invalidated if there is a mass jettison after the coast (multiple upper stages after the coast). I have no idea how it seems to be kinda sorta working correctly right now in the fixed-final-time and multiple-unguided-upper-stages approaches.

It also isn't possible to use the |pv(t_ci)| = |pv(t_cf)| condition across a coast that initiates with a mass jettison, that condition is only valid for a coast that begins and end at an optimum time. That arises because even if the hamiltonian isn't zero, it is a fixed value across a coast so H(t_ci) = H(t_cf) that means that H0(t_ci) + T(t_ci)S(t_ci) = H0(t_cf) + T(t_cf)S(t_cf). You have to use S(t_ci) = S(t_cf) across an optimally starting and ending coast, along with the the mass and mass costate being constant across a coast to arrive at the expression |pv(t_ci)| = |pv(t_cf)|. That isn't valid if the start of the coast isn't optimal and happens after a stage burns out (although it doesn't require that the hamiltonian is zero so is useful for fixed-final-time and later-jettisoned-coasts).

The solution to all of this is to use S(t_cf) = 0 directly. This, however, means integrating the mass costate. For stability it necessary to integrate this backwards from Pm(t_f) = 1 condition at the end of the burn (although fixed-final-time maximum-energy burns probably need to use Pm(t_f) = 0 as the terminal value since Pm(t_f) = 1 converts the problem into a J = -m_f optimization problem, but that is advantageous for certain problems like stage-and-a-half anyway). The backwards integration probably needs to be done with RK using actual dense output and storing that for each phase since the integration of Pm depends on the instantaneous value of Pv. Then the phases replay backwards and RK is used to integrate the mass costate. Analytical expressions for Pv could be used in the analytical computations. It is possible that acceptable results could be found by using analytical expressions for Pv in the backwards integration in the RK case? That would certainly be easier and would be work that needs to get done anyway.

One paper that utilizes the mass coast and maximum-final-mass problem (although doesn't handle stage jettison):

• Murillo, Oscar, and Ping Lu. “Fast Ascent Trajectory Optimization for Hypersonic Air-Breathing Vehicles.” In AIAA Guidance, Navigation, and Control Conference. Toronto, Ontario, Canada: American Institute of Aeronautics and Astronautics, 2010. https://doi.org/10.2514/6.2010-8173.

At the mass jettison events there will be a discontinuous jump in the mass costate. I don't know the expression for this yet, and some papers indicate you need to include a new adjoint variable to account for the jump (and use continuity of the hamiltonian across the jump as a new constraint?). Maybe a closed form solution can be found by some algebra? It should be possible to validate that the answer is correct by taking existing solutions found without integrating the mass costate and validating that they produce the correct answers when using the mass costate and the jump conditions and that the new value of the hamiltonian is constant across the mass jettison events.

One paper that mentions jump discontinuities in the mass costate:

• Colasurdo, Guido, Dario Pastrone, and Lorenzo Casalino. “Optimization of Rocket Ascent Trajectories Using an Indirect Procedure.” In Guidance, Navigation, and Control Conference. Baltimore,MD,U.S.A.: American Institute of Aeronautics and Astronautics, 1995. https://doi.org/10.2514/6.1995-3323.

Furthermore, the calculation of the optimum coast length should be moved out of being an optimizer variable and computed on every function evaluation of the optimizer through a root finding problem. There are a series of papers that go into this in depth, but unfortunately only for the case of an optimal coast, not for the case of a coast-after-jettison. To start with it may be sufficient to get much better results by doing a search for every 10% step between coast_min and coast_max times for changes in the switching condition. The right solution will likely result in having to replicate some of the papers here for the mass jettison case:

• Xu, Yunjun. “Enhancement in Optimal Multiple-Burn Trajectory Computation by Switching Function Analysis.” Journal of Spacecraft and Rockets 44, no. 1 (January 2007): 264–72. https://doi.org/10.2514/1.25082.
• Jamison, Brian Roger, and Victoria Coverstone. “Analytical Study of the Primer Vector and Orbit Transfer Switching Function.” Journal of Guidance, Control, and Dynamics 33, no. 1 (January 2010): 235–45. https://doi.org/10.2514/1.41126.
• Pan, Binfeng, Ping Lu, and Zheng Chen. “Coast Arcs in Optimal Multiburn Orbital Transfers.” Journal of Guidance, Control, and Dynamics 35, no. 2 (March 2012): 451–61. https://doi.org/10.2514/1.54655.
• Pan, Binfeng, Ping Lu, and Zheng Chen. “Three-Dimensional Closed-Form Costate Solutions in Optimal Coast.” Acta Astronautica 77 (August 2012): 156–66. https://doi.org/10.1016/j.actaastro.2012.04.009.

And finally, a general good enhancement to the overall PVG architecture would be to validate the solution results against running it through the SQP solver, without the transversality conditions and with an explicit cost function and determining that the results agree. If the PVG-solved value is optimal then the SQP solver should quickly exit with the same answer and stability shouldn't be an issue.

Pretty sure that due to time not being a variable part of the 'constraint' on the mass jettison that the Hamiltonian is constant across the mass jettison. That means that H+ - H- = 0 (eqn 48 in Colasurdo). So there's no jump condition in the Hamiltonian. There is a jump condition in the mass costate which can be assigned a variable mu so that Pm+ = Pm- + mu. That latter equation replaces the continuity condition on the mass costate across the stage jettison. The condition that the Hamiltonian be constant is the extra condition due to the introduction of mu. Due to continuity of the other states and costates we have that H0+ - H0- = 0 automatically so that T+S+ - T-S- = 0 becomes our condition for mass jettison between two burn arcs. For the case of a coast after a mass jettison then T+ = 0 and T- != 0 so that recovers that S- = 0 is the constraint.

That suggests that the mass costate for the booster stages before a coast just doesn't matter since you can introduce a mass costate jump whose only constraint is that it makes the switching function zero there. The initial mass costate for the booster stages before a coast could be set to zero and the continuity condition on the mass costate could be switched to set the terminal mass costate of that stage to be zero with the mass costate for the next stage left free. Actual integration of the mass costate could start with the [first] coast stage.

Then the question of backwards integration happens, but the initial mass costate can be selected to be any number > 0 which will serve to scale the overall costate of the problem (similar to the Costate.magnitude = 1 constraint which much be replaced due to needing to use the H(tf) = 0 constraint). It may be possible to select a mass costate guessed value which works well with the guess of the rest of the costate for the given rocket and target, or to force the mass costate to some constant value and scale the rest of the costate guess to get it within the radius of convergence. Since we don't need the mass costate for the initial coastless + infinite ISP guess generation that solution might be useful to bootstrap the initial mass costate guess close enough to be within the radius of convergence.

So with unguided upper stages there's a condition for the termination of the coast-before-upper-unguided-stages, which I cannot figure out at all. For a given problem, I can calculate the optimal length of the coast in an outer loop. What it looks like I need then is a constraint on H0 == <small nonzero value> at the start or end of the coast. And I failed to figure out what that value would be or how to go about it some other way.

The only way forward I think is through converting the problem to a hybrid problem, dropping the transversality condition on the coast length and feeding the appropriate metric (J = -m_f or whatever) into an SQP solver or something of that nature instead of using a simple root finder like I have been using (Levenburg-Marquardt). Experiments with that don't work terribly well, though, and I suspect I need better Jacobian information than just finite differencing and need to care about smoothness of the problem a lot more.