# Optimal Control Approach to [PERSON]'s Problem and Gibbs' Method
## 5 Numerical Simulations  
### Lobatto Pseudospectral Method  
Generally, optimal control problems are very difficult to find an analytical solution for and, thus, numerical methods are widely used. Many different numerical methods can be found in literature [11; 12] and, typically, pseudospectral methods are preferred for their computational efficiency and accuracy. In this paper, LPM is used because LPM can calculate the control values at both end points, which is a required property in [PERSON]'s problem solutions.  
In LPM, the continuous-time optimal control problem is transformed to a nonlinear programming (NLP) problem by discretizing the state and control variables of the dynamic equations at Legendre-Gauss-Lobatto (LGL) points; the cost function is approximated using a Gauss quadrature.  
\begin{table}
\begin{tabular}{c c} \hline \hline  & \(\mathcal{E}_{set}\gets EnergyComputate\left(\overrightarrow{r}_{1}, \overrightarrow{r}_{2},\overrightarrow{r}_{3}\right)\) \\ \hline
1 & \(r_{1}=\sqrt{\overrightarrow{r}_{1}\cdot\overrightarrow{r}_{1}}\), \(r_{2}=\sqrt{\overrightarrow{r}_{2}\cdot\overrightarrow{r}_{2}}\), \(r_{3}=\sqrt{\overrightarrow{r}_{3}\cdot\overrightarrow{r}_{3}}\) \\
2 & \(\Delta\theta_{12}=\text{cos}^{-1}\frac{r_{1}r_{2}}{r_{1}r_{3}}\), \(\Delta\theta_{13}=\text{cos}^{-1}\frac{r_{1}r_{3}}{r_{1}r_{3}}\), \(\Delta\theta_{23}=\Delta\theta_{13}-\Delta\theta_{12}\) \\
3 & \(p=\frac{\sin\Delta\theta_{13}-sin\Delta\theta_{12}-sin\Delta\theta_{23}}{ \left(\frac{sin\Delta\theta_{13}}{r_{2}}-\frac{sin\Delta\theta_{23}}{r_{3}} \right)}\) \\
4 & \(e\sin\theta_{1}=\frac{1+\left(\frac{p}{r_{1}}-1\right)\cos\Delta\theta_{12}- \frac{p}{r_{2}}}{sin\Delta\theta_{12}}\) \\
5 & \(e\cos\theta_{1}=\frac{p}{r_{1}}-1\) \\
6 & \(\theta_{1}=\tan^{-1}\left(\frac{e\sin\theta_{1}}{e\cos\theta_{1}}\right)\) \\
7 & \(e=\frac{\frac{p}{r_{1}}-1}{\alpha\cos\theta_{1}}\) \\
8 & \(\mathcal{E}_{set}=-\frac{\mu\left(1-e^{2}\right)}{2p}\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Orbital energy computation algorithm.  
Since a detailed explanation of LPM can be found easily in literature [11; 12], the final results are summarized here as follows.  
First, the original time interval where \(t\in[t_{1},t_{2}]\) is transformed to the time interval where \(\tau\in[-1,1]\) by the following equation:  
\[\tau=\frac{2t}{t_{2}-t_{1}}-\frac{t_{2}+t_{1}}{t_{2}-t_{1}}. \tag{36}\]  
The state and control variables of Equation (2) are approximated in terms of \(N\) Lagrange interpolating polynomials as:  
\[\overset{\rightarrow}{r}(\tau)\approx\mathbf{R}(\tau)=\sum_{i=1}^{N}\mathbf{R}_{i}L _{i}(\tau) \tag{37}\]  
\[\overset{\rightarrow}{v}(\tau)\approx\mathbf{V}(\tau)=\sum_{i=1}^{N}\mathbf{V}_{i}L _{i}(\tau)\]  
where \(\mathbf{R}_{i}=\mathbf{R}(\tau_{i})\), \(\mathbf{V}_{i}=\mathbf{V}(\tau_{i})\), and \(L_{i}(\tau)\) is defined as:  
\[L_{i}(\tau)=\sum_{\begin{array}{c}k=1\\ k\
eq i\end{array}}^{N}\frac{\tau-\tau_{k}}{\tau_{i}-\tau_{k}} \tag{38}\]  
Next, the dynamic Equation (2) is discretized at the LGL points as:  
\[\sum_{i=1}^{N}D_{ji}\mathbf{R}_{i}=\frac{t_{2}-t_{1}}{2}\mathbf{V}_{j},\ \ \ \ j=1,\ldots,N \tag{39}\]  
where \(D_{ji}=\dot{L}_{i}\big{(}\tau_{j}\big{)}\) is an element of the Lobatto pseudospectral differential matrix. The cost function of Equation (5) is approximated using a Gauss quadrature as:  
\[J=\phi(\mathbf{R}_{N},t_{2})+\frac{t_{2}-t_{1}}{2}\sum_{i=1}^{N}w_{i}\Bigg{(}\frac {\mathbf{V}_{i}^{T}\mathbf{V}_{i}}{2}+\frac{\mu}{\sqrt{\mathbf{R}_{i}^{T}\mathbf{R}_{i}}} \Bigg{)} \tag{40}\]  
where \(w_{i}\) are the LGL weights. Finally, the discretized constraints of Equation (6) are also expressed, as follows:  
\[\psi(\mathbf{R}_{1},t_{1},\mathbf{R}_{N},t_{2})=\mathbf{0} \tag{41}\]  
The cost function of Equation (40), along with the algebraic constraints of Equations (39) and (41), formulates the converted NLP problem; the solution of this NLP is the approximate solution to the original optimal control problem.