How to numerically solve an optimal control problem with an indirect method (in MATLAB)?

This question was also posted at the Space Exploration StackExchange, but I thought it might have a better chance of receiving an answer here.

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 a 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 makes a transfer 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, 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?

Thanks a lot!

In any case, you would need a relatively good initial guess at the adjoint variables, which is not easy at all. I suggest that you first try to solve the problem as a direct optimization problem by minimizing $J$ while satisfying the end-point constraints. Then you can try to recover the adjoint variables from the suboptimal control $\bar\beta(t)$ and use them as the initial guess.