Here The basic theory of UPG and AAPDG will be explained.
UPG basically use a linear gravity assumption, which makes it easy to compute final state.
The equations of the system with linear gravity could be written as:
$$
\dot{\boldsymbol{r}} = \boldsymbol{v}\
\dot{\boldsymbol{v}} = A\boldsymbol{r} + \boldsymbol{a_T}\
\dot{m} = -\frac{T}{v_{ex}}
$$
where $A=-\mu/R_0^3$is the approximated gravity coefficient,
The powered descent goal can be written as such
$$
\boldsymbol{r}f = \boldsymbol{r}{target}\
\boldsymbol{v}f = \boldsymbol{v}{target}\
$$
and we want to use as little fuels as possible, that is, to minimize
$$
J_1 = \int_0^{t_f} \frac{T}{v_{ex}} \mathrm{d}t
$$
This is called the pinpoint landing problem. However, this target can be difficult to solve, and therefore we can use a less strict contrain, which only requires the craft to land somewhere near the target. The new contraints are:
$$
|\boldsymbol{r}f| = |\boldsymbol{r}{target}|\
\boldsymbol{v}f = \boldsymbol{v}{target}\
$$
and to minimize
$$
J_2 = \kappa (\boldsymbol{r}f - \boldsymbol{r}{target})^T(\boldsymbol{r}f - \boldsymbol{r}{target}) + \int_0^{t_f} \frac{T}{v_{ex}} \mathrm{d}t
$$
where
To solver the above mentioned landing problems, one needs to calculate the final state of the craft when time reaches
It is known that with the optimal control theory, we have Hamiltonian
To solve the landing problems, we need to find the root of the system, with 7 variables $\boldsymbol{z}= \begin{pmatrix} \boldsymbol{p}_r\ \boldsymbol{p}_v\ t_f \end{pmatrix}$. The target contraints have given 6 equations for pinpoint landing and 4 constraints for bolza landing, so we still need to find 1 more constraint for pinpoint and 3 more for bolza.
The seventh constraint comes from the transversality conditions, which is
For bolza landing, the last 2 equations comes from the transversality conditions associated with the terminal constraints,
$$
\boldsymbol{\lambda}(t_f)=-\kappa\frac{\partial \phi}{\partial \boldsymbol{x}_f}+\left[\frac{\partial \boldsymbol{s(\boldsymbol{x}f)}}{\partial \boldsymbol{x}f}\right]^T\boldsymbol{\nu}
$$
where $\phi = [\boldsymbol{r}(t_f)-\boldsymbol{r}{target}]^T[\boldsymbol{r}(t_f)-\boldsymbol{r}{target}]$,
To solve the root finding problem, one can use the well developed algorithms, such as powell's method, lm, etc. It is found that powell's method works well. Note that the final states are analytically calculated, one can write the jacobians with close form to save time for numerically evaluate jacobians. In this part, the Jacobians of the target function are given.
For convenience, we will breifly summarize how the target function computed.\
Given the target state $\boldsymbol{r}^$ and $\boldsymbol{v}^
Firstly, we calculate the state at time
Then we do the same thing for
$$ \boldsymbol{x}_2 = \boldsymbol{\Phi}(t_2 - t_1) \boldsymbol{x}_1 + \boldsymbol{\Gamma}(t_2) , \boldsymbol{I}(t_2,t_1)\
\boldsymbol{x}_f = \boldsymbol{\Phi}(t_f - t_2) \boldsymbol{x}_2 + \boldsymbol{\Gamma}(t_f) , \boldsymbol{I}(t_f,t_2) $$
After we get the final state, we can calculate the 7 constraints:
$$ s_1=\boldsymbol{r}_f^T \boldsymbol{r}_f - \boldsymbol{r}^{T} \boldsymbol{r}^\[0.4em]
s_2=\boldsymbol{v}_{f1} - \boldsymbol{v}^*_1\[0.4em]
s_3=\boldsymbol{v}_{f2} - \boldsymbol{v}^*_2\[0.4em]
s_4=\boldsymbol{v}_{f3} - \boldsymbol{v}^*_3\[0.4em]
s_5=\boldsymbol{p}_{rf}^T \boldsymbol{v}_f -
\boldsymbol{p}_{vf}^T \boldsymbol{r}_f +
\frac{T_{max}}{m_f , g_0} | \boldsymbol{p}_{vf} | -\boldsymbol{1}\[0.4em]
s_6= \boldsymbol{y}1^T \left[ \boldsymbol{p}{rf} + 2\kappa (\boldsymbol{r} _f-\boldsymbol{r}^*) \right]\[0.4em]
s_7=\boldsymbol{y}2^T \left[ \boldsymbol{p}{rf} + 2\kappa (\boldsymbol{r}_f-\boldsymbol{r}^*) \right]\ $$ here the * represent the target positon and velocity.
Now that we have obtained the target function
\frac{\partial \boldsymbol{x}_2}{\partial \boldsymbol{\lambda}} =\boldsymbol{\Phi}(t_2 - t_1) \frac{\partial \boldsymbol{x}_1}{\partial \boldsymbol{\lambda}} + \boldsymbol{\Gamma}(t_2)\frac{\partial \boldsymbol{I}_2}{\partial \boldsymbol{\lambda}}\[0.4em]
\frac{\partial \boldsymbol{x}_f}{\partial \boldsymbol{\lambda}} =
\boldsymbol{\Phi}(t_f - t_2) \frac{\partial \boldsymbol{x}_2} {\partial \boldsymbol{\lambda}} +
\boldsymbol{\Gamma}(t_f)\frac{\partial \boldsymbol{I}_f}{\partial \boldsymbol{\lambda}}
$$
Here
As we use the numerical method to evaluate the thrust integral, it can be written as $$\boldsymbol{i}(\tau)=\begin{pmatrix} {\boldsymbol{1}_p}_v\cos(\tau)\frac{T}{m(\tau)g_0}\ {\boldsymbol{1}_p}_v\sin(\tau)\frac{T}{m(\tau)g_0} \end{pmatrix}$$
$$\begin{aligned}
\boldsymbol{I} &= \int_{t_1}^{t_2} \boldsymbol{i}(\tau)\mathrm{d}\tau\
&=\frac{t_2-t_1}{N}\sum_{j=0}^4 a_j \boldsymbol{i}(t_1+ j\delta)
\end{aligned}$$
Where ${\boldsymbol{1}p}v=\frac{p_v(\tau)}{|p_v(\tau)|}$ is the unit vector of $p_v$, and we have
$$\delta = \frac{t_2-t_1}{4}, N=90, a_0=7, a_1=32, a_2=12, a_3=32, a_4=7$$
for Milne’s rule.
Taking derivatives of this expression is quite simple, remember that
$$\boldsymbol{\lambda}(t)=\boldsymbol{\Phi}(t-t_0)\boldsymbol{\lambda}(t_0)$$
The Jacobian of $I$ can be written as
$$
\frac{\partial \boldsymbol{I}}{\partial \boldsymbol{\lambda}}=
\frac{t_2-t_1}{N}\sum{j=0}^4 a_j
\frac{\partial \boldsymbol{i}(t_1+ j\delta)}{\partial \boldsymbol{\lambda}}
$$
And the Jacobian of $i$ can be written as:
$$
\frac{\partial \boldsymbol{i}(t_1+ j\delta)}{\partial \boldsymbol{\lambda}} =
\frac{T}{m(\tau)g_0}
\begin{pmatrix}
\cos(t_1+ j\delta)I{3\times3}\
\sin(t_1+ j\delta)I_{3\times3}
\end{pmatrix}
\frac{\partial {\boldsymbol{1}_p}_v}
{\partial \boldsymbol{\lambda}}
$$
Now the key part is to calculate the Jacobian
\frac{\partial {{\boldsymbol{1}_p}_v}} {{\partial \boldsymbol{p}_v}} &=\frac{1} {|\boldsymbol{p}_v|} \left( 1-{\boldsymbol{1}_p}_v{\boldsymbol{1}_p}_v^T \right) \equiv \boldsymbol{K} \end{aligned} $$ and $$ \begin{aligned} \frac{\partial \boldsymbol{p}_v} {\partial \boldsymbol{\lambda}} &= \begin{pmatrix} \frac{\partial \boldsymbol{p}_v} {\partial {\boldsymbol{p}_r}_0} & \frac{\partial \boldsymbol{p}_v} {\partial {\boldsymbol{p}v}0} \end{pmatrix}\ &= \begin{pmatrix} -\sin(\delta)\boldsymbol{I}{3\times3} & \cos(\delta)\boldsymbol{I}{3\times3} \end{pmatrix}\equiv \boldsymbol{\Lambda}(\delta)\
\frac{\partial \boldsymbol{p}_r}{\partial \boldsymbol{\lambda}} &=
\begin{pmatrix}
\frac{\partial \boldsymbol{p}_r}
{\partial {\boldsymbol{p}_r}_0} &
\frac{\partial \boldsymbol{p}_r}
{\partial {\boldsymbol{p}_v}_0}
\end{pmatrix}\\
&=
\begin{pmatrix}
\cos(\delta)\boldsymbol{I}_{3\times3} &
\sin(\delta)\boldsymbol{I}_{3\times3}
\end{pmatrix}
\end{aligned}
$$
where
Hence, the jacobian of
\frac{t_2-t_1} {N} \sum_{j=0}^4 a_j \frac{T} {m(\tau)g_0} \begin{pmatrix} \cos(t_1+ j\delta)\boldsymbol{I}{3\times3}\ \sin(t_1+ j\delta)\boldsymbol{I}{3\times3} \end{pmatrix} \boldsymbol{K}(t_1+ j\delta)\boldsymbol{\Lambda}(j\delta) $$
Apart from the Jacobians of the thrust integrals, the state's derivatives respect to
\frac{\partial \boldsymbol{\Gamma}(t_f)} {\partial t_f}= \begin{pmatrix} \cos(t_f)\boldsymbol{I}{3\times3} & \sin(t_f)\boldsymbol{I}{3\times3}\ -\sin(t_f)\boldsymbol{I}{3\times3} & \cos(t_f)\boldsymbol{I}{3\times3} \end{pmatrix} = \boldsymbol{\Phi}(t_f)\[0.4em]
\frac{\partial \boldsymbol{I}(t_f,t_2)}
{\partial t_f}=\boldsymbol{i}(t_f)
$$
Thus, the state's derivatives respect to
After we obtained the Jacobians of the , the jacobian of the target function can be easily calculated.
We write target function
\frac{\partial s_2}
{\partial \boldsymbol{\lambda}}&
\frac{\partial s_2}
{\partial t_f}\\
\vdots & \vdots\\
\frac{\partial s_7}
{\partial \boldsymbol{\lambda}}&
\frac{\partial s_7}
{\partial t_f}
\end{pmatrix} $$ Now we will carefully analyze the 7 equations.
let
$$
\begin{pmatrix}
\frac{\partial \boldsymbol{r}_f}
{\partial \bullet}\[0.5em]
\frac{\partial \boldsymbol{v}_f}
{\partial \bullet}
\end{pmatrix}=
\frac{\partial \boldsymbol{x}_f}
{\partial \bullet}
$$
For
$$\begin{aligned}
\frac{\partial s_1}
{\partial \boldsymbol{\lambda}} &=
2({r_f}_1 \frac{\partial {r_f}_1}
{\partial \boldsymbol{\lambda}}+
{r_f}_2 \frac{\partial {r_f}_2}
{\partial \boldsymbol{\lambda}}+
{r_f}_2 \frac{\partial {r_f}_2}
{\partial \boldsymbol{\lambda}})\
&=2r_f^T \frac{\partial \boldsymbol{r}_f}
{\partial \boldsymbol{\lambda}}\
\frac{\partial s_1}
{\partial t_f}&=
2r_f^T \frac{\partial r_f}
{\partial t_f}
\end{aligned}
$$
For
$$ \begin{pmatrix} \displaystyle \frac{\partial s_2} {\partial \boldsymbol{\lambda}}& \displaystyle \frac{\partial s_2} {\partial t_f}\[1.3em]
\displaystyle \frac{\partial s_3}
{\partial \boldsymbol{\lambda}}&
\displaystyle \frac{\partial s_3}
{\partial t_f}\\[1.3em]
\displaystyle \frac{\partial s_4}
{\partial \boldsymbol{\lambda}}&
\displaystyle \frac{\partial s_4}
{\partial t_f}
\end{pmatrix} = \boldsymbol{I}_{3\times3} \begin{pmatrix} \displaystyle \frac{\partial v_f} {\partial \boldsymbol{\lambda}} & \displaystyle \frac{\partial v_f} {\partial t_f} \end{pmatrix} $$ The last three equations' derivate can be obtained with similar process, here we gave the results without detailed analysis. $$ \frac{\partial s_5} {\partial \boldsymbol{\lambda}}= v_f^T\frac{\partial {\boldsymbol{p}_r}_f} {\partial \boldsymbol{\lambda}}+ {p_r}_f^T\frac{\partial \boldsymbol{v}_f} {\partial \boldsymbol{\lambda}}- r_f^T\frac{\partial {\boldsymbol{p}_v}_f} {\partial \boldsymbol{\lambda}}- {p_v}_f^T\frac{\partial \boldsymbol{r}f} {\partial \boldsymbol{\lambda}}+ \frac{T{max}} {m_f g_0} \frac{{\boldsymbol{p}_v}_f^T} {|{\boldsymbol{p}_v}_f|} \frac{\partial {\boldsymbol{p}_v}_f} {\partial \boldsymbol{\lambda}}\[0.4em]
\frac{\partial s_5} {\partial t_f}= v_f^T \frac{\partial {\boldsymbol{p}_r}_f} {\partial t_f}+ {p_r}_f^T\frac{\partial \boldsymbol{v}_f} {\partial t_f}- r_f^T \frac{\partial {\boldsymbol{p}_v}_f} {\partial t_f}- {p_v}_f^T\frac{\partial \boldsymbol{r}f} {\partial t_f}+ \frac{T{max}} {m_f g_0} \frac{{\boldsymbol{p}_v}_f^T} {|\boldsymbol{p}_v|} \frac{\partial {\boldsymbol{p}v}f} {\partial t_f}+ \frac{T{max}^2} {m_f^2 g_0 v{ex}} |{\boldsymbol{p}_v}_f|\sqrt{\frac{r_0}{g_0}}\[0.4em]
\frac{\partial s_{6,7}} {\partial \boldsymbol{\lambda}}= \left{ [{\boldsymbol{p}_r}_f+2\kappa(\boldsymbol{r}_f-\boldsymbol{r}^*)]^T \boldsymbol{A}_i \frac{\partial \boldsymbol{r}_f} {\partial \boldsymbol{\lambda}}+\boldsymbol{r}_f^T \boldsymbol{A}_i^T ( \frac{\partial {\boldsymbol{p}_r}_f} {\partial \boldsymbol{\lambda}}+ 2\kappa\frac{\partial \boldsymbol{r}_f} {\partial \boldsymbol{\lambda}} ) \right}\[0.4em]
\frac{\partial s_{6,7}}
{\partial t_f}=
\left{
[{\boldsymbol{p}_r}_f+2\kappa(\boldsymbol{r}_f-\boldsymbol{r}^*)]^T \boldsymbol{A}_i
\frac{\partial \boldsymbol{r}_f}
{\partial t_f}+\boldsymbol{r}_f^T \boldsymbol{A}_i^T
(
\frac{\partial {\boldsymbol{p}_r}_f}
{\partial t_f}+
2\kappa\frac{\partial \boldsymbol{r}_f}
{\partial t_f}
)
\right}
$$
where
Now the root finding problem is ready to be solved.
The augmented APDG algorithm is capable to perform landing without suffering solved fail to converge problems.
If we assume the gravity is contant, then the state equations will be
$$ \dot{\boldsymbol{r}} = \boldsymbol{v}\ \dot{\boldsymbol{v}} = \boldsymbol{g} + \boldsymbol{a_T}\ $$ we want to arrive the target postion with target speed in the given time, that is,
$$ \boldsymbol{r}(t_f) = \boldsymbol{r}^\ \boldsymbol{v}(t_f) = \boldsymbol{v}^\ $$ In E-guidance, we assume the acceleration is linear with time,
$$
\boldsymbol{a}_T=\boldsymbol{c_0}+\boldsymbol{c_1}t
$$
coefficients
$$
\boldsymbol{a}T=
-\frac{2}
{t{go}}
[\boldsymbol{v}^-\boldsymbol{v}]+
\frac{6}{t_{go}^2}[\boldsymbol{r}^-\boldsymbol{r}-\boldsymbol{v}t_{go}]-\boldsymbol{g}
$$
where
$$ \boldsymbol{a}T= -\frac{6} {t{go}} [\boldsymbol{v}^-\boldsymbol{v}]+ \frac{12}{t_{go}^2}[\boldsymbol{r}^-\boldsymbol{r}-\boldsymbol{v}t_{go}]+\boldsymbol{a}_T^* $$ It is known that E-guidance is more propellant optimal than APDG, but the latter gives ability to conform the final acceleration condition.
It can be observed that E-guidance and APDG have a similar structure, with different coefficients. If we assume a implicit guidance law
$$
\boldsymbol{a}T= \boldsymbol{a}d
-\frac{k_v}
{t{go}}
[\boldsymbol{v}^*-\boldsymbol{v}]-
\frac{k_r}{t{go}^2}[\boldsymbol{r}^*-\boldsymbol{r}]
$$
and assume the acceleration to be linear, the coefficient of jerk can be eliminate by setting
$$
\boldsymbol{a}T=
\frac{2}
{t{go}}
(1-\frac{k_r}{3})
[\boldsymbol{v}^-\boldsymbol{v}]+
\frac{k_r}{t_{go}^2}[\boldsymbol{r}^-\boldsymbol{r}-\boldsymbol{v}t_{go}]+\frac{1}{6}(k_r-6)\boldsymbol{a}_T^*+\frac{1}{6}(k_r-12)\boldsymbol{g}
$$
when
However, the final acceleration will be different from
$$
\boldsymbol{a}_T(t_f)=(\frac{k_r}{6}-1)\boldsymbol{a}_T^+(\frac{k_r}{6}-2)\boldsymbol{g}
$$
typically we want the final acceleration be $\boldsymbol{a}_T^=-\kappa\boldsymbol{g}$, where
For further discussion see the original paper for details.