This question is a follow-up to this question and is linked to this question, where if needed more details can be found on the parameters used in the following.
For an in-plane non-impulsive orbital maneuver, I'd like to find the thrust-direction history $\beta(t)$ to minimize the time needed to transfer the spacecraft from a specified initial state to a specified terminal state, subject to the following state equations:
$$\dot{\textbf{x}}(\textbf{x},\beta,t)=\begin{bmatrix}\dot{r} & \dot{V}_r & \dot{V}_{\theta}\end{bmatrix}$$ $$\dot{r}=v_r$$ $$\dot{v}_r=\frac{v_{\theta}^2}{r}-\frac{\mu}{r^2}+\frac{T\sin{\beta}}{m_0-|\dot{m}|t}$$ $$\dot{v}_{\theta}=-\frac{v_{r}v_{\theta}}{r}+\frac{T\cos{\beta}}{m_0-|\dot{m}|t}$$
where $\mu$, $T$, $m_0$ and $\dot{m}$ are constants, with state variables:
$$\textbf{x}(t)=\begin{bmatrix}r & V_r & V_{\theta}\end{bmatrix}$$
and initial and terminal conditions:
$$\begin{align*} \textbf{x}(t_0) &=\begin{bmatrix}1745100 & 0 & 2027\end{bmatrix}\\ \textbf{x}(t_f) &=\begin{bmatrix}1765100 & 0 & 0\end{bmatrix}\\ \end{align*}$$
i.e. we deal with a two-point boundary value problem where a spacecraft transfers from an initial position to a higher position, while decreasing its initial velocity to zero terminal velocity. The problem is cast as an optimal control problem for which I have formulated the following minimum-time cost functional:
$$J=\min_{\beta}\int_{t_0}^{t_f}L(\textbf{x},\beta,t)dt=\min_{\beta}\int_{t_0}^{t_f}(1)dt$$
where $L$ denotes the Lagrangian. The Hamiltonian is, therefore,
$$\begin{align*} H(\textbf{x},\beta,t) &= L(\textbf{x},\beta,t)+\lambda^T\dot{\textbf{x}}(\textbf{x},\beta,t) \\ &= 1+\lambda_{r}\dot{r}+\lambda_{v_r}\dot{v}_r+\lambda_{v_{\theta}}\dot{v_{\theta}} \\ &= 1+\lambda_{r}\cdot v_r+\lambda_{v_r}\cdot (\frac{v_{\theta}^2}{r}-\frac{\mu}{r^2}+\frac{T\sin{\beta}}{m_0-|\dot{m}|t})+\lambda_{v_{\theta}}\cdot (-\frac{v_{r}v_{\theta}}{r}+\frac{T\cos{\beta}}{m_0-|\dot{m}|t}) \end{align*}$$
where $\mathbf{\lambda}(t)$ is the Lagrange multiplier vector, whose elements are the costate variables. The $\lambda$-dynamics that capture the behavior of the Lagrange multipliers, also known as the costate equations, are given by:
$$\begin{align*} \dot{\lambda}_{r} =-\frac{\delta H}{\delta r} &= -\lambda_{v_r} \cdot (-\frac{v_{\theta}^2}{r^2}+\frac{2\mu}{r^3}) - \lambda_{v_{\theta}} \cdot (\frac{v_{r}v_{\theta}}{r^2})\\ \dot{\lambda}_{v_r} =-\frac{\delta H}{\delta v_r} &= -\lambda_{r} + \lambda_{v_{\theta}} \cdot (\frac{v_{\theta}}{r})\\ \dot{\lambda}_{v_{\theta}} =-\frac{\delta H}{\delta v_{\theta}} &= -\lambda_{v_r} \cdot (\frac{2v_{\theta}}{r}) + \lambda_{v_{\theta}} \cdot (\frac{v_r}{r})\\ \end{align*}$$
The partial of $H$ with respect to $\beta$ must equal zero, so
$$\frac{\delta H}{\delta \beta}=\lambda_{v_r} \cdot (\frac{T\cos{\beta}}{m_0-|\dot{m}|t}) - \lambda_{v_{\theta}} \cdot (\frac{T\sin{\beta}}{m_0-|\dot{m}|t})=0$$
leading to the following control law:
$$\tan{\beta}=\frac{-\lambda_{v_r}}{-\lambda_{v_{\theta}}}$$
or,
$$\sin{\beta}=\frac{-\lambda_{v_r}}{\sqrt{\lambda_{v_r}^2+\lambda_{v_{\theta}}^2}} \qquad \cos{\beta}=\frac{-\lambda_{v_{\theta}}}{\sqrt{\lambda_{v_r}^2+\lambda_{v_{\theta}}^2}}$$
Note: the signs in the ratio are negative because we are looking for a minimum.
This is the first time I'm working with an optimal control problem, and I have some trouble understanding how to proceed from here and solve the problem with MATLAB. Is there perhaps someone out here who knows how to do this and could explain what the next steps would be?
