Consider two orbits $x(t),\space y(t)$ representing the origin and destination for some spaceflight of interest. These could be, for example, cycloids describing LEO and another orbit circling, say, the moon or Mars.
Suppose I want to leave at $x(t_i)$ and arrive at $y(t_f)$, departing and joining those orbits 'smoothly', so their derivatives $\dot{x}(t_t)$ and $\dot{y}(t_i)$ are also imposing boundary conditions on my trajectory.
To accomplish this, I allow for rocket engines to be fired, most notably while leaving one orbit and decelerating into the other, but also at any time mid-flight. On the other hand, I restrict the total spent fuel along the trip, or the total "delta-v" changes, or some equivalent.
What would be the best way to solve for this trajectory, including the firing pattern of the engine?
What I thought about was devising some adjusted Lagrangian that includes the impulse caused by the engines and adding to it a Lagrange multiplier imposing fuel consumption. This is then used to derive new equations of motion. For example, I would consider an action such as: $$S=\int_{t_i}^{t_f}L[q,\dot{q},F_e,t]dt+\lambda\left(J-\int_{t_i}^{t_f}|F_e(t)|dt \right)$$ Where $q(t)$ is the desired orbit, $F_e(t)$ is the force exerted by the engines and unknown a priori, and $J$ the total impulse supplied by them along the journey$^\dagger$. Now derive EOMs and solve for $q,F_e$.
Eventually, I'd like to explore if this 'formalism' produces solutions (trajectories) that can be inverted, solving either for $J$ (allowing minimization of fuel) or $t_f,t_i$ (minimizing flight time, picking launch window).
Where I am struggling is how to modify said Lagrangian. Should $F_e$ be considered as another generalized coordinate? If so, then a kinetic term for it is needed. Should it be treated as another term in the potential energy? Maybe it should be restated as a force originating from enforcing some constrained orbit? Maybe it's actually another 'Lagrangianesque' function of $q, \dot{q}$ and perhaps $\ddot{q}$ as well?
$^\dagger$Admittedly, the instantaneous thrust is not constrained, but it's a start.
EDIT
Following the reference supplied in the answer below, I managed to make some progress with the problem. This was done by "inverting roles", so now the "Lagrangian" is the engine force and a Lagrange multiplier enforces the Newtonian equation of motion:
$$S=\int_{t_0}^{t_f}dt\left( F^2 + 2\lambda(t)[m\ddot{x}+m \nabla \phi-F] \right)$$
Now, taking the differential of $S$ with regards to $x, F$ and $\lambda$, and doing the routine integration by parts on the variation $\delta\ddot{x}$, we obtain the following equations of motion:
$$(1) \space F=\lambda$$
$$(2) \space m\ddot{x} = -m\nabla\phi +F$$
$$(3) \space \ddot{\lambda} = - H_{\phi}*\lambda$$
Where $H_{\phi}$ denotes the Hessian matrix of $\phi$, i.e. $(H_{\phi})_{ij}=\partial_i\partial_j\phi$.
These equations can know be solved with given boundary conditions $x(t_0),x(t_f),\dot{x}(t_0),\dot{x}(t_f)$, which proved to be pretty straight-forward to do in MATLAB.
The next thing on the agenda is making $t_0, t_f$ free variable subject to optimization as well, so we can find the best launch window. This is proving to be a bit more difficult. Introducing these extra degrees of freedom requires more Lagrange multipliers and produces algebraic equations for each of them and the multipliers. The major difficulty is that the equations of motion above can't be easily integrated and then solved for the time $t$, so they can be evaluated at the launch and landing times. Any suggestions?
