# Rough or crumpled: Phases in kinetic growth with surface relaxation

Sudip Mukherjee<EMAIL_ADDRESS>Barasat Government College, 10, KNC Road,
Gupta Colony, Barasat, Kolkata 700124, West Bengal, India  Abhik Basu
abhik.123<EMAIL_ADDRESS>Theory Division, Saha Institute of
Nuclear Physics, 1/AF Bidhannagar, Calcutta 700064, West Bengal, India

###### Abstract

We show that generic kinetic growth processes with surface relaxations can
exhibit a hitherto unexplored crumpled phase with short-range orientational
order at dimensions $d<4$. A sufficiently strong spatially non-local part of
the chemical potential associated with the particle current above a threshold
in the system can trigger this crumpling. The system can also be in a
perturbatively accessible rough phase with long range orientational order but
short range positional order at $d<4$ with known scaling exponents.
Intriguingly, in $d>4$ we argue that there is no crumpling transition;
instead, there is a roughening transition from a smooth to a rough phase for
large enough non-local particle chemical potential. Experimental and
theoretical implications of these results are discussed.

Nonequilibrium fluctuations are far more difficult to characterise
theoretically due to the combined effects of the nonlinear interactions and
the lack of an Fluctuation-Dissipation-Theorem chaikin . The absence of a
general theoretical framework for nonequilibrium systems has prompted
physicists to work on simple nonequilibrium models, which are amenable to
well-controlled calculations, but still produce novel results with wider
ramifications. A prominent example is the Kardar-Parisi-Zhang equation (KPZ)
equation, originally proposed as a simple nonlinear model of nonconserved
surface growth without any overhangs kpz ; stanley , that is reinvented as
paradigmatic nonequilibrium model undergoing nonequilibrium roughening
transition between a smooth and a perturbatively inaccessible rough phase at
dimensions $d>2$ kpz ; stanley ; kpz-rough . A conserved version of the KPZ
(CKPZ) equation driven by a conserved noise ckpz-basic  has been proposed that
unsurprisingly belongs to a different universality class with no roughening
transition at all. A subsequent generalisation of the CKPZ equation having a
non-local part of the chemical potential $\mu$ associated with the particle
current $\bf J$ mike  reveals a more complex phase behaviour including new
growth phases in a parameter regime for $d>1$. A dynamical model that is
closely related to the CKPZ equation is a nonlinear model, originally proposed
for ideal Molecular Beam Epitaxy das-sarma ; stanley  (hereafter the Lai-Das
Sarma or LDS equation) and also subsequently used as a model for tumor growth
tumor . It has the same conservation law structure as the CKPZ equation ckpz-
basic , but is driven by a nonconserved noise stanley ; das-sarma . Due to the
difference in the noise statistics, the CKPZ and LDS equations belong to two
different universality classes. This naturally begs the question whether a
generalisation of the LDS equation similar to the generalised CKPZ equation
shows a similar transition, or something else, or nothing at all.

In this Letter, we investigate the universal scaling in a generalised LDS
model with $\mu$ having both local and non-local parts, analogous to the
generalised CKPZ equation mike . This reduces to the LDS equation das-sarma ,
when the non-local part of $\mu$ vanishes. We show that for weaker non-local
$\mu$, the model belongs to the LDS universality class: At $d<4$, it only has
a rough phase whose exponents can be obtained in systematic renormalised
perturbation expansions; the scaling exponents are identical to those obtained
in Ref. das-sarma . For larger non-local $\mu$, the surface however crumples
at $d<4$, as soon as the system size exceeds a finite threshold, with a
concomitant loss of orientational long range order (LRO). This forms a
hitherto unstudied nonequilibrium analogue of membrane crumpling, an
intriguing phenomenon that is well-studied in the statistical mechanics of
membranes in thermal equilibrium tethered ; david-guitter ; john-tethered ;
john-tethered1 , but not in systems out of equilibrium. For $d>4$, the model
shows a roughening transition from a smooth to rough surfaces, controlled by
the non-local part of $\mu$, with the attendant loss of positional LRO.

We now derive these results. We start with the equation of motion

$\frac{\partial h}{\partial
t}=-\nu\nabla^{4}h-\frac{\lambda}{2}\nabla^{2}({\bm{\nabla}}h)^{2}-\lambda_{1}{\bm{\nabla}}\cdot[({\bm{\nabla}}h)\nabla^{2}h]+\eta,$
(1)

where $\nu>0$ is a damping coefficient, $\eta({\bf x},t)$ is a Gaussian noise,
and $\lambda,\,\lambda_{1}$ are the nonlinear coupling constants whose signs
are arbitrary; $\lambda_{1}=0$ gives the LDS equation das-sarma . Noise
$\eta({\bf x},t)$ is a zero-mean, Gaussian white noise with a variance

$\langle\eta({\bf x},t)\eta(0,0)\rangle=2D\delta({\bf x})\delta(t).$ (2)

Equation (1) has the form $\partial h/\partial t=-{\bm{\nabla}}\cdot{\bf
J}+\eta$, where

${\bf
J}=\nu{\bm{\nabla}}\nabla^{2}h+\frac{\lambda}{2}{\bm{\nabla}}({\bm{\nabla}}h)^{2}+\lambda_{1}({\bm{\nabla}}h)\nabla^{2}h.$
(3)

The first term of the rhs of (1), or equivalently, of (3), is a
“curvature”-dependent part of the particle current at the linear order in $h$
das-sarma . It contributes to (3) only if there is a gradient in the local
mean curvature. Further, the two nonlinear terms with coefficients $\lambda$
and $\lambda_{1}$ are the symmetry-permitted lowest order nonlinear terms,
which are of nonequilibrium origin, and hence cannot be obtained from a free
energy. Writing ${\bf J}=-{\bm{\nabla}}\mu$ we get

$\mu=-\nu\nabla^{2}h-\frac{\lambda}{2}({\bm{\nabla}}h)^{2}-\lambda_{1}\nabla^{-2}{\bm{\nabla}}\cdot[({\bm{\nabla}}h)\nabla^{2}h],$
(4)

giving $\nabla^{2}\mu=-{\bm{\nabla}}\cdot{\bf J}$. Thus, the
$\lambda_{1}$-term acts as a non-local chemical potential, as opposed to the
local chemical potential term with coefficient $\lambda$. The
$\lambda_{1}$-term, originally introduced in Ref. mike  to generalise the CKPZ
equation, makes a contribution to $\bf J$ in (3) that is bilinear in the local
curvature and tilt. In fact, its effect can been interpreted as “blind
geodesic jumping”, in which particles move a small fixed geodesic distance
along the surface in a random direction mike . As illustrated in Ref. mike ,
this potentially leads to a current contribution as in the $\lambda_{1}$-term.
This $\lambda_{1}$-term is responsible for a roughening transition in the
generalised CKPZ equation, absent in the original the CKPZ equation, where
$\lambda_{1}=0$ identically mike .

Universality of the statistical steady states of (1) are described by the
universal spatio-temporal scaling of the time-dependent correlation function
$C$ of $h$:

$C(r,t)\equiv\langle[h({\bf x},t)-h(0,0)]^{2}\rangle\sim
r^{2\chi}f_{h}({r^{z}}/{t}),$ (5)

where, $r=|{\bf x}|$, $\chi_{h}$ and $z$, respectively, are the roughness and
dynamic exponents.

Renormalisation group (RG) analysis on the LDS equation shows that
fluctuations modify the linear theory values of the scaling exponents at
$d<4$; it gives a rough phase with $\chi=\epsilon/3$ and $z=4-\epsilon/3$ to
${\cal O}(\epsilon)$ where $\epsilon\equiv 4-d$ das-sarma . At $d>4$,
fluctuation-corrections are irrelevant and the linear theory values of the
scaling exponents are obtained with $z=4$ and $\chi=(4-d)/2$.

Is the $\lambda_{1}$-term a relevant perturbation on the LDS equation? It has
the same number of fields and gradients as the $\lambda$-term, and hence, both
the terms are expected to be naïvely equally relevant in the RG sense. The
simplest way to ascertain that is by checking their scaling dimensions:
Rescale ${\bf x}\rightarrow b{\bf x},\,t\rightarrow b^{z}t,\,h\rightarrow
b^{\chi}h$ with the linear theory values of $\chi=2-d/2$ and $z=4$, such that
$\nu$ and $D$ do not scale. Under this rescaling, both
$(\lambda,\,\lambda_{1})\rightarrow b^{4-d}(\lambda,\,\lambda_{1})$, showing
their equal relevance in the RG sense with 4 as their common critical
dimension. In fact, there are no symmetry arguments to discard the
$\lambda_{1}$-term, since both of them have the same symmetry being invariant
under $h\rightarrow h+const$. This naturally calls for inclusion of the
$\lambda_{1}$-term in a complete RG study, that will allow us to
systematically explore the “new steady states” due to a non-zero
$\lambda_{1}$, not found in the LDS equation. What is the nature of these new
steady states? The nonlinear terms preclude exact enumeration of the scaling
exponents. We therefore take a perturbative approach. Conservation law and the
invariance under a constant shift of the base plane ensure that the
fluctuations of $h$ are long lived; the life-time of the fluctuations of size
$q^{-1}$ diverges as wavevector $q\rightarrow 0$. As a result, perturbative
corrections to the model parameters diverge in the infra-red limit. These
divergences are systematically handled within the dynamic RG framework stanley
; fns ; halpin ; janssen . The RG procedure is conveniently implemented by
using a path integral description, equivalent to and constructed from Eq. (1)
together with (2), in terms of the field $h({\bf x},t)$, and its dynamic
conjugate field $\hat{h}({\bf x},t)$ janssen .

We perform the one-loop Wilson momentum shell dynamic RG procedure fns ;
stanley ; halpin . Since the $\lambda_{1}$-term in Eq. (1) can be written as a
combination of $\nabla^{2}({\bm{\nabla}}h)^{2}$ and
$\nabla_{i}\nabla_{j}(\nabla_{i}h\,\nabla_{j}h)$, fluctuations do not generate
a $\nabla^{2}h$-term in Eq. (1). The resulting fluctuation-corrections to the
model parameters are represented by the one-loop Feynman diagrams; see
Supplemental Material (SM) sm . There are no relevant fluctuation-corrections
to the noise amplitude $D$, which is an exact statement. This is due to the
fact that the noise in (1) is non-conserved, whereas the equation of motion
(1) has the form of a conservation law. This means the corrections to $D$ is
at ${\cal O}(q^{2})$, or higher, whereas the bare noise amplitude is ${\cal
O}(q^{0})$, ruling out any relevant (in the RG sense) corrections to it.
Furthermore, there are no relevant one-loop corrections to $\lambda$ or
$\lambda_{1}$ janssen-prl . There are however diverging one-loop corrections
to $\nu$. By dimensional analysis we identify two dimensionless effective
coupling constants

$g=\frac{\lambda^{2}D}{\nu^{3}}K_{d}\Lambda^{d-4},\;\gamma=\frac{\lambda_{1}}{\lambda}.$
(6)

Here, $K_{d}$ is the surface area of a $d$-dimensional hypersphere of unit
radius. Further $g$ has a critical dimension $d_{c}=4$ as in Ref. das-sarma .
By following the standard steps of RG outlined in stanley , we obtain the
following differential RG recursion relation to the linear order in
$\epsilon\equiv 4-d$:

$\displaystyle\frac{d\nu}{dl}$ $\displaystyle=$
$\displaystyle\nu[z-4+g\Delta(\gamma)],$ (7)
$\displaystyle\frac{d(\lambda,\,\lambda_{1})}{dl}$ $\displaystyle=$
$\displaystyle(\lambda,\,\lambda_{1})[\chi+z-4],$ (8)
$\displaystyle\frac{dD}{dl}$ $\displaystyle=$ $\displaystyle D[z-d-2\chi],$
(9) $\displaystyle\frac{dg}{dl}$ $\displaystyle=$ $\displaystyle
g[\epsilon-3g\Delta(\gamma)],$ (10)

where,
$\Delta(\gamma)=\left(-\frac{3\gamma^{2}}{8}+\gamma+\frac{1}{8}\right)$,
$b\equiv\exp(l)$ is a length-scale; $\xi=b^{\chi}$, with $\chi$ being the
roughness exponent. With $\lambda_{1}=0$, (7) and (10) reduce to those in the
LDS equation das-sarma . The surface is rough or smooth, if $\chi>0$ or $<0$.
Since $\lambda$ and $\lambda_{1}$ can have either sign independently of each
other, the dimensionless ratio $\gamma$ can also be positive or negative. The
nonrenormalisation of both $\lambda$ and $\lambda_{1}$ at the one-loop order
means their ratio $\gamma\equiv\lambda_{1}/\lambda$ too does not renormalise,
and is marginal: $d\gamma/dl=0$.

Depending upon the signs of $\epsilon$ and $\Delta$, there are four distinct
cases obtained from $dg/dl=0$ at the RG fixed point. We first focus on the
cases with $\epsilon>0$ (i.e., $d<4$), which includes the physically relevant
dimension 2.

(i) If $\epsilon>0$ ($d<4$) and $\Delta>0$, then

$g=0,\;g=\frac{\epsilon}{3\Delta}=g^{*}(\gamma)$ (11)

as the fixed points of $g$.Out of these two fixed points, $g=g^{*}(\gamma)$ is
stable. It is, in fact, a fixed line in the $g-\gamma$ plane, parametrised by
$\gamma$. On this fixed line, by using $d\nu/dl=0$, we get

$\displaystyle z$ $\displaystyle=$ $\displaystyle
4-g\Delta=4-\frac{\epsilon}{3},\;\chi=\frac{z-d}{2}=\frac{\epsilon}{3}.$ (12)

These are identical to their counterparts in the LDS equation (i.e.,
$\gamma=0$) das-sarma , and form the LDS universality class. Thus, at this
stable fixed point, valid when $\Delta>0$ together with $\epsilon>0$, model
equation (1) belongs to the LDS universality class das-sarma  for all
$\gamma$.

(ii) Now consider $\epsilon>0$ ($d<4$) and $\Delta<0$. This gives

$\frac{dg}{dl}=g\left[\epsilon+3g|\Delta|\right].$ (13)

In this case, $g=0$ is the only fixed point, which is unstable. We argue below
that this represents crumpling of the surface.

We note that the above two cases are distinguished by the sign of $\Delta$.
Naturally, $\Delta=0$ gives the separatrix between the two kinds of behaviour,
giving two solutions $\gamma_{\pm}$ for $\gamma$, or $\lambda_{1}$ in terms of
$\lambda$:

$\displaystyle\gamma_{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\bigg{[}4\pm\sqrt{19}\bigg{]}>(<)0,$ (14)
$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{3}\bigg{[}4\pm\sqrt{19}\bigg{]},$ (15)

giving the boundaries of stability in the $\lambda-\lambda_{1}$ plane; see
Fig. 1(a) for a schematic phase diagram and discussions below. For
$\gamma>\gamma_{+}$, or $\gamma<\gamma_{-}$, $\Delta<0$, corresponding to case
(ii), giving instability. In contrast, when $\gamma$ belongs to the other
range, i.e., $\gamma_{+}>\gamma>\gamma_{-}$, $\Delta>0$, we get
$g=g^{*}(\gamma)=\epsilon/(3\Delta)$ as the stable fixed point. As noted
earlier, the stable case (i) belongs to the LDS universality class das-sarma .
The unstable case (ii) is new, and is due to the nonlinear term with
coefficient $\lambda_{1}$ in Eq. (1). Since $\gamma$ is marginal at the one-
loop order, we obtain a fixed line in the $g-\gamma$ plane in the stable case.
In fact from (11), if $\gamma=0$, unsurprisingly the LDS fixed point is
recovered: $g^{*}(\gamma=0)=8\epsilon/3$. For $\gamma_{1}>\gamma\geq 0$, we
note that initially for very small $\gamma>0$, $g^{*}(\gamma)$ decreases as
$\gamma$ rises, but eventually starts to rise as $\gamma$ crosses a threshold
$\gamma_{+}^{*}=4/3$, which maximises $\Delta$, and diverges as
$\gamma\rightarrow\gamma_{+}$ from below. Thus $g^{*}$ has a nonmonotonic
dependence on $\gamma$ when $\gamma_{+}>\gamma\geq 0$. We find that
$g^{*}(\gamma=\gamma_{+}^{*})=19\epsilon/24$, which is the smallest value of
$g^{*}$ possible. In contrast, in the range $\gamma_{-}<\gamma<0$, $g^{*}$
rises monotonically as $\gamma$ decreases, eventually diverging as
$\gamma\rightarrow\gamma_{-}$ from above. Since $\gamma$ is marginal in the
one-loop approximation, within the range $\gamma_{-}<\gamma<\gamma_{+}$ the RG
flow lines run parallel to the $g$-axis towards the “fixed line” described
above. Outside this regions, the flow lines flow parallel to the $g$-axis
towards infinite $g$, indicating breakdown of the perturbation theory. See
Fig. 1(b) for a schematic RG flow diagram in the $g-\gamma$ plane.

We can solve Eq. (13) to get

$g(l)=\frac{\epsilon}{3|\Delta|}\frac{A\exp(\epsilon l)}{1-A\exp(\epsilon
l)}.$ (16)

The constant of integration $A$ can be evaluated by using the “initial”
condition: at $l=0,\,g(0)=g_{0}$. This gives
$A=({g_{0}}/{\epsilon}){3|\Delta|+g_{0}})$. If $\Delta\rightarrow 0$ from
below, we get $A\approx g_{0}(\frac{3|\Delta|}{\epsilon})$, corresponding to
$g(l)\approx g_{0}e^{\epsilon l},$ giving the initial growth of $g(l)$ before
the nonlinear effects become important. For any finite $\Delta<0$, $g(l)$
diverges at a finite $l=l_{c}$ if $1-A\exp(\epsilon l_{c})\rightarrow 0_{+}$,
giving $l_{c}=\frac{1}{\epsilon}\ln A$ as the critical “RG time” in which
$g(l)$ diverges.

We now use (16) to obtain the scale-dependent renormalised $\nu(l)$.
Substituting $g(l)$ in (7), we find

$\nu(l)=\nu_{0}\exp\bigg{[}-\frac{\epsilon}{3}\int dl\frac{A\exp(\epsilon
l)}{1-A\exp(\epsilon l)}+(z-4)l\bigg{]},$ (17)

which clearly shows that $\nu(l)$ vanishes as $l\rightarrow l_{c}\sim{\cal
O}(1)$ from below, or equivalently, as the length-scale dependent $\nu(L)$
vanishes as $L\rightarrow\xi$ from below, where $\xi=a_{0}\exp(l_{c})$;
$a_{0}$ is a small-scale cutoff. We thus identify $\xi$ as the persistence
length, beyond which the system crumples. Of course, the perturbation theory
actually breaks down as soon as $g(l)$ becomes ${\cal O}(1)$.

We now argue that case (ii) represents a crumpled phase, quite distinct from a
rough surface. A rough surface is characterised by (a) a positive $\chi>0$,
which implies $\langle h^{2}({\bf x},t)\rangle$ growing with the system size
$L$, which in turn implies a lack of positional LRO, together with (b) a
finite (i.e., $L$-independent in the limit of large $L$)
$\langle({\bm{\nabla}}h)^{2}\rangle$ meaning orientational LRO. For instance,
in the LDS equation, $\chi=\epsilon/3>0$. Therefore, $\langle h^{2}({\bf
x},t)\rangle\sim L^{2\epsilon/3}$, growing indefinitely with $L$. On the other
hand, $\langle({\bm{\nabla}}h)^{2}\rangle$ remains $L$-independent. In
contrast, a crumpled surface not only lacks positional LRO, it also lacks
orientational LRO. This is reflected in the divergence of
$\langle({\bm{\nabla}}h)^{2}\rangle$ in the limit of large $L$. It is clear
that the LDS rough surfaces are not crumpled, since
$\langle({\bm{\nabla}}h)^{2}\rangle$ remains $L$-independent, which is similar
to the equilibrium tethered membranes in their low temperature ($T$)
orientationally ordered phase tethered . In contrast, finite but sufficiently
large two-dimensional (2D) equilibrium fluid or lipid membranes necessarily
crumple at any non-zero temperature chaikin ; peliti . Thus a rough (but non-
crumpled) and a crumpled surface are distinguished by the presence or absence
of long-range orientational order. While the phenomenon of crumpling should be
rather generic, independent of the conditions of equilibrium, a nonequilibrium
example of crumpling in a fluctuating surface has not yet been found.

We further note that in the crumpled phase of the present study, the
fluctuation-corrections to $\nu$ are negative, meaning $\nu(l)$ vanishes in a
finite RG time $l_{c}=\ln(\xi/a_{0})$; see (17) above. Indeed, at this stage
we can make a formal connection with the crumpling of equilibrium statistical
mechanics of membranes. Notice that $D/\nu_{\text{eff}}$ has the direct
correspondence with $T/\kappa_{\text{eff}}$ of the membrane, where
$\kappa_{\text{eff}}$ is an “effective” bend modulus of the membrane, and
$\nu_{\text{eff}}$ is the fluctuation-corrected, effective $\nu$ of the
surface, calculated in the bare perturbation theory. In the case of fluid
peliti ; tirtha-mem1  or asymmetric tethered john-tethered ; john-tethered1
membranes, $\kappa_{\text{eff}}$ vanishes, and consequently both both $\langle
h^{2}({\bf x},t)\rangle$ and $\langle({\bm{\nabla}}h)^{2}\rangle$ diverge as
soon as the membrane size exceeds a finite threshold, a telltale signature of
crumpling. Likewise here, $\nu_{\text{eff}}$ vanishes and along with it both
$\langle h^{2}({\bf x},t)\rangle$, and $\langle({\bm{\nabla}}h)^{2}\rangle$
diverge, as soon as the system size exceeds a finite threshold, giving a
nonequilibrium crumpling. Clearly, the surface morphology in the crumpled
phase is fundamentally different from that in the rough phase in the LDS
equation.

Thus in the $\lambda-\lambda_{1}$ plane, for $d<4$ (15) gives the boundaries
between regions with a rough phase belonging to the LDS universality class,
and a crumpled phase; see the schematic phase diagram in Fig. 1(a). We further
show the RG flow lines ($d<4$) in the $g-\gamma$ plane in Fig. 1(b), which run
parallel to the $g$ axis. For $\gamma_{-}<\gamma<\gamma_{+}$, flow lines
approach a fixed line, whose each point corresponds to the LDS universality
class. Outside this window, the flow lines run to infinity along the
$g$-direction, indicating crumpling of the surface.

We have found that $g^{*}\rightarrow\infty$ as
$\gamma\rightarrow\gamma_{+}(\gamma_{-})$ from below (above). However, long
before $g$ diverges, our perturbation theory breaks down. In fact, it loses
validity as soon as $g^{*}\sim{\cal O}(1)$. While we cannot obviously follow
the RG flows all the way to infinity, we can speculate about the nature of the
phases in the region of the parameter space where the RG flow lines appear to
run away to infinity. For this, we are guided by the general understanding
that for large enough noise any system, equilibrium or nonequilibrium, should
undergo a phase transition from a low-noise (low-$T$ in equilibrium systems)
“ordered phase” to a high-noise (high-$T$ in equilibrium systems) “disordered
phase”. Accordingly, the physical expectation in the present model is that for
large enough $g$ (which here means large enough noise or $D$, generalization
of $T$ in equilibrium), the system at $d<4$ undergoes a phase transition from
a low-noise orientationally ordered rough phase to a phase with short-range
order (SRO) only (implying no orientational LRO), i.e., the crumpled phase.
Thus, upon increasing $g$, the rough phase ($d<4$) should be unstable, and the
surface should eventually crumple. This should hold even in the original LDS
equation ($\gamma=0$), for sufficiently $g$: we expect an unstable “crumpling
fixed point” beyond which the RG flow runs away to infinity signalling loss of
orientational LRO. Accordingly, there must be an unstable critical point on
the $g$-axis, controlling this transition to this putative crumpled phase. If
we now consider the full RG flows for a surface in the two dimensional
parameter space ($g,\,\gamma$) and connect this putative flow with our flows
for small $g$ and $\gamma$ in the simplest possible way (i.e., without
introducing any other new fixed points), we are then led to Fig. 1(c). This is
basically an “Occam’s razor”-style argument: Fig. 1(c) has the simplest flow
topology that naturally reduces to the known flow trajectories for small
$g,\,\gamma$ (as shown in Fig. 1(b)). At the same time, it gives the putative
global flow lines, allowing for a transition to a presumed crumpled phase.

Figure 1: (colour online)(a) Schematic phase diagram in the
$\lambda-\lambda_{1}$ plane ($d<4$), showing the regions corresponding to a
rough (stripes) and a crumpled (checkerboard) phase. (b) Schematic phase
diagram and the RG flow lines in the $g-\gamma$ plane ($d<4$). The small
filled circle on the $g$-axis is the original LDS fixed point $g=8\epsilon/3$.
The red curved line is the fixed line $g=\epsilon/(3\Delta(\gamma))$. (c)
Conjectured “Occam’s razor” global RG flows in the $g-\gamma$-plane ($d<4$).
Arrows indicate the flow directions. The green and blue broken lines give the
boundary between the stable (rough) and the crumpled phases. The small filled
red square on the $g$-axis is the putative unstable fixed point not accessible
in our perturbative RG (see text).

We now consider the cases with $d>4$ (i.e., $\epsilon<0$), which although are
not physically accessible, provide interesting theoretical insights.

(iii) Consider now $\epsilon<0$ ($d>4$) and $\Delta<0$. This gives

$\frac{dg}{dl}=g\left[-|\epsilon|+3g|\Delta|\right],$ (18)

giving $g^{*}=0$ as the stable fixed point for a smooth phase with positional
LRO, and $g^{*}=|\epsilon|/(3|\Delta(\gamma)|)$ as the unstable fixed line,
parametrised by $\gamma$. This is reminiscent of a roughening transition of
the KPZ equation stanley . We also note that the fluctuation-corrections to
$\nu$ remain finite even for $L\rightarrow\infty$ for $d>4$, which rules out
vanishing of $\nu_{\text{eff}}$, in turn precluding crumpling as above.
Instead, we obtain a roughening transition between a smooth phase ($g^{*}=0$)
and a perturbatively inaccessible rough phase, akin to the roughening
transition in the KPZ equation at $d>2$ stanley .

(iv) Finally, for $\epsilon<0$ ($d>4$) and $\Delta>0$

$\frac{dg}{dl}=g\left[-|\epsilon|-3g\Delta\right],$ (19)

giving $g^{*}=0$ is the only fixed point, that is stable.

The RG flow for cases (iii) and (iv) are shown in Fig. 2(a). We could again
apply an Occam’s razor style argument as above. Unlike $d<4$, we now expect
loss of positional LRO of the smooth phase upon increasing $g$, resulting into
transitions to a rough phase with no positional LRO (but with orientational
LRO): we assume a putative unstable roughening fixed point in the $g$-axis for
high enough $g$, and arrive at the putative global RG flow diagram shown in
Fig. 2(b). For $d>4$, Fig. 1(a) gives the phase boundaries between regions
having only smooth phases with no phase transitions and regions with
roughening transitions.

Figure 2: (colour online) (a) Schematic phase diagram and the RG flow lines in
the $g-\gamma$ plane ($d>4$). The red curved lines are the unstable fixed
lines given by $g=\epsilon/3|\Delta|,\,\Delta<0$. The shaded region has only a
smooth phase; elsewhere there is a roughening transition, (b) Conjectured
“Occam’s razor” global RG flows in the $g-\gamma$-plane ($d>4$). Arrows
indicate the flow directions. The red curved lines are the unstable fixed
lines that reduces to $g=\epsilon/3|\Delta|,\,\Delta<0$, in the lowest order
perturbation theory. (c) Phases in the $\gamma-\epsilon$ plane. Loss of order
across the transitions is shown (see text).

For $d>4$, (15) give the boundaries between smooth phase ($\chi<0$) and
roughening transition to a perturbatively inaccessible rough phase akin to the
rough phase of the KPZ equation for $d>2$.

The four distinct cases depending upon $\gamma$ and $\epsilon$ are
schematically shown in Fig. 2(c).

To summarise, by suitably generalising the LDS equation, we have uncovered an
intriguing crumpling instability in $d<4$ in addition to the well-known
perturbatively-accessible rough phase which is controlled by the LDS fixed
point. The crumpling is controlled by the relative strength of a spatially
non-local chemical potential vis-á-vis the corresponding local part. We
further show that for $d>4$, instead of a crumpling transition, the system
undergoes a roughening transition for a sufficiently large non-local chemical
potential, which is akin to the roughening transition of the KPZ equation at
$d>2$. Thus, in a direct analogy with the significance of the KPZ equation as
a paradigmatic nonequilibrium model with a roughening transition, our work
establishes the generalised LDS equation as a simple nonequilibrium model that
shows both crumpling and roughening transitions. We note that in Eq. (1) a
$\nabla^{2}h$-term can in principle be added without violating the
conservation law. Such a term remains unaffected by fluctuations. Our theory
applies when the coefficient of this term vanishes. Our work can be further
extended and complemented by applying various nonperturbative and numerical
methods to study the crumpled and rough phases in this model. We conclude with
a cautionary remark that in order to distinguish between a rough (having
orientational LRO) and a crumpled surface (with SRO), one must measure the
variance $\langle({\bm{\nabla}}h)^{2}\rangle$ of the orientation fluctuations,
which for a rough surface is finite in thermodynamic limit, but diverges for a
crumpled surface.

Acknowledgement:- The authors thank J. Toner for many useful discussions. S.M.
thanks the SERB, DST (India) for partial financial support through the TARE
scheme [file no.: TAR/2021/000170] (2022). A.B. thanks the SERB, DST (India)
for partial financial support through the MATRICS scheme [file no.:
MTR/2020/000406] (2021).

## References

  * (1) P. M. Chaikin and T. C. Lubensky, Principles of condensed matter physics, Vol. 1 (Cambridge University Press, Cambridge, 2000).
  * (2) M. Kardar, G. Parisi and Y-C. Zhang, Dynamic Scaling of Growing Interfaces, Phys. Rev, Lett., 56, 889 (1986).
  * (3) A-L Barabási and H. E. Stanley, Fractal concepts in surface growth (Cambridge university press, 1995).
  * (4) T. Halpin-Healy and K-C. Zhang, Kinetic roughening phenom- ena, stochastic growth, directed polymers and all that, Phys. Rep. 254, 215 (1995); A. Basu and E. Frey, “Novel universality classes of cou- pled driven diffusive systems,” Phys. Rev. E 69, 015101 (2004); A. Basu and E. Frey, Scaling and universality in coupled driven diffusive models, J. Stat. Mech.: Theory Exp. 2009, P08013 (2009).
  * (5) T. Sun, H. Guo and M. Grant, Dynamics of driven interfaces with a conservation law, Phys. Rev. A 40, R6763 (1989).
  * (6) F. Caballero et al, Strong Coupling in Conserved Surface Roughening: A New Universality Class?, Phys. Rev. Lett. 121, 020601 (2018).
  * (7) Z. -W. Lai and S. Das Sarma, Kinetic Growth with Surface Relaxation: Continuum versus Atomistic Models, Phys. Rev. Lett. 66, 2348 (1991).
  * (8) A. Brú, J. M. Pastor, I. Fernaud, I. Brú, S. Melle, and C. Berenguer, Phys. Rev. Lett. 81, 4008 (1998); A. Brú, S. Albertos, J. L. Subiza, J. L. García-Asenjo, and I. Brú, Biophys. J. 85, 2948 (2003).
  * (9) M. Paczuski, M. Kardar, and D. R. Nelson, Landau Theory of the Crumpling Transition, Phys. Rev. Lett. 60, 2638 (1988).
  * (10) T. Banerjee, N. Sarkar, J. Toner, and A. Basu, Rolled up or Crumpled: Phases of Asymmetric Tethered Membranes, Phys. Rev. Lett. 122, 218002 (2019).
  * (11) T. Banerjee, N. Sarkar, J. Toner and A. Basu, Statistical mechanics of asymmetric tethered membranes: Spiral and crumpled phases, Phys. Rev. E 99, 053004 (2019).
  * (12) F. David and E. Guitter, Crumpling Transition in Elastic Membranes: Renormalization Group Treatment, Europhys. Lett. 5, 709 (1988).
  * (13) D. Forster, D. R Nelson, and M. J Stephen, Large-distance and long-time properties of a randomly stirred fluid, Phys. Rev. A 16, 732 (1977).
  * (14) P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. 49, 435 (1977).
  * (15) R Bausch, H. Janssen, and H. Wagner, Renormalized field theory of critical dynamics, Z. Phys. B 24, 113 (1976); U. Täuber, Critical Dynamics (Cambridge University Press, Cambridge, 2014).
  * (16) See Supplemental Material for details of the RG calculations including the Feynman diagrams.
  * (17) At the two-loop order, $\lambda$ does renormalise; see, H. Janssen, On Critical Exponents and the Renormalization of the Coupling Constant in Growth Models with Surface Diffusion, Phys. Rev. Lett. 78, 1082 (1997). In the absence of a symmetry argument, presumably not only $\lambda_{1}$ too renormalises, $\lambda$ and $\lambda_{1}$ may mix under RG. This provides a technical justification to include the $\lambda_{1}$-term in the theory.
  * (18) L. Peliti and S. Leibler, Effects of Thermal Fluctuations on Systems with Small Surface Tension, Phys. Rev. Lett. 54, 1690 (1985).
  * (19) T. Banerjee and A. Basu, Thermal fluctuations and stiffening of symmetric heterogeneous fluid membranes, Phys. Rev. E 91, 012119 (2015).

##  Appendix A Action functional

Model equation (1) of the main text can be converted to a path integral over
configurations given by $h({\bf x},t)$ and its dynamic conjugate $\hat{h}({\bf
x},t)$, which can be used to conveniently calculate the various correlation
and vertex functions in the problem. The generating functional is given by

${\cal Z}=\int{\cal D}h{\cal D}\hat{h}\exp(-S),$ (20)

where $S$ is the action functional given by

$S=\int d^{d}rdt\bigg{[}-D\hat{h}({\bf r},t)\hat{h}({\bf r},t)\bigg{]}+\int
d^{d}xdt\,\hat{h}\bigg{[}\partial_{t}h+\nabla^{2}\\{\nu\nabla^{2}h+\frac{\lambda}{2}({\bm{\nabla}}h)^{2}\\}+\lambda_{1}{\bm{\nabla}}\cdot\\{(\nabla^{2}h){\bm{\nabla}}h\\}\bigg{]}.$
(21)

We can now write down the “bare” two-point correlation functions in the
problem by using (20) and the harmonic or Gaussian part of (21); i.e., by
setting $\lambda=0=\lambda_{1}$. We obtain

$\displaystyle\langle|\hat{h}({\bf q},\omega)|^{2}\rangle=0,$ (22)
$\displaystyle\langle\hat{h}({\bf-q},-\omega)h({\bf
q},\omega)\rangle=\frac{1}{-i\omega+\nu q^{2}}.$ (23)
$\displaystyle\langle|h({\bf
q},\omega)|^{2}\rangle=\frac{2D}{\omega^{2}+\nu^{2}q^{4}}.$ (24)

##  Appendix B Renormalization group analysis

We outline here the basic steps of the RG calculations performed in the main
text. The RG is done by tracing over the short wavelength Fourier modes of the
fields $h({\bf x},t)$ and $\hat{h}({\bf x},t)$ (see Refs. [14-16] of the main
text). In particular, we follow the usual approach of initially restricting
the wavevectors to be within a bounded spherical Brillouin zone: $|{\bf
k}|<\Lambda$. However, the precise value of the upper cutoff $\Lambda$ has no
effect on our final results. The fields $h({\bf x},t)$ and $\hat{h}({\bf
x},t)$ are separated into the high and low wave vector parts: $h({\bf
x},t)=h^{<}({\bf x},t)+h^{>}({\bf x},t)$, $\hat{h}({\bf x},t)=\hat{h}^{<}({\bf
x},t)+\hat{h}^{>}({\bf x},t)$, where $h^{>}({\bf x},t)$ and $\hat{h}^{>}({\bf
x},t)$ have support in the large wave vector (short wavelength) range $\Lambda
e^{-l}<|{\bf k}|<\Lambda$, while $h^{<}({\bf x},t)$ and $\hat{h}^{<}({\bf
x},t)$ have support in the small wave vector (long wavelength) range $|{\bf
k}|<e^{-l}\Lambda$; $b\equiv e^{l}>1$. We then integrate out $h^{>}({\bf
x},t)$ and $\hat{h}^{>}({\bf x},t)$ in the anhamornic coupling $\lambda$ and
$\lambda_{1}$; as usual, this resulting perturbation theory of $h^{<}({\bf
x},t)$, $\hat{h}^{<}({\bf x},t)$ can be represented by Feynman graphs, with
the order of perturbation theory reflected by the number of loops in the
graphs we consider. We restrict ourselves here up to the one-loop order. After
this perturbative step, we rescale lengths, with ${\bf x}\rightarrow{\bf
x}^{\prime}e^{l}$, which restores the UV cutoff back to $\Lambda$, together
with rescaling of time $t\rightarrow t^{\prime}e^{zl}$ (equivalently in the
Fourier space, the momentum and frequency are rescaled as ${\bf
k}\rightarrow{\bf k}^{\prime}/b$ and $\omega\rightarrow\omega^{\prime}/b^{z}$,
respectively), where $z$ is the dynamic exponent. This is then followed by
rescaling the long wave length part of the fields that we define in the
Fourier space for calculational convenience. We scale $h^{<}({\bf
q},\omega)=\xi h(b{\bf q},b^{z}\omega),\,\hat{h}^{<}({\bf
q},\omega)=\hat{\xi}\hat{h}(b{\bf q},b^{z}\omega^{\prime})$.

Next, we set $\xi\hat{\xi}=b^{d+2z}$ by demanding the term in action (21)
$\int d^{d}{\bf k}d\omega\,\hat{h}h\,\omega$ does not scale with the scaling
of momenta and frequencies. Accordingly the model parameters scale as given
below:

$\displaystyle\nu^{\prime}=\nu^{<}b^{z-4},\,D^{\prime}_{h}=D_{h}^{<}b^{z-d-2\chi_{h}},$
(25)
$\displaystyle(\lambda^{\prime},\lambda_{1}^{\prime})=(\lambda^{<},\lambda_{1}^{<})b^{z-2+\chi}.$
(26)

###  B.1 One-loop Feynman diagrams for $\nu$

Figure 3: One-loop Feynman diagrams which contribute to the fluctuation-
correction of $\nu$.

###  B.2 One-loop vertex corrections

Figure 4: Effective vertex $\tilde{\lambda}({\bf
k,\,q})=-\frac{1}{2}[q^{2}{\bf k\cdot(k-q)}+({\bf k-q})^{2}{\bf k\cdot q}]$.

Figure 5: One-loop Feynman diagrams which contribute to the fluctuation-
corrections of $\lambda$.

Sum of these diagrams in Fig. (5) vanishes.

Figure 6: One-loop Feynman diagrams which contribute to the fluctuation-
corrections of $\tilde{\lambda}$.

Sum of these diagrams in Fig. (6) vanishes.

##  Appendix C Perturbative corrections to the model parameters

Combining all the diagrams in Fig. (3), and evaluating them, we obtain the
fluctuation-correction diffusivity $\nu^{<}$:

$\displaystyle\nu^{<}$ $\displaystyle=$
$\displaystyle\nu\bigg{[}1+\frac{{\lambda}^{2}D}{{\nu}^{3}}\frac{6-d}{4d}K_{d}\int_{\Lambda/b}^{\Lambda}dq\,q^{d-5}$
(27)
$\displaystyle-\frac{\lambda_{1}^{2}D}{{\nu}^{3}}K_{d}\int_{\Lambda/b}^{\Lambda}dq\,q^{d-5}$
(28)
$\displaystyle+\frac{\lambda_{1}^{2}D}{{\nu}^{3}d}\left(\frac{6-d}{4}-\frac{6}{d+2}+3\right)K_{d}\int_{\Lambda/b}^{\Lambda}dq\,q^{d-5}$
(29)
$\displaystyle+\frac{{\lambda}{\lambda_{1}}D}{{\nu}^{3}}K_{d}\int_{\Lambda/b}^{\Lambda}dq\,q^{d-5}\bigg{]},$
(30) $\displaystyle\lambda^{<}$ $\displaystyle=$ $\displaystyle\lambda,$ (31)
$\displaystyle\lambda_{1}^{<}$ $\displaystyle=$ $\displaystyle\lambda_{1},$
(32) $\displaystyle D^{<}$ $\displaystyle=$ $\displaystyle D.$ (33)

We now rescale space, time and fields so as to raise the upper wavevector
cutoff back to $\Lambda$. This gives the following differential RG flow
equations for $\nu$ and $g$:

$\displaystyle\frac{d\nu}{dl}$ $\displaystyle=$
$\displaystyle\nu\bigg{[}z-4+\bigg{\\{}\frac{{\lambda}^{2}D}{{\nu}^{3}}\frac{6-d}{4d}-\frac{\lambda_{1}^{2}D}{{\nu}^{3}}$
(34)
$\displaystyle+\frac{\lambda_{1}^{2}D}{{\nu}^{3}d}\left(\frac{6-d}{4}-\frac{6}{d+2}+3\right)$
(35)
$\displaystyle+\frac{{\lambda}{\lambda_{1}}D}{{\nu}^{3}}\bigg{\\}}K_{d}{\Lambda}^{d-4}\bigg{]},$
(36) $\displaystyle\frac{dg}{dl}$ $\displaystyle=$ $\displaystyle
g\bigg{[}4-d-3\bigg{\\{}\frac{{\lambda}^{2}D}{{\nu}^{3}}\frac{6-d}{4d}-\frac{\lambda_{1}^{2}D}{{\nu}^{3}}$
(37)
$\displaystyle+\frac{\lambda_{1}^{2}D}{{\nu}^{3}d}\left(\frac{6-d}{4}-\frac{6}{d+2}+3\right)$
(38)
$\displaystyle+\frac{{\lambda}{\lambda_{1}}D}{{\nu}^{3}}\bigg{\\}}K_{d}{\Lambda}^{d-4}\bigg{]}.$
(39)

This clearly shows that $d=4$ is the critical dimension of $g$. Near $d=4$,
(39) reduces to (10) of the main text.

The Feynman diagrams in Fig. 3 can be used to calculate $\nu_{\text{eff}}$,
whence the wavevector loop integrals are done from $\Lambda$ all the way to
$2\pi/L$. Clearly for $d<4$, the wavevector loop integrals depend sensitively
on $L$, the system size. In fact, they diverge with $L$, a fact that opens up
the intriguing possibility of crumpling with increasing $L$. In contrast,
these wavevector integrals remain finite for $d>4$, even if
$L\rightarrow\infty$ limit is taken.


# DEEP REINFORCEMENT LEARNING FOR SAFE LANDING SITE SELECTION WITH CONCURRENT
CONSIDERATION OF DIVERT MANEUVERS 111This paper is an updated version of Paper
AAS 20-583 presented at the AAS/AIAA Astrodynamics Specialist Conference,
Online, in 2020.

Keidai Iiyama Master Student, Department of Aeronautics and Astronautics, The
University of Tokyo, 7-3-1, Hongo,Bunkyo-ku,Tokyo, Japan. Kento Tomita Ph.D.
Student, School of Aerospace Engineering, Georgia Institute of Technology,
Atlanta, GA, 30332. Bhavi A. Jagatia Master Student, School of Aerospace
Engineering, Georgia Institute of Technology, Atlanta, GA, 30332. Tatsuwaki
Nakagawa Bachelor Student, School of Aerospace Engineering, Georgia Institute
of Technology, Atlanta, GA, 30332. and Koki Ho Assistant Professor, School of
Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA, 30332.

###### Abstract

This research proposes a new integrated framework for identifying safe landing
locations and planning in-flight divert maneuvers. The state-of-the-art
algorithms for landing zone selection utilize local terrain features such as
slopes and roughness to judge the safety and priority of the landing point.
However, when there are additional chances of observation and diverting in the
future, these algorithms are not able to evaluate the safety of the decision
itself to target the selected landing point considering the overall descent
trajectory. In response to this challenge, we propose a reinforcement learning
framework that optimizes a landing site selection strategy concurrently with a
guidance and control strategy to the target landing site. The trained agent
could evaluate and select landing sites with explicit consideration of the
terrain features, quality of future observations, and control to achieve a
safe and efficient landing trajectory at a system-level. The proposed
framework was able to achieve 94.8 $\%$ of successful landing in highly
challenging landing sites where over 80$\%$ of the area around the initial
target lading point is hazardous, by effectively updating the target landing
site and feedback control gain during descent.

##  1 Introduction

Figure 1: Landing site selection and control

Figure 2: Observability and trajectory

On-board hazard detection and avoidance (HDA) capabilities are essential to
enable new mission concepts that involve planetary surface operations. With a
quick assessment of the perceived terrain data (e.g., DEM, visible spectrum
map, or a combination thereof) from optical and/or LIDAR sensors, the HDA
technology creates a map of probability of safety for prioritizing candidate
landing zones. NASA has been actively developing the technology for Precise
landing and Hazard Avoidance (PL& HA) [1, 2, 3, 4, 5, 6], and Safe and Precise
Landing Integrated Capabilities Evolution (SPLICE) project is currently
underway to develop next generation PL&HA technologies [7, 8, 9]. The HDA
algorithm being developed for SPLICE directly leverages the HDA algorithms
from the previous Autonomous Landing Hazard Avoidance Technology (ALHAT)
project. ALHAT is capable of quickly assessing the DEM on-board and in real
time during the descent, and assigns a probability of a safe landing to each
pixel on the map. [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. ALHAT’s HDA
capability along with the terrain relative navigation and hazard detection
functions was successfully demonstrated during the hardware-in-the-loop
testing on the Morpheus Vertical Testbed [16, 21]. However, a critical caveat
of ALHAT is that it only selects a list of discrete landing points based on
local static information (e.g., slopes, roughness, etc.). To consider the
feasibility of the divert maneuver, a landing selection algorithm that
calculates approximate landing footprint has been presented [22, 23]. To
reflect predetermined scientific values of each landing site in the landing
site selection process, landing site selection methods that leverages Bayesian
networks has been proposed [6]. Cui, et al. (2017) proposed a method that
calculates synthetic landing area assessment criteria based on terrain safety,
fuel consumption, and touchdown performance during descent [24].

In these previously proposed landing site assessment methods, future changes
in the target landing site are not considered. However, since the field of
view and the quality of the observation data changes depending on the
spacecraft state, the best landing site changes each time new observation data
is obtained. This effect could not be ignored, especially when observability
is limited or safe landing sites are sparse. When multiple chances of
observation and divert maneuver are considered, the following two aspects
should be considered.

Figure 3: General Concept

The first aspect is the coupling between the target landing site selection and
control maneuver planning. Existing landing site selection algorithms assume
that the lander is guided by a simple predetermined control law, while the
control law is designed to guide the lander to the decided landing target.
However, when additional divert maneuvers are available in the future, the
landing site selection policy is dependent on the control policy, and vice-
versa. Therefore, the landing site and control plan has to be decided
simultaneously in order to achieve an overall optimal trajectory. Fig.2 shows
an example of this coupling. Suppose landing site A and landing site B both
have uncertainty in the safety level that could be judged from the
observation. If there is little residual fuel, it could be safer to take the
fuel optimal path to the closer landing site B (Trajectory B). However, when
there is sufficient residual fuel, it could be safer to take a longer path to
landing site B (Trajectory A), if additional diverting to landing site A could
happen in case landing site B turns out to be unsafe as the lander approaches
the ground.

The second aspect is the observability of the terrain during descent. As
illustrated in Figure 2, the change in the control policy and the target
landing site affects the quality and number of the obtained observation data
during descent. Therefore, the controller has to be aware that sufficient
information to judge a good landing site could be obtained by following the
planned trajectory. In general landing trajectory design, glide slope
constraints are applied as a path constraint to ensure observations from
higher elevations. However, the actual quality and quantity of observation
information could not be explained only by slope angles. To tackle this
problem, Crane (2013) developed an on-line information-seeking trajectory
modification method that selects a trajectory that minimizes the weighted sum
of the estimated entropy and field of view of the image obtained in the
future, fuel consumption, and the coverage of the entire field [25]. However,
the research focused on modifying the nominal trajectory for a fixed landing
site, and updating the landing site successively during descent was not
considered in the research. The impact of landing site selection on the
observation data quality should also be considered, since changing the landing
site affects the observation through changes in the trajectory and the
observation target,

In order to tackle these problems, this paper proposes a learning-based method
that selects landing site and design guidance trajectory successively during
the HDA phase. The agent in the proposed method seeks to maximize the total
probability of landing to a safe terrain landing site while minimizing total
fuel consumption, landing error to the target, and final velocity. Model-free
reinforcement learning techniques are leveraged to learn a policy that
concurrently selects landing site and guidance strategy to the target, by
interacting with the simulator environment and learning how to maximize the
total probability of successful landing with minimal fuel consumption. The
general concept is shown in Fig.3.

The outline is as follows. First, we explain the assumptions used in modeling
the HDA phase and show that the sequence of action choices can be formulated
as a partially observable Markov decision process (POMDP). Next, the method to
optimize landing site selection policies and the controller is introduced.
Finally, the performance of the obtained controller is analyzed, and
qualitative interception is given.

##  2 Problem Formulation

###  2.1 Dynamics

A soft lunar landing scenario is assumed in this paper. In this paper, the 3
degrees of freedom (3-DOF) problem is considered. The equation of motion
governing the dynamics of the problem are given as follows.

$\begin{split}\boldsymbol{v}&=\boldsymbol{\dot{r}}\\\
\boldsymbol{a}&=\frac{\boldsymbol{T}}{m}+\boldsymbol{g}\\\
\dot{m}&=-\frac{||\boldsymbol{T}||}{g_{ref}I_{sp}}\end{split}$ (1)

where $r=[r_{x}\ r_{y}\ r_{z}]$ is the position in the target centered
orthonormal frame with the z-axis pointing upward, and $\boldsymbol{g}=[0\ 0\
-1.62]^{T}$ is the gravity. In this paper, gravity is considered to be
constant during the entire mission, and the effect of planetary rotation is
ignored. In addition, limitation in the thrust range is applied as follows.

$0\leq|\boldsymbol{T}|\leq T_{max}$ (2)

The following glide slope constraint are applied in most planetary landing
problems

$\theta_{g}=\arctan{\biggl{(}\frac{\sqrt{r_{x}^{2}+r_{y}^{2}}}{r_{z}}\biggr{)}}<\theta_{gmax}$
(3)

This constraint is applied in order to avoid the lander from hitting the
terrain at low altitudes, and to ensure that the inclination to the target is
kept over a certain amount for navigation and hazard detection purposes.

###  2.2 Terrain Generation

Figure 4: Terrain generation process

Figure 5: Example of the generated terrain

Observation data consists of DEM generated by HD LIDAR, and spacecraft state
information (position, velocity, mass) with noise. In order to generate LIDAR
DEMs, we need a true DEM as a reference. Lunar Reconnaissance Orbiter (LROC)
database provides Digital Terrain Model (DTM) of the lunar surface [26].
However, the resolution of the DTM ($\approx$ 50m/pixel) in the database is
not sufficient for the LIDAR DEM model generation. Therefore, we made an
original terrain generator leveraging crater distribution models proposed by
Pike, and boulder distribution models proposed by Bernard [27][28]. While the
model was originally developed for Mars, studies of rock density on the moon
show that the model can be substituted for the lunar topography by adjusting
the parameters. The generated terrain has a size of 1000m x 1000m and a
resolution of 1m. The generation process of the terrain is described in Fig.5.
An example of the generated terrain is shown in Fig.5.

###  2.3 Safety Map Generation

In order to assess if the landing has succeeded in the simulator, the safety
of each landing point in the generated terrain map has to be assessed. We
evaluated the safety by applying the safety assessment algorithm developed in
the ALHAT project to the entire generated terrain. The algorithm for assessing
the deterministic safety value $V_{D}\in\\{0,1\\}$ is shown in Algorithm 1.

Algorithm 1 Deterministic and Stochastic safey map generation

DEM $D(m\times n)$

deterministic safety value map $V_{D}(m\times n)\in\\{0,1\\}$, stochastic
safety value map $V_{P}(m\times n)\in[0,1]$

Initialize:

$S\leftarrow$ [0, …, 0] $\triangleright$ list of slope for each orientation
$R\leftarrow$ [0, …, 0] $\triangleright$ list of roughness for each
orientation
$P\leftarrow$ [0, …, 0] $\triangleright$ list of safety probability for each
orientation
$O\leftarrow$ [0, $\frac{1}{n_{o}}\pi,\ldots,\frac{n_{o}-1}{n_{o}}\pi$]
$\triangleright$ possible lander footpad orientation

for $row=1,\cdots,m$ do $\triangleright$ for each pixel in the DEM

for $col=1,\cdots,n$ do

Calculate pad placement and contact

for $o_{i}=1,\cdots,n_{o}$ do $\triangleright$ calculate worst case slope for
all orientations

$S[o_{i}]=\textsc{Get\\_slope}(row,col,D)$

end for

$s_{max}\leftarrow\textsc{Max}(S)$

if $s_{max}>safe\ slope\ threshold$ then $\triangleright$ slope not safe

$V[row][col]\leftarrow 0$

else

for $o_{i}=1,\cdots,n_{o}$ do

$R[o_{i}]=\textsc{Get\\_roughness}(row,col,D)$

$P[o_{i}]=\textsc{Roughness\\_To\\_Safety\\_Probability}(row,col,D,R)$

end for

$V_{P}[row][col]\leftarrow\textsc{Min}(P)$

$r_{max}\leftarrow\textsc{Max}(R)$

if $r_{max}>safe\ roughness\ threshold$ then $\triangleright$ roughness not
safe

$V_{D}[row][col]\leftarrow 0$

else$\triangleright$ slope and roughness safe

$V_{D}[row][col]\leftarrow 1$

end if

end if

end for

end for

return $V_{D}.V_{P}$

###  2.4 Observation Data Generation

The control agent utilizes two observation data for control: estimated lander
state (position, velocity, and mass) and 2d map of the terrain. For the lander
state, we assumed that we have perfect information without errors.
Incorporating navigation errors and hazard relative navigation algorithms in
the framework is important future work. For the 2d map of the terrain, we
considered two options. The first option is to directly pass the LIDAR DEM,
while the second option is to pass the stochastic safety map obtained by
applying Algorithm1 on the obtained LIDAR DEM. The first option requires the
control agent to assess the safety of each landing site from scratch, while in
the second option, the control agent has access to more direct information
about the safety of each landing point based on roughness and slope values. In
this paper, we chose the second approach as we failed to design a controller
using the first approach.

The simulation of the LIDAR DEM during descent requires complex calculations.
To simulate a LIDAR DEM, detector pattern formulation, ray interception with
the terrain, addition of the range bias, transformation to the point clouds,
and transformation to a digital elevation map has to be conducted even if we
assume that the LIDAR position is fixed to a known position. In lunar descent
cases, errors from vehicle state knowledge, LIDAR misalignment, map assembly
errors when considering vehicle motions, and system latency should also have
to be considered. Since the accurate modeling of the entire procedure is
extremely complicated, as an initial concept study, we chose to generate
pseudo observation safety map at each observation timing by cutting out the
portion of the stochastic safety map of the entire field, and adding noise to
it. The stochastic safety map could be obtained by directly applying
Algorithm1 to the entire terrain DEM. By this way, we are able to greatly save
the computation time required for simulation since we can avoid generating
DEMs and running Algorithm1 at each observation timing.

When generating the ”observation safety map”, three main relationship between
the lander’s state (or trajectory) and the LIDAR DEM (and the generated safety
map) were modeled. The first effect is the field of view (FOV) of the LIDAR
DEM. For simplicity, we assumed the obtained DEM is square, and its two axes
is always aligned with the $x,y$ coordinate of the map. The length of one side
of the square $w_{x}=w_{y}$ is calculated using the following equation.

$w_{x}=r_{s}\tan{\phi_{l}}$ (4)

where $r_{s}$ is the slant range (distance between the Lander and the target
position), and $\phi_{l}$ is the field of view of the LIDAR which was set to
$11.4$ [deg].(Fig. 6) In reality, the DEM field of view is not rectangle, and
the field of view differs between the horizontal and vertical direction of the
lander, but we believe this assumption is fair enough for an initial concept
study.

Figure 6: Illustration of the field of view of the LIDAR DEM

The second effect is the effect of the slant range. As the slant range
increases, errors in the range measurement in the LIDAR DEM increases. This
effect was simulated by adding each pixel a Gaussian noise with an error
proportional to the distance to each pixel. In addition, as the slant range
gets longer, the sampling distance increases, small boulders in the terrain
will be overlooked. This effect was modeled by fixing the LIDAR DEM size to
64x64. The third effect is the effect of slant angle. Slant angle is the angle
between the horizontal plane and beam direction. As the slant angle decreases,
holes in the DEM appear behind surface hazards. Therefore, regions where
safety value could not be calculated appear in the safety map which leads to
fewer safe regions in the map [29]. In this paper, this effect was simply
modeled by assigning unsafe labels to pixels with slant angles between the
lander smaller than 70 degrees.

In our settings, it is assumed that LIDAR DEMs could be obtained every 5
seconds, and the DEM size was fixed to 64x64. The interval time is based on
the maximum processing time of the ALHAT algorithm to process DEM and generate
a safety map. The entire observation data generation process is described in
Algorithm 2.

Algorithm 2 Observation Data Generation at Each Timestep

true probabilistic safety value map of the entire terrain $V_{p}(1000\times
1000)$

true position of the lander in a surface fixed frame
$X_{L}=[x_{l},y_{1},z_{l}]$

noisy probabilistic safety value map $O_{P}(64\times 64)\in[0,1]$

// calculate the FOV of the DEM and return 2D array of center position of each
pixel in the

// observed DEM

$X_{D},Y_{D},w_{x},w_{y}=\textsc{Calc\\_FOV}(X_{L},\phi_{l})$ $\triangleright$
$\phi_{l}$: LIDAR FOV, Eq.4

for $row=1,\cdots,m$ do $\triangleright$ for each pixel in the DEM

for $col=1,\cdots,n$ do

$x_{d},y_{d}=X_{d}[row][col],Y_{d}[row][col]$

$\theta=\textsc{Slant\\_Angle}(x_{d},y_{d},X_{L})$

if $\theta>20^{\circ}$ then $\triangleright$ slant angle to the pixel is over
thershold

$O_{P}[row][col]=0$

else

// calculate the safety probability of the pixel by finding nearest neighbor
pixel in the

// whole $V_{P}$ map

$v_{p}=\textsc{Nearest\\_Neighbor}(x_{d},y_{d},V_{P})$

$r=\textsc{Slant\\_Range}(x_{d},y_{d},X_{L})$ $\triangleright$ [m]

$O_{P}[row][col]=\textsc{Clip}(v_{p}+\mathcal{N}(0,\frac{r}{500}*0.05),0,1)$
$\triangleright$ Add noise to the safety value

end if

end for

end for

return $O_{P},w_{x},w_{y}$

###  2.5 Modeling the Problem as POMDP

The HDA sequence is modeled as a POMDP
$P=<\mathcal{S},\mathcal{A},T,R,\Omega,O>$ as follows.

  * •

State space $s\in\mathcal{S}$: Spacecraft state (position, velocity, mass),
true DEM of the entire field

  * •

Action space $a\in\mathcal{A}$: Lander thrust output, target landing position
(determines next target position for LIDAR measurement)

  * •

Reward space $r\in R$: Consumed fuel, the safety of the final landing point,
final velocity

  * •

State transition function $T:S\times A\rightarrow\Pi(S)$: Stochastic dynamics
of the spacecraft. $T(s,a,s^{\prime})$ is the probability of moving to state
$s^{\prime}\in\mathcal{S}$ when the agent at state $s\in\mathcal{S}$ takes
action $a\in\mathcal{A}$.

  * •

Observation space $\Omega$: LIDAR DEM, predicted spacecraft state

  * •

Observation function $O:S\times A\rightarrow\Pi(\Omega)$: LIDAR DEM generation
models, spacecraft state observation model

Obtaining an optimal policy (a policy that could maximize expected reward) in
POMDP is difficult due to the partial observability of the problem. In
general, an agent requires the entire history of the observation and action
pairs $h_{t}=\\{a_{0},r_{0},o_{1},\ldots,a_{t-1},r_{t-1},o_{t}\\}$. POMDP
could be converted into a belief Markov decision process (belief MDP) by
introducing belief state $b$. Belief state is a conditional probability
function for $s\in S$ given the history $h_{t}$ . When state transition
function $T$ and observation function $O$ is known, an optimal policy of the
POMDP could be obtained by approximating the belief MDP and applying a value
iteration method.

In this paper, we seek to obtain the optimal policy of the above POMDP without
modeling the $T$ and $O$ function, but by learning from transition data
collected by interacting with the simulator. There are two approaches in
model-free approach for POMDPs. The first approach is the memory-less
approach, which learns Markov policy by simply considering the most recent
observation $o$ as a state. In POMDP, the Bellman equation is not strictly
satisfied, so deterministic policies leveraging only current state information
is not guaranteed to be optimal. However, when the partial observability of
the problem is weak, memory-less approach may be sufficient. The second
approach is the memory-based approach, which learns history-dependent policy
that uses the entire history data. This is usually achieved by using Recurrent
Neural Networks (RNNs) that could store history information as a single state.
In this paper, both memory-less and memory-based approaches are tested.

##  3 Guidance and Control

###  3.1 Observation Data Interpretation

Since $64\times 64$ map data is used as part of the observation, policies that
take in full image data taken as an input have too large dimension to
optimize. Therefore, the whole history data is transferred into a low-
dimension internal state value by leveraging an auto-encoder. The role of the
auto-encoder is to extract an abstract, compressed representation of the LIDAR
DEM data as a latent vector $z$, to avoid directly passing down large input
data to the reinforcement learning agent. The auto-encoder was trained
separately with the reinforcement learning agent, as proposed in the ”World
Models” paper by Ha (2018). The dataset for training was collected through
5000 random rollouts, and the auto-encoder was trained to minimize the
difference between the observed safety map and the reconstructed safety map
produced by the decoder. The following Fig.7 shows examples of the training
data and the reconstructed image.

Figure 7: Original observation data (Up) and reconstructed observation
data(Down) by the Autoencoder

###  3.2 General Concept of Control

For control, the most direct approach is to train the agent to output a three-
dimensional thrust vector that guides the lander along a robust and fuel-
efficient trajectory to the selected landing point. However, learning this
feed-forward control policy from scratch requires high computational cost for
learning. Several research have tackled this problem by shaping the reward in
an effective way[30], but we took a different approach in order to relate the
control policy with target landing points. We adopted Zero-Effort-Miss/Zero
Effort-Velocity (ZEM-ZEV) feedback guidance algorithm as a baseline guidance
law [31], and designed the agent that controls the target landing point
$r_{f}$, and adaptive hyper-parameters ($K_{R},K_{V},t_{go}$) in the ZEM-ZEV
guidance algorithm. The idea of learning adaptive hyper-parameters in ZEM-ZEV
feedback guidance algorithm with reinforcement learning has already been
proposed in previous research, to design an ZEM-ZEV feedback controller that
avoids slant angle constraint violation during descent. [32] The difference in
our paper is that the target position also changes during the descent.

###  3.3 Controller

ZEM/ZEV feedback guidance algorithm calculates the optimal acceleration using
the ZEM and ZEV, which represents the distance between the final target
position and velocity and the projected final position and velocity when no
additional control is added from current time $t$. The optimal acceleration is
calculated as

$\boldsymbol{a}=\frac{K_{R}}{t_{go}^{2}}\boldsymbol{ZEM}-\frac{K_{V}}{t_{go}}\boldsymbol{ZEV}$
(5)

where $K_{R}$, $K_{V}$ are control gains, and $t_{go}$ is time-to-go. When
there are no limitations in the thrust magnitude and no constraints in the
trajectory, it is proved that energy-optimal trajectory could be obtained by
setting the control gain as $K_{R}=6,K_{V}=2$. Energy optimal time-to-go
$t_{go}$ could be obtained by solving the following equation.

$t_{go}^{4}\boldsymbol{g^{T}}\boldsymbol{g}-2t_{go}^{2}(\boldsymbol{v^{T}}\boldsymbol{v}+\boldsymbol{v_{f}^{T}}\boldsymbol{v}+\boldsymbol{v_{f}^{T}}\boldsymbol{v_{f}})+12t_{go}{(\boldsymbol{r_{f}}-\boldsymbol{r})}^{T}(\boldsymbol{v}+\boldsymbol{v_{f}})-18{(\boldsymbol{r_{f}}-\boldsymbol{r})}^{T}(\boldsymbol{r_{f}}-\boldsymbol{r})=0$
(6)

When there is constraint $|\boldsymbol{T}|=m|\boldsymbol{a}|\leq T_{max}$ in
thrust magnitude, saturated acceleration is used for control as follows.

$\displaystyle\boldsymbol{a}=\begin{cases}\boldsymbol{\bar{a}}&|\boldsymbol{\bar{a}}|\leq
T_{max}/m\\\
\boldsymbol{\bar{a}}T_{max}/m|\boldsymbol{\bar{a}}|&|\boldsymbol{\bar{a}}|>T_{max}/m\end{cases}$
(7)
$\displaystyle\boldsymbol{\bar{a}}=\frac{K_{R}}{t_{go}^{2}}\boldsymbol{ZEM}-\frac{K_{V}}{t_{go}}\boldsymbol{ZEV}$
(8)

The goal of the research is to train a reinforcement learning agent that
decides $K_{R},K_{V},t_{go}$ and target landing position $r_{f}$ at each time
step.

##  4 Training

###  4.1 Reinforcement Learning

In standard reinforcement learning problems, the agent interacts with the
fully-observed environment $E$. At each time step $t$, the agent outputs an
action $a_{t}$, and the environment returns the next observation
$o_{t+1}=s_{t+1}$ and a step reward $r_{t+1}=g(s_{t},a_{t})$. The agent
decides its action based on its policy $\pi(s_{t})$, which maps the state to a
probability distribution over possible actions. The environment follows the
transition dynamics $p(s_{t+1}|s_{t},a_{t})$ and initial state distribution
$p_{s0}$. The entire process could be modeled as a Markov decision process
$\mathcal{M}=\\{\mathcal{S},\mathcal{A},p_{s0},p,g\\}$. The goal of the agent
is to learn a policy that maximizes the expectation of discounted future
reward $\mathbb{E}[R_{t}|\mathcal{M},\pi]$ where
$R_{t}=\Sigma_{i=t}^{T}\gamma^{T}g(s_{i},a_{i})$. Reinforcement learning
algorithms utilize the following action-value function (or Q-function), which
describes the expected return after taking action $a_{t}$ in state $s_{t}$,
then following policy $\pi$.

$Q^{\pi}(s_{t},a_{t})=\mathbb{E}[R_{t}|s_{t},a_{t}]$ (9)

For the Q function, it is known that the following recursive equation called
the Bellman equation holds.

$Q^{\pi}(s_{t},a_{t})=\mathbb{E}_{r_{t},s_{t+1}\sim
E}[g(s_{t},a_{t})+\gamma\mathbb{E}_{a_{t+1}\sim\pi}[Q^{\pi}(s_{t+1},a_{t+1})]]$
(10)

###  4.2 Deep Deterministic Policy Gradient (DDPG)

The Deep Deterministic Policy Gradient (DDPG) is an actor-critic, off-policy,
model-free reinforcement learning algorithm for continuous action spaces. [33]
In actor-critic methods, the critic estimates the action-value function
$Q(s,a)$, while the actor produces the action given the current state based on
its policy $\pi(a|s)$. As for the actor, we consider a parameterized
deterministic policy $a=\pi_{\phi}(s)$ with parameter $\phi$. $\phi$ could be
trained using the following deterministic policy gradient theorem and the
estimated Q-values by the critic. [34]

$\nabla_{\phi}J(\phi)=\mathbb{E}_{s\sim\rho_{\pi}}[\nabla_{a}Q(s,a|\theta)|_{s=s_{t},a=\mu(s_{t})}\nabla_{\phi}\mu(s|\phi)|_{s=s_{t}}]$
(11)

For the critic, in order to handle continuous state and action spaces,
$Q(s,a)$ is also approximated using a function approximater parameterized by
$\theta$. The critic uses the collected data to learn parameters that better
approximate $Q(s,a)$. This could be achieved by minimizing the following loss
$L(\theta)$

$L(\theta)=\mathbb{E}_{s_{t}\sim\rho^{\beta},a_{t}\sim\beta,r_{t}\sim
E}[Q(s_{t},a_{t}|\theta)-y_{t})^{2}]$ (12)

where $\beta$ is an arbitrary stochastic behavior policy, and $\rho^{\beta}$
is a state visitation distribution using policy $\beta$. $y_{t}$ is the target
value

$y_{t}=g(s_{t},a_{t})+\gamma Q(s_{t+1},a_{t+1}|\theta)$ (13)

In DDPG, deterministic policy is considered. For deterministic policy
$a=\mu(s)$, the inner expectation in Eq.10 could be eliminated as follows.

$Q^{\mu}(s_{t},a_{t})=\mathbb{E}_{r_{t},s_{t+1}\sim E}[g(s_{t},a_{t})+\gamma
Q^{\mu}(s_{t+1},\mu(s_{t+1}))]$ (14)

Therefore, when calculating the loss in Eq.12, transitions obtained from
different stochastic policy could be used. The collected transition data
$(s_{t},a_{t},r_{t},s_{t+1})$ is stored in a replay buffer, and then minibatch
is sampled from the replay-buffer for loss calculation in Eq.12. This enables
the samples for training to be distributed identically and independently which
is the basic assumption for neural network training. When exploring the
environment for collecting transition data, noise is added to the actor
policy.

$\bar{\mu}(s_{t})=\mu(s_{t},\theta_{t})+\mathcal{N}$ (15)

In addition, target networks are introduced in order to stabilize the training
process and achieve greater convergence. When calculating Eq.13 with collected
transition data $(s_{i},a_{i},r_{i},s_{i+1})$, target critic with parameter
$\theta^{\prime}$, and target policy with parameter $\phi^{\prime}$ are used.

$y_{i}=r_{i}+\gamma
Q^{\prime}(s_{i+1},\mu^{\prime}(s_{i+1}|\phi^{\prime})|\theta^{\prime})$ (16)

Target networks are updated in order to slowly track the learned networks, as
follows

$\displaystyle\theta^{\prime}$
$\displaystyle\leftarrow\tau\theta+(1-\tau)\theta^{\prime}$ (17)
$\displaystyle\phi^{\prime}$
$\displaystyle\leftarrow\tau\phi+(1-\tau)\phi^{\prime}$ (18)

where $\tau<<1$.

The key advantage of the DDPG method is that it enables to train deep
reinforcement learning for continuous action space off-policy by assuming a
deterministic policy. This greatly improves the sample efficiency because
transition data obtained by different policies in the past could be re-used
for training the agent.

###  4.3 Twin Delayed Deep Deterministic Policy Gradient (TD3)

While DDPG has achieved significant performance in continuous control
problems, the critical drawback is that it was sensitive to hyper parameter
settings. This is because the overestimation of the Q-values makes the policy
fall into local optima. The Twin Delayed Deep Deterministic Policy Gradient
(TD3) algorithm made three major changes in the DDPG algorithm to overcome
this drawback. [35]. The first change (clipped double Q-learning) is to use
two seperate Q value approximators (critics), and use the smaller Q value of
the two networks to calculate the target value $y_{t}$. The second change is
target policy smoothing. For action in the target Q value, clipped noise is
added to the actual action. This prevents the Q function from having sharp
peaks by returning similar Q values for similar actions. With the clipped
double Q-learning and target policy smoothing, the target value from
transition data $(s_{i},a_{i},r_{i},s_{i+1})$ could be calculated as follows.

$y_{i}=r_{i}+\gamma\min_{i=1,2}Q^{\prime}(s_{i+1},\mu^{\prime}(s_{i+1}|\phi^{\prime})+\epsilon|\theta_{i}^{\prime}),\quad\epsilon\sim
clip(\mathcal{N}(0,\sigma),-c,c)$ (19)

The last change is the delayed updates of the policy and the target networks.
In order to stabilize the learning process, updates of the target networks and
policy networks are carried out less frequently (example: once per every $N$
steps) than the critic and networks.

###  4.4 TD3 with memory

When the system is partially observed, optimal policy and action-value
function both become function of the entire observation-action history
$h_{t}$. Recurrent Deterministic Policy Gradient (RDPG) algorithm uses
recurrent neural networks (RNN) in the policy and the action-value networks to
preserve (limited) information of the history as part of the state. [36] In
the RDPG setting, $\mu(s)$ and $Q(s,a)$ could be replaced with $\mu(h)$ and
$Q(h,a)$. Thus, the policy update could be obtained as follows.

$\nabla_{\phi}J(\phi)=\mathbb{E}_{\tau\sim\nu_{\mu}}\biggl{[}\sum_{t}\gamma^{t-1}\nabla_{a}Q(h,a|\theta)|_{h=h_{t},a=\mu(h_{t})}\nabla_{\phi}\mu(h|\phi)|_{h=h_{t}}\biggr{]}$
(20)

The critic and actor networks are where
$\tau=(s_{0},o_{0},a_{0},s_{1},o_{2},a_{2},\ldots)$ are drawn from the
trajectory distribution $\nu_{\mu}$ generated by the current policy $\mu$ . In
this paper, the RDPG algorithm was implemented by adding LSTM networks in the
critic and actor networks. LSTM is a RNN architecture that is composed of a
cell, an input gate, an output gate, and a forget gate. The critic and actor
agent inputs the observation and action history $h_{t}$ to their LSTM, and
uses the output of the output gate $\tilde{h_{t}}$ to calculate their outputs.
The network structure for the critic and the actor is shown in Fig.8,9. This
architecture is based on the works of Peng (2018) of robotic control. The
upper network in the figure uses the current observation (and action for
critic), while the bottom network utilizes past information stored in the
LSTM. The network was trained with the TD3 algorithm, as shown in Algorithm 3.
The difference from the original paper is that in our implementation, we
created multiple replay buffers to store transition sequences of different
episode length $T$. In our environment, the episode length varied from 4-10
time steps. We stored transition data with different length to separate
replay-buffers so that we can create mini-batches with same episode length
when training.

Figure 8: Architecture of the policy network Figure 9: Architecture of the
critic network Algorithm 3 TD3 with LSTM

Initialize Critic Network $Q_{\theta_{1}},Q_{\theta_{2}}$, and actor network
$\mu_{\phi}$ with random parameters
$\theta_{1}^{\prime},\theta_{2}^{\prime},\mu^{\prime}$

Initialize target networks
$\theta^{\prime}_{1}\leftarrow\theta_{1},\theta^{\prime}_{2}\leftarrow\theta_{2},\phi^{\prime}\leftarrow\phi$

Initialize list of replay buffer list $\mathcal{B}[]$

for episodes = 1, M do

initialize empty history $h_{0}$

$a_{0}\leftarrow\textsc{Random}(\mathcal{A}),\quad t\leftarrow 0$

while not done do

$t\leftarrow t+1$

receive observation $o_{t}$, reward $r_{t-1}$

$h_{t}\leftarrow h_{t-1},a_{t-1},o_{t}$ $\triangleright$ Append observation
and previous action to history

Select action
$a_{t}=\mu(h_{t}|\phi)+\epsilon,\quad\epsilon\sim\mathcal{N}(0,\sigma)$

If done $d_{t}\leftarrow 1$, else $d_{t}\leftarrow 0$

end while

Store the sequence
$(a_{0},o_{1},r_{0},d_{0},\ldots,a_{t-1},o_{t},r_{t},d_{t})$ in
$\mathcal{B}[t]$

$T\leftarrow\textsc{Random}(t_{max})$ $\triangleright$ Select the sequence
length to sample from

Sample a mini-batch of $N$ episodes with length $T$:

$(a_{0}^{i},o_{1}^{i},r_{0}^{i},d_{1}^{i},\ldots,a_{T-1}^{i},o_{T}^{i},r_{T-1}^{i},d_{T}^{i})_{i=1,\ldots,N}$
from $\mathcal{B}[T]$

Construct $T\times N$ set of sequences
$h_{t}^{i}=(a_{0}^{i},o_{1}^{i},r_{0}^{i},d_{1}^{i}\ldots,a_{t-1}^{i},o_{t}^{i},r_{t-1}^{i},d_{t}^{i})$

Compute target values for each sample episode for $t=0,\ldots,T-1$

$\quad\quad
y_{t}^{i}=r_{t}^{i}+(1-d_{t+1}^{i})\gamma\min_{k=1,2}Q^{\prime}(h_{t+1}^{i},\mu^{\prime}(h_{t+1}^{i}|\phi^{\prime})+\epsilon|\theta_{k}^{\prime}),\quad\epsilon\sim
clip(\mathcal{N}(0,\sigma),-c,c)$

Update critic by minimizing the loss

$\quad\quad\theta_{k}\leftarrow
argmin_{\theta_{k}}\frac{1}{NT}\sum_{i}\sum_{t}(y_{t}^{i}-Q(h_{t}^{i},a_{t}^{i}))^{2}$

Update actor parameter $\phi$ by the deterministic policy gradient

$\quad\quad\nabla_{\phi}J(\phi)=\frac{1}{NT}\sum_{i}\sum_{t}\nabla_{a}Q(h,a|\theta)|_{h=h_{t}^{i},a=\mu(h_{t}^{i})}\nabla_{\phi}\mu(h|\phi)|_{h=h_{t}^{i}}$

end for

###  4.5 Overall Framework and Training Process

The observation $o_{t}$ which is used as an input of the agent is summarized
in Table.1 and the output of the agent $a_{t}$ is summarized in Table.2. For
the training agent, we tested two architectures: a simple TD3 algorithm for
memory-less approach, and a TD3 with LSTM approach for memory-based approach.
For the architecture of TD3 agent without memory, we used the same
architecture as the upper network of the policy and critic networks shown in
Fig.8, 9. In order to stabilize the training, all action parameters were
scaled to [0,1] during training. In the actor, the target point at timestep
$t$ ($=x_{f}^{(t)},y_{f}^{(t)}$) is not outputted directly. Instead 2
variables $\alpha_{x},\alpha_{y}\in[-0.25,0.25]$ is outputted from the actor,
and calculated the next target point using the following equation.

$x_{f}^{(t)}=x_{f}^{(t-1)}+\alpha_{x}w_{x},\quad
y_{f}^{(t)}=y_{f}^{(t-1)}+\alpha_{y}w_{y}$ (21)

In this way, the next target landing point is always selected to be within the
FOV of the latest observed DEM. This prevents control failures caused by large
deviations of the target between observations. The range of $\alpha_{x}$ and
$\alpha_{y}$ was set to $\alpha_{x},\alpha_{y}\in[-0.25,0.25]$. The calculated
new target point will be the center of the FOV of the DEM in the next step.
Therefore, $\alpha_{x}$ and $\alpha_{y}$ are important output that has impacts
on both guidance and observation. In addition, the actor does not directly
output $t_{go}$ of the next step, but instead output the degrade of the time-
to-go until the next observation timing which will be conducted 5 seconds
later ($=\delta t_{go}$). The initial $t_{go}$ is calculated using the initial
lander position and target landing site and Eq.6.

Table 1: Elements of the observed state $o_{t}$ content | symbol | size
---|---|---
(true) lander position | $r_{t}$ | 3
(true) lander velocity | $v_{t}$ | 3
(true) lander mass | $m_{t}$ | 1
width of the observed DEM range in x,y direction | $w_{x}^{(t)},w_{y}^{(t)}$ | 2
current target position | $x_{f}^{(t)},y_{f}^{(t)},z_{f}^{(t)}$ | 3
latent vector of the generated safety map encoded through autoencoder | $z_{t}$ | 32
total |  | 44
Table 2: Elements of the action $a_{t}$ content | symbol | size | range
---|---|---|---
gain of the ZEM-ZEV control | $K_{R}$ | 1 | [5, 7]
gain of the ZEM-ZEV control | $K_{V}$ | 1 | [1, 3]
degradation of time-to-go from previous timestep | $\delta t_{go}$ | 1 | [4.25,5.75]
target landing point position within the DEM FOV | $\alpha_{x},\alpha_{y}$ | 2 | [-0.5, 0.5]
total |  | 5 |

Reward setting is the key in reinforcement learning. The reward at time step
$t$ was implemented as follows.

$\displaystyle\begin{aligned}
r_{t}&=\alpha_{m}\frac{m_{t-1}-m_{t}}{m_{0}-m_{dry}}+\alpha_{f}(U(z_{max}-z_{t})+U(m_{dry}-m_{t}))+\alpha_{v}d_{t}(|\boldsymbol{v_{t}}|)\\\
&+\alpha_{r}d_{t}(|\boldsymbol{r_{t}}-\boldsymbol{r_{f}^{(t-1)}}|)+\alpha_{s}d_{t}V_{D}(r_{t})\end{aligned}$
(22)

where

$\displaystyle\begin{aligned} U(x)&=\begin{cases}1&\text{if $x\geq 0$}\\\
0&\text{otherwise}\end{cases}\\\ d_{t}&=\begin{cases}1&\text{if episode
done}\\\ 0&\text{otherwise}\end{cases}\\\
V_{D}(r_{t})&=\begin{cases}1&\text{if $r_{t}$ is a safe landing point}\\\
-1&\text{if $r_{t}$ is not a safe landing point}\end{cases}\\\
z_{max}&=50\quad\text{(maximum height of the terrain)}\end{aligned}$ (23)

The term with $\alpha_{m}$ is a penalty for fuel consumption. The term with
$\alpha_{f}$ is a penalty term for breaking the two critical path constraints:
hitting the ground and running out of fuel. The term with $\alpha_{v}$ is a
penalty for velocity norm at final episode to ensure soft landing, and the
term with $\alpha_{r}$ is a penalty for not landing to the final target
landing point. These two terms are required because ZEM-ZEV control with
adaptive gain and $t_{go}$ parameters are not guaranteed to achieve a pinpoint
soft landing to the target. Setting the weights
$\alpha_{m},\alpha_{v},\alpha_{r},\alpha_{s}$ is a difficult problem which
depends not only on the mission designer focus but the stability of the
training process. It is natural to assume that $\alpha_{f},\alpha_{v}$ and
$\alpha_{s}$ should take relatively large values because hitting the ground at
a high rate of speed and landing in a danger zone could both lead to the
immediate loss of the lander. We used
$\alpha_{m}=1,\alpha_{f}=-10,\alpha_{v}=0.1,\alpha_{r}=0.01,\alpha_{s}=1$ in
our research. Note that slope angle constraints are not considered as penalty
in the reward settings. This is because slope angle constraints are not
directly related to the primary objective of the planetary landing, but rather
conservative constraints to assure preferable geometry for observation. By
interacting with the simulator that simulates the observation process and the
dynamics, we expect the agent to learn to take a trajectory with adequate
slope angle in order to achieve safe landing. The overall framework is
summarized in Fig 10.

Figure 10: Overall framework

##  5 Simulation

###  5.1 Initial Condition and Hyper-parameter Settings

The initial state and other specs of the lander are summarized in Table.3.
Parameters that have range are randomly initialized within the range every
time the episode starts. The true DTM is picked randomly from the list of DTMs
generated before the training.

Table 3: Simulation condition content | symbol | value
---|---|---
Initial altitude (downrange) [m] | $z_{0}$ | [900, 1000]
Initial crossrange distance [m] | $\sqrt{x_{0}^{2}+y_{0}^{2}}$ | [200 250]
Initial downrange velocity [m/s] | $|v_{z}|$ | [20, 35]
Initial crossrange velocity [m/s] | $|v_{z}|$ | [5, 10]
Initial mass of the lander [kg] | $m_{0}$ | 1200
Altitude of the final target point [m] | $z_{f}$ | 50
Dry mass of the lander [kg] | $m_{dry}$ | 1150
Specific Impulse of the lander [s] | $I_{sp}$ | 325
Maximum thrust of the lander [N] | $T_{max}$ | 12000
Thrust magnitude error ratio | $|T-T_{plan}|/T_{plan}$ | 0.05
Number of different DTMs used for training | $N_{D-training}$ | 70
Number of different DTMs used for variation | $N_{D-variation}$ | 18

###  5.2 Training Process

Fig.11 illustrates the transition of the training loss of the critic and actor
network, reward, and the ratio of landing to safe landing site during
training, for both the agent without LSTM and with LSTM. Values in the graph
represent the moving average of 5000 episodes. When validating the performance
during training, DTM datasets that are different from training datasets were
used. It required around 250k iteration of training for both agents to reach
around the maximum performance.

Contrary to our expectations, the agent without LSTM showed better performance
compared with the agent with LSTM. The training of agent with LSTM was much
more sensitive to its hyper-parameters, and we did not manage to achieve
higher performance than TD3 without memory. Another reason why agent without
LSTM could achieve high performance is likely because the partial
observability of our problem is not strong, thanks to the access to the
stochastic safety map with relatively low errors. Adding larger errors into
the observation data such as navigation errors might require the use of memory
to cope with strong partial observability.

Figure 11: Training log. It should be noted that the x-axis of the above two
figures represent the number of training updates, while the x-axis of the
bottom two figures represent the number of training episodes.

###  5.3 Comparison with other strategies

In order to assess the performance of the trained agent, the performance of
the two agents were compared with two other strategies. The first strategy is
the ”Fixed Control Policy”. This policy ignores the $K_{R},K_{V},t_{go}$
output from the agent, and uses the agent only to select the target point. The
gains are fixed to $K_{R}=6,K_{V}=2$, and $t_{go}$ was calculated by Eq.6.
This policy was introduced to see how the agent’s adjustment of gain and
flight time affected the trajectory, final velocity, and landing errors to the
target. The second strategy is the ”Single Divert Policy”. This strategy makes
no use of trained agents. This policy targets the initial target point until
the altitude reaches 500m, and then selects the pixel with maximum safety
value in the stochastic safety map obtained at 500m altitude as the final
target point. The gains are fixed to $K_{R}=6,K_{V}=2$, and $t_{go}$ was
calculated by Eq.6. This policy was introduced to represent a simplified
version of existing landing techniques.

The performance of the four policies was tested by landing simulation of 1000
episodes with random initial conditions. In order to assess the performance in
cases that the agent needs to change its initial landing site, the initial
target position was chosen so that at least 80% of the area within 100 meter
radius is hazardous. Fig.12 shows the histogram of the reward, landing error
to the final target point[m], final velocity magnitude[m/s], maximum
thrust[N], and minimum slant angle [deg]. The mean value is summarized in
Table.4.

The trained memory-less agent achieved maximum ratio ($94.8\%$) of landing to
safe landing site. The fixed control policy had a lower safe landing ratio to
the agent policy, due to its limitation in changing the thrust magnitude
flexibly when long divert maneuver is required. Fig.13 shows an example of
such a situation. In this case, the initial target point (the center of the
red square) is inside the hazardous region (black area of the figure), and the
target landing point has to be shifted toward upside of the figure in order to
achieve a safe landing. While the trained agent reduces the $z$ axis velocity
immediately as shown in Fig.13-(b) in order to maintain the altitude and wide
field of view for searching safe landing sites, the fixed controller does not
aggressively thrust the propellant as shown in Fig.13-(d). Therefore the field
of view of the fixed gain controller is kept limited, and the fixed controller
fails in finding a wide safe area to land.

The trained agent with memory had the lowest probability of landing to safe
landing site, and the largest landing position error and velocity error among
the four policies. On the contrary, the average fuel consumption was the
lowest among the four policies. This implies that the agent was trapped in a
local optima policy that increases reward by reducing fuel consumption rather
than improving other metrics.

The single observation policy fails in selecting safe landing sites when safe
landing sites could not be obtained in the single observation due to
limitations in the field of view or slant angles. The rate of failure may be
reduced by selecting optimal altitude for the divert decision or implementing
better landing site selection algorithms from safety maps. However, it should
be noted that in reality, the risk of single observation policy is likely to
increase than our simulator due to additional error sources (navigation
errors, attitude constraints, etc) which are not simulated in our research.
These errors degrade the quality of the DEM (or safety map) and the accuracy
of lander guidance to the target point.

Figure 12: Evaluation of performance. ’Agent’ is the policy of the trained memory-less agent, ’LSTM’ is the policy of the trained memory-based agent, ’Fixed’ is the policy with fixed ZEM-ZEV parameters, and ’Single’ is the policy that conducts single divert maneuver and selects the landing site by directly using the stochastic safety map. Table 4: Comparison of performance between methods (1000 episodes. Mean value) Criteria | Memory-Less Agent | Memory-Based | Fixed Control | Single Divert
---|---|---|---|---
Safe Landing Ratio [%] | 94.8 % | 88.3 % | 93.2 % | 89.7 %
Distance to final target [m] | 2.492 | 4.547 | 3.786 | 0.005
Final Velocity [m/s] | 0.637 | 1.520 | 1.189 | 0.557
Propellant Consumption [%] | 69.88 | 63.70 | 67.84 | 68.89
Minimum Slant Angle [deg] | 72.429 | 72.570 | 72.795 | 70.85
Maximum Thrust [N] | 5818 | 6863 | 4553 | 4556

2

(a) Trajectory by agent controller (White:Safe Black:Unsafe)

(b) Control history by the agent controller (Up:Velocity Down:Thrust)

(c) Trajectory by fixed controller (White:Safe Black:Unsafe)

(d) Control history by the fixed controller ((Up:Velocity Down:Thrust)

Figure 13: Comparison of the obtained agent and a ZEM-ZEV controller with
fixed gain. Red line in figure(a),(c) shows the lander trajectory, while the
square shows the field of view (FOV) of the LIDAR DEMs. Figure 14: Changes in
the control gain $K_{R},K_{V}$ while control in Fig.13-(a),(b)

While the trained agent showed high potential in adaptly changing its target
position and control outputs to find safe landing sites, they also have
limitations in reducing velocities and landing error to the final target point
decided by the agent. Since convergence to the target position and velocity is
not guaranteed when control gain in ZEM-ZEV control is changed, this
disadvantage is inevitable without further refinements. In addition, fuel
consumption of the memory-less trained agent was larger compared with the two
agents for comparison, due to the large thrust magnitude during descent and
frequent changes of the target point. Addressing the trade-off between safety
and fuel consumption from the point of the frequency of the target point
changes is an interesting future work.

##  6 Conclusion

In this paper, we proposed a new learning-based framework for Hazard Detection
and Avoidance (HDA) phase which successively updates the target landing site
and control parameters simultaneously after each observation, in order to cope
with the coupling between observation, guidance, and control. We modeled the
HDA sequence as a POMDP, and a reinforcement learning agent that interprets
the obtained map using auto-encoder and outputs control parameters for ZEM-ZEV
feedback control law was developed to find an optimal policy. The agent was
trained by interacting with the simulator, and the trained agent was able to
achieve over 90% probability of successful landing at difficult landing sites
where over 80% of the terrain was hazardous around the initial target point,
by gradually updating the target landing point towards safe regions.

In order to incorporate a more realistic and higher level of uncertainty to
the environment, our future work will mainly focus on the refinement of the
simulator model. This includes the incorporation of accurate LIDAR DEM
generation models, expansion of dynamics to 6-DOF, and incorporation of lander
state measurement errors. As the partial observability of the environment
increases by incorporating various error sources, a memory-based agent might
be required for sufficient performance. We are also planning to test other
baseline control policies and observation data interpretation architectures to
improve the agent’s performance and stability.

##  7 Acknowledgements

This material is partially based upon work supported by the National
Aeronautics and Space Administration under Grant No.80NSSC20K0064 through the
NASA Early Career Faculty Program.

## References

  * [1] A. E. Johnson, A. R. Klumpp, J. B. Collier, and A. A. Wolf, “Lidar-based hazard avoidance for safe landing on Mars,” Journal of guidance, control, and dynamics, Vol. 25, No. 6, 2002, pp. 1091–1099.
  * [2] N. Serrano and H. Seraji, “Landing site selection using fuzzy rule-based reasoning,” Proceedings 2007 IEEE International Conference on Robotics and Automation, 2007, pp. 4899–4904.
  * [3] A. Huertas, Y. Cheng, and L. H. Matthies, “Real-time hazard detection for landers,,” NASA Science Technology Conference, 2007.
  * [4] L. Matthies, A. Huertas, Y. Cheng, and A. Johnson, “Stereo vision and shadow analysis for landing hazard detection,” 2008 IEEE International Conference on Robotics and Automation, 2008, pp. 2735–2742.
  * [5] Y. Cheng, A. E. Johnson, L. H. Mattheis, and A. A. Wolf, “Passive imaging based hazard avoidance for spacecraft safe landing,” I-SIRAS 2001: 6th International Symposium on Artifical Intelligence, Robotics and Automation in Space, 2001, pp. 1–14.
  * [6] N. Serrano, “A bayesian framework for landing site selection during autonomous spacecraft descent,” 2006 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2006, pp. 5112–5117.
  * [7] C. I. Restrepo, R. Lovelace, R. R. Sostaric, J. M. Carson, and N. G. Spaceflight, “NASA SPLICE Project : Developing the Next Generation Hazard Detection System,” 2nd RPI Space Imaging Workshop, 2019, pp. 2–3.
  * [8] A. D. Cianciolo, S. Striepe, J. Carson, R. Sostaric, D. Woffinden, C. Karlgaard, R. Lugo, R. Powell, and J. Tynis, “Defining navigation requirements for future precision lander missions,” AIAA Scitech 2019 Forum, No. January, 2019, pp. 1–18, 10.2514/6.2019-0661.
  * [9] J. M. Carson, M. M. Munk, R. R. Sostaric, J. N. Estes, F. Amzajerdian, J. Bryan Blair, D. K. Rutishauser, C. I. Restrepo, A. D. Cianciolo, G. T. Chen, and T. Tse, “The splice project: Continuing nasa development of gn&c technologies for safe and precise landing,” AIAA Scitech 2019 Forum, No. January, 2019, pp. 1–9, 10.2514/6.2019-0660.
  * [10] T. Brady, E. Robertson, S. P. C. App, and D. Zimpfer, “Hazard detection methods for lunar landing,” 2009 IEE Aerospace Conference, 2009, pp. 1–18.
  * [11] C. D. App and T. B. Smith, “Autonomous precision landing and hazard detection and avoidance technology (ALHAT),” 2007 IEEE Aerospace Conference, 2007, pp. 1–7.
  * [12] T. Ivanov, A. Huertas, and J. M. Carson, “Probabilistic Hazard Detection for Autonomous Safe Landing,” AIAA Guidance, Navigation, and Control (GNC) Conference, 2013, p. 5019.
  * [13] T. Brady and J. Schwartz, “ALHAT system architecture and operational concept,” 2007 IEEE Aerospace Conference, 2007, pp. 1–13.
  * [14] S. A. Striepe, C. D. Epp, and E. A. Robertson, “Autonomous precision landing and hazard avoidance technology (ALHAT) project status as of May 2010,” International Planetary Probe Workshop 2010 (IPPW-7), 2010.
  * [15] D. Rutishauser, C. Epp, and E. Robertson, “Free-Flight Terrestrial Rocket Lander Demonstration for NASA’s Autonomous Landing and Hazard Avoidance Technology (ALHAT) System,” AIAA SPACE 2012 Conference & Exposition, 2012, p. 5239.
  * [16] C. Epp, E. Robertson, and J. M. Carson, “Real-time hazard detection and avoidance demonstration for a planetary lander,” AIAA SPACE 2014 Conference and Exposition, 2014, p. 4312.
  * [17] A. E. Johnson, A. Huertas, R. A. Werner, and J. F. Montgomery, “Analysis of on-board hazard detection and avoidance for safe lunar landing,” 2008 IEEE Aerospace Conference, 2008, pp. 1–9.
  * [18] A. Huertas, A. E. Johnson, R. A. Werner, and R. A. Maddock, “Performance evaluation of hazard detection and avoidance algorithms for safe Lunar landings,” 2010 IEEE Aerospace Conference, 2010, pp. 1–20.
  * [19] B. E. Cohanim and B. K. Collins, “Landing point designation algorithm for lunar landing,” Journal of Spacecraft and Rockets, Vol. 46, No. 4, 2009, pp. 858–864.
  * [20] S. Paschall and T. Brady, “Demonstration of a safe & precise planetary landing system on-board a terrestrial rocket,” 2012 IEEE Aerospace Conference, 2012, pp. 1–8.
  * [21] N. Trawny, A. Huertas, M. E. Luna, C. Y. Villalpando, K. Martin, J. M. Carson, A. E. Johnson, C. Restrepo, and V. E. Roback, “Flight testing a Real-Time Hazard Detection System for Safe Lunar Landing on the Rocket-Powered Morpheus Vehicle,” AIAA Guidance, Navigation, and Control (GNC) Conference, 2015\.
  * [22] N. Serrano, “A Bayesian framework for landing site selection during autonomous spacecraft descent,” IEEE International Conference on Intelligent Robots and Systems, 2006, pp. 5112–5117, 10.1109/IROS.2006.282603.
  * [23] S. R. Ploen, H. Seraji, and C. E. Kinney, “Determination of spacecraft landing footprint for safe planetary landing,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 45, No. 1, 2009, pp. 3–16, 10.1109/TAES.2009.4805259.
  * [24] P. Cui, D. Ge, and A. Gao, “Optimal landing site selection based on safety index during planetary descent,” Acta Astronautica, Vol. 132, No. July 2016, 2017, pp. 326–336, 10.1016/j.actaastro.2016.10.040.
  * [25] E. S. Crane and S. M. Rock, “Guidance augmentation for reducing uncertainty in vision-based hazard mapping during lunar landing,” IEEE Aerospace Conference Proceedings, 2013, 10.1109/AERO.2013.6496946.
  * [26] M. Barker, E. Mazarico, G. Neumann, M. Zuber, J. Haruyama, and D. Smith, “A new lunar digital elevation model from the Lunar Orbiter Laser Altimeter and SELENE Terrain Camera,” Icarus, Vol. 273, 2016, pp. 346 – 355, https://doi.org/10.1016/j.icarus.2015.07.039.
  * [27] R. J. Pike, “Size-dependence in the shape of fresh impact craters on the moon.,” Impact and Explosion Cratering: Planetary and Terrestrial Implications (D. J. Roddy, R. O. Pepin, and R. B. Merrill, eds.), Jan. 1977, pp. 489–509.
  * [28] D. E. Bernard and M. P. Golombek, “Crater and rock hazard modeling for Mars landing,” AIAA Space 2001 Conference and Exposition, Albuqurque, NM, 2001, 10.2514/6.2001-4697.
  * [29] C. I. Restrepo, P.-T. Chen, R. R. Sostaric, and J. M. Carson, “Next-Generation NASA Hazard Detection System Development,” 2020, pp. 1–10, 10.2514/6.2020-0368.
  * [30] B. Gaudet, R. Linares, and R. Furfaro, “Deep Reinforcement Learning for Six Degree-of-Freedom Planetary Landing,” Advances in Space Research, 01 2020, 10.1016/j.asr.2019.12.030.
  * [31] B. Ebrahimi, M. Bahrami, and J. Roshanian, “Optimal sliding-mode guidance with terminal velocity constraint for fixed-interval propulsive maneuvers,” Acta Astronautica, Vol. 62, No. 10, 2008, pp. 556 – 562, https://doi.org/10.1016/j.actaastro.2008.02.002.
  * [32] R. Furfaro, A. Scorsoglio, R. Linares, and M. Massari, “Adaptive Generalized ZEM-ZEV Feedback Guidance for Planetary Landing via a Deep Reinforcement Learning Approach,” 2020, 10.1016/j.actaastro.2020.02.051.
  * [33] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning.,” ICLR (Y. Bengio and Y. LeCun, eds.), 2016.
  * [34] D. Silver, G. Lever, N. Heess, T. Degris, D. Wierstra, and M. Riedmiller, “Deterministic Policy Gradient Algorithms,” Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32, ICML’14, JMLR.org, 2014, p. I–387–I–395.
  * [35] S. Fujimoto, H. Hoof, and D. Meger, “Addressing Function Approximation Error in Actor-Critic Methods,” International Conference on Machine Learning, 2018, pp. 1582–1591.
  * [36] N. M. O. Heess, J. J. Hunt, T. P. Lillicrap, and D. Silver, “Memory-based control with recurrent neural networks,” ArXiv, Vol. abs/1512.04455, 2015\.