# Cumulative Residual Extropy of Minimum Ranked Set Sampling with Unequal
Samples

Mohammad Reza Kazemi  Saeid Tahmasebi  Camilla Calì  Maria Longobardi
<EMAIL_ADDRESS>Department of Statistics, Faculty of Science, Fasa
University, Fasa, Iran  Department of Statistics, Persian Gulf University ,
Bushehr, Iran  Department of Biology, University of Napoli Federico II,
Napoli, Italy

###### Abstract

Recently, an alternative measure of uncertainty called cumulative residual
extropy (CREX) was proposed by Jahanshahi et al. [2019]. In this paper, we
consider uncertainty measures of minimum ranked set sampling procedure with
unequal samples (MinRSSU) in terms of CREX and its dynamic version and we
compare the uncertainty and information content of CREX based on MinRSSU and
simple random sampling (SRS) designs. Also, using simulation, we study on new
estimators of CREX for MinRSSU and SRS designs in terms of bias and mean
square error. Finally, we provide a new discrimination measure of disparity
between the distribution of MinRSSU and parental data SRS.

###### keywords:

Cumulative residual extropy, Discrimination measure, Minimum ranked set
sampling, Stochastic ordering.

###### MSC:

[2010] 62B10 , 60E15 , 62D05 , 94A17

††journal: Journal of LaTeX Templates

##  1 Introduction

Ranked set sampling (RSS) design is a cost-effective sampling for situations
where taking actual measurements on units is expensive but ranking units is
easy. For the first time, based on the RSS sampling design, McIntyre [1952]
provided a more efficient estimator of the population mean comparing to the
simple random sampling (SRS) counterpart. To learn more about this concept,
the readers can refer to Patil et al. [1999]. There are many available studies
that have developed and generalized the method of sampling used in RSS scheme
and they efficiently estimate the population parameter comparing to the SRS
scheme. Recently, Qiu & Eftekharian [2020] studied information content of
minimum ranked set sampling procedure with unequal samples (MinRSSU) as useful
modification of RSS procedure in terms of extropy. In the MinRSSU, we draw $m$
simple random samples, where the size of the $i$th samples is $i$,
$i=1,...,m$. The one-cycle MinRSSU involves an initial ranking of $m$ samples
of size $m$ as follows:

$\begin{array}[]{ccccccc}1:&\underline{\boldsymbol{X_{(1:1)1}}}&&&&\rightarrow&\tilde{X}_{1}=X_{(1:1)1}\\\
2:&\underline{\boldsymbol{X_{(1:2)2}}}&\
X_{(2:2)2}&&&\rightarrow&\tilde{X}_{2}=X_{(1:2)2}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\\
m:&\underline{\boldsymbol{X_{(1:m)m}}}&X_{(2:m)m}&\cdots&X_{(m:m)m}&\rightarrow&\tilde{X}_{m}=X_{(1:m)m}\\\
\end{array}$

where $X_{(i:i)j}$ denotes the $i$th order statistic from the $j$th SRS of
size $i$. The resulting sample is called one-cycle MinRSSU of size $m$ and
denoted by $\boldsymbol{X}^{(m)}_{MinRSSU}=\\{\tilde{X}_{i},i=1,\dots,m\\}$.
The parameter $m$ should be kept small because the ranking should not be
difficult in this sense and the ranking may be done for example by using an
easily measurable covariate, then it is not difficult to identify the minimum
of ranked individuals in each subset. Note that $\tilde{X}_{i}$ has the same
distribution as $X_{(1)i}$ which is the smallest order statistic in a set of
size $i$ with probability density function (pdf)
$f_{(1)i}(x)=if(x)[1-F(x)]^{i-1}$ and survival function
$\bar{F}_{(1)i}(x)=[1-F(x)]^{i}=\bar{F}^{i}(x)$, where $f(.)$, $F(.)$ and
$\bar{F}(.)$ are the underlying pdf, cumulative distribution function (cdf)
and survival function. In reliability theory, $\tilde{X}_{i}$ measures the
lifetime of a series system.
Several authors have worked on measures of information for RSS and its
variants. Jozani & Ahmadi [2014] explored the notions of information content
of RSS data and compared them with their counterparts in SRS data. Tahmasebi
et al. [2016] obtained some results of residual (past) entropy for ranked set
samples. Eskandarzadeh et al. [2016] studied information measures for record
ranked set sampling. Eskandarzadeh et al. [2018] considered information
measures of maximum ranked set sampling procedure with unequal samples in
terms of Shannon entropy, Rényi entropy and Kullback-Leibler information,
instead Tahmasebi et al. [2020] in terms of Tsallis entropy. More recently,
Qiu & Eftekharian [2020] studied information content of MinRSSU in terms of
extropy and Raqab & Qiu [2019] considered the problems of uncertainty and
information content of RSS data based on extropy measure and the related
monotonic properties and stochastic comparisons.
Let $X$ denotes a continuous random variable with pdf $f$. Frank et al. [2015]
introduced a new measure termed by extropy associated with $X$ as

$J(X)=-\frac{1}{2}\int_{-\infty}^{+\infty}[f(x)]^{2}dx=-\frac{1}{2}\int_{0}^{1}f(F^{-1}(u))du,$
(1.1)

where $F^{-1}(.)$ is the quantile function of $X$. Qiu [2017] explored some
characterization results, monotone properties, and lower bounds of extropy of
order statistics and record values. Also, Qiu & Eftekharian [2020] and Raqab &
Qiu [2019] considered the information measure of extropy $J(X)$ based on
MinRSSU and RSS schemes, respectively and compared the results with their
counterpart under SRS design.
This paper is organized as follows: Section 2 deals with the results of
cumulative residual extropy (CREX) for MinRSSU data by comparing to its
counterpart under SRS data. In Section 3, new estimators are proposed for CREX
in SRS and MinRSSU designs using empirical approach. Also, by using simulation
study, the behavior of estimators of CREX in MinRSSU and SRS are compared in
terms of bias and mean square error. Furthermore, we show that how MinRSSU
scheme can efficiently reduce the uncertainty measure comparing to SRS design.
Section 4 provides a new discrimination measure of disparity between the
distribution of MinRSSU and parental data SRS. Section 5 concludes the paper.

##  2 Cumulative residual extropy of MinRSSU

Let $X$ denotes the lifetime of a system with survival function $\bar{F}$.
Recently, a new measure of information is proposed by Jahanshahi et al. [2019]
with substituting the function $\bar{F}$ in extropy formula (1.1). This new
measure is called CREX and defined as

$\displaystyle{\mathcal{\xi
J}}(X)=-\frac{1}{2}\int_{0}^{+\infty}\bar{F}^{2}(x)dx.$ (2.1)

Note that $-\infty<{\mathcal{\xi J}}(X)\leq 0$. If the CREX of $X$ is less
than that of another random variable, say $Y$, i.e. ${\mathcal{\xi
J}}(X)\leq{\mathcal{\xi J}}(Y)$, then $X$ has less uncertainty than $Y$. Now
let ${\mathcal{\xi J}}(X)<+\infty$. Then, for the MinRSSU and SRS designs, we
have

$\displaystyle{\mathcal{\xi J}}(\boldsymbol{X}^{(m)}_{MinRSSU})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\prod\limits_{i=1}^{m}[-2{\mathcal{\xi
J}}(X_{(1:i)})]=-\frac{1}{2}\prod\limits_{i=1}^{m}\int_{0}^{+\infty}\bar{F}^{2i}(x)dx$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\prod\limits_{i=1}^{m}\int_{0}^{1}\frac{(1-u)^{2i}}{f(F^{-1}(u))}du$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\prod\limits_{i=1}^{m}\mathbb{E}\left[\frac{(1-U)^{2i}}{f(F^{-1}(U))}\right],$
(2.3)

and

$\mathcal{\xi
J}(\boldsymbol{X}^{(m)}_{SRS})=-\frac{1}{2}\left[\int_{0}^{+\infty}\bar{F}^{2}(x)dx\right]^{m}=-\frac{1}{2}[-2{\mathcal{\xi
J}}(X)]^{m}.$ (2.4)

To compare the above measures, let us consider the following examples.

######  Example 2.1.

If $U\sim Uniform(0,1)$, then

${\mathcal{\xi
J}}(\boldsymbol{U}^{(m)}_{MRSSU})=-\frac{1}{2}\prod\limits_{i=1}^{m}\frac{1}{2i+1}=-\frac{1}{2}\left(\frac{\sqrt{\pi}}{2^{m}\Gamma(m+\frac{1}{2})}\right)<{\mathcal{\xi
J}}(\boldsymbol{U}^{(m)}_{SRS})=-0.5\left(\frac{1}{3}\right)^{m}.$ (2.5)

######  Example 2.2.

If $Z$ is exponentially distributed with mean $\frac{1}{\lambda}$. Then, we
have

${\mathcal{\xi
J}}(\boldsymbol{Z}^{(m)}_{MRSSU})=-\frac{1}{2}\prod\limits_{i=1}^{m}\frac{1}{2i\lambda}<{\mathcal{\xi
J}}(\boldsymbol{Z}^{(m)}_{SRS})=-\frac{1}{2}\left(\frac{1}{2\lambda}\right)^{m}.$
(2.6)

######  Example 2.3.

Let $X$ is finite range distribution with
$\bar{F}(x)=(1-ax)^{b},\;0<x<\frac{1}{a},\;a>0,\;b>0$. Then, we have

${\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MRSSU})=-\frac{1}{2}\prod\limits_{i=1}^{m}\frac{1}{a(1+2ib)}<{\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{SRS})=-\frac{1}{2}\left(\frac{1}{a(1+2b)}\right)^{m}.$
(2.7)

######  Theorem 2.1.

Let $\boldsymbol{X}^{(m)}_{MinRSSU}$ be the MinRSSU from population $X$ with
pdf $f$ and cdf $F$. Then, ${\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})\leq{\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{SRS})$ for $m>1$.

###### Proof.

Since $\bar{F}^{2}(x)\geq\bar{F}^{2i}(x)$ for $i\geq 1$, we have

$\left(\int_{0}^{+\infty}\bar{F}^{2}(x)dx\right)^{m}\leq\prod_{i=1}^{m}\int_{0}^{+\infty}\bar{F}^{2i}(x)dx.$

The proof follows by recalling (2) and (2.4). ∎

######  Remark 2.1.

If $f(F^{-1}(u))\geq 1$, $0<u<1$, then ${\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})$ is increasing in $m\geq 1$.

###### Proof.

From (2), we get

$\displaystyle\frac{{\mathcal{\xi
J}}(\boldsymbol{X}^{(m+1)}_{MinRSSU})}{{\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})}=\int_{0}^{1}\frac{(1-u)^{2m+2}}{f(F^{-1}(u))}du\leq\frac{1}{2m+3}\leq
1.$

The result follows readily, since the extropy is negative . ∎

In the following, we provide some results on the cumulative residual extropy
of $\boldsymbol{X}^{(m)}_{MinRSSU}$ in terms of stochastic ordering
properties. Now, we state important properties of ${\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})$ using the stochastic ordering. For that
we present the following definitions:

######  Definition 2.2.

(Shaked and Shanthikumar, 2007) Let $X$ and $Y$ be two non-negative random
variables with pdfs $f$ and $g$, cdfs $F$ and $G$, and hazard functions
$\lambda_{X}(x)=\frac{f(x)}{\bar{F}(x)}$ and
$\lambda_{Y}(y)=\frac{g(y)}{\bar{G}(y)}$, respectively. Then
1\. $X$ is said to be smaller than $Y$ in the usual stochastic order (denoted
by $X\leq_{st}Y$) if $P(X\geq x)\leq P(Y\geq x)$ for all $x\in\mathbb{R}$.
2\. $X$ is smaller than $Y$ in the hazard rate order (denoted by
$X\leq_{hr}Y$) if $\lambda_{X}(x)\geq\lambda_{Y}(x)$ for all $x$.
3\. $X$ is smaller than $Y$ in the dispersive order (denoted by
$X\leq_{disp}Y$) if $f(F^{-1}(u))\geq g(G^{-1}(u))$ for all $u\in(0,1)$, where
$F^{-1}$ and $G^{-1}$ are right continuous inverses of $F$ and $G$,
respectively.
4\. $X$ is said to have decreasing failure rate (DFR) if $\lambda_{X}(x)$ is
decreasing in $x$.
5\. $X$ is smaller than $Y$ in the convex transform order (denoted by
$X\leq_{c}Y$) if $G^{-1}F(x)$ is a convex function on the support of $X$.
6\. $X$ is smaller than $Y$ in the star order (denoted by $X\leq_{*}Y$) if
$\frac{G^{-1}F(x)}{x}$ is increasing in $x\geq 0$.
7\. $X$ is smaller than $Y$ in the superadditive order (denoted by
$X\leq_{su}Y$) if $G^{-1}(F(t+u))\geq G^{-1}(F(t))+G^{-1}(F(u))$ for $t\geq
0,u\geq 0$. 8\. $X$ is said to have an increasing reversed hazard rate (IRHR)
if $\tilde{\lambda}_{X}(x)=\frac{f(x)}{F(x)}$ is increasing in $x$.

######  Theorem 2.3.

If $X\leq_{st}Y$, then ${\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})\geq{\mathcal{\xi
J}}(\boldsymbol{Y}^{(m)}_{MinRSSU}),\;m>1.$

###### Proof.

By the assumption of the stochastic order,
$\bar{F}^{2i}(x)\leq\bar{G}^{2i}(x)$ for all $x\geq 0$. Now using (2), for
$m>1$, we get the desired result. ∎

######  Theorem 2.4.

Let $X$ and $Y$ be two non-negative random variable. If $X\leq_{disp}Y$, then
${\mathcal{\xi J}}(\boldsymbol{X}^{(m)}_{MinRSSU})\geq{\mathcal{\xi
J}}(\boldsymbol{Y}^{(m)}_{MinRSSU})$ for $m>1$.

###### Proof.

By the assumption of the dispersive order , $f(F^{-1}(u))\geq g(G^{-1}(u))$
for all $u\in(0,1)$. Using (2), for $m>1$ the result follows. ∎

######  Theorem 2.5.

If $X\leq_{hr}Y$, and $X$ or $Y$ is DFR, then ${\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})\geq{\mathcal{\xi
J}}(\boldsymbol{Y}^{(m)}_{MinRSSU})$ for $m>1$.

###### Proof.

If $X\leq_{hr}Y$, and $X$ or $Y$ is DFR, then $X\leq_{disp}Y$, due to Bagai &
Kochar [1986]. Thus, from Theorem (2.4) the desired result follows. ∎

######  Theorem 2.6.

Let $X$ and $Y$ be two non-negative random variable with pdf’s $f$ and $g$,
respectively, such that $f(0)\geq g(0)>0$ . If
$X\leq_{su}Y(X\leq_{*}Y\;or\;X\leq_{c}Y)$ , then ${\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})\geq{\mathcal{\xi
J}}(\boldsymbol{Y}^{(m)}_{MinRSSU})$ for $m>1$.

###### Proof.

If $X\leq_{su}Y(X\leq_{*}Y\;or\;X\leq_{c}Y)$, then $X\leq_{disp}Y$, due to
Ahmed et al. [1986]. So, from Theorem (2.4) the desired result follows. ∎

######  Proposition 2.7.

Let $\boldsymbol{X}^{(m)}_{MinRSSU}$ and $\boldsymbol{X}^{(m)}_{RSS}$ be
MinRSSU and RSS data from distribution $X$ with DFR ageing property,
respectively. Then for $m>1$ we have

${\mathcal{\xi J}}(\boldsymbol{X}^{(m)}_{MinRSSU})\geq{\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{SRS}).$

######  Remark 2.2.

Let $\varphi$ be a non-negative function such that the derivative
$\varphi^{\prime}(x)\geq 1$ for all $x$. Then $X\leq_{disp}\varphi(X)$. Thus,

${\mathcal{\xi J}}(\boldsymbol{X}^{(m)}_{MinRSSU})\geq{\mathcal{\xi
J}}(\varphi(\boldsymbol{X})^{(m)}_{MinRSSU}).$

If $X\leq_{disp}Y$, Theorem 3.B.26 in Shaked & Shanthikumar [2007] claims that
$X_{i:m}\leq_{disp}Y_{i:m},\;\;i=1,2,...,m$. Thus, according to Theorems 4.5
and 4.6 in Qiu [2017], we obtain the following Proposition .

######  Proposition 2.8.

Let $\boldsymbol{X}^{(m)}_{MinRSSU}$ be a sample from MinRSSU design.

(i) If $X$ is DFR ageing property, then ${\mathcal{\xi J}}(X_{(1)m})$ is
increasing in $m\geq 1$.

(ii) If $X$ is IRHR property, then ${\mathcal{\xi J}}(X_{(1)m})$ is decreasing
in $m\geq 1$.

######  Proposition 2.9.

Let $\boldsymbol{Y}^{(m)}_{MinRSSU}=a\boldsymbol{X}^{(m)}_{MinRSSU}+b$ with
$a>0$ and $b\geq 0$. Then, ${\mathcal{\xi
J}}(\boldsymbol{Y}^{(m)}_{MinRSSU})=a{\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})$.

######  Proposition 2.10.

Let $X$ be a symmetric random variable with respect to the finite mean
$\mu=\mathbb{E}(X)$. Then

$\displaystyle{\mathcal{\xi
J}}(\boldsymbol{X}^{(m)}_{MinRSSU})={\mathcal{CJ}}(\boldsymbol{X}^{(m)}_{MinRSSU}),$

where ${\mathcal{CJ}}(X)=-\frac{1}{2}\int_{0}^{+\infty}[F_{X}(x)]^{2}dx$ is
the cumulative extropy (see Jahanshahi et al. [2019]).

Let $X$ be the random lifetime of a system, recall that $X_{[t]}=[X-t\mid
X\leq t]$ describes the residual lifetime of a system. For all $t\geq 0$ the
mean residual lifetime is given by

$\displaystyle\mu(t)=E[X-t\mid X\geq
t]=\frac{1}{\bar{F}(t)}\int_{t}^{+\infty}\bar{F}(x)dx.$

Now, we can define a generalized measure of cumulative residual extropy as

$\displaystyle\mathcal{\xi
J}(X;t)=-\frac{1}{2}\int_{t}^{+\infty}\left[\frac{\bar{F}(x)}{\bar{F}(t)}\right]^{2}dx.$
(2.8)

Note that $\mathcal{\xi J}(X;t)\geq-\mu(t)/(2\bar{F}(t))$. Moreover , we have

${\mathcal{\xi J}}(\boldsymbol{X}^{(m)}_{SRS};t)=-\frac{1}{2}[-2{\mathcal{\xi
J}}(X,t)]^{m}.$ (2.9)

Under the MinRSSU design, it is clear that

$\displaystyle{\mathcal{\xi J}}(\boldsymbol{X}^{(m)}_{MinRSSU};t)$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\prod\limits_{i=1}^{m}[-2{\mathcal{\xi
J}}(X_{(i:i)};t)]=-\frac{1}{2}\prod\limits_{i=1}^{m}\int_{t}^{+\infty}\left[\frac{\bar{F}(x)}{\bar{F}(t)}\right]^{2i}dx$
(2.10) $\displaystyle=$
$\displaystyle-\frac{1}{2}\prod\limits_{i=1}^{m}\mathbb{E}\left[\frac{U^{2i}\bar{F}(t)}{f(F^{-1}(1-U\bar{F}(t)))}\right].$

######  Theorem 2.11.

Let $X$ be a random lifetime variable with cdf $F(\cdot)$, Then, for $m>1$

$\displaystyle\mathcal{\xi
J}(\boldsymbol{X}^{(m)}_{MinRSSU};t)\leq\mathcal{\xi
J}(\boldsymbol{X}^{(m)}_{SRS};t).$ (2.11)

###### Proof.

The proof is similar to Theorem 2.1. ∎

######  Remark 2.3.

If $f(F^{-1}(1-u\bar{F}(t)))\geq 1$, $0<u<1$, then $\mathcal{\xi
J}(\boldsymbol{X}^{(m)}_{MinRSSU};t)$ is increasing in $m\geq 1$.

##  3 Results on empirical measure of CREX

This section focuses on the estimation of the ${\mathcal{\xi J}}(X)$ based on
the SRS and MinRSSU schemes. Let

$X_{(1)}\leq X_{(2)}\leq...\leq X_{(n)}$ be the order statistics of the random
sample $X_{1},X_{2},...,X_{n}$ from cdf $F.$ Then the empirical measure of $F$
is defined as

$\hat{F}_{n}(x)=\left\\{\begin{array}[]{ll}0,&\ x<X_{(1)},\\\
\frac{k}{n},&X_{(k)}\leq x\leq X_{(k+1)},\ \ \ \ k=1,2,...,n-1\\\ 1,&\
x>X_{(n)}.\end{array}\right.$

Thus, the empirical measure of ${\mathcal{\xi J}}(X)$ is obtained by replacing
the distribution function F by the empirical distribution function
$\hat{F}_{n}$ as

$\displaystyle V_{n}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int\hat{\bar{F}}_{n}^{2}(x)dx=-\frac{1}{2}\sum_{k=1}^{n-1}\int_{X_{(k)}}^{X_{(k+1)}}\left(1-\frac{k}{n}\right)^{2}dx$
(3.1) $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{k=1}^{n-1}U_{k+1}\left(1-\frac{k}{n}\right)^{2},$

where $U_{k+1}=X_{(k+1)}-X_{(k)},k=1,...,n-1$. Jahanshahi et al. [2019] showed
that $V_{n}$ almost surely converges to the CREX of $X$, i.e.

$V_{n}\overset{a.s.}{\rightarrow}{\mathcal{\xi J}}(X),\;\;\text{as
}\;n\rightarrow+\infty.$

The problem of estimation of ${\mathcal{\xi J}}(X)$ based on MinSSU scheme can
be deduced in the same line of the estimation based on SRS design $V_{n}$. In
this part, we assume that instead of one-cycle MinRSSU, the process is
repeated $l$ cycles to have a sample of size $n=ml$. In this case, the
resulting MinRSSU is denoted by $\left\\{X_{\left(1:i\right)j},\text{
}i=1,...,m;\text{ }j=1,...,l\right\\}$, where $X_{\left(1:i\right)j}$ is the
lowest order statistic from the $i$th sample in the $j$th cycle. Let
$Y_{\left(1\right)},...,Y_{\left(n\right)}$ be the ordered values of the
MinRSSU design $\left\\{X_{\left(1:i\right)j},\text{ }i=1,...,m;\text{
}j=1,...,l\right\\}.$ Then, the natural estimation of ${\mathcal{\xi J}}(X)$
based on the MinRSSU, can be obtained as

$R_{n}=-\frac{1}{2}\sum_{k=1}^{n-1}Z_{k+1}\left(1-\frac{k}{n}\right)^{2},$

where $Z_{k+1}=Y_{(k+1)}-Y_{(k)},k=1,2,...,n-1$. Our preliminary computations
and simulations showed that $R_{n}$ has some deficiencies to be unbiased and
have low mean square error (MSE) to estimate the ${\mathcal{\xi J}}(X).$ We
observed that the term $\left(1-\frac{k}{n}\right)$ should be slightly
modified so that the estimator has optimal properties to estimate
${\mathcal{\xi J}}(X).$ We propose to modify this term with
$\left(1-\frac{k}{n+m+w}\right),$ where $w$ is a number that resulting
estimator has optimally low bias and MSE. The resulting estimator of the
${\mathcal{\xi J}}(X)$ has the following form

$R_{m,n}=-\frac{1}{2}\sum_{k=1}^{n-1}Z_{k+1}\left(1-\frac{k}{n+m+w}\right)^{2}.$
(3.2)

###  3.1 A new estimation

In the previous section, the estimator $V_{n}$ of ${\mathcal{\xi J}}(X)$ is a
linear function of sample spacing $U_{k+1}=X_{(k+1)}-X_{(k)},k=1,...,n-1$. The
asymptotic distribution of this linear function of sample spacing can be
found, for example in Di Crescenzo & Longobardi [2009] and Tahmasebi [2019]
only for exponential and standard uniform distributions. So we provide another
estimator for ${\mathcal{\xi J}}(X)$ which is a linear function of order
statistics.

######  Proposition 3.1.

Let $X$ be an absolutely continuous non-negative random variable with survival
function $\bar{F}$, then

${\mathcal{\xi J}}(X)=-\int^{+\infty}_{0}x\bar{F}(x)d{F}(x).$ (3.3)

###### Proof.

By (2.1) and Fubini’s theorem, we obtain

$\displaystyle-\int^{+\infty}_{0}x\bar{F}(x)d{F}(x)$ $\displaystyle=$
$\displaystyle-\int^{+\infty}_{0}\left(\int^{x}_{0}dt\right)\bar{F}(x)f(x)dx$
$\displaystyle=$
$\displaystyle-\int^{+\infty}_{0}\left(\int^{+\infty}_{t}\bar{F}(x)f(x)dx\right)dt=-\frac{1}{2}\int_{0}^{+\infty}\bar{F}^{2}(t)dt.$

Hence, the proof is completed ∎

The new estimator can be obtained replacing $\bar{F}(x)$ with
$\hat{\bar{F}}_{n}(x)$ in (3.3), so $\widehat{\mathcal{\xi J}}(X)$ has the
following form

$\widehat{\mathcal{\xi
J}}(X)=-\int_{0}^{+\infty}x\hat{\bar{F}}_{n}\left(x\right)d\hat{F}_{n}\left(x\right)=-\frac{1}{n}\sum_{i=1}^{n}\left(1-\frac{i}{n}\right)X_{\left(i\right)}.$
(3.4)

Let $J\left(x\right)=\left(1-x\right).$ Then (3.4) has the form

$\widehat{{\mathcal{\xi
J}}}(X)=-\frac{1}{n}\sum_{i=1}^{n}J\left(\frac{i}{n}\right)X_{\left(i\right)},$

which is a linear function of order statistics. The natural estimation of
${\mathcal{\xi J}}(X)$ based on the MinRSSU, can be obtained as

$\widehat{{\mathcal{\xi J}}}(Y):=\widehat{{\mathcal{\xi
J}}}(X_{MinRSSU})=-\frac{1}{n}\sum_{i=1}^{n}J\left(\frac{i}{n}\right)Y_{\left(i\right)}.$

Stigler [1974] showed that asymptotic distribution of such a linear
combination is normal distribution. The results of Stigler [1974] also hold if
the independent observations are not identically distributed. These properties
help us to obtain the asymptotic distribution of $\widehat{{\mathcal{\xi
J}}}(X)$ for both SRS (observations are independent and identical) and MinRSSU
(observations are only independent) designs.

######  Theorem 3.2.

Assume that $E\left(X^{2}\right)<+\infty$. Then

$\displaystyle\sqrt{n}\left(\widehat{{\mathcal{\xi J}}}(X)-{\mathcal{\xi
J}}(X)\right)\overset{d}{\rightarrow}N\left(0,\sigma^{2}\left(J,F\right)\right),$
$\displaystyle\sqrt{n}\left(\widehat{{\mathcal{\xi J}}}(Y)-{\mathcal{\xi
J}}(X)\right)\overset{d}{\rightarrow}N\left(0,\sigma_{MinRSSU}^{2}(J,\tilde{F},K)\right),$

where

$\sigma^{2}\left(J,F\right)=\int_{0}^{+\infty}\int_{0}^{+\infty}J(F(x))J(F(y))[F(min(x,y))-F(x)F(y)]dxdy,$
$\sigma_{MinRSSU}^{2}(J,\tilde{F},K)=\int_{0}^{+\infty}\int_{0}^{+\infty}J(\tilde{F}(x))J(\tilde{F}(y))K(x,y)dxdy,$

and $\tilde{F}(.)$ and $K(x,y)$ are given as

$\displaystyle\tilde{F}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{i=1}^{m}F_{(1)i}(x),$ $\displaystyle K(x,y)$
$\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{i=1}^{m}[F_{(1)i}(min(x,y))-F_{(1)i}(x)F_{(1)i}(y)].$

###### Proof.

For $\tilde{F}(x)$ and $K(x,y)$, we have

$\displaystyle\tilde{F}(x)$ $\displaystyle=$
$\displaystyle\lim_{n\to+\infty}\frac{l}{n}\sum_{i=1}^{m}F_{(1)i}(x)=\frac{1}{m}\sum_{i=1}^{m}F_{(1)i}(x),$
$\displaystyle K(x,y)$ $\displaystyle=$
$\displaystyle\lim_{n\to+\infty}\frac{l}{n}\sum_{i=1}^{m}[F_{(1)i}(min(x,y))-F_{(1)i}(x)F_{(1)i}(y)]$
$\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{i=1}^{m}[F_{(1)i}(min(x,y))-F_{(1)i}(x)F_{(1)i}(y)],$

where as before $l$ is the size of the cycle of MinRSSU design with $n=ml$.
The rest of the proof is done by using the results of Stigler [1974] for both
SRS and MinRSSU designs. ∎

As we formerly stated for providing the estimator $R_{m,n}$, some adjusted
forms of estimator $\widehat{{\mathcal{\xi J}}}(Y)$ can be used since the
estimator $\widehat{{\mathcal{\xi J}}}(Y)$ has some deficiencies to estimate
CREX which need to be fixed. We again observed that the the term
$\left(1-\frac{i}{n}\right)$ needs to be slightly adjusted so that the
resulting estimator has optimal properties to estimate CREX. This new term has
the following form

$1-\frac{i}{n+\psi(m,w)},$ (3.5)

where the function $\psi(.,.)$ is a challenging factor that can be
specifically determined for each given distribution and consequently the
resulting estimator has optimal low bias and MSE. In this case, the resulting
estimator has the following form

$\widehat{{\mathcal{\xi
J}}}_{m,n}(Y)=-\frac{1}{n}\sum_{i=1}^{n}J\left(\frac{i}{n+\psi(m,w)}\right)Y_{\left(i\right)}.$
(3.6)

In the next section, we determine the form of the function $\psi(.,.)$ for
exponential, uniform and beta distributions and we show that the optimal
choice of this function can reduce the bias and MSE in the estimate of CREX.

###  3.2 Simulation study

In advance, we explain the role of the parameter $w$ and function $\psi(m,w)$
for which $R_{m,n}$ and $\widehat{{\mathcal{\xi J}}}_{m,n}(X_{MinRSSU})$ have
optimally low bias and MSE. For this purpose, we examine some distributions to
obtain the function $\psi(m,w)$ and optimal value of $w$ in $R_{m,n}$ and
$\widehat{{\mathcal{\xi J}}}_{m,n}(Y)$. In Tables 1-3, for four estimators
$R_{n}$, $R_{m,n}$, $\widehat{{\mathcal{\xi J}}}(Y)$ and
$\widehat{{\mathcal{\xi J}}}_{m,n}(Y)$, we compute the bias and MSE to
estimate the parameter ${\mathcal{\xi J}}(X)$ for some different values of
$w$. Here, the exponential ($Exp\left(\lambda\right)$), uniform
($Unif\left(0,b\right)$) and beta ($Beta\left(\alpha,1\right),\alpha>1$)
distributions are considered. For each configuration, the simulation study was
carried out with 5000 repetitions. The number of cycle and size of the sample
in each cycle are taken as $l=2,3$ and $m=2,...,5$, respectively. We compute
the bias and root of MSE (RMSE) of each estimator of parameter ${\mathcal{\xi
J}}(X)$. In Tables 1-3, it can be seen that results of biases and RMSEs of
estimator $\widehat{{\mathcal{\xi J}}}(Y)$ are not comparable to those of
$\widehat{{\mathcal{\xi J}}}_{m,n}(Y)$ for different values of $w$. We
intuitively obtain the form of the function $\psi(m,w)$ for different
distributions. We found that the function $\psi$ has the form $5m-4k_{m}+w$
with $k_{2}=3$,…, $k_{5}=0$, $3m-(2k_{m}+1)+w$ with $k_{2}=-1$,…, $k_{5}=2$
and $m-w$ for exponential, uniform and beta distributions, respectively. It is
observed that changing the value of $l$ has no effect on the whole results
verified from biases and RMSEs. In all tables, we see that choosing the proper
function $\psi(m,w)$ and parameter $w$ as the challenging factors can improve
the efficiency of estimators $R_{m,n}$ and $\widehat{{\mathcal{\xi
J}}}_{m,n}(Y)$ against $R_{n}$ and $\widehat{{\mathcal{\xi J}}}(Y)$,
respectively, in estimating the parameter $\mathcal{\xi J}$.

Table 1: The biases and MSEs of the different estimators: Exponential
distribution

|  |  Results based on $R_{m,n}$ |  |  |  |  Results based on $\widehat{\mathcal{\xi J}}_{m,n}(Y)$
---|---|---|---|---|---|---
|  | $l=2$ |  |  | $l=3$ |  |  |  |  | $l=2$ |  |  | $l=3$ |
$m$ | $w$ | Bias | RMSE |  | Bias | RMSE |  | $m$ | $w$ | Bias | RMSE |  | Bias | RMSE
2 | -2 | 0.321 | 0.407 |  | 0.370 | 0.477 |  | 2 | -11 | 0.632 | 0.650 |  | 0.536 | 0.562
| -1 | 0.131 | 0.354 |  | 0.133 | 0.460 |  |  | -10 | 0.402 | 0.483 |  | 0.354 | 0.427
| 0 | -0.051 | 0.423 |  | -0.069 | 0.580 |  |  | -9 | 0.263 | 0.435 |  | 0.223 | 0.370
| 1 | -0.216 | 0.556 |  | -0.234 | 0.733 |  |  | -8 | 0.171 | 0.434 |  | 0.125 | 0.361
| $R_{n}$ | 0.321 | 0.407 |  | 0.370 | 0.477 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.402 | 0.483 |  | 0.354 | 0.427
3 | -1 | 0.124 | 0.367 |  | 0.164 | 0.304 |  | 3 | -7 | 0.033 | 0.374 |  | 0.006 | 0.307
| 0 | -0.017 | 0.429 |  | 0.040 | 0.315 |  |  | -6 | 0.012 | 0.382 |  | -0.020 | 0.317
| 1 | 0.124 | 0.367 |  | 0.164 | 0.304 |  |  | -5 | -0.006 | 0.391 |  | -0.044 | 0.328
| 2 | -0.253 | 0.632 |  | -0.191 | 0.472 |  |  | -4 | -0.022 | 0.399 |  | -0.065 | 0.340
| $R_{n}$ | 0.437 | 0.478 |  | 0.411 | 0.444 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.452 | 0.489 |  | 0.430 | 0.458
4 | 0 | 0.024 | 0.383 |  | 0.118 | 0.291 |  | 4 | -3 | 0.029 | 0.329 |  | -0.007 | 0.269
| 1 | -0.077 | 0.452 |  | 0.029 | 0.313 |  |  | -2 | 0.020 | 0.332 |  | -0.019 | 0.274
| 2 | -0.170 | 0.532 |  | -0.057 | 0.363 |  |  | -1 | 0.012 | 0.334 |  | -0.030 | 0.278
| 3 | -0.256 | 0.616 |  | -0.139 | 0.427 |  |  | 0 | 0.004 | 0.337 |  | -0.040 | 0.283
| $R_{n}$ | 0.486 | 0.509 |  | 0.474 | 0.490 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.509 | 0.527 |  | 0.489 | 0.503
5 | 1 | -0.004 | 0.390 |  | 0.114 | 0.281 |  | 5 | 1 | 0.072 | 0.277 |  | 0.050 | 0.232
| 2 | -0.082 | 0.448 |  | 0.046 | 0.296 |  |  | 2 | 0.067 | 0.277 |  | 0.043 | 0.233
| 3 | -0.156 | 0.512 |  | -0.021 | 0.329 |  |  | 3 | 0.063 | 0.278 |  | 0.036 | 0.233
| 4 | -0.224 | 0.578 |  | -0.085 | 0.373 |  |  | 4 | 0.058 | 0.278 |  | 0.030 | 0.234
| $R_{n}$ | 0.533 | 0.545 |  | 0.524 | 0.532 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.554 | 0.563 |  | 0.540 | 0.547

Table 2: The biases and MSEs of the different estimators: Uniform distribution

|  |  Results based on $R_{m,n}$ |  |  |  |  Results based on $\widehat{\mathcal{\xi J}}_{m,n}(Y)$
---|---|---|---|---|---|---
|  | $l=2$ |  |  | $l=3$ |  |  |  |  | $l=2$ |  |  | $l=3$ |
$m$ | $w$ | Bias | RMSE |  | Bias | RMSE |  | $m$ | $w$ | Bias | RMSE |  | Bias | RMSE
2 | -2 | 0.360 | 0.392 |  | 0.298 | 0.333 |  | 2 | -4 | 0.098 | 0.275 |  | 0.058 | 0.227
| -1 | 0.249 | 0.305 |  | 0.212 | 0.264 |  |  | -3 | 0.045 | 0.281 |  | -0.001 | 0.236
| 0 | 0.158 | 0.253 |  | 0.134 | 0.214 |  |  | -2 | 0.005 | 0.293 |  | -0.049 | 0.254
| 1 | 0.085 | 0.233 |  | 0.066 | 0.186 |  |  | -1 | -0.026 | 0.306 |  | -0.088 | 0.275
| $R_{n}$ | 0.360 | 0.392 |  | 0.298 | 0.333 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.284 | 0.340 |  | 0.238 | 0.293
3 | -1 | 0.191 | 0.243 |  | 0.201 | 0.237 |  | 3 | -2 | -0.006 | 0.247 |  | -0.023 | 0.200
| 0 | 0.125 | 0.204 |  | 0.147 | 0.196 |  |  | -1 | -0.026 | 0.255 |  | -0.047 | 0.209
| 1 | 0.067 | 0.183 |  | 0.096 | 0.165 |  |  | 0 | -0.043 | 0.263 |  | -0.068 | 0.220
| 2 | 0.016 | 0.181 |  | 0.055 | 0.148 |  |  | 1 | -0.058 | 0.270 |  | -0.087 | 0.231
| $R_{n}$ | 0.349 | 0.372 |  | 0.318 | 0.339 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.318 | 0.350 |  | 0.290 | 0.317
4 | 0 | 0.14 | 0.195 |  | 0.175 | 0.204 |  | 4 | 0 | -0.012 | 0.205 |  | -0.008 | 0.164
| 1 | 0.092 | 0.171 |  | 0.135 | 0.173 |  |  | 1 | -0.023 | 0.209 |  | -0.021 | 0.168
| 2 | 0.049 | 0.159 |  | 0.098 | 0.148 |  |  | 2 | -0.034 | 0.213 |  | -0.034 | 0.172
| 3 | 0.018 | 0.159 |  | 0.062 | 0.131 |  |  | 3 | -0.043 | 0.217 |  | -0.046 | 0.177
| $R_{n}$ | 0.372 | 0.387 |  | 0.349 | 0.360 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.352 | 0.369 |  | 0.337 | 0.350
5 | 1 | 0.116 | 0.17 |  | 0.175 | 0.198 |  | 5 | 2 | 0.020 | 0.181 |  | 0.028 | 0.146
| 2 | 0.079 | 0.153 |  | 0.144 | 0.172 |  |  | 3 | 0.013 | 0.183 |  | 0.019 | 0.146
| 3 | 0.044 | 0.144 |  | 0.114 | 0.150 |  |  | 4 | 0.006 | 0.184 |  | 0.010 | 0.147
| 4 | 0.012 | 0.144 |  | 0.085 | 0.132 |  |  | 5 | -0.001 | 0.186 |  | 0.002 | 0.148
| $R_{n}$ | 0.392 | 0.401 |  | 0.381 | 0.388 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.387 | 0.398 |  | 0.375 | 0.383

Table 3: The biases and MSEs of the different estimators: Beta distribution

|  |  Results based on $R_{m,n}$ |  |  |  |  Results based on $\widehat{\mathcal{\xi J}}_{m,n}(Y)$
---|---|---|---|---|---|---
|  | $l=2$ |  |  | $l=3$ |  |  |  |  | $l=2$ |  |  | $l=3$ |
$m$ | $w$ | Bias | RMSE |  | Bias | RMSE |  | $m$ | $w$ | Bias | RMSE |  | Bias | RMSE
2 | -2 | 0.296 | 0.298 |  | 0.276 | 0.279 |  | 2 | -3 | 0.202 | 0.204 |  | 0.146 | 0.149
| -1 | 0.280 | 0.283 |  | 0.264 | 0.267 |  |  | -2 | 0.108 | 0.116 |  | 0.081 | 0.089
| 0 | 0.267 | 0.271 |  | 0.253 | 0.257 |  |  | -1 | 0.053 | 0.070 |  | 0.035 | 0.055
| 1 | 0.257 | 0.261 |  | 0.244 | 0.248 |  |  | 0 | 0.015 | 0.053 |  | 0.000 | 0.045
| $R_{n}$ | 0.296 | 0.298 |  | 0.276 | 0.279 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.108 | 0.116 |  | 0.081 | 0.089
3 | -1 | 0.251 | 0.255 |  | 0.238 | 0.242 |  | 3 | -3 | 0.096 | 0.102 |  | 0.079 | 0.085
| 0 | 0.241 | 0.245 |  | 0.230 | 0.234 |  |  | -2 | 0.052 | 0.064 |  | 0.045 | 0.056
| 1 | 0.233 | 0.238 |  | 0.222 | 0.227 |  |  | -1 | 0.019 | 0.045 |  | 0.017 | 0.039
| 2 | 0.225 | 0.231 |  | 0.216 | 0.221 |  |  | 0 | -0.007 | 0.043 |  | -0.006 | 0.037
| $R_{n}$ | 0.275 | 0.278 |  | 0.256 | 0.259 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.096 | 0.102 |  | 0.079 | 0.085
4 | 0 | 0.227 | 0.231 |  | 0.218 | 0.222 |  | 4 | -3 | 0.059 | 0.068 |  | 0.056 | 0.063
| 1 | 0.220 | 0.225 |  | 0.212 | 0.216 |  |  | -2 | 0.031 | 0.047 |  | 0.033 | 0.044
| 2 | 0.214 | 0.219 |  | 0.206 | 0.211 |  |  | -1 | 0.007 | 0.038 |  | 0.013 | 0.033
| 3 | 0.208 | 0.213 |  | 0.201 | 0.206 |  |  | 0 | -0.012 | 0.041 |  | -0.004 | 0.031
| $R_{n}$ | 0.264 | 0.267 |  | 0.246 | 0.249 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.095 | 0.100 |  | 0.082 | 0.087
5 | 1 | 0.217 | 0.221 |  | 0.210 | 0.214 |  | 5 | -3 | 0.043 | 0.053 |  | 0.047 | 0.054
| 2 | 0.211 | 0.216 |  | 0.205 | 0.209 |  |  | -2 | 0.022 | 0.039 |  | 0.030 | 0.040
| 3 | 0.206 | 0.210 |  | 0.200 | 0.204 |  |  | -1 | 0.004 | 0.034 |  | 0.015 | 0.031
| 4 | 0.201 | 0.205 |  | 0.196 | 0.201 |  |  | 0 | -0.012 | 0.036 |  | 0.002 | 0.028
| $R_{n}$ | 0.257 | 0.259 |  | 0.239 | 0.242 |  |  | $\widehat{{\mathcal{\xi J}}}(Y)$ | 0.097 | 0.101 |  | 0.088 | 0.091

##  4 Discrimination information

This section considers a new discrimination measure of disparity between the
distribution of MinRSSU and parental data SRS. Raqab & Qiu [2019] defined the
discrimination information between the density function of the $i$th order
statistic $f_{(i)m}$ and the underlying density function $f$ as

$D_{m}\left(f_{(i)m}:f\right)=\frac{1}{2}\int_{-\infty}^{+\infty}f_{(i)m}\left(x\right)\left(f_{(i)m}\left(x\right)-f\left(x\right)\right)dx.$
(4.1)

Analogously to (4.1), we define the discrimination information between the
survival function of the smallest ordered statistic $\bar{F}_{(1)i}$ and the
underlying survival function $\bar{F}$ as

$D\left(\bar{F}_{(1)i}:\bar{F}\right)=-\frac{1}{2}\int_{-\infty}^{+\infty}\bar{F}_{(1)i}\left(x\right)\left(\bar{F}_{(1)i}\left(x\right)-\bar{F}\left(x\right)\right)dx.$
(4.2)

The discrimination information in (4.2) may be rewritten in a simpler way. It
can be shown that

$D\left(\bar{F}_{(1)i}:\bar{F}\right)=-\frac{1}{2}\left[E(X_{(1)2i})-E(X_{(1)i+1})\right],$
(4.3)

where $X_{(1)j}$ is the smallest order statistic in a random sample of size
$j$.

######  Example 4.1.

Let $U\sim Uniform(0,1)$. We know that the order statistics from the standard
uniform distribution follow the beta distribution. Then discrimination
information based on (4.2) is

$D\left(\bar{F}_{(1)i}:\bar{F}\right)=-\frac{1}{2}\left[\frac{1}{2i+1}-\frac{1}{i+2}\right]=\frac{i-1}{2\left(2i+1\right)\left(i+2\right)}.$

In the following theorem, we obtain the discrimination information $D$ between
MinRSSU and SRS designs.

######  Theorem 4.1.

For $\boldsymbol{X}_{MinRSSU}^{(m)}$ and $\boldsymbol{X}_{SRS}^{(m)}$, we have

$D\left(\boldsymbol{X}_{MinRSSU}^{(m)}:\boldsymbol{X}_{SRS}^{(m)}\right)=-\frac{1}{2}\left(\prod\limits_{i=1}^{m}E\left(X_{\left(1\right)2i}\right)-\prod\limits_{i=1}^{m}E\left(X_{\left(1\right)i+1}\right)\right).$

###### Proof.

From (2.1) and (4.2) we have

$\displaystyle
D\left(\boldsymbol{X}_{MinRSSU}^{(m)}:\boldsymbol{X}_{SRS}^{(m)}\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\left(\prod\limits_{i=1}^{m}\int_{-\infty}^{+\infty}\bar{F}^{2i}\left(x_{1}\right)-\prod\limits_{i=1}^{m}\int_{-\infty}^{+\infty}\bar{F}^{i+1}\left(x_{1}\right)\right)$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\left(\prod\limits_{i=1}^{m}E\left(X_{\left(1\right)2i}\right)-\prod\limits_{i=1}^{m}E\left(X_{\left(1\right)i+1}\right)\right).$

The proof is completed. ∎

######  Example 4.2.

Let $U\sim Uniform(0,1)$. Using the results of Example (4.1), the
discrimination information $D$ for the MinRSSU and SRS designs of the same
size $m$ is

$D\left(\boldsymbol{U}_{MinRSSU}^{(m)}:\boldsymbol{U}_{SRS}^{(m)}\right)=-\frac{1}{2}\left(\prod\limits_{i=1}^{m}\frac{1}{2i+1}-\prod\limits_{i=1}^{m}\frac{1}{i+2}\right).$

##  5 Conclusion

This paper has introduced the uncertainty measure of the cumulative residual
extropy based on the MinRSSU and SRS data. Several results of the CREX measure
including stochastic orders were obtained for MinRSSU and SRS data. Also, we
provided two estimators of CREX measure for both SRS and MinRSSU data.
Furthermore, it was shown that MinRSSU scheme can efficiently reduce the
uncertainty measure of CREX. Also, by providing a discrimination measure, we
derived the distance size between MinRSSU and SRS data.

## Acknowledgements

C. Calì and M. Longobardi are partially supported by the GNAMPA research group
of INDAM (Istituto Nazionale di Alta Matematica) and MIUR-PRIN 2017, Project
”Stochastic Models for Complex Systems” (No. 2017JFFHSH).

## References

## References

  * Ahmed et al. [1986] Ahmed, A. N., Alzaid, A., Bartoszewicz, J., & Kochar, S. C. (1986).  Dispersive and superadditive ordering.  Advances in Applied Probability, 18, 1019–1022.
  * Bagai & Kochar [1986] Bagai, I., & Kochar, S. C. (1986).  On tail-ordering and comparison of failure rates.  Communications in Statistics - Theory and Methods, 15, 1377–1388.
  * Di Crescenzo & Longobardi [2009] Di Crescenzo, A., & Longobardi, M. (2009).  On cumulative entropies.  Journal of Statistical Planning and Inference, 139, 4072–4087.
  * Eskandarzadeh et al. [2018] Eskandarzadeh, M., Di Crescenzo, A., & Tahmasebi, S. (2018).  Measures of information for maximum ranked set sampling with unequal samples.  Communications in Statistics - Theory and Methods, 47, 4692–4709.
  * Eskandarzadeh et al. [2016] Eskandarzadeh, M., Tahmasebi, S., & Afshari, M. (2016).  Information measures for record ranked set samples.  Ciência e Natura, 38, 554–563.
  * Frank et al. [2015] Frank, L., Sanfilippo, G., & Agrò, G. (2015).  Extropy: Complementary dual of entropy.  Statistical Science, 30, 40–58.
  * Jahanshahi et al. [2019] Jahanshahi, S. M. A., Zarei, H., & Khammar, A. H. (2019).  On cumulative residual extropy.  Probability in the Engineering and Informational Sciences, 264, 1–21. doi:10.1017/S0269964819000196.
  * Jozani & Ahmadi [2014] Jozani, M. J., & Ahmadi, J. (2014).  On uncertainty and information properties of ranked set samples.  Information Sciences, 264, 291–301.
  * McIntyre [1952] McIntyre, G. A. (1952).  A method for unbiased selective sampling, using ranked sets.  Australian Journal of Agricultural Research, 3, 385–390.
  * Patil et al. [1999] Patil, G., Sinha, A., & Taillie, C. (1999).  Ranked set sampling: a bibliography.  Environmental and Ecological Statistics, 6, 91–98.
  * Qiu [2017] Qiu, G. (2017).  The extropy of order statistics and record values.  Statistics & Probability Letters, 120, 52–60.
  * Qiu & Eftekharian [2020] Qiu, G., & Eftekharian, A. (2020).  Extropy information of maximum and minimum ranked set sampling with unequal samples.  Communications in Statistics - Theory and Methods, . doi:10.1080/03610926.2019.1678640.
  * Raqab & Qiu [2019] Raqab, M., & Qiu, G. (2019).  On extropy properties of ranked set sampling.  Statistics, 53, 210–226.
  * Shaked & Shanthikumar [2007] Shaked, M., & Shanthikumar, J. G. (2007).  Stochastic Orders.  Springer, New York.
  * Stigler [1974] Stigler, S. (1974).  Linear functions of order statistics with smooth weight functions.  The Annals of Statistics, 2, 676–693.
  * Tahmasebi [2019] Tahmasebi, S. (2019).  Weighted extensions of generalized cumulative residual entropy and their applications.  Communications in Statistics - Theory and Methods, . doi:10.1080/03610926.2019.1615094.
  * Tahmasebi et al. [2016] Tahmasebi, S., Jafari, A. A., & Eskandarzadeh, M. (2016).  Some results on residual entropy of ranked set samples.  Statistics & Probability Letters, 112, 137–145.
  * Tahmasebi et al. [2020] Tahmasebi, S., Longobardi, M., Kazemi, M. R., & Alizadeh, M. (2020).  Cumulative tsallis entropy for maximum ranked set sampling with unequal samples.  Physica A. Statistical Mechanics and its Applications, . doi:10.1016/j.physa.2020.124763.


# Neural representation of a time optimal, constant acceleration rendezvous

Dario Izzo111Scientific coordinator, Advanced Concepts Team, ESA-ESTEC and
Sebastien Origer222Scientist, Advanced Concepts Team, ESA-ESTEC European Space
Technology and Research Centre, 2201 AZ Noordwijk, The Netherlands

###### Abstract

We train neural models to represent both the optimal policy (i.e. the optimal
thrust direction) and the value function (i.e. the time of flight) for a time
optimal, constant acceleration low-thrust rendezvous. In both cases we develop
and make use of the data augmentation technique we call backward generation of
optimal examples. We are thus able to produce and work with large dataset and
to fully exploit the benefit of employing a deep learning framework. We
achieve, in all cases, accuracies resulting in successful rendezvous
(simulated following the learned policy) and time of flight predictions (using
the learned value function). We find that residuals as small as a few m/s,
thus well within the possibility of a spacecraft navigation $\Delta V$ budget,
are achievable for the velocity at rendezvous. We also find that, on average,
the absolute error to predict the optimal time of flight to rendezvous from
any orbit in the asteroid belt to an Earth-like orbit is small (less than 4%)
and thus also of interest for practical uses, for example, during preliminary
mission design phases.

## Introduction

Deep neural models for optimal spacecraft guidance and control have been
recently proposed [1, 2, 3, 4, 5, 6] and studied as a possible alternative to
more classical architectures. Previously, researchers had proposed the use of
“neurocontrollers” [7] as part of the guidance and control system, but mostly
using a neural model to track some precomputed reference guidance signal. More
recent works succeeded in learning directly the solution to Bellman’s
equations [8] hence synthesising an optimal state feedback for the underlying
nonlinear equations of motion and thus allowing the use of networks for
guidance and control simultaneously. These new types of networks, which we
call Guidance and Control Networks (G&CNETs), have a modest computational cost
during inference, with clear applications for on-board, autonomous guidance
and control [2, 1, 9]. Training them, on the other hand, is mostly done via
supervised learning and requires the generation and use of a great number of
optimal trajectories to be used as labels to learn from. The generation of a
suitable dataset is key to the success in training any neural model to achieve
the precision requested by real mission design constraints. To this end, an
indirect method based on Pontryagin’s approach [10, 11] can be used to
generate (and later learn from) and label optimal decisions. As a consequence,
the solution to many Two Points Boundary Value Problems (TPBVPs) needs to be
computed. Even in the optimistic case of good co-states guesses being
available [12, 3] and homotopy methods employed, these numerical procedures
are time consuming and limit the amount of reference trajectories one can
learn from [6, 3]. A similar conclusion can be made also when direct methods
are used to generate optimal examples. The development of neural models
representing the optimal policy of some interplanetary transfer, such as
G&CNETs, is not the only case where the availability of a large amount of
optimal solutions is needed. Attempts to learn directly machine learning
models (mostly deep neural models) for the propellant consumption, the $\Delta
V$, the arrival mass or the time of flight have also been performed. These can
all be seen, from an optimal control point of view, as attempts to learn with
different degrees of approximation the value function of some underlying
fundamental optimal control problem. Complex low-thrust interplanetary
trajectories between objects of the main belt [13, 14], transfers to Earth co-
orbiting asteroids [15] as well as orbiting debris scenarios[16] have been
considered as applications in the past. A critic to these approaches for value
function learning has been recently worded as: “… one has to balance the
advantage of computation speed against a remarkable time for data sets
generating and network training.”[16].

In this work we utilise a constant acceleration interplanetary rendezvous
problem, to show the use of the data augmentation technique called Backward
Generation of Optimal Examples [2] as a methodology to reduce the time
required to generate datasets by orders of magnitude. The method allows to
generate optimal samples avoiding entirely the solution of TPBVPs by simply
integrating backward in time the augmented system of differential equations
derived from Pontryagin’s theory. We then also show how, contrary to a common
sentiment, accuracies compatible with real world requirements can be reached
both when learning the optimal policy and the value function.

The paper is structured as follows. We begin giving the necessary details and
formalism on the optimal control problem considered throughout the paper. In
doing so we also follow recent practical guidelines [11] for solving low-
thrust problems via the maximum principle, which reveal to have added specific
advantages for the augmentation technique. We then show how to apply the
backward generation of optimal examples to augment a single nominal optimal
trajectory found and discuss the results obtained training several G&CNETs for
different initial orbits. In the following section we move on to consider the
value function problem and discuss the results one can obtain when augmenting
each trajectory in a dataset aimed at representing the time of transfer from
orbits in the asteroid belt to a target rendezvous.

## Methods

### The optimal control problem

Consider a time optimal, constant acceleration rendezvous with a body in a
perfectly circular orbit. This problem was recently brought to the attention
of a wide international community during the 11th edition of the Global
Trajectory Optimisation Competition (GTOC11) where it represented low-thrust
transfers of objects from the belt. Let us introduce a rotating frame
$\mathcal{F}=\left[\hat{\mathbf{i}},\hat{\mathbf{j}},\hat{\mathbf{k}}\right]$
having angular velocity
$\boldsymbol{\omega}=\sqrt{\frac{\mu}{R^{3}}}\hat{\mathbf{k}}$ so that the
target body position is, in $\mathcal{F}$, stationary and indicated with
$R\hat{\mathbf{i}}$. The dynamics is thus described by the following simple
set of ordinary differential equations:

$\left\\{\begin{array}[]{l}\dot{\mathbf{r}}=\mathbf{v}\\\
\dot{\mathbf{v}}=-\frac{\mu}{r^{3}}\mathbf{r}-2\boldsymbol{\omega}\times\mathbf{v}-\boldsymbol{\omega}\times\boldsymbol{\omega}\times\mathbf{r}+\Gamma\mathbf{\hat{t}}(t)\end{array}\right.$
(1)

where $\mu$ is the gravitational constant of the Sun and $\Gamma$ a constant
acceleration assumed to be acting along some direction (to be controlled)
described by the unit vector $\mathbf{\hat{t}}$. We will consider, for the
rest of this work, $\Gamma=0.1$ mm/s2 which is a value compatible with
successful interplanetary low-thrust missions such as Dawn [17], Bepi-Colombo
[18] or Smart-1 [19]. We consider the problem of finding $t_{f}$ and a (piece-
wise continuous) function $\mathbf{\hat{t}}(t)$ with $t\in[t_{0},t_{f}]$ so
that, under the dynamics defined in Eq.(1), the state is steered from any
initial state $\mathbf{r}_{0},\mathbf{v}_{0}$ to the target state
$\mathbf{r}_{t}=R\hat{\mathbf{i}},\mathbf{v}_{t}=\mathbf{0}$ in a minimal
amount of time, thus minimizing the following cost function:
$J=t_{f}-t_{0}=\int^{t_{f}}_{t_{0}}dt$. Following Pontryagin’s Minimum
Principle [10] and some of the practical guidelines suggested by Jiang and
Baoyin [11] we introduce the following Hamiltonian:

$\mathcal{H}(\mathbf{r},\mathbf{v},\boldsymbol{\lambda}_{r},\boldsymbol{\lambda}_{v},\mathbf{\hat{t}})=\boldsymbol{\lambda}_{r}\cdot\mathbf{v}+\boldsymbol{\lambda}_{v}\cdot\bigg{(}-\frac{\mu}{r^{3}}\mathbf{r}-2\boldsymbol{\omega}\times\mathbf{v}-\boldsymbol{\omega}\times\boldsymbol{\omega}\times\mathbf{r}+\Gamma\mathbf{\hat{t}}\bigg{)}+\lambda_{J}$
(2)

where we have introduced the co-states functions $\mathbf{\lambda_{r}}$ and
$\mathbf{\lambda_{v}}$ and an additional $\lambda_{J}$ constant coefficient
that multiplies the actual cost function considered
$J=\lambda_{J}(t_{f}-t_{0})$. The additional coefficient $\lambda_{J}$ results
in increased numerical stability and, as we shall see, offers a useful extra
degree of freedom for data augmentation.

Since the Hamiltonian must be minimized by the optimal control (i.e.
$\mathbf{\hat{t}}^{*}$) during an optimal trajectory, it is straight forward
to derive the classical necessary condition:

$\mathbf{\hat{t}}^{*}=-\frac{\boldsymbol{\lambda}_{\mathbf{v}}}{\lambda_{v}}$

The augmented dynamics is then derived from the Hamiltonian performing the
derivatives $\dot{\mathbf{x}}=\frac{\partial\mathcal{H}}{\partial\mathbf{x}}$,
$\dot{\boldsymbol{\lambda}}=-\frac{\partial\mathcal{H}}{\partial\boldsymbol{\lambda}}$:

$\left\\{\begin{array}[]{l}\dot{\mathbf{r}}=\mathbf{v}\\\
\dot{\mathbf{v}}=-\frac{\mu}{r^{3}}\mathbf{r}-2\boldsymbol{\omega}\times\mathbf{v}-\boldsymbol{\omega}\times\boldsymbol{\omega}\times\mathbf{r}-\Gamma\frac{\boldsymbol{\lambda}_{\mathbf{v}}}{\lambda_{v}}\\\
\dot{\boldsymbol{\lambda}}_{\mathbf{r}}=\mu\left(\frac{\boldsymbol{\lambda}_{\boldsymbol{v}}}{r^{3}}-3(\boldsymbol{\lambda}_{\boldsymbol{v}}\cdot\mathbf{r})\frac{\mathbf{r}}{r^{5}}\right)-\boldsymbol{\omega}\times\boldsymbol{\omega}\times\boldsymbol{\lambda}_{v}\\\
\dot{\boldsymbol{\lambda}}_{\mathbf{v}}=-\boldsymbol{\lambda}_{\mathbf{r}}+2\boldsymbol{\omega}\times\boldsymbol{\lambda}_{v}\end{array}\right.$
(3)

According to Pontryagin’s theory, the optimal transfer will necessarily be a
solution to the above differential equations with the added condition
$\mathcal{H}=0$ as we consider a free time problem. The inverse is also true:
any solution to the above system of differential equations, respecting the
added condition on the Hamiltonian, is an optimal trajectory and can thus be
used to learn from, regardless of where it actually starts from. A point of
view that, as we shall see, is at the basis of the backward generation of
optimal examples here used to form the database of state action pairs.

Figure 1: Nominal trajectory used for policy learning

### Finding one optimal trajectory

In order to find one solution to the optimal control problem outlined in the
previous section, we introduce the following shooting function:

$\phi(\boldsymbol{\lambda}_{r_{0}},\boldsymbol{\lambda}_{v_{0}},\lambda_{J},t_{f})=\left\\{\mathbf{r}(t_{f})-\mathbf{r}_{t},\mathbf{v}(t_{f})-\mathbf{v}_{t},\mathcal{H},||\boldsymbol{\lambda}||-1\right\\}$
(4)

where $\mathbf{r}(t_{f})$ and $\mathbf{v}(t_{f})$ (and $\mathcal{H}$) are
obtained by numerically propagating Eq.(3) from the assigned initial
conditions
$\mathbf{r}_{0},\mathbf{v}_{0},\boldsymbol{\lambda}_{r_{0}},\boldsymbol{\lambda}_{v_{0}}$
for a time $t_{f}$. The shooting function is a system of eight non-linear
relations between the variables
$\boldsymbol{\lambda}_{r_{0}},\boldsymbol{\lambda}_{v_{0}},\lambda_{J},t_{f}$:
its root is here found numerically applying the sequential quadratic
programming solver SNOPT [20] and defines entirely an optimal trajectory. Note
that the condition on the initial co-state norm is not strictly necessary in
this case, but it helps with the numerical stability of the root finder.

As an example, let us here consider an orbit having, at epoch $t_{0}$, the
following osculating orbital Keplerian elements in the inertial frame:
$a=2.687$ AU, $e=0.23$, $i=0.116$ rad., $\Omega=3.137$ rad., $\omega=4.453$
rad. and eccentric anomaly $E=3.01$ rad. Without loss of generality we assume
that, at $t_{0}$, the rotating frame is coincident to the inertial frame, so
that a rotation matrix between the two frames will take the simple form:

$\mathbf{R}(t)=\left[\begin{array}[]{ccc}\cos{\omega(t-t_{0})}&\sin{\omega(t-t_{0})}&0\\\
-\sin{\omega(t-t_{0})}&\cos{\omega(t-t_{0})}&0\\\ 0&0&1\end{array}\right]$

Given the Keplerian osculating elements at $t$ in the inertial frame we thus
find the initial conditions $\mathbf{r}_{0},\mathbf{v}_{0}$ in the rotating
frame $\mathcal{F}$ and find the root of Eq.(4) when the target orbit is at
$R=1.3$ AU. The resulting trajectory is shown in Figure 1 and has a total time
of flight of $t_{f}^{*}=4.62$ years. Being a solution to the shooting equation
is a necessary condition for optimality but not a sufficient one. As a
consequence, there may be multiple roots for the shooting equation
corresponding to different transfer strategies. This corresponds to the
existence of multiple local minima and is a problem hard to tackle rigorously
within Pontryagin’s theory. Typically these multiple roots can be found
restarting the root finder from different initial guesses for the co-states,
an approach that, although not rigorous, is often taken and deemed as
unavoidable. For this particular nominal trajectory we verified experimentally
that no other root exists via multiple restarts, but as we shall see later
when assembling a large dataset of nominal trajectories, it is challenging to
guarantee mathematically that this is the case.

### Augmenting one optimal trajectory

To create a dataset of optimal trajectories one typically would solve multiple
times the shooting function in Eq.(4) using different values of
$\mathbf{r}_{0}$ and $\mathbf{v}_{0}$. In this way, a significant amount of
computational resources are needed. In more difficult cases, where a homotopy
technique is necessary (see [21] for an example) and local minima and
convergence are issues, it may just be impossible to compute the number of
optimal trajectories required to train an accurate model. A far more efficient
method is to use the Backward Generation of Optimal Examples [2]. This
technique can be seen as a data augmentation technique where each sample in a
dataset of optimal trajectories can be used, at low computational cost, to
generate many more different samples. In order to apply this technique to the
particular case at hand, we indicate the final values of the co-states
$\boldsymbol{\lambda}_{r}$, $\boldsymbol{\lambda}_{v}$ relative to one optimal
trajectory with $\boldsymbol{\lambda}_{f}$. We thus introduce the update rule:

$\boldsymbol{\lambda}_{f}^{+}=\boldsymbol{\lambda}_{f}+\boldsymbol{\lambda}_{f}\cdot\boldsymbol{\Delta}$
(5)

Where $\boldsymbol{\Delta}$ is a perturbation vector chosen as to satisfy all
of the conditions dictated by Pontryagin’s principle, in our case only
$\mathcal{H}_{f}=0$. In this work we will often use for each element of the
perturbation vector a uniformly distributed number
$\mathcal{U}(-\delta,\delta)$. The benefit of having introduced the
normalisation factor $\lambda_{J}$ appears now clearly. In fact, setting
Eq.(2) to zero, we see how we are allowed to choose all of the perturbations
on the co-states freely since we can always set $\lambda_{J}^{new}$ as follow:

$\lambda_{J}^{+}=-\boldsymbol{\lambda}_{\mathbf{r}_{f}}^{+}\cdot\mathbf{v}_{f}-\boldsymbol{\lambda}_{\mathbf{v}_{f}}^{+}\cdot\bigg{(}-\frac{\mu}{r_{f}^{3}}\mathbf{r}_{f}-2\boldsymbol{\omega}\times\mathbf{v}_{f}-\boldsymbol{\omega}\times\boldsymbol{\omega}\times\mathbf{r}_{f}-\Gamma\frac{\boldsymbol{\lambda}_{\mathbf{v}_{f}}^{+}}{{\lambda}_{v_{f}}^{+}}\bigg{)}$

and thus guarantee, for any arbitrary perturbation vector, that
$\mathcal{H}_{f}=0$. Note that with respect to the case treated in [2] we have
here simplified the overall complexity of the dataset generation further by
avoiding the solution of a nonlinear equation representing the transversality
condition on the Hamiltonian. Once a new valid value
$\boldsymbol{\lambda}_{f}^{+}$ for the final co-states is available we
propagate backward Eq.(3) and thus obtain an optimal trajectory to be added to
the dataset. In this work we use a Taylor technique [22] to perform a highly
accurate numerical propagation and we are able to generate 10,000 optimal
trajectories in the order of seconds: orders of magnitude less than directly
solving the corresponding shooting function. Of course, it must be noted, the
exact initial conditions of all these optimal trajectories are not chosen, but
a consequence of the new final co-states values.

## Neural representation of the optimal policy

In this section we train a neural model $\mathcal{N}$ to represent the optimal
policy $\mathbf{\hat{t}}^{*}$ as a function of the spacecraft state
$\mathbf{r},\mathbf{v}$. The resulting neural state feedback
$\mathbf{\hat{t}}^{*}=\mathcal{N}(\mathbf{r},\mathbf{v})$, called a G&CNET, is
then used in Eq.(1) to simulate the spacecraft dynamics.

### The dataset

Figure 2: Bundle of 200 optimal trajectories obtained applying the backward
generation of optimal examples. Co-state perturbation is 1 ‰(black), 8%
(yellow).

We generate a dataset of optimal state-action pairs augmenting one single
nominal trajectory via the Backward Generation of Optimal Examples. The
generic single entry in our dataset thus contains the values
$[X_{i},Y_{i}]=(x_{i},y_{i},z_{i},{v_{x}}_{i},{v_{y}}_{i},{v_{z}}_{i}),({t_{x}}_{i}^{*},{t_{y}}_{i}^{*},{t_{z}}_{i}^{*})$.
These are computed generating 200,000 optimal trajectories sampled in 100
points equally distanced in time. The resulting dataset thus contains
20,000,000 optimal state-action pairs. Half of the dataset is generated
perturbing the final co-states $\lambda_{f}^{+}$ using $\delta=0.001$ and
back-propagating the augmented equations of motion for a time $(1+c)t^{*}_{f}$
with the constant $c$ uniformly distributed $c\in[0,1]$. The resulting bundle,
shown in Figure 2, grows very thin as the spacecraft approaches the final
rendezvous. As a result, small deviations from the nominal trajectory may
lead, during a real transfer, the spacecraft state to be outside of the
training data and the resulting error to grow to unacceptable levels. To
overcome this issue, the remaining half of the dataset is generated using a
larger perturbation $\delta=0.08$. The resulting dataset is visualised in
Figure 2, where the different colours correspond to the two halves of the
dataset. We note how the larger perturbation leads to trajectories that may
also depart significantly from the nominal trajectory. Interestingly, during
our experiments, we found these trajectories important to a successful
training procedure, likely as they thicken the bundle during the second half
of the interplanetary transfer. The CPU time required to assemble the entire
dataset, on a standard modern CPU (Intel® Core™ i5-6600K CPU, 3.50 GHz), is
less than one minute. We also generate a completely independent test set
containing 4,000 whole trajectories computed using $\delta=0.0005$ and back-
propagating the augmented equations of motion for a time $(1+c)t^{*}_{f}$. The
smaller perturbation vector used to generate the test set is aimed at making
sure to only probe the model within the bounds of the training data. It is
here worth mentioning, though, that we also used and generated a test set with
the same perturbation vector as that used to generate the training set and
obtained, essentially, similar results. All numerical integrations are
performed using Heyoka [22], a Taylor based method, and setting the numerical
tolerance to machine level, i.e. $\epsilon=1\cdot 10^{-16}$.

$\vdots$$\vdots$$\vdots$$\vdots$$x$$y$$z$$v_{x}$$v_{y}$$v_{z}$$t_{x}$$t_{y}$$t_{z}$(Softplus)
Hiddenlayer(Softplus) Hiddenlayer(Softplus) Hiddenlayer(Softplus)
Hiddenlayer(Linear) Outputlayer Figure 3: Feed forward neural network
architecture

### The neural model

We consider learning the optimal thrust from the spacecraft states as a
supervised learning task and we make use of a simple feed forward neural
network, as introduced in [1, 2] for some related tasks. A softplus activation
function is used for the input and the hidden layers, allowing to obtain a
continuous representation of the optimal controls differentiable to high order
[9]. We use four hidden layers of 700 neurons each. For the output layer we
use a linear activation function, avoiding possible saturation issues. The
network inputs are the states ($x$, $y$, $z$, $v_{x}$, $v_{y}$, $v_{z}$) and
the outputs are the thrust direction ($t_{x}$, $t_{y}$, $t_{z}$), see Figure
3. We will also use the notation
$\mathcal{N}(\mathbf{r},\mathbf{v})=\hat{\mathbf{t}}_{nn}$ to indicate the
network predictions.

As a loss function we make use of the cosine similarity between the estimated
thrust direction vector $\mathbf{\hat{t}}_{nn}$ and the ground truth
$\mathbf{\hat{t}}^{*}$. In particular, the cosine similarity is defined as:

$\text{cosine\\_similarity}\left(\mathbf{\hat{t}}^{*},\mathbf{\hat{t}}_{nn}\right)=\frac{\mathbf{\hat{t}}^{*}\cdot\mathbf{\hat{t}}_{nn}}{\mathbf{|\hat{t}}^{*}||\mathbf{\hat{t}}_{nn}|}$

and the loss used:

$\mathcal{L}\left(\mathbf{\hat{t}}_{nn},\mathbf{\hat{t}}^{*}\right)=1-\text{cosine\\_similarity}\left(\mathbf{\hat{t}}_{nn},\mathbf{\hat{t}}^{*}\right)$

This formulation of the loss function allows the network to focus on
predicting the thrust direction only, ignoring any effect of the predicted
norm. Together with the last linear layer used, this allows for an efficient
training. We use $80\%$ of the dataset for training, and the remaining $20\%$
for validation. The network weights are learned using the Adam optimizer [23]
with a starting learning rate $l=0.5\cdot 10^{-4}$ and no weight decay. A
learning rate schedule was imposed by reducing the rate by a factor $f=0.9$
after $p=10$ consecutive epochs without any improvement on the loss as
computed on a random batch of the validation data. Numerical experiments were
performed with other losses as well as alternative representations of the
thrust direction, revealing, in our case, this setup to be the most accurate.

### Results and discussion

Figure 4: Loss during a typical training (left). Networks predictions vs.
ground truth during one optimal trajectory in the test set (centre and right).

We train several neural models following the methodology outlined in the
previous sections. In Figure 4 we report the loss decrease during a typical
training, as well as an example of the predictions made by the neural network
on a previously unseen optimal trajectory from the test set. In particular,
the network prediction is transformed into spherical coordinates
$\theta=\arccos{t_{z}}$ and $\varphi=\arctan{\frac{t_{y}}{t_{x}}}$ and
compared to the optimal ground truth. The final loss value decreases below
$10^{-3}$ and corresponds, as evaluated on the test set, to an average error
of $0.59$ degrees in the prediction of the optimal direction. In order to
assess if this error is small enough to allow the network to successfully
steer the spacecraft optimally, we simulate the spacecraft dynamics in Eq.(1)
when $\hat{\mathbf{i}}(t)=\mathcal{N}(\mathbf{r}(t),\mathbf{v}(t))$, that is
when the thrust direction is instantaneously given by the network predictions.
We consider all initial conditions in the test set, numerically integrating
the equations of motion from there. We record the system state achieved at the
time $t_{f}^{*}$ as well as at the time when the target semi-major axis
$a_{t}$ is acquired exactly, at a time past $t_{f}^{*}$. The results, for the
nominal trajectory used as test case in the previous sections, are given in
Table 1 (case N) and show how the spacecraft is able to rendezvous with the
target body with good accuracy when the low-thrust profile follows the trained
neural model. The network, acting as a state-feedback, is able to cancel out,
rather than accumulate, the prediction errors possibly made during the
transfer. To confirm quantitatively the method accuracy and study its
reliability, we repeat the very same experiment starting from 9 different
nominal trajectories, identified with the capital letters A,B,..I, having all
the same orbital elements, but a different initial eccentric anomaly and hence
a difference phasing to the target. In all cases the trained G&CNET is able to
steer the spacecraft to the target rendezvous reliably and with similar
accuracies. When integrating until $t_{f}^{*}$, the spacecraft is left within
$~{}500,000$ Km and with a residual relative velocity of the order of $~{}100$
[m/s]. In the better cases (trajectories B, F) residuals of less than $~{}20$
m/s are reached on the final velocity and of the order of 100,000 km on the
position. These residuals are compatible with an operational use of the
network for automated on-board guidance and control. It has to be noted here
that the numbers appearing in table 1 have been assembled using the very same
dataset generation and training parameters, e.g. $\delta$, $c$, $l$. Since the
sensitivity of the training accuracy to these parameters is rather
significant, it is expected that during a real design phase one would be able
to tune these parameters to the specific case (and dynamics) as to obtain a
G&CNET fulfilling operational constraints to the requested accuracies.

Table 1: Mean errors at target when propagating the neural dynamics. The numerical integration stops either at $t^{*}_{f}$ or when the spacecraft achieves the target orbital energy (i.e. semi major axis $a_{t}$). Case |  $E_{0}$, rad |  |  $||\boldsymbol{\delta}\mathbf{r}||$, AU |  $||\boldsymbol{\delta}\mathbf{v}||$, km/s |  $\delta a$, AU | $\delta e$ |  $\delta i$, deg |  $\delta t$, years
---|---|---|---|---|---|---|---|---
N | $3.01$ | stop at $t_{f}^{*}$ | $1.15\cdot 10^{-3}$ | $7.65\cdot 10^{-2}$ | $7.60\cdot 10^{-4}$ | $2.88\cdot 10^{-3}$ | $1.57\cdot 10^{-2}$ | 0
stop at $a_{t}$ | $1.16\cdot 10^{-3}$ | $7.46\cdot 10^{-2}$ | 0 | $2.77\cdot 10^{-3}$ | $1.60\cdot 10^{-2}$ | $2.52\cdot 10^{-3}$
A | $0.0$ | stop at $t_{f}^{*}$ | $1.48\cdot 10^{-3}$ | $4.65\cdot 10^{-2}$ | $3.16\cdot 10^{-5}$ | $1.55\cdot 10^{-3}$ | $1.13\cdot 10^{-2}$ | 0
stop at $a_{t}$ | $1.49\cdot 10^{-3}$ | $4.68\cdot 10^{-2}$ | 0 | $1.56\cdot 10^{-3}$ | $1.14\cdot 10^{-2}$ | $4.87\cdot 10^{-5}$
B | $0.698$ | stop at $t_{f}^{*}$ | $5.09\cdot 10^{-4}$ | $1.68\cdot 10^{-2}$ | $1.39\cdot 10^{-4}$ | $5.63\cdot 10^{-4}$ | $1.79\cdot 10^{-2}$ | 0
stop at $a_{t}$ | $4.94\cdot 10^{-4}$ | $1.66\cdot 10^{-2}$ | 0 | $5.57\cdot 10^{-4}$ | $1.79\cdot 10^{-2}$ | $4.98\cdot 10^{-4}$
C | $1.396$ | stop at $t_{f}^{*}$ | $2.73\cdot 10^{-3}$ | $1.58\cdot 10^{-1}$ | $2.27\cdot 10^{-3}$ | $5.93\cdot 10^{-3}$ | $7.98\cdot 10^{-2}$ | 0
stop at $a_{t}$ | $2.76\cdot 10^{-3}$ | $1.48\cdot 10^{-1}$ | 0 | $5.43\cdot 10^{-3}$ | $8.23\cdot 10^{-2}$ | $7.60\cdot 10^{-3}$
D | $2.094$ | stop at $t_{f}^{*}$ | $3.50\cdot 10^{-3}$ | $1.50\cdot 10^{-1}$ | $2.20\cdot 10^{-3}$ | $5.87\cdot 10^{-3}$ | $2.37\cdot 10^{-2}$ | 0
stop at $a_{t}$ | $3.41\cdot 10^{-3}$ | $1.44\cdot 10^{-1}$ | 0 | $5.40\cdot 10^{-3}$ | $2.37\cdot 10^{-2}$ | $7.11\cdot 10^{-3}$
E | $2.793$ | stop at $t_{f}^{*}$ | $1.81\cdot 10^{-3}$ | $8.94\cdot 10^{-2}$ | $8.77\cdot 10^{-4}$ | $3.39\cdot 10^{-3}$ | $1.04\cdot 10^{-2}$ | 0
stop at $a_{t}$ | $1.84\cdot 10^{-3}$ | $8.79\cdot 10^{-2}$ | 0 | $3.30\cdot 10^{-3}$ | $1.08\cdot 10^{-2}$ | $2.96\cdot 10^{-3}$
F | $3.491$ | stop at $t_{f}^{*}$ | $1.21\cdot 10^{-3}$ | $1.56\cdot 10^{-2}$ | $2.09\cdot 10^{-4}$ | $5.57\cdot 10^{-4}$ | $6.40\cdot 10^{-3}$ | 0
stop at $a_{t}$ | $1.15\cdot 10^{-3}$ | $1.49\cdot 10^{-2}$ | 0 | $5.0\cdot 10^{-4}$ | $6.12\cdot 10^{-3}$ | $7.21\cdot 10^{-4}$
G | $4.189$ | stop at $t_{f}^{*}$ | $6.46\cdot 10^{-4}$ | $3.54\cdot 10^{-2}$ | $4.90\cdot 10^{-4}$ | $1.17\cdot 10^{-3}$ | $7.02\cdot 10^{-3}$ | 0
stop at $a_{t}$ | $6.31\cdot 10^{-4}$ | $3.18\cdot 10^{-2}$ | 0 | $9.69\cdot 10^{-4}$ | $5.14\cdot 10^{-3}$ | $1.51\cdot 10^{-3}$
H | $4.887$ | stop at $t_{f}^{*}$ | $1.28\cdot 10^{-3}$ | $5.22\cdot 10^{-2}$ | $2.49\cdot 10^{-4}$ | $2.0\cdot 10^{-3}$ | $8.48\cdot 10^{-3}$ | 0
stop at $a_{t}$ | $1.28\cdot 10^{-3}$ | $5.19\cdot 10^{-2}$ | 0 | $1.98\cdot 10^{-3}$ | $8.52\cdot 10^{-3}$ | $7.99\cdot 10^{-4}$
I | $5.585$ | stop at $t_{f}^{*}$ | $3.18\cdot 10^{-3}$ | $1.41\cdot 10^{-1}$ | $6.02\cdot 10^{-4}$ | $3.78\cdot 10^{-3}$ | $3.64\cdot 10^{-2}$ | 0
stop at $a_{t}$ | $3.20\cdot 10^{-3}$ | $1.42\cdot 10^{-1}$ | 0 | $3.81\cdot 10^{-3}$ | $3.77\cdot 10^{-2}$ | $2.11\cdot 10^{-3}$

## Neural representation of the value function

Figure 5: Initial conditions in the value function learning dataset.
Astronomical units used.

In this section we train a neural model $\mathcal{N}$ to represent the value
function, or cost-to-go for our problem, that is
$v(\mathbf{r},\mathbf{v})=J=t_{f}-t_{0}$: the optimal time of flight as a
function of the state. Note that, in theory, if we were able to learn the
value function exactly, the optimal policy could also be derived from
Hamilton-Jacobi-Bellman (HJB) equations as:

$\mathbf{t}^{*}=\operatorname*{arg\,min}_{|\mathbf{t}|=1}(\nabla_{\mathbf{r}}v\cdot\dot{\mathbf{r}}+\nabla_{\mathbf{v}}v\cdot\dot{\mathbf{v}})\rightarrow\mathbf{t}^{*}=-\frac{\nabla_{\mathbf{v}}v}{|\nabla_{\mathbf{v}}v|}$

this fact, in a different context, was recently successfully used by Izzo and
Ozturk who used the HJB equations to directly inform the network training [2].
The technique, developed independently but related to the emerging field of
physically informed neural networks [24], made use of a training set of
optimal trajectories defined in an ample bundle surrounding one nominal
solution. In general, when the aim is to represent the optimal policy (i.e.
the optimal low-thrust profile) in proximity to one nominal trajectory, for
example for guidance and control purposes, the approach pursued in the
previous section (i.e. a G&CNET) is already able to produce very accurate
results and while alternatives exist, they will not here be pursued. Instead,
we here seek to approximate the value function for the purpose of predicting
the optimal rendezvous time-of-flight from any initial starting condition in a
wide area of the entire state space, as opposed to targeting only a bundle
around one nominal trajectory. In particular we focus here in a scenario where
any orbit in an asteroid belt has to rendezvous to a target object placed in a
circular, zero-inclined orbit.

### The dataset

We generate a dataset of optimal state-time-of-flight pairs by solving the
shooting function in Eq.(4) for 4,000 different values of
$\mathbf{r}_{0},\mathbf{v}_{0}$. We generate these values by randomly sampling
the osculating Keplerian elements in the inertial frame and then converting
them into the Cartesian state in $\mathcal{F}$. The sampling bounds are chosen
as to cover a good portion of the inner and middle asteroid belt with
$a\in[2,3]$ AU, $e\in[0,0.25]$, $i\in[0,5]$ deg. and
$\Omega,\omega,E\in[0,360]$ deg. In Figure 5 the initial conditions are
visualised as well as the target position for the rendezvous. The generic
single entry in our dataset thus contains the values
$[X_{i},Y_{i}]=[(x_{i},y_{i},z_{i},{v_{x}}_{i},{v_{y}}_{i},{v_{z}}_{i}),(t_{f_{i}}^{*})]$.
In this setup, the problem of local minima must be addressed with particular
care as to avoid inserting in the dataset wrongful and deceiving information
as well as avoiding database biases. To this aim, we try to solve each of the
4,000 shooting functions in $N_{trials}=120$ separate experiments ,
corresponding to the number of CPUs available to us. In each case, a maximum
of $N_{res}=30$ restarts are tried to find one solution starting from randomly
sampled co-states in $[-1,1]$ and times-of-flight in $[1.5,3.]$ years. Out of
the 120 guesses, if at least $N_{sol}=50$ result in a solution we then find,
among them, the one corresponding to the smallest time of flight and insert it
in the dataset, else we generate a new initial condition for the state and try
again. We report that roughly 8% of the initial conditions failed to meet this
criteria (i.e. for those cases $N_{sol}<50$) and were thus substituted. Note
that this kind of filter is likely to create a bias on the dataset towards
simpler problems, as it discards the instances leading to a poor convergence
of the root finder which may be also correlated to difficult cases for the
regression, but the small fraction of discarded conditions was deemed to be
acceptable for our purposes. We also report that in the creation of the
dataset roughly 20% of cases resulted in multiple locally optimal solutions
(mostly two, of which we inserted in the dataset the one with smallest
$t_{f}^{*}$) while the rest converged consistently (all $N_{trials}$ times) to
the very same transfer. This procedure, although computationally expensive
(taking less than half a day of computations on our supercomputer), allows to
compile a dataset that includes correctly difficult cases for a low-thrust
transfer such as those, for example, close to the hyper-surface of the state
space separating trajectories requiring a different number of revolutions.
These are in fact the conditions that often result in local minima or mis-
convergence and that, if excluded, would result in a biased (and easier)
dataset. If these conditions were instead included without caring for the
possibility of local minima, the dataset would contain deceiving information
for the model training. Eventually, we use 75% trajectories for training and
the remaining as a test set.

We then produce also a second dataset, augmenting all the 3,000 training
trajectories using the backward generation of optimal examples. For each
nominal trajectory we produce 32 new optimal examples perturbing the final co-
states by $\delta=0.0001$ and integrating back in time up to the optimal
$t_{f}^{*}$. The resulting augmented dataset is created in less than 10
seconds on a standard CPU and contains 96,000 optimal transfers from the belt
to the target.

We release openly both dataset produced as well as the exact train/test split
used in this manuscript on the CERN owned platform Zenodo 333The dataset can
be downloaded here: https://doi.org/10.5281/zenodo.6365326 and was assigned
the following d.o.i. 10.5281/zenodo.6365326 .

### The neural model

We consider the task of learning the value function (i.e. the optimal time of
flight) for the 1,000 trajectories of the test set as a supervised learning
task, essentially a nonlinear regression problem. We use a simple feed-forward
neural network with rectified linear units (ReLu) and 4 hidden layers with 100
neurons each. From the Cartesian elements $X$ contained in the dataset, we
compute the modified equinoctial elements [25]. In order to avoid issues
related to the discontinuity of its representation, instead of the true
longitude $L$, we use $\cos L,\sin L$ and we indicate the new set of
attributes, now with an increased dimension, with $X_{eq}$. As a last step we
scale all attributes uniformly in the range $[-1,1]$ via the simple
transformation:

$\begin{array}[]{l}X_{in}=2(X_{eq}-X_{min})/(X_{max}-X_{min})-1\\\
\end{array}$ (6)

For the optimizer we use the same setup (Adam [23]) as in the experiments on
the optimal policy, changing only the initial learning rate to $l=0.5\cdot
10^{-3}$. Note that this regression problem, having a relatively small size,
can also be tackled using different machine learning approaches (see for
example [15]) potentially offering improvements over the neural model here
described. The neural approach here used, on the other hand, is of particular
interest when coupled with the data augmentation technique employed, the
backward generation of optimal examples, which creates a much larger training
set to exploit.

### Results and discussion

Figure 6: Error on the training/test set.

The results for the best models trained over both the original and the
augmented dataset are visualised in Figure 6. In general we note how, in
relative terms, the error on the value function is now much larger than the
one we obtained in a previous section when we trained the network to represent
the optimal policy (even though still smaller than 4% on average). This is
expected and consequential to two different reasons. Firstly, we are tackling
a much more difficult problem in that the diversity of the trajectories used
is vastly larger, covering the entire asteroid belt, not only the
neighbourhood of one nominal trajectory. Secondly, we are here learning from a
much smaller and less dense dataset, a consequence of the fact that many
nominal trajectories have here to be computed, not only one. Nevertheless, the
accuracy of $25.42$ days is greatly higher than, for example, the $468.56$
days one obtains using the Edelbaum approximation [26]. This shows how the
neural model developed is capable of capturing the complex phenomenon of
phasing as well as effects connected to the orbits eccentricity and spatial
orientation, aspects entirely neglected by the Edelbaum approach. A second
fact worth noting is the advantage offered by the use of data augmentation.
While all the augmented trajectories are in some neighbourhood of one entry of
the original dataset, and thus add only a limited amount of new information,
their presence allows for a significant improvement on the prediction (26%
decrease in MAE). It is worth mentioning here that a larger dataset, such as
that produced by the data augmentation, also means higher training times for
the model, something that would likely be accepted in a real application given
the substantial gain introduced and the fact that the database augmentation
cost, in terms of computational time is essentially, negligible. On the right
plot of Figure 6 we see the generalization capabilities of the model on the
test set. These are trajectories whose initial elements are still sampled
uniformly in the same bounds as the training set, but that have never been
encountered while training the model. It is interesting to note how very large
errors are made by the model on a few trajectories, mostly in the test set,
but also in the training set. We thus performed an investigation on these
outliers as to understand their origin. We find that in the test set there are
13 data points whose error on the predicted time of flight is larger than 300
days. For these 13 cases we find that the prediction of the trained neural
model can be explained by one of three reasons.

  * •

The predicted time of flight is actually correct and the label in the dataset,
instead, correspond to a local minimum which slipped through our dataset
creation pipeline. This actually happened for 6 of the 13 outliers: the
majority of cases. This is particularly interesting as it hints at a possible
use of the neural model to find or suggest possible cases where some full
optimization model converges to a local minima.

  * •

The predicted time of flight corresponds, within a few days, to a local
minima. This actually happened for 2 of the 13 outliers.

  * •

The predicted time of flight does not represent any real optimal trajectory.
This actually happened for 5 of the 13 outliers. While these cases are always
bound to exist in the learned mode, removing erroneous labels in the training
set and increasing the dataset size is likely beneficial to alleviate this
issue.

In Figure 7 we visualise one of the troublesome transfers, erroneously
inserted in the dataset, but correctly predicted by the trained model,
together with the correct globally optimal strategy. For this case, our root
finder (applied to the shooting function) converged to the global solution
only after many restarts, confirming this transfer as a particularly difficult
case to find for any dataset creation procedure. One could of course try to
mitigate this problem when creating a dataset by employing more restarts, but
no number guarantees the absence of these occurrences. On the other hand, as
here shown, the trained model does provide hints as to where they may occur,
something that can be exploited _a posteriori_.

Figure 7: An outlier in the neural model predictions over the test set. The
globally optimal transfer is shown on the left, while the one consistently
found by our numerical approach (and thus ending up in the dataset) is shown
on the right.

## Conclusion

The time optimal low-thrust problem for a constant acceleration rendezvous has
been studied in the context of neural modelling. Both the optimal policy and
value function representation are studied. We find that a G&CNET, a deep feed
forward neural network trained on optimal examples, is able to represent to
great accuracy the optimal policy. Our results suggest that the on-board use
of similar neural models for the spacecraft guidance and control is able to
fulfil well possible operational requirements on the velocity and position
residuals, thus constituting a valid alternative to currently proposed
methods. We also show how the dataset of optimal trajectories employed to
backpropagate the model errors can be created in seconds rather than days, if
the data augmentation technique called backward generation of optimal examples
is used. In the case of value function learning we find that a neural model is
able to represent the value function in a vast region of the state space
corresponding to the entire asteroid belt and with errors that are below 4% on
average and mostly concentrated in a few outliers related to the presence of
local minima solutions and, ultimately, of very complex globally optimal
transfer strategies that compete with much simpler locally optimal ones. We
show that employing the data augmentation technique also in this case allows
to significantly improve the predictions made by the neural model, leading to
a 26% decrease in MAE. While the dynamical model used in this paper as a test
case is not appropriate for many cases of practical interest, it served its
purpose of resulting in quick optimisation procedures and thus allowing to
refine the overall work logic and results. We hope that our work will
facilitate, in the future, the training of neural models representing more
complex cases.

## References

  * Sánchez-Sánchez and Izzo [2018] Sánchez-Sánchez, C., and Izzo, D., “Real-time optimal control via deep neural networks: study on landing problems,” _Journal of Guidance, Control, and Dynamics_ , Vol. 41, No. 5, 2018, pp. 1122–1135.
  * Izzo and Öztürk [2021] Izzo, D., and Öztürk, E., “Real-Time Guidance for Low-Thrust Transfers Using Deep Neural Networks,” _Journal of Guidance, Control, and Dynamics_ , Vol. 44, No. 2, 2021, pp. 315–327.  10.2514/1.G005254.
  * Li et al. [2019a] Li, H., Baoyin, H., and Topputo, F., “Neural networks in time-optimal low-thrust interplanetary transfers,” _IEEE Access_ , Vol. 7, 2019a, pp. 156413–156419.
  * Federici et al. [2021] Federici, L., Benedikter, B., and Zavoli, A., “Deep Learning Techniques for Autonomous Spacecraft Guidance During Proximity Operations,” _Journal of Spacecraft and Rockets_ , 2021, pp. 1–12.
  * Hovell and Ulrich [2020] Hovell, K., and Ulrich, S., “On Deep Reinforcement Learning for Spacecraft Guidance,” _AIAA Scitech 2020 Forum_ , 2020, p. 1600.
  * Cheng et al. [2018] Cheng, L., Wang, Z., Jiang, F., and Zhou, C., “Real-time optimal control for spacecraft orbit transfer via multiscale deep neural networks,” _IEEE Transactions on Aerospace and Electronic Systems_ , Vol. 55, No. 5, 2018, pp. 2436–2450.
  * Mohammadzaheri et al. [2012] Mohammadzaheri, M., Chen, L., and Grainger, S., “A critical review of the most popular types of neuro control,” _Asian Journal of Control_ , Vol. 14, No. 1, 2012, pp. 1–11.
  * Bardi et al. [1997] Bardi, M., Dolcetta, I. C., et al., _Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations_ , Vol. 12, Springer, 1997.
  * Tailor and Izzo [2019] Tailor, D., and Izzo, D., “Learning the optimal state-feedback via supervised imitation learning,” _Astrodynamics_ , Vol. 3, No. 4, 2019, pp. 361–374.
  * Pontryagin [1987] Pontryagin, L., _Mathematical Theory of Optimal Processes_ , CRC Press, 1987\.
  * Jiang et al. [2012] Jiang, F., Baoyin, H., and Li, J., “Practical techniques for low-thrust trajectory optimization with homotopic approach,” _Journal of guidance, control, and dynamics_ , Vol. 35, No. 1, 2012, pp. 245–258.
  * Li et al. [2019b] Li, H., Topputo, F., and Baoyin, H., “Autonomous Time-Optimal Many-Revolution Orbit Raising for Electric Propulsion GEO Satellites via Neural Networks,” _arXiv preprint arXiv:1909.08768_ , 2019b.
  * Hennes et al. [2016] Hennes, D., Izzo, D., and Landau, D., “Fast approximators for optimal low-thrust hops between main belt asteroids,” _2016 IEEE Symposium Series on Computational Intelligence (SSCI)_ , IEEE, 2016, pp. 1–7.
  * Li et al. [2020] Li, H., Chen, S., Izzo, D., and Baoyin, H., “Deep networks as approximators of optimal low-thrust and multi-impulse cost in multitarget missions,” _Acta Astronautica_ , Vol. 166, 2020, pp. 469–481.
  * Mereta et al. [2017] Mereta, A., Izzo, D., and Wittig, A., “Machine learning of optimal low-thrust transfers between near-earth objects,” _International Conference on Hybrid Artificial Intelligence Systems_ , Springer, 2017, pp. 543–553.
  * Shen and Casalino [2021] Shen, H.-X., and Casalino, L., “Simple $\Delta$ V Approximation for Optimization of Debris-to-Debris Transfers,” _Journal of Spacecraft and Rockets_ , Vol. 58, No. 2, 2021, pp. 575–580.
  * Russell and Raymond [2011] Russell, C., and Raymond, C., “The dawn mission to Vesta and Ceres,” _The dawn mission to minor planets 4 vesta and 1 ceres_ , Springer, 2011, pp. 3–23.
  * Benkhoff et al. [2010] Benkhoff, J., Van Casteren, J., Hayakawa, H., Fujimoto, M., Laakso, H., Novara, M., Ferri, P., Middleton, H. R., and Ziethe, R., “BepiColombo—Comprehensive exploration of Mercury: Mission overview and science goals,” _Planetary and Space Science_ , Vol. 58, No. 1-2, 2010, pp. 2–20.
  * Foing et al. [2006] Foing, B., Racca, G. D., Marini, A., Evrard, E., Stagnaro, L., Almeida, M., Koschny, D., Frew, D., Zender, J., Heather, J., et al., “SMART-1 mission to the Moon: status, first results and goals,” _Advances in Space Research_ , Vol. 37, No. 1, 2006, pp. 6–13.
  * Gill et al. [2005] Gill, P. E., Murray, W., and Saunders, M. A., “SNOPT: An SQP algorithm for large-scale constrained optimization,” _SIAM review_ , Vol. 47, No. 1, 2005, pp. 99–131.
  * Bertrand and Epenoy [2002] Bertrand, R., and Epenoy, R., “New smoothing techniques for solving bang–bang optimal control problems—numerical results and statistical interpretation,” _Optimal Control Applications and Methods_ , Vol. 23, No. 4, 2002, pp. 171–197.
  * Biscani and Izzo [2021] Biscani, F., and Izzo, D., “Revisiting high-order Taylor methods for astrodynamics and celestial mechanics,” _Monthly Notices of the Royal Astronomical Society_ , Vol. 504, No. 2, 2021, pp. 2614–2628.
  * Kingma and Ba [2014] Kingma, D. P., and Ba, J., “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * Raissi et al. [2017] Raissi, M., Perdikaris, P., and Karniadakis, G. E., “Physics informed deep learning (part I): Data-driven solutions of nonlinear partial differential equations,” _arXiv preprint arXiv:1711.10561_ , 2017.
  * Walker et al. [1985] Walker, M., Ireland, B., and Owens, J., “A set modified equinoctial orbit elements,” _Celestial mechanics_ , Vol. 36, No. 4, 1985, pp. 409–419.
  * Edelbaum [1965] Edelbaum, T. N., “Optimum power-limited orbit transfer in strong gravity fields,” _AIAA Journal_ , Vol. 3, No. 5, 1965, pp. 921–925.