# Reinforced random walks under memory lapses

Manuel González-Navarrete and Ranghely Hernández

###### Abstract

We introduce a one-dimensional random walk, which at each step performs a
reinforced dynamics with probability $\theta$ and with probability $1-\theta$,
the random walk performs a step independent of the past. We analyse its
asymptotic behaviour, showing a law of large numbers and characterizing the
diffusive and superdiffusive regions. We prove central limit theorems and law
of iterated logarithm based on the martingale approach.

##  1 Introduction

In the area of mathematical physics, the study of random dynamics that combine
simplicity in their definition and complexity in their behaviors has drawn
vital relevance. In this sense, the so-called elephant random walk (ERW) has
generated a lot of interest in recent years. The ERW was proposed by Schutz
and Trimper [22] and can be formulated as follows. Let the sequence
$\\{X_{1},X_{2},\ldots\\}$ where $X_{i}\in\\{-1,+1\\}$, for all $i\geq 1$. The
position of the elephant at instant $n$ is given by
$S_{n}=\sum_{i=1}^{n}X_{i},$ and $S_{0}=0$. It is assumed that the elephant
completely remembers its past and performs the $(n+1)$-step by choosing
$t\in\\{1,....n\\}$ uniformly at random, then the elephant decides to repeat
the step $X_{t}$ with probability $p$ or to do the opposite with probability
$1-p$.

The main characteristic of ERW is the existence of a phase transition from a
diffusive to a super-diffusive behavior ([1, 11]). This refers to the fact
that the second moment of the position is scaled linearly (diffusive) or as a
power law (super-diffusive). This fact has motivated the development of new
models and results that generalize the ERW, for instance, Bercu and Laulin [4,
5], Bertoin [7, 8, 9], González-Navarrete and Lambert [13, 15] and Gut and
Stadtmüller [16, 17].

This work is motivated by the recent formulation in [16], which defines an ERW
with delays. In particular, we propose the addition of memory lapses to the
dynamics of the process, this kind of evolution was introduced by [14] and
also used in [2, 10]. In this sense, at each step, with probability $\theta$
the elephant performs its movement by following the dynamics in [16],
otherwise the elephant performs a step independent of its history.

We aim to characterize the asymptotic behaviors of the proposed model. In
particular, it seeks to analyse the influence of the parameters $p\in[0,1]$
and $\theta\in[0,1)$. Specifically, we prove a law of large numbers and
characterize the diffusive and superdiffusive regions, showing convergences in
distribution and laws of iterated logarithm. The main technique we used is the
martingale theory, based on the works [3, 12, 18, 19, 20].

The rest of the paper is organized as follows. In Section 2 we include the
mathematical formulation of the model and state the main results. Section 3
includes some auxiliary results and the proofs of the theorems.

##  2 The model and main results

Gut and Stadtmüller [16] proposed an ERW with memory defined as follows. Let
the first step given by

$X_{1}=\left\\{\begin{array}[]{ll}+1&\ \mbox{, with probability $p$,}\\\ -1&\
\mbox{, with probability $q$,}\\\ 0&\ \mbox{, with probability $r$.}\\\
\end{array}\right.$ (2.1)

where $p+q+r=1$. The next steps are performed by the rule

$X_{n+1}=\left\\{\begin{array}[]{lll}X_{t}&,\ \mbox{with probability $p$,}\\\
-X_{t}&,\ \mbox{with probability $q$,}\\\ 0&,\ \mbox{with probability
$r$,}\end{array}\right.$ (2.2)

where $t$ is uniformly distributed on $\\{1,2,\ldots\,n\\}$. The sequence
$\\{X_{i}\\}$ generates a one-dimensional random walk $S_{n}$ given by:

$\begin{array}[]{ccccc}S_{0}:=0&\text{ and
}&S_{n}=\displaystyle\sum_{i=1}^{n}X_{i}&\forall~{}n=1,2,\ldots\end{array}$
(2.3)

In this sense, [16] proved for $0\leq p-q<1/2$ that

$\frac{S_{n}}{\sqrt{n}}\overset{d}{\rightarrow}(p+q)N\left(0,\frac{p+q}{1-2(p-q)}\right)+r\delta_{0}(x),\
\text{ as }n\to\infty$ (2.4)

where $\overset{d}{\rightarrow}$ means convergence in distribution. Moreover,
if $p-q=1/2$, then

$\frac{S_{n}}{\sqrt{n\cdot\log
n}}\overset{d}{\rightarrow}(p+q)N\left(0,p+q\right)+r\delta_{0}(x),\ \text{ as
}n\to\infty.$ (2.5)

we use the standard $\delta_{0}(x)$ to denote the distribution function with a
jump of height one at $0$. The results in (2.4) and (2.5) are given by the
fact that, if $X_{1}=0$, then $S_{n}=0$, $\forall~{}n\geq 1$, which occurs
with probability $r$. As usual, in the region $1/2<p-q<1$, there exists a
random variable with positive variance, such that
$\lim\limits_{n\to\infty}\frac{S_{n}}{n^{p-q}}=W$ almost surely as
$n\to\infty$.

Based on (2.2) we propose a model with random tendency as described in [13].
At each step $n\geq 1$, we draw a Bernoulli random variable $Y_{n}$ with
parameter $\theta$, such that for $Y_{n}=0$ the random walk performs a step
independent of its past, otherwise if $Y_{n}=1$ the elephant follows the
dynamics proposed in (2.2). The conditional probabilities can be formulated as
follows. Let $\mathcal{F}_{n}$ the $\sigma$-algebra generated by
$\\{X_{1},X_{2},\ldots,X_{n}\\}$, for all $n=1,2,\ldots$

$\begin{array}[]{c}{\mathbb{P}}(X_{n+1}=+1|\mathcal{F}_{n})=\left(\dfrac{n_{+}}{n}\cdot
p+\dfrac{n_{-}}{n}\cdot q\right)\cdot\theta+(1-\theta)\cdot p\\\\[14.22636pt]
{\mathbb{P}}(X_{n+1}=-1|\mathcal{F}_{n})=\left(\dfrac{n_{-}}{n}\cdot
p+\dfrac{n_{+}}{n}\cdot q\right)\cdot\theta+(1-\theta)\cdot q\\\\[14.22636pt]
{\mathbb{P}}(X_{n+1}=0|\mathcal{F}_{n})=\left(\dfrac{n_{0}}{n}\cdot(p+q)+r\right)\cdot\theta+(1-\theta)\cdot
r=\dfrac{n_{0}}{n}\cdot(p+q)\cdot\theta+r\end{array}$ (2.6)

where $\theta\in[0,1)$, $n_{+}=\\#\\{i\in\\{1,2,\ldots,n\\}:X_{i}=+1\\}$, is
the number of $+1$ steps until the $n$-step, and $n_{-}$ and $n_{0}$ represent
the number of steps in the directions $-1$ and $0$ until the $n$-step,
respectively.

###  2.1 Main results

We use the following notation

$\alpha=(p-q)\cdot\theta,\ \omega=(p-q)(1-\theta),\ \tau=(1-\theta)(p+q),\
\gamma=(p+q)\theta\text{ and
}\phi=\frac{\tau}{1-\gamma}-\left(\frac{\omega}{1-\alpha}\right)^{2}$ (2.7)

and let denote the sequence of constants

$\begin{array}[]{cccc}a_{1}:=1&\text{ and
}a_{n}:=\displaystyle\prod_{k=1}^{n-1}\left(1+\dfrac{\alpha}{k}\right)=\dfrac{\Gamma(n+\alpha)}{\Gamma(n)\cdot\Gamma(\alpha+1)}&\text{
for }&n=2,\ldots\end{array}$ (2.8)

The main results are

######  Theorem 2.1.

Let the RW given by (2.6), for all $\alpha\in[0,1)$

$\begin{array}[]{lll}\lim\limits_{n\to\infty}\frac{S_{n}-{\mathbb{E}}(S_{n})}{n}=0&\text{
a.s}\end{array}$ (2.9)

and

$\begin{array}[]{ll}\lim\limits_{n\to\infty}\frac{S_{n}}{n}=\frac{\omega}{1-\alpha}&a.s.\end{array}$
(2.10)

######  Theorem 2.2.

Let the RW given by (2.6):

  * i)

If $0\leq\alpha<1/2$, then

$\frac{S_{n}-{\mathbb{E}}[S_{n}]}{\sqrt{\frac{\phi}{1-2\alpha}n}}\overset{d}{\rightarrow}N(0,1),\text{
as }\ n\to\infty$ (2.11)

and

$\limsup_{n\to\infty}\pm\frac{S_{n}-{\mathbb{E}}[S_{n}]}{\sqrt{2\frac{\phi}{1-2\alpha}n\log\log\left(\frac{\phi\Gamma(\alpha+1)^{2}}{1-2\alpha}n^{1-2\alpha}\right)}}=1\
\text{ a.s }$ (2.12)

  * ii)

If $\alpha=1/2$, then

$\frac{S_{n}-{\mathbb{E}}[S_{n}]}{\sqrt{\phi\cdot n\cdot\log
n}}\overset{d}{\rightarrow}N(0,1),\text{ as }\ n\to\infty$ (2.13)

and

$\limsup_{n\to\infty}\pm\frac{S_{n}-{\mathbb{E}}[S_{n}]}{\sqrt{2\phi\cdot
n\cdot\log n\cdot\log\log(\phi\Gamma(3/2)\log n)}}=1\ \text{ a.s }$ (2.14)

######  Theorem 2.3.

If $1/2<\alpha<1$, then there exists a random variable $W$ with positive
variance such that

$\lim\limits_{n\to\infty}\frac{S_{n}-{\mathbb{E}}[S_{n}]}{a_{n}}=W\
\text{a.s.}\ \text{and in }L^{2},$ (2.15)

Moreover,

$\frac{S_{n}-{\mathbb{E}}[S_{n}]-W\cdot
a_{n}}{\sqrt{\frac{\phi}{2\alpha-1}n}}\overset{d}{\rightarrow}N(0,1),\text{ as
}\ n\to\infty$ (2.16)

and

$\limsup_{n\to\infty}\pm\frac{S_{n}-{\mathbb{E}}[S_{n}]-W\cdot
a_{n}}{\sqrt{2\frac{\phi}{2\alpha-1}n\log\left|\log\left(\frac{\phi\Gamma(\alpha+1)^{2}}{2\alpha-1}n^{1-2\alpha}\right)\right|}}=1\
\text{ a.s }$ (2.17)

We remark that the distribution and higher moments of $W$ are of interest. In
the case of original ERW we refer [3, 6] for details.

##  3 Proofs

First, we need to introduce some definitions and auxiliary results

######  Definition 3.1.

Let define the sequence

$\begin{array}[]{ccc}M_{n}=\dfrac{S_{n}-{\mathbb{E}}(S_{n})}{a_{n}}&\text{ for
}&n=1,2,\ldots\end{array}$ (3.1)

where $a_{n}$ as defined in (2.8).

We will prove that $M_{n}$ is a martingale (see Lemma 3.4) and the asymptotic
behaviour of the martingale is strictly related with the sum:

$v_{n}=\displaystyle\sum_{k=1}^{n}\frac{1}{a_{k}^{2}}$ (3.2)

Now, note that by Stirling formula for the gamma function

$\begin{array}[]{cccc}a_{n}\thicksim\frac{n^{\alpha}}{\Gamma(\alpha+1)}&\text{
as }&n\rightarrow\infty\end{array}$ (3.3)

Therefore, in the diffusive region, where $0\leq\alpha<1/2$, we will prove in
Lemma 3.9 that $v_{n}$ diverges. In particular, by (2.8) we have:

$v_{n}=\displaystyle\sum_{k=1}^{n}\left(\frac{\Gamma(k)\Gamma(\alpha+1)}{\Gamma(k+\alpha)}\right)^{2}\sim(\Gamma(\alpha+1))^{2}\displaystyle\sum_{k=1}^{n}\frac{1}{k^{2\alpha}}$
(3.4)

as $n\rightarrow\infty$. Then, by the p-series we know that
$v_{n}\sim(\Gamma(\alpha+1))^{2}\cdot\frac{n^{1-2\alpha}}{1-2\alpha}$. Then,

$\lim\limits_{n\to\infty}\frac{v_{n}}{n^{1-2\alpha}}=\frac{(\Gamma(\alpha+1))^{2}}{1-2\alpha}$
(3.5)

In the critical region, where $\alpha=1/2$, we have
$v_{n}\sim(\Gamma(3/2))^{2}\log n$ and

$\begin{array}[]{ll}v_{n}&\sim(\Gamma(3/2))^{2}\displaystyle\sum_{k=1}^{n}\frac{1}{k}\end{array}$
(3.6)

Then $v_{n}$ diverges with velocity $\log n$ (see Lemma 3.9), by the harmonic
series, we obtain

$\lim\limits_{n\to\infty}\frac{v_{n}}{\log n}=\frac{\pi}{4}$ (3.7)

Finally, in the superdiffusive region, if $1/2<\alpha\leq 1$, $v_{n}$
converges to the finite value

$\begin{array}[]{ll}\displaystyle\lim_{n\to\infty}v_{n}&=\displaystyle\sum_{k=0}^{\infty}\left(\frac{\Gamma(k+1)\Gamma(\alpha+1)}{\Gamma(k+\alpha+1)}\right)^{2}=\displaystyle\sum_{k=0}^{\infty}\frac{(1)_{k}(1)_{k}(1)_{k}}{(\alpha+1)_{k}(\alpha+1)_{k}k!}\\\\[8.5359pt]
&={}_{3}F_{2}\left(\begin{array}[]{c|c}1,1,1&\\\
(\alpha+1),(\alpha+1)&\end{array}1\right)\end{array}$ (3.8)

where for $a\in\mathbb{R}$, $(a)_{k}=a(a+1)\cdots(a+k)$ for $k\geq 0$ and
$(a)_{0}=1$ represent the Pochhammer symbol and ${}_{3}F_{2}$ the
hypergeometric generalized function.

We also will use the following definitions.

######  Definition 3.2.

The predictable quadratic variation associated to $(M_{n})$ is given by
$\langle M\rangle_{0}=0$ and, for all $n\geq 1$

$\langle M\rangle_{n}=\displaystyle\sum_{k=1}^{n}{\mathbb{E}}[\Delta
M^{2}_{k}|\mathcal{F}_{k-1}].$ (3.9)

######  Definition 3.3.

Let denote $Z_{n}=\displaystyle\sum_{i=1}^{n}(X_{i})^{2}$,
$\forall~{}n=1,2,\ldots$, and
$\\{M^{*}_{n}\\}=\frac{Z_{n}-{\mathbb{E}}[Z_{n}]}{b_{n}}$ where $b_{n}$ is
given by $b_{1}=1$ and
$b_{n}:=\displaystyle\prod_{k=1}^{n-1}\left(1+\dfrac{\gamma}{k}\right)$ for
all $n\geq 2$, with $\gamma$ as in (2.7).

Now, we prove some lemmas useful in the proofs of the Theorems.

######  Lemma 3.4.

The sequences $\\{M_{n}\\}$ and $\\{M^{*}_{n}\\}$ are square integrable
martingales.

Proof First we prove that $\\{M_{n}\\}$ is a martingale, in fact from (2.6)

$\begin{array}[]{lll}{\mathbb{E}}(X_{n+1}|S_{n})&=&1\cdot{\mathbb{P}}(X_{n+1}=1|\mathcal{F}_{n})-1\cdot{\mathbb{P}}(X_{n+1}=-1|\mathcal{F}_{n})+0\cdot{\mathbb{P}}(X_{n+1}=0|\mathcal{F}_{n})\\\\[14.22636pt]
&=&\dfrac{(p-q)\cdot\theta}{n}\cdot S_{n}+(1-\theta)\cdot(p-q)\end{array}$
(3.10)

and

${\mathbb{E}}(S_{n+1}|S_{n})=\lambda_{n}S_{n}+\omega$ (3.11)

where $\lambda_{n}=1+\dfrac{\alpha}{n}$ with $\alpha=(p-q)\cdot\theta$ and
$\omega=(1-\theta)\cdot(p-q)$.

$\begin{array}[]{ll}{\mathbb{E}}(M_{n+1}|\mathcal{F}_{n})&={\mathbb{E}}\left[\frac{S_{n}+X_{n+1}-{\mathbb{E}}(S_{n}+X_{n+1})}{a_{n+1}}|\mathcal{F}_{n}\right]\\\\[8.5359pt]
&=\frac{(n+\alpha)(S_{n}-{\mathbb{E}}(S_{n}))}{n\cdot\left(1+\frac{\alpha}{n}\right)a_{n}}=M_{n}\end{array}$
(3.12)

Analogously, note that ${\mathbb{E}}(M^{*}_{n+1}|\mathcal{F}_{n})$ equals

$\begin{array}[]{ll}{\mathbb{E}}\left[\frac{Z_{n+1}-{\mathbb{E}}(Z_{n+1})}{b_{n+1}}|\mathcal{F}_{n}\right]&={\mathbb{E}}\left[\frac{Z_{n}+X^{2}_{n+1}-{\mathbb{E}}(Z_{n}+X^{2}_{n+1})}{b_{n+1}}|\mathcal{F}_{n}\right]\\\\[8.5359pt]
&=\frac{nZ_{n}+\gamma
Z_{n}-n{\mathbb{E}}(Z_{n})-\gamma{\mathbb{E}}(Z_{n})}{nb_{n+1}}=M^{*}_{n}\end{array}$
(3.13)

In addition, for any $n\geq 1$, $X_{n}$ are random variables with values on
$\\{+1,0,-1\\}$. Then, $|S_{n}|,|Z_{n}|\leq n$, which implies that $(M_{n})$
are $(M^{*}_{n})$ are locally square integrable martingales.

######  Lemma 3.5.

The martingale $(M_{n})$ can be written in the additive form

$M_{n}=\displaystyle\sum_{k=1}^{n-1}\Delta
M_{k}=\displaystyle\sum_{k=1}^{n-1}(M_{k}-M_{k-1})=\displaystyle\sum_{k=1}^{n-1}\frac{X_{k}-{\mathbb{E}}(X_{k}|\mathcal{F}_{k-1})}{a_{k}}$
(3.14)

Proof First we calculate the expectation

$\begin{array}[]{c}{\mathbb{E}}(S_{n}-\lambda_{n-1}S_{n-1})={\mathbb{E}}\left[S_{n-1}+X_{n}-\left(1+\frac{\alpha}{n-1}\right)S_{n-1}\right]\\\\[11.38092pt]
={\mathbb{E}}\left[{\mathbb{E}}\left[X_{n}|S_{n-1}\right]-\frac{\alpha}{n-1}S_{n-1}\right]=\omega\end{array}$
(3.15)

Now, the martingale difference $\Delta M_{n}=M_{n}-M_{n-1}$, using (3.10) and
(3.15):

$\begin{array}[]{ll}\Delta
M_{n}&=\frac{S_{n}-{\mathbb{E}}(S_{n})}{a_{n}}-\frac{S_{n-1}-{\mathbb{E}}(S_{n-1})}{a_{n-1}}=\frac{S_{n}-\lambda_{n-1}S_{n-1}-({\mathbb{E}}[S_{n}-\lambda_{n-1}S_{n-1}])}{a_{n}}\\\\[11.38092pt]
&=\frac{X_{n}-{\mathbb{E}}(X_{n}|\mathcal{F}_{n-1})}{a_{n}}\end{array}$

######  Remark 3.6.

Given that $X_{n}$ assumes values on $\\{-1,0,1\\}$ it is easy to see that

$\frac{|\Delta M_{n}|}{a_{n}}\leq\frac{2}{a_{n}}$ (3.16)

######  Lemma 3.7.

For $n=1,2,\ldots,$

${\mathbb{E}}[S_{n}]=\beta a_{n}+\omega\cdot
a_{n}\displaystyle\sum_{l=1}^{n-1}\frac{1}{a_{l+1}}$ (3.17)

and

${\mathbb{E}}[Z_{n}]=\psi b_{n}+\tau\cdot
b_{n}\displaystyle\sum_{l=1}^{n-1}\frac{1}{b_{l+1}}$

Proof
Note that: $x_{1}={\mathbb{E}}[S_{1}]={\mathbb{E}}[X_{1}]=p-q=\beta$, and
$x_{n+1}=E[S_{n+1}]={\mathbb{E}}[{\mathbb{E}}[S_{n+1}|S_{n}]]=\lambda_{n}{\mathbb{E}}(S_{n})+\omega$
by (3.11). Now, using $x_{1}=f_{0}=\beta$, $f_{n}=\omega,g_{n}=\lambda_{n}$
and $x_{n}={\mathbb{E}}(S_{n})$ for $n=1,2,\ldots,$ we have by Lemma A1 from
[21] that:

${\mathbb{E}}[S_{n}]=\beta\cdot\displaystyle\prod_{k=1}^{n-1}\lambda_{n}+\omega\displaystyle\sum_{l=1}^{n-1}\displaystyle\prod_{k=l+1}^{n-1}\lambda_{k}\\\
=\beta a_{n}+\omega\cdot a_{n}\displaystyle\sum_{l=1}^{n-1}\frac{1}{a_{l+1}}.$

In a similar form we obtain ${\mathbb{E}}[Z_{n}]$.

######  Lemma 3.8.

Let $a_{n}$ defined in (2.8), then
$\displaystyle\sum_{l=1}^{n-1}\frac{1}{a_{l+1}}\sim\frac{\Gamma(\alpha+1)n^{1-\alpha}}{(1-\alpha)}$.

Proof Note that

$\displaystyle\sum_{l=1}^{n-1}\frac{1}{a_{l+1}}=\displaystyle\sum_{l=2}^{n}\frac{1}{a_{l}}=\displaystyle\sum_{l=1}^{n}\frac{\Gamma(l)\Gamma(\alpha+1)}{\Gamma(l+\alpha)}-\frac{1}{a_{1}}=\Gamma(\alpha+1)\displaystyle\sum_{l=1}^{n}\frac{\Gamma(l)}{\Gamma(l+\alpha)}-1$

By using Lemma B.1 in [3], that is
$\sum_{l=1}^{n}\frac{\Gamma(l)}{\Gamma(l+\alpha)}=\frac{\Gamma(n+1)}{(\alpha-1)\Gamma(n+\alpha)}\cdot\left(\frac{\Gamma(n+\alpha)\Gamma(1)}{\Gamma(n+1)\Gamma(\alpha)}-1\right)$,
the equation (2.8) and taking $a=0$, $b=\alpha$, with $\alpha\neq 1$, it
holds:

$\displaystyle\sum_{l=1}^{n-1}\frac{1}{a_{l+1}}=\frac{n}{(1-\alpha)a_{n}}+\frac{1}{\alpha-1}\\\
\sim\frac{\Gamma(\alpha+1)n^{1-\alpha}}{(1-\alpha)}$

######  Lemma 3.9.

The series $\sum\frac{1}{a_{n}^{2}}$ converges, if and only if,
$\alpha>\frac{1}{2}$

Proof Let $C_{n}=\frac{1}{a_{n}^{2}}$ by the Raabe criteria, we need to
analyse $\lim\limits_{n\to\infty}n\left(1-\frac{C_{n+1}}{C_{n}}\right)$. Note
that,

$n\left(1-\frac{C_{n+1}}{C_{n}}\right)=n\left(1-\frac{1}{\lambda_{n}^{2}}\right)=n\left(\frac{\frac{2\alpha}{n}+\frac{\alpha^{2}}{n^{2}}}{1+\frac{2\alpha}{n}+\frac{\alpha^{2}}{n^{2}}}\right)=2\alpha$

The series converges if $\alpha>1/2$, and diverges if $\alpha<1/2$. In the
case $\alpha=1/2$, from (2.8) it follows from (3.3) that
$a_{n}\thicksim\frac{n^{1/2}}{\sqrt{\pi/4}}$. Since
$\frac{1}{a_{n}^{2}}=\frac{\pi/4}{n}\geq\frac{1}{2n}$, and the series
$\sum\frac{1}{a_{n}^{2}}$ diverges.

Finally, we state the following lemma without its proof. The details can be
obtained in Lemma 3.10 in [21], and the proof is based on Theorem 1 (b) from
[19].

######  Lemma 3.10.

Suppose that $\\{M_{n}\\}$ is a square-integrable martingale with mean $0$.
Let $\Delta M_{k}=M_{k}-M_{k-1}$, for $k=1,2,\ldots$, where $M_{0}=0$. If

$\displaystyle\sum_{k=1}^{\infty}{\mathbb{E}}[(\Delta M_{k})^{2}]<+\infty$

holds in addition, then we have the following: let
$r_{n}^{2}=\displaystyle\sum_{k=n}^{\infty}{\mathbb{E}}[(\Delta M_{k})^{2}]$

  * (i)

The limit $M_{\infty}:=\sum_{k=1}^{\infty}\Delta M_{k}$ exists almost surely
and $M_{n}\overset{L^{2}}{\rightarrow}M_{\infty}$

  * (ii)

Assume that

    * a)

$\displaystyle\frac{1}{r_{n}^{2}}\sum_{k=n}^{\infty}(\Delta M_{k})^{2}\to 1$
as $n\to\infty$ in probability, and

    * b)

$\displaystyle\lim_{n\to\infty}\frac{1}{r_{n}^{2}}{\mathbb{E}}\left[\sup_{k\geq
n}\ (\Delta M_{k})^{2}\right]=0$. Then, we have

$\displaystyle\frac{M_{\infty}-M_{n}}{r_{n+1}}=\frac{\sum_{k=n+1}^{\infty}\Delta
M_{k}}{r_{n+1}}\overset{d}{\rightarrow}N(0,1).$

  * (iii)

Assume that the following three conditions hold

    * a’)

$\displaystyle\frac{1}{r_{n}^{2}}\sum_{k=n}^{\infty}(\Delta M_{k})^{2}\to 1$
as $n\to\infty$ a.s

    * c)

$\displaystyle\sum_{k=1}^{\infty}\frac{1}{r_{k}}{\mathbb{E}}[|\Delta
M_{k}|:|\Delta M_{k}|>\varepsilon r_{k}]<+\infty$ for any $\varepsilon>0$, and

    * d)

$\displaystyle\sum_{k=1}^{\infty}\frac{1}{r_{k}^{4}}{\mathbb{E}}[(\Delta
M_{k})^{4}:|\Delta M_{k}|>\delta r_{k}]<+\infty$ for some $\delta>0$.

Then,
$\displaystyle\limsup_{n\to\infty}\pm\frac{M_{\infty}-M_{n}}{\sqrt{2\cdot
r^{2}_{n+1}\log|\log r^{2}_{n+1}|}}=1$ a.s.

###  3.1 Proof of Theorem 1

First, note that

$\frac{a_{n}}{n}=\frac{1}{n}\displaystyle\prod_{k=1}^{n-1}\left(1+\dfrac{\alpha}{k}\right)=\displaystyle\prod_{k=1}^{n-1}\left(\dfrac{k+\alpha}{k+1}\right)=\displaystyle\prod_{k=1}^{n-1}\left(1-\dfrac{1-\alpha}{k+1}\right).$
(3.18)

Moreover $-1\leq\alpha\leq 1$, then $0\leq\dfrac{k+\alpha}{k+1}\leq 1$. In
addition, $\displaystyle\sum_{k=1}^{\infty}\frac{1-\alpha}{k+1}=\infty$
implies that $\lim\frac{a_{n}}{n}=0$. It follows that $(a_{n})_{n\geq 1}$ is
non increasing. Let define $N_{j}=\frac{a_{j}}{j}\cdot\Delta M_{j}$. By (3.16)
we conclude that $(N_{j})_{j\geq 1}$ is a sequence of martingales such that

$\displaystyle\sum_{j=1}^{\infty}{\mathbb{E}}[N_{j}^{2}|\mathcal{F}_{j-1}]\leq\displaystyle\sum_{j=1}^{\infty}\frac{2^{2}}{j^{2}}<\infty$

From Theorem 2.17 in [18], $\displaystyle\sum_{j=1}^{\infty}N_{j}$ converges
almost surely. Given that $\frac{n}{a_{n}}\rightarrow\infty$, by a direct
application of Kronecker’s lemma,

$\frac{a_{n}}{n}\displaystyle\sum_{j=1}^{n}\Delta
M_{j}=\frac{a_{n}}{n}M_{n}\overset{a.s}{\rightarrow}0$

Then, we showed (2.9). On the other hand, using (3.17), and $\alpha<1$, we
have that the first term is $o(n)$ as $n\rightarrow\infty$, and by (3.3) and
Lemma 3.8

$\omega\cdot
a_{n}\displaystyle\sum_{l=1}^{n-1}\frac{1}{a_{l+1}}\sim\omega\frac{n^{\alpha}}{\Gamma(\alpha+1)}\Gamma(\alpha+1)\frac{n^{1-\alpha}}{1-\alpha}=\frac{\omega}{1-\alpha}n$

###  3.2 Proof of Theorem 2

Let denote $\xi_{n}=X_{n}-{\mathbb{E}}(X_{n}|\mathcal{F}_{n-1})$, from (3.9)
we have ${\mathbb{E}}[\xi_{n}|\mathcal{F}_{n-1}]=0$ and

${\mathbb{E}}[(\xi_{n+1})^{2}|\mathcal{F}_{n}]={\mathbb{E}}[X_{n+1}^{2}|\mathcal{F}_{n}]-({\mathbb{E}}(X_{n+1}|\mathcal{F}_{n}))^{2},$
(3.19)

where

${\mathbb{E}}[(X_{n+1})^{2}|\mathcal{F}_{n}]=\dfrac{n_{+}+n_{-}}{n}\cdot(p+q)\cdot\theta+(1-\theta)\cdot(p+q).$
(3.20)

Let denote $Z_{n}=n_{+}+n_{-}$, $\gamma=(p+q)\theta$ and
$\tau=(1-\theta)\cdot(p+q)$. Then, using (3.10)

${\mathbb{E}}[(\xi_{n+1})^{2}|\mathcal{F}_{n}]=\frac{\gamma}{n}Z_{n}+\tau-\left(\frac{\alpha}{n}S_{n}+\omega\right)^{2}.$
(3.21)

Similarly,

$\begin{array}[]{l}{\mathbb{E}}[(\xi_{n+1})^{4}|\mathcal{F}_{n}]={\mathbb{E}}\left[\left(X_{n+1}-{\mathbb{E}}(X_{n+1}|\mathcal{F}_{n})\right)^{4}|\mathcal{F}_{n}\right]\\\\[14.22636pt]
={\mathbb{E}}\left[X^{4}_{n+1}-4X^{3}_{n+1}{\mathbb{E}}(X_{n+1}|\mathcal{F}_{n})+6X^{2}_{n+1}({\mathbb{E}}(X_{n+1}|\mathcal{F}_{n}))^{2}\right.\\\\[14.22636pt]
\left.-4X_{n+1}\left({\mathbb{E}}(X_{n+1}|\mathcal{F}_{n})\right)^{3}+\left({\mathbb{E}}(X_{n+1}|\mathcal{F}_{n})\right)^{4}|\mathcal{F}_{n}\right]\\\\[14.22636pt]
={\mathbb{E}}(X^{2}_{n+1}|\mathcal{F}_{n})-4({\mathbb{E}}(X_{n+1}|\mathcal{F}_{n}))^{2}+6{\mathbb{E}}(X^{2}_{n+1}|\mathcal{F}_{n})\left({\mathbb{E}}(X_{n+1}|\mathcal{F}_{n})\right)^{2}\\\\[14.22636pt]
-3\left({\mathbb{E}}(X_{n+1}|\mathcal{F}_{n})\right)^{4},\end{array}$ (3.22)

For any $n\geq 1$, we have that $X_{n}$ is a random variable which takes
values on $\\{+1,0,-1\\}$. In consequence,
${\mathbb{E}}[(\xi_{n+1})^{2}|\mathcal{F}_{n}]$ and
${\mathbb{E}}[(\xi_{n+1})^{4}|\mathcal{F}_{n}]$, are bounded and by (3.21),
(3.22) we have the almost sure limits

$\begin{array}[]{lll}\displaystyle\sup_{n\geq
0}~{}{\mathbb{E}}[\xi^{2}_{n+1}|\mathcal{F}_{n}]\leq 2&\text{ and
}&\displaystyle\sup_{n\geq
0}~{}{\mathbb{E}}[\xi^{4}_{n+1}|\mathcal{F}_{n}]\leq 16\\\ \end{array}$ (3.23)

We deduce from (3.9) and (3.21) that

$\begin{array}[]{ll}\langle
M\rangle_{n}&=\displaystyle\sum_{k=1}^{n}\frac{1}{a_{k}^{2}}\left[\frac{\gamma}{k}Z_{k}+\tau-\left(\frac{\alpha}{k}S_{k}+\omega\right)^{2}\right]\\\\[14.22636pt]
\end{array}$ (3.24)

Now, note that using (2.6) and Theorem 2.1 we obtain
${\mathbb{E}}(X_{n+1}|\mathcal{F}_{n})=\frac{\alpha}{n}S_{n}+\omega\rightarrow\alpha\frac{\omega}{1-\alpha}+\omega=\frac{\omega}{1-\alpha}$
Moreover, by (3.20) we have that
${\mathbb{E}}[X_{n+1}^{2}|\mathcal{F}_{n}]=\frac{\gamma}{n}Z_{n}+\tau$, in the
Lemma 3.4 we constructed a martingale for $Z_{n}$, and in similar form than
Theorem 2.1 it is possible to prove that:

$\lim\limits_{n\to\infty}\frac{Z_{n}-{\mathbb{E}}(Z_{n})}{n}=0\text{ and
}\lim\limits_{n\to\infty}\frac{{\mathbb{E}}(Z_{n})}{n}=\frac{\tau}{1-\gamma}$
(3.25)

Then,

${\mathbb{E}}[X_{n+1}^{2}|\mathcal{F}_{n}]=\frac{\gamma}{n}Z_{n}+\tau\sim\frac{\tau}{1-\gamma}$
(3.26)

Finally, from (3.19)

${\mathbb{E}}[(\xi_{n+1})^{2}|\mathcal{F}_{n}]\sim\frac{\tau}{1-\gamma}-\left(\frac{\omega}{1-\alpha}\right)^{2}=\phi$
(3.27)

as defined in (2.7). Given that $\Delta M_{k}=\frac{\xi_{k}}{a_{k}}$, the
expectation of the square of the martingale difference satisfies,

$\begin{array}[]{ll}{\mathbb{E}}[(\Delta
M_{k})^{2}|\mathcal{F}_{k-1}]&\sim\frac{\phi}{a_{k}^{2}}\end{array}$ (3.28)

Diffusive region, $0\leq\alpha<1/2$

We use central limit theorem for martingales as stated in Corollary 2.1.10 in
[12]. First, note that from (3.2) and using Lemma 3.9, we have that $v_{n}$
diverges in this region. Moreover, by (3.9) and (3.28) we obtain that

$\begin{array}[]{cc}\lim\limits_{n\to\infty}\frac{\langle
M\rangle_{n}}{v_{n}}=\phi&a.s\\\ \end{array}$ (3.29)

Now, we need to prove that $(M_{n})$ satisfies the Lindeberg’s condition, that
is, for any $\varepsilon>0$,

$\frac{1}{v_{n}}\displaystyle\sum_{k=1}^{n}{\mathbb{E}}\left[|\Delta
M_{k+1}|^{2}\mathbb{I}_{|\Delta
M_{k+1}|\geq\varepsilon\sqrt{v_{n}}}|\mathcal{F}_{n}\right]\overset{P}{\rightarrow}0,$
(3.30)

it can be seen, for instance, in Theorem 2.1.9 from [12]. Then using (3.23),
we have:

$\begin{array}[]{l}\displaystyle\frac{1}{v_{n}}\displaystyle\sum_{k=1}^{n}{\mathbb{E}}\left[|\Delta
M_{k+1}|^{2}\mathbb{I}_{|\Delta
M_{k+1}|\geq\varepsilon\sqrt{v_{n}}}|\mathcal{F}_{n}\right]\leq\frac{1}{\varepsilon^{2}v_{n}^{2}}\displaystyle\sum_{k=1}^{n}{\mathbb{E}}\left[|\Delta
M_{k+1}|^{4}|\mathcal{F}_{n}\right]\\\\[14.22636pt]
\leq\displaystyle\sup_{1\leq k\leq
n}{\mathbb{E}}\left[\xi^{4}_{k+1}|\mathcal{F}_{n}\right]\frac{1}{\varepsilon^{2}v_{n}^{2}}\displaystyle\sum_{k=1}^{n}\frac{1}{a_{k}^{4}}\leq\frac{16}{\varepsilon^{2}v_{n}^{2}}\displaystyle\sum_{k=1}^{n}\frac{1}{a_{k}^{4}}\\\\[14.22636pt]
\end{array}$

Now, note that $a_{n}^{-4}v_{n}^{-2}$ is equivalent to $(1-2\alpha)n^{-2}$ and
it is well known that $\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\frac{\pi^{2}}{6}$.
Therefore

$\displaystyle\sum_{n=1}^{\infty}\frac{1/a_{n}^{4}}{v_{n}^{2}}<\infty$ (3.31)

By the Kronecker’s lemma (see [12]), we obtain

$\lim\limits_{n\to\infty}\frac{1}{v_{n}^{2}}\displaystyle\sum_{k=1}^{n}\frac{1}{a_{k}^{4}}=0$
(3.32)

In this sense, by central limit theorem for martingales, we have

$\lim\limits_{n\to\infty}\frac{M_{n}}{\sqrt{v_{n}}}\overset{d}{\rightarrow}N(0,\phi)$
(3.33)

Since $M_{n}=\frac{S_{n}-{\mathbb{E}}(S_{n})}{a_{n}}$, by (2.8) and (3.4), the
equation (2.11) holds. Finally by (3.29) and given that $v_{n}/v_{n+1}\to 1$,
as $n\to\infty$, we use Lemma 3.5 from [20] to obtain,

$\displaystyle\limsup_{n\to\infty}\pm\frac{M_{n}}{\sqrt{2\phi
v_{n}\log\log(\phi v_{n})}}=1\ \text{ a.s. }$ (3.34)

Critical region, $\alpha=1/2$

In this region we proceed in a similar way than previous case. Note that
$v_{n}$ diverges and using Lemma 3.9, and (3.24) we obtain (3.29). Moreover
(3.31) holds since $a_{n}^{-4}v_{n}^{-2}$ is equivalent to $(n\log n)^{-2}$ by
(3.6) and it is well known that $\sum_{n=1}^{\infty}\frac{1}{(n\log
n)^{2}}<\infty$. We realized in the diffusive region that (3.31) is a
sufficient condition to have $(M_{n})$ satisfying the Lindeberg’s condition.
In consequence, we use the central limit theorem for martingales (3.33),
equation (3.1) and (3.6) to obtain

$\lim\limits_{n\to\infty}\frac{M_{n}}{\sqrt{v_{n}}}=\frac{S_{n}-{\mathbb{E}}(S_{n})}{\sqrt{n\log
n}}$

Then (2.13) holds. Moreover, to prove (2.14) we use Lemma 3.5 from [20], we
obtain the analogous of (3.34), given (3.6) we complete the proof.

###  3.3 Proof of Theorem 3

In superdiffusive region, $1/2<\alpha\leq 1$, we use Lemma 3.10. Note that by
(3.28) and the bounded convergence theorem,
$\displaystyle\sum_{k=1}^{\infty}{\mathbb{E}}[(\Delta
M_{k})^{2}]\sim\phi\Gamma(\alpha+1)^{2}\displaystyle\sum_{k=1}^{\infty}\frac{1}{k^{2\alpha}}$.
Then, if $\alpha>1/2$, we have
$\displaystyle\sum_{k=1}^{\infty}{\mathbb{E}}[(\Delta M_{k})^{2}]<\infty$.
Therefore, from item (i) of the cited Lemma, the following limit exists with
probability one,

$W:=\displaystyle\sum_{k=1}^{\infty}\Delta
M_{k}=\lim\limits_{n\to\infty}\frac{S_{n}-{\mathbb{E}}[S_{n}]}{a_{n}},$

and $M_{n}\overset{L^{2}}{\rightarrow}W$. Then ${\mathbb{E}}(W)=0$ and
${\mathbb{E}}(W^{2})=\sum_{k=1}^{\infty}{\mathbb{E}}[(\Delta M_{k})^{2}]>0$.

To prove (2.16) we use item ii) of Lemma 3.10, for that, we need to prove for
$\alpha>1/2$ that items a) and b) hold.
A sufficient condition for a) is:

$\frac{1}{r_{n}^{2}}\sum_{k=n}^{\infty}{\mathbb{E}}\left[(\Delta
M_{k})^{2}|\mathcal{F}_{k-1}\right]\rightarrow 1~{}\text{ as
}~{}n\to\infty~{}\text{ in probability}$ (3.35)

by Corollary 1 item b‘) in [19], where
$r_{n}^{2}=\sum_{k=n}^{\infty}{\mathbb{E}}[(\Delta M_{k})^{2}]$ and given
(3.28), we have:

$\begin{array}[]{ll}\displaystyle\sum_{k=n}^{\infty}{\mathbb{E}}[(\Delta
M_{k})^{2}|\mathcal{F}_{k-1}]&=\displaystyle\sum_{k=n}^{\infty}\frac{\phi}{a_{k}^{2}}\sim\phi\Gamma(\alpha+1)^{2}\displaystyle\sum_{k=n}^{\infty}\frac{1}{k^{2\alpha}}\\\\[11.38092pt]
&\sim\displaystyle\frac{\phi\Gamma(\alpha+1)^{2}}{(2\alpha-1)n^{2\alpha-1}}\sim\frac{\phi
n}{(2\alpha-1)a_{n}^{2}}~{}a.s\end{array}$

Using the bounded convergence theorem, we obtain that:

$\sum_{k=n}^{\infty}{\mathbb{E}}[(\Delta
M_{k})^{2}]=\sum_{k=n}^{\infty}\frac{\phi}{a_{k}^{2}}\sim\frac{\phi
n}{(2\alpha-1)a_{n}^{2}}~{}a.s$

Then, as $n\to\infty$, we have:

$\displaystyle\lim_{n\to\infty}\frac{1}{r_{n}^{2}}\sum_{k=n}^{\infty}{\mathbb{E}}\left[(\Delta
M_{k})^{2}|\mathcal{F}_{k-1}\right]=1~{}a.s$

This shows that conditions a) and a’) in Lemma 3.10 are satisfied. Moreover

$r_{n}^{4}\sim\frac{\Gamma(\alpha+1)^{4}\phi^{2}}{(2\alpha-1)^{2}n^{4\alpha-2}}\text{
and }(\Delta
M_{n})^{4}\leq\frac{16}{(a_{n}^{4})}\sim\frac{16\Gamma(\alpha+1)^{4}}{n^{4\alpha}}.$

Now, the sufficient condition to have b) is the Remark A1 in [21], in this
sense, using Theorem 2.1.9 from [12] and (3.23):

$\begin{array}[]{ll}\frac{1}{r_{n}^{2}}\displaystyle\sum_{k=n}^{\infty}{\mathbb{E}}\left[(\Delta
M_{k+1})^{2}:|\Delta M_{k+1}|\geq\varepsilon
r_{k}\right]&\leq\frac{1}{\varepsilon^{2}r_{n}^{4}}\displaystyle\sum_{k=n}^{\infty}{\mathbb{E}}\left[|\Delta
M_{k+1}|^{4}\right]\\\\[11.38092pt]
&\leq\displaystyle\frac{16\Gamma(\alpha+1)^{4}}{\varepsilon^{2}r_{n}^{4}}\sum_{k=n}^{\infty}\frac{1}{a_{k}^{4}}\leq\frac{C_{1}}{r_{n}^{4}}\displaystyle\sum_{k=n}^{\infty}\frac{1}{a_{k}^{4}}\end{array}$

At this point, we need to prove that

$\lim\limits_{n\to\infty}\frac{1}{r_{n}^{4}}\displaystyle\sum_{k=n}^{\infty}\frac{1}{a_{k}^{4}}=0$
(3.36)

which clearly holds since $a_{n}^{-4}r_{n}^{-4}$ is equivalent to
$\frac{1}{a_{n}^{2}\phi}$ and by Lemma 3.9, in the superdiffusive case
$\displaystyle\sum_{n=1}^{\infty}\frac{1}{a_{n}^{2}}<\infty$. Therefore,
$\displaystyle\sum_{n=1}^{\infty}\frac{1/a_{n}^{4}}{r_{n}^{4}}<\infty$ and by
Lemma 1 item ii) in [19], (3.36) holds. Then, the proof finishes by using the
Theorem 1 (b) in [19]

$\frac{W-M_{n}}{r_{n+1}}\overset{d}{\rightarrow}N(0,1)$

Finally, for $\varepsilon>0$ we have that:

$\frac{1}{r_{k}}{\mathbb{E}}\left[|\Delta M_{k+1}|:|\Delta
M_{k+1}|\geq\varepsilon
r_{k}\right]\leq\frac{1}{r_{k}}\frac{1}{\varepsilon^{3}r_{n}^{3}}{\mathbb{E}}\left[|\Delta
M_{k+1}|^{4}\right]\leq C_{3}k^{4\alpha-2}k^{4\alpha}\leq\frac{C_{3}}{k^{2}}$

which implies that c) holds. In the case of d), it is enough to note that
$\sum_{n=1}^{\infty}\frac{1}{r_{k}^{4}}{\mathbb{E}}[(\Delta
M_{k})^{4}]<\infty$. Then, (2.17) holds.

## Acknowledgements

The authors thank the anonymous referees for their valuable comments. Also
several corrections by Christian Caamaño and Rodrigo Lambert. This work was
partially supported by Fondecyt nº 11200500

## References

  * [1] Baur, E. and Bertoin, J. (2016)  Elephant random walks and their connection to Pólya-type urns. Phys. Rev. E 94, 052134.
  * [2] Baur, E. (2020)  On a class of random walks with reinforced memory. J. Stat. Phys.
  * [3] Bercu, B. (2018) A martingale approach for the elephant random walk. J. Phys. A: Math. Theor. 51(1).
  * [4] Bercu, B. and Laulin, L. (2019)  On the multi-dimensional elephant random walk. J. Stat. Phys. 175, 1146-1163.
  * [5] Bercu, B. and Laulin, L. (2021) On the center of mass of the elephant random walk. Stochastic Process. Appl. 133, 111-128.
  * [6] Bercu, B., Chabanol, M.-L. and Ruch, J.-J. (2019) Hypergeometric identities arising from the elephant random walk. J. Math. Anal. Appl. 480(1), 123360, 12.
  * [7] Bertoin, J. (2020) Noise reinforcement for Levy processes. Ann. Inst. H. Poincaré Probab. Statist. 56(3), 2236-2252.
  * [8] Bertoin, J. (2020)  Universality of Noise Reinforced Brownian Motions. In and Out of Equilibrium 3: Celebrating Vladas Sidoravicius, pp 147-161.
  * [9] Bertoin, J. (2020)  How linear reinforcement affects Donsker’s Theorem for empirical processes. Probab. Theory Relat. Fields 178, 1173-1192.
  * [10] Businger, S. (2018) The Shark Random Swim. J. Stat. Phys. 172(3), 701-717.
  * [11] Coletti, C., Gava, R. and Schütz, G. (2017) Central limit theorem for the elephant random walk. J. Math. Phys. 56, 05330.
  * [12] Duflo, M. (1997) Random Iterative Models. Applications of Mathematics, vol. 34. Springer, Berlin.
  * [13] González-Navarrete, M. (2020) Multidimensional walks with random tendency. J. Stat. Phys. 181(3): 1138-1148.
  * [14] González-Navarrete, M. and Lambert, R. (2018) Non-Markovian random walks with memory lapses. J. Math. Phys. 59, 113301.
  * [15] González-Navarrete, M. and Lambert, R. (2019) The diffusion of opposite opinions in a randomly biased environment. J. Math. Phys. 60, 113301.
  * [16] Gut, A. and Stadtmüller, U. Elephant random walks with delays. Preprint arXiv:1906.04930v1, 2019.
  * [17] Gut, A. and Stadtmüller, U. (2021) The number of zeroes in Elephant random walks with delays. Statist. Probab. Lett. 174, 109112.
  * [18] Hall, P. and Heyde, C.C. (1980) Martingale Limit Theory and Its Application. Academic Press, New York.
  * [19] Heyde, C.C. (1977)  On central limit and iterated logarithm supplements to the martingale convergence theorem. J. Appl. Probab. 14(4), 758-775.
  * [20] James, B., James, K. and Qi, Y. (2008)  Limit theorems for correlated Bernoulli random variables. Stat. Probab. Lett. 78, 2339–2345.
  * [21] Kubota, N. and Takei, M. (2019).  Gaussian fluctuation for superdiffusive elephant random walks. J. Stat. Phys. 177, 1157-1171.
  * [22] Schütz, G. and Trimper, S. (2004) Elephants can always remember: Exact long-range memory effects in a non-Markovian random walk. Phys. Rev. E 70, 045101.

Departamento de Estadística, Universidad del Bío-Bío. Avda. Collao 1202, CP
4051381, Concepción, Chile. E-mail address<EMAIL_ADDRESS>


# Online Trajectory Optimization for Persistent Monitoring Problems in
Partitioned Environments

Systems Engineering
Boston University
Massachusetts, USA  Jonas Hall1, Christos G. Cassandras1,2, Sean B.
Andersson1,3 This work was supported in part by the NSF through
ECCS-19316001Division of Systems Engineering, Boston University,
USA2Department of Electrical and Computer Engineering3Division of Mechanical
Engineering, Boston University, USA{hallj, cgc<EMAIL_ADDRESS>
###### Abstract

We consider the problem of using an autonomous agent to persistently monitor a
collection of dynamic targets distributed in an environment. We generalize
existing work by allowing the agent’s dynamics to vary throughout the
environment, leading to a hybrid dynamical system. This introduces an
additional layer of complexity towards the planning portion of the problem: we
must not only identify in which order to visit the points of interest, but
also in which order to traverse the regions. We design an offline high-level
sequence planner together with an online trajectory optimizer realizing the
computed visiting sequence. We provide numerical experiments to illustrate the
performance of our approach.

##  I Introduction

Persistent Monitoring (PM) describes a broad class of problems in which an
agent moves through an environment to collect information about specific
targets over time. It is applicable across a wide range of applications, such
as ocean monitoring [1], forest fire surveillance [2], tracking of individual
biological macromolecules [3], and data harvesting [4]. As a specific
motivating application, consider a disaster scenario. An efficient response to
a catastrophic event such as an earthquake, hurricane, or tsunami requires the
persistent and simultaneous state estimation of many locations in order to
make time-sensitive resource distribution decisions. The affected area may
consist of various types of terrain, each having its own characteristics. For
instance, some parts may be urban and thus contain obstacles, others may be
coastal regions with the presence of strong winds. In this paper, we propose
an online trajectory optimization scheme to minimize the average estimation
error of the site states for piecewise continuous agent dynamics. These hybrid
dynamics have the capability of capturing the various terrains the agent must
move through; in turn, this requires the extension of standard trajectory
optimizers for PM in order to exploit the local structures.

In this paper, we assume that there are multiple sites to be monitored, each
of which is located in a different region. We model the state of each site
with a stochastic Linear Time-Invariant (LTI) system and assume the agent has
the ability to interact with each site in its vicinity through noisy
measurements of this state. To obtain good estimates of the state of site, the
agent must navigate to the site and remain nearby to collect measurements over
time. However, since the agent is tasked with monitoring all the sites, it
cannot spend too much time at each individual site since its knowledge of
other sites not being monitored decays. Thus, we want to optimize both the
time spent at each location and the sequence of site visits.

The optimization of the visiting sequence is often approached by abstracting
the problem to a graph structure and using Traveling Salesperson Problem (TSP)
inspired solutions [5] in an offline preparation phase. The generation of such
a graph requires the estimation of travel times between sites, which is itself
a challenging problem in our setting of hybrid agent dynamics. One conceivable
approach to tackle this problem is the utilization of nonsmooth optimal
control solvers [6, 7]. However, since they are numerical solvers, they are
prone to getting stuck at local minima. Furthermore, they cannot exploit any
structure that is specific to each terrain. Since the visiting sequence can be
found offline, we do not require solutions in real-time and introduce an
algorithm inspired by Rapidly-exploring Random Tree (RRT*) [8]. The proposed
algorithm makes the assumption that regional path planners exist for each
terrain type, and then generates a tree of how to reach a site by sampling
nodes on the boundaries of the region partitioning. The usage of RRT* is not
novel within PM, e.g., the authors in [9] have proposed a rapidly-exploring
random cycle generation algorithm to tackle the entire trajectory planning
problem. The search space of all cycles, however, can become extremely large
for PM problems. Here, we only use it for the abstract graph generation to
plan an appropriate visiting sequence. After determining a visiting sequence,
we complete the trajectory by solving local monitoring Optimal Control
Problems (OCP) for each site visit. Finally, we optimize the monitoring
durations in an online bilevel optimization scheme.

The bilevel optimization approach builds upon our prior work in [10]. The
present paper differs from that effort in many key aspects, including a more
realistic model of target uncertainty (described now in terms of variance of
the state estimates), a more complex environment, and a focus on minimizing
the average uncertainty rather than the period of the trajectory subject to
constraints on the uncertainty reduction in each cycle. The contributions of
this paper are as follows:

  * •

We propose a novel PM formulation by introducing piecewise continuous agent
dynamics.

  * •

We introduce the Rapidly-exploring Random Boundary Tree (RRBT*) algorithm for
global path planning in the introduced context.

  * •

We introduce an online optimization scheme to minimize the steady-state
average estimation error.

The remainder of this paper is organized as follows. In Sec. II we introduce
the problem formulation. Sec. III introduces the RRBT* algorithm, which lays
the basis for computing the interregional trajectory segments. We then
generate a TSP solution from the graph generated by successive utilization of
the RRBT* algorithm. Using the computed visiting sequence, we discuss the
trajectory decomposition into local control segments in Sec. IV. Sec. IV-C
then provides the online bilevel optimization of the global cost function by
alternating between updates of the local control segments and the global
optimization parameters. The efficacy of the method is demonstrated in Sec. V
and we conclude the article in Sec. VI while providing open questions for
future work.

##  II Problem Formulation

As depicted in Fig. 1, we consider a partitioned environment in which a single
agent moves to carry out its PM mission. The partition is defined by the
agent’s dynamics: these dynamics are smooth within any given region but may
change abruptly when moving between regions. Each region is assumed to contain
at most one target; each target is associated with a set of states that evolve
according to stochastic LTI dynamics. The agent is equipped with a sensor that
can take noisy measurements of a target’s internal states whenever it is
within the respective region. In this section we introduce the necessary
notation and formally define our PM problem.

For a vector $x\in\mathbb{R}^{n}$ we denote its Euclidean norm by $\|x\|$. For
a given vector $x$ and a nonnegative scalar $\alpha\geq 0$ we denote the set
of all vectors $v\in\mathbb{R}^{n}$ that satisfy $\|x-v\|\leq\alpha$ by
$B_{\alpha}(x)$. For a set $S\subset\mathbb{R}^{n}$ we denote its interior by
$\stackrel{{\scriptstyle\circ}}{{S}}$, it’s closure by $\bar{S}$, and its
boundary by $\partial
S=\bar{S}~{}\setminus\stackrel{{\scriptstyle\circ}}{{S}}$. By
$\mathcal{B}(\mathbb{R}^{n})$ we denote the set of all subsets of
$\mathbb{R}^{n}$. Given two paths
$u_{1}:[0,T_{1}]\to\mathbb{R}^{n},u_{2}:[0,T_{2}]\to\mathbb{R}^{n}$, we define
their concatenation $u_{1}\circ u_{2}:[0,T_{1}+T_{2}]\to\mathbb{R}^{n}$ as
$(u_{2}\circ u_{1})(t)=\alpha(t)u_{1}(t)+(1-\alpha(t))u_{2}(t)$, where
$\alpha(t)\in\\{0,1\\}$ is equal to $1$ iff $t\in[0,T_{1}]$.

Figure 1: Illustrating the partitioned mission space based on the piecewise
smooth (uncontrolled) agent dynamics together with a comparison of periodic
trajectories.

Let us assume that the compact and connected mission space
$S\subset\mathbb{R}^{2}$ is partitioned into $P$ polyhedra,

$S_{i}=\\{x\in\mathbb{R}^{2}\mid g_{j}^{\top}x\leq
b_{j}\mathrm{~{}for~{}}j=1,2,\dots,P_{i}\\},$ (1)

for $i=1,2,\dots,P$, where $P_{i}$ denotes the number of inequalities required
to describe the set $S_{i}$. We assume that

$\stackrel{{\scriptstyle\circ}}{{S}}_{i}\cap\stackrel{{\scriptstyle\circ}}{{S}}_{j}=\emptyset\mathrm{~{}for~{}}i\neq
j\quad\mathrm{~{}and~{}}\quad\bigcup_{i=1}^{P}S_{i}=S,$

i.e., we assume that any two sets in the the partition have disjoint interiors
and that the union of all sets provides the entire space. The partitioning is
illustrated in Fig. 1.

Now, let $a(t)\in\mathbb{R}^{2}$ denote the agent’s position at time $t$. To
capture the changing dynamics of the agent, let us define the piecewise smooth
vector field

$f(x)=f_{i}(x)\mathrm{~{}if~{}}x\in\ \stackrel{{\scriptstyle\circ}}{{S}}_{i},$
(2)

where each $f_{i}$ is assumed to be smooth. To extend the dynamics on the
switching set $\mathcal{S}_{f}=\cup_{i}\partial S_{i}$ of the partition, we
use the notion of a Filippov set-valued map:

$\begin{split}\mathcal{F}[f](x)&=\bar{co}\left\\{\lim_{i\to\infty}f(x_{i})~{}:~{}x_{i}\to
x,~{}x_{i}\notin S_{f}\right\\}.\end{split}$

Note that in the interior of a region $i$ the set-valued map is the singleton
set $\\{f_{i}(x)\\}$ due to the continuity of $f$, which shows its consistency
with the original vector field (2). Using this notion, we define the agent’s
dynamics as

$\dot{a}(t)\in\mathcal{F}[f](a(t))+u(t),$

where $u:[0,\infty)\to B_{1}(0)$ is assumed to be an integrable control
function. For more details on Filippov differential inclusions, we refer the
reader to [11, 12].

We further consider $M$ targets located at $p_{i}\in S_{j_{i}}$ in some region
$j_{i}$. We assume that each region contains at most one target, i.e.,
$j_{i_{1}}\neq j_{i_{2}}$ for $i_{1}\neq
i_{2}\in\mathcal{T}=\\{1,2,\dots,M\\}$. This assumption is not restrictive
since we could introduce an additional partition with identical dynamics
separating the targets. Each target contains an internal state
$\phi_{i}(t)\in\mathbb{R}^{m_{i}}$ following the stochastic LTI dynamics

$\dot{\phi}_{i}(t)=A_{i}\phi_{i}(t)+\omega_{i}(t),$ (3)

where $\omega_{i}$ is a zero-mean white noise process with covariance
$Q_{i}\succ 0$.

The agent is tasked to maintain estimates of all targets’ internal states. The
only form of interaction between the agent and target $i$ is described by the
sensor’s measurement model

$z_{i}(a(t),\phi_{i}(t))=\gamma_{i}(a(t))H_{i}\phi_{i}(t)+\nu_{i}(t).$ (4)

Here $\nu_{i}$ is a zero mean white noise process with covariance matrix
$R_{i}$ and $\gamma_{i}:\mathbb{R}^{2}\to[0,1]$ captures the quality of the
measurement depending on the position of the agent. The precise form of this
function is not crucial, however, we assume it to (i) be differentiable with
respect to the agent position within the region $S_{i}$ and (ii) have compact
support $K_{i}\subseteq S_{i}$. We can then write

$\gamma_{i}(a(t))=\begin{cases}0,&\mathrm{~{}if~{}}a(t)\notin\
\stackrel{{\scriptstyle\circ}}{{S}}_{i},\\\
\mathring{\gamma}_{i}(a(t))&\mathrm{~{}if~{}}a(t)\in\
\stackrel{{\scriptstyle\circ}}{{S}}_{i},\end{cases}$ (5)

where the $\mathring{\gamma}_{i}:\mathbb{R}^{2}\to[0,1]$ is a differentiable
function. Fig. 1 illustrates different potential forms of this function, e.g.,
the region of target 3 has a nonconvex quality mapping; whereas the other
regions have Gaussian quality mappings.

From previous results [13] we know that the optimal unbiased estimator of
$\phi_{i}$ using the measurements from (4) is the Kalman-Bucy filter. The
covariance matrix of this estimator is given by

$\begin{split}\dot{\Omega}_{i}(t)&=A_{i}\Omega_{i}(t)+\Omega_{i}(t)A_{i}^{\top}+Q_{i}-\Omega_{i}(t)G_{i}(t)\Omega_{i}(t)\\\
&=\vcentcolon f_{\Omega_{i}}(t),\end{split}$ (6)

where

$G_{i}(t)=\gamma_{i}^{2}(a(t))H_{i}^{\top}R_{i}^{-1}H_{i}.$

The final ingredient required to formulate the introduced problem as an OCP is
the choice of a cost function. Some examples of objectives typical to
persistent monitoring applications are the minimization of the (i) maximum
revisit time [5]; (ii) worst case estimation error [14, 15]; or (iii) average
estimation error [16]. In this paper we choose to minimize the average
estimation error. Combining all ingredients allows us to formulate the
following OCP. ulim_T →∞ 1T∫_0^T ∑_i=1^M tr(Ω_i(t)) dt a(0)= ¯a Ω(0)= ¯Ω ˙a(t)
∈F[f](a(t)) + u(t) ˙Ω_i(t)= f_Ω_i(t) for all i ∥ u(t) ∥≤1, where
$\bar{a},~{}\bar{\Omega}$ fix the initial states.

######  Assumptions.

Throughout this paper we assume that

  * •

the pair $(A_{i},H_{i})$ is observable for all $i\in\mathcal{T}$; and

  * •

$Q_{i}$ and $\Omega_{i}(0)$ are positive definite for all $i\in\mathcal{T}$.

These assumptions ensure that the covariance matrix $\Omega_{i}(t)$ is
positive definite for all $t\geq 0$ and that the covariance matrix
trajectories eventually become periodic under any periodic control law [17].
It is well known that periodic trajectories have the capability of
approximating the optimal cost of average infinite horizon problems
arbitrarily well [18], and we thus restrict our attention to periodic control
laws. This allows the reformulation of the infinite horizon OCP (1) to a
finite horizon problem by optimizing the steady-state average cost

$\begin{split}J=\frac{1}{T}\int_{0}^{T}\sum_{i=1}^{M}\mathrm{tr}(\Omega_{i}(t))~{}\textup{d}t,\end{split}$
(7)

over a cycle of period $T$. The following section introduces the offline
global planner, which we utilize in order to determine a periodic visiting
sequence.

##  III Global Planning

As illustrated by the yellow dashed curve in Fig. 1, an optimal trajectory is
comprised of (i) the sequence in which the targets are visited; (ii) the
points on the partition set boundaries where the agent changes its mode; (iii)
the trajectories connecting the points from (ii) of one target region to the
next; and (iv) the local monitoring segments within each target region. Note
that since the agent does not sense any of the targets along trajectory
segments in (iii), these trajectories must be time-optimal for the overall
trajectory to be optimal. Optimizing (i)-(iii) can be done offline; in this
section we develop heuristics for this. Sec. IV then introduces the
optimization of (iv), which is done online.

The goal of this section is to find a periodic target visiting sequence
$(i_{1},i_{2},\dots,i_{K})$, where $K$ is the length of that sequence,
together with the time optimal switching segments between the target regions
$i_{k}$ to $i_{k+1}$. We do this by abstracting the mission space to a graph,
where each node represents a target and each directed edge denotes the minimum
time required to travel from one node to the other. Given such a graph we seek
a minimum time cycle that visits all targets, i.e., we want to solve a TSP
problem, which is NP-hard in and of itself. Such a path is illustrated by the
(directed) red line in Fig. 1. Depending on the number of targets we can use
an exact solver or resort to heuristics.

Figure 2: Illustrating the used global path planner based on an RRBT*
algorithm and we use it to find a global path from the initial node (green
circle) to the root node (yellow star).

In order to generate the graph introduced above, we must compute the minimum
time required to travel between any pair of targets. To do this, we created an
algorithm inspired by RRT* [8] in order to generate a global search tree
capturing minimum-time paths to each target $i_{\ell}$ from anywhere in the
mission space. The main difference from the formulation in [8] is the presence
of the hybrid dynamics. We exploit the partition structure by assuming the
existence of a local controller for each region $i$ with the capability of
computing a (time optimal) control law that steers the agent from $x$ to $y$
within region $S_{i}$. We denote the time required to do so by

$\delta_{i}(x,y):S_{i}\times S_{i}\to\mathbb{R}_{\geq 0}\cup\\{\infty\\}.$ (8)

Note that while we assume $\delta$ to be time optimal, the rest of our
approach does not strictly require it; the resulting paths would, of course,
be suboptimal. With this abstraction, we can modify the original RRT*
algorithm to produce a search tree on which all nodes, apart from the root
node, lie on the boundary sets of the partition. Below we describe the RRBT*
algorithm in more detail. A brief pseudo-code overview is provided in
Algorithm 1. In Fig. 2 we illustrate an example of a search tree generated for
the root node (yellow) together with a particular path from another node
(green).

Given a point $x$ in the mission space, we initialize the tree of the RRBT*
algorithm by creating a root node at $x$, and we initially _activate_ all
regions which contain the root node. Active regions are all regions that have
a representative within the tree; this is a necessary condition that
guarantees a newly sampled node from that region can be connected to the tree
by using the respective local controller (8). Note that only one region is
initially activated if the root lies in the interior of a partition set.
Otherwise multiple regions are activated, since each set of the partition is
closed (1).

Let us now walk through a main iterate. The sampling phase aims to generate a
new node to be included into the tree. We do this by first uniformly sampling
a region $r$ from the active regions. We then uniformly sample a boundary
point $b$ from that region $r$, which provides the new node. The neighborhood
of a boundary point $b$, denoted $N(b)$, is defined to be the set of all nodes
in the tree that share a region with the sampled boundary point. Note that
$N(b)$ is not empty by design. Given a potential parent $p\in N(b)$, we can
compute an upper bound for the cost-to-root $c(b,x)$ by passing through $p$,
i.e., $c(b,x)\leq\delta_{r}(b,p)+c(p,x)$. We then connect the newly sampled
node to the parent that minimizes this cost-to-root metric. We repeat this
procedure for a predefined number of iterations.

We remark here that if all points in the space are reachable from $x$ then all
regions will eventually become active if we sample long enough, and we can
then connect each target to any other target, which is how we plan to use this
algorithm. To see this, note that each region only has a finite number
boundary segments, each of which has a positive length. Thus, the probability
of sampling a point from a specific boundary segment has positive probability.
If we sample sufficiently often from one region, every boundary segment will
eventually have a representative and all neighbors will become active. Now let
$r_{A}$ be some active region and $r_{1}$ an arbitrary region. Since the
regions are given by a partitioning of a connected space, there is a finite
sequence of regions $(r_{1},r_{2},\dots,r_{A})$, where each $r_{i}$ is a
neighbor of $r_{i\pm 1}$. Since all neighboring regions of $r_{A}$ will
eventually become active, $r_{A-1}$ will also become active, and thus also
$r_{A-2}$, and so forth, until eventually $r_{1}$ will become active.

Once we have generated a search tree from a point $x$, we can find a global
trajectory from any $y$ that lies in an active region by performing the
connection phase of Algorithm 1 with $b=y$. We use this mechanism to generate
a tree for each target and then compute the target-to-target distance for each
combination.

1

2# initialization phase

3 initialize tree with root at $x$

4 activate all regions containing $x$

5[0.5em]

6# main loop

7 while  _max iter not exceeded_ do

8 # sample phase

9 sample $r$ from all active regions

10 sample boundary point $b$ from $r$

11 [0.5em]

12 # connection phase

13 get neighborhood $N(b)$

14 find parent $p\in N(b)$ minimizing cost-to-root

15 connect $b$ to $p$

16 [0.5em]

17 # activation phase

18 activate any new region containing $b$

Algorithm 1 RRBT* (from root node $x$)

##  IV Optimizing the Average Cycle Cost

The algorithm in the previous section generates a periodic visiting sequence
$(i_{1},i_{2},\dots,i_{K})$ together with an initial trajectory of some period
$T$. Using this initialization, we decompose the control law
$u:[0,T]\to\mathcal{U}$ into its local segments

$u=(u^{\mathrm{s}}_{K-1}\circ
u^{\mathrm{m}}_{K-1})\circ\dots(u^{\mathrm{s}}_{1}\circ u^{\mathrm{m}}_{1}),$

where $u^{\mathrm{s}}_{k}$ describes the $k$th switching control law from
target region $i_{k}$ to $i_{k+1}$, which is obtained using the optimal path
generated using the RRBT* algorithm, and $u^{\mathrm{m}}_{k}$ describes the
$k$th monitoring control law. Fig. 3 depicts this decomposition. The next
section focuses on the optimization of the monitoring control laws.

###  IV-A The Local Monitoring Control Problem

The boundary conditions of the monitoring segment are given by the start and
end points of the next and previous switching segments, respectively. We
denote the agent’s entry point of the $k$th monitoring segment by
$a^{\varphi}_{k}\in\partial S_{i_{k}}$, the exit point by
$a^{\psi}_{k}\in\partial S_{i_{k}}$, the associated duration by
$\tau_{k}\in\mathbb{R}_{\geq 0}$, and the initial values of the estimator’s
covariance matrices by $\bar{\Omega}$. We then define the $k$th monitoring
OCP: τ, u^m_k∫_0^τ ∑_i=1^M tr(Ω_i(t)) dtM^∗_k(τ_k) = τ= τ_k a(0) = a^φ_k
a(τ_k) = a^ψ_k Ω(0)= ¯Ω ˙a(t) = f_i_k(x) + u^m_k(t) ∥ u^m_k(t) ∥≤1 ˙Ω_i(t)=
f_Ω_i(t) for all i a(t) ∈S_i_k. Let us remark here that the dynamics (IV-A)
and (IV-A) are smooth, as long as the resulting trajectory remains within the
region $i_{k}$ for all $t\in(0,\tau_{k})$, which is enforced by the linear
path constraints (IV-A). Note that the local monitoring OCP only affects the
estimator’s covariance matrix corresponding to the current target $i_{k}$
since the sensing quality function for all other targets are zero (see (5)).
In view of optimizing (IV-A), removing the states of the targets that are not
monitored from (IV-A) and their associated cost terms leads to the same
optimizer. However, in view of the decomposition described in the subsection
below, they must be included. Whether these terms are included before or after
the optimization has no effect. We include the terms here in order to be
consistent with the presentation later on.

Fixing the local monitoring durations $\tau_{k}$ allows us to solve the OCP in
(IV-A). Analytical solutions may exist for some very particular scenarios, but
in order to achieve flexibility with respect to $f_{i_{k}}$ and the choice of
the sensing quality function we solve this OCP using numerical optimization
techniques [19]. Crucially, in order for the local monitoring problem (IV-A)
to be consistent, we require that the exit point $a^{\psi}_{k}$ is reachable
for the agent from the entry point $a^{\varphi}_{k}$ within $\tau_{k}$ units
of time, i.e., we must ensure

$\tau_{k}\geq\delta_{i_{k}}(a^{\varphi}_{k},a^{\psi}_{k}),$ (9)

where we recall that $\delta_{i_{k}}(a^{\varphi}_{k},a^{\psi}_{k})$ is the
(minimum) time required to travel from $a^{\varphi}_{k}$ to $a^{\psi}_{k}$ in
region $i_{k}$. This constraint will be enforced in the bilevel optimization
problem (IV-B) described below.

Figure 3: Illustrating the decomposition into switching and monitoring
segments.

###  IV-B Global Cost Decomposition

The global cost function (7) can be decomposed into the local cost
contributions of the monitoring and switching periods (see again Fig. 3 for an
illustration of the decomposition). Recall that the cost of the $k$th
monitoring period is reflected by $M_{k}^{\ast}(\tau_{k})$ in (IV-A). The cost
of the switching period between the $k$th and $(k+1)$st monitoring period is
denoted by

$S_{k}=\int_{t^{\mathrm{sw}}_{k}}^{t^{\mathrm{sw}}_{k}+\Delta_{k}}\sum_{i=1}^{M}\mathrm{tr}(\Omega_{i}(t))~{}\textup{d}t,$

where $\Delta_{k}=\delta(a^{\psi}_{k},a^{\varphi}_{k+1})$ is the (fixed)
duration of the $k$th switching period and $t^{\mathrm{sw}}_{k}$ denotes the
starting time of the $k$th switching period. Note that the dynamics of the
estimator’s covariance matrices $\Omega_{i}$ are all smooth during the
switching segment since no measurements are taken during this interval.
Although the agent may traverse multiple regions when transitioning from one
target region to the next, none of these traversed regions, by design, contain
a target.

With this decomposition, we establish the global bilevel optimization problem
τ1T(τ)∑_k=1^K ( M^∗_k(τ_k) + S_k ) τ_k≥δ_i_k(a^φ_k, a^ψ_k) for k = 1, …, K,
which couples the local segments. The cycle duration in dependence of the
monitoring durations is given by

$T(\tau)=\sum_{k=1}^{K}\left(\tau_{k}+\Delta_{k}\right).$

###  IV-C Optimizing the Monitoring Durations

The gradient of the cost function (IV-B) with respect to the parameter
$\tau_{k}$ is given by

$\frac{\frac{dM^{\ast}_{k}(\tau_{k})}{d\tau_{k}}\cdot
T(\tau)-\sum_{k=1}^{K}\bigl{(}M^{\ast}_{k}(\tau_{k})+S_{k}(u^{\mathrm{s}}_{k})\bigr{)}}{T(\tau)^{2}}.$
(10)

Note that the gradient $\frac{dM_{k}^{\ast}}{d\tau_{k}}$ is the _sensitivity_
of the $k$th optimal monitoring cost with respect to the parameter $\tau_{k}$.
Assuming that strong duality holds in the monitoring OCP (IV-A), we can
utilize the shadow price equation [20]

$\frac{dM_{k}^{\ast}}{d\tau_{k}}=-\lambda_{k},$

where $\lambda_{k}\in\mathbb{R}$ is the dual variable of the equality
constraint (IV-A). This allows us to compute the gradient of the global cost
function with respect to the local monitoring durations. We then solve (IV-B)
in an alternating fashion, where we fix the monitoring durations and then
simulate one (or multiple) cycle(s). Note that simulating a cycle includes re-
solving each monitoring OCP with the updated initial covariance matrices. We
utilize the gradient (10) and update the monitoring durations based on a
simple projected gradient method with diminishing step size. We compare two
optimization methods in the next section:

  1. (i)

Simulating to steady state before updating $\tau$.

  2. (ii)

Updating $\tau$ after every cycle $(i_{1},i_{2},\dots,i_{K})$.

The first method is motivated from the fact that the computed gradient (10)
describes the steepest ascending direction of the global cycle cost under a
given configuration consisting of the monitoring durations and the estimator
covariance matrices at the beginning of the cycle. If the cycle is not in
steady state, the computed gradient may be inaccurate. On the other hand, the
initial configuration, and thus the initial steady-state covariance matrices
may be far from the optimal, such that a precise gradient may be unnecessary
to obtain. As we decrease the amount of change in the configuration, the cycle
will eventually converge to a steady-state cycle, where the gradient will
become precise. This motivates the second variant. Although convergence in
bilevel optimization problems is difficult to establish, the next section
supports that both variants work well in terms of performance.

##  V Numerical Results

We discretize the local monitoring OCP (IV-A) via a multiple shooting method
using the RK4 integrator [19]. We model and solve this via CasADi [21]
utilizing the interior-point solver IPOPT [22]. Optimization of the monitoring
durations is done using a projected gradient descent with diminishing step
size. All simulations were conducted on a hardware featuring an Intel i5
processor running at 1.60GHz with 16GB of RAM. In the following results all
regions have randomly generated constant dynamics with norm bounded by one.
All targets either have a Gaussian sensing model or a nonconvex sensing model,
the exact form of which is omitted here for space reasons. All specific
details on problem set up can be found in the provided repository
github.com/hallfjonas/hytoperm.

We first consider the scenario illustrated in Fig. 1, which contains four
targets and ten regions. Fig. 4 compares the two variants (i) and (ii). Note
that in this particular scenario, variant (ii) completes its optimization
procedure before variant (i) even converges to a steady-state cycle. This is
due to the fact that this variant is capable of modifying the monitoring
periods more rapidly and as a result the steady-state is reached more quickly.
For both variants we see that the objective function is flat near the local
optimum, since the cost does not change for the final few iterations while the
monitoring durations are still updated. Fig. 5 depicts the controls for the
optimal cycle together with the mean estimation errors for variant (i). The
results show that throughout the cycle the control constraints $\|u(t)\|\leq
1$ are respected and that the mean estimation errors are indeed periodic. Both
variants converge to the same cycle, which is identical to the one shown in
Fig. 1.

To show that the method extends to more complex scenarios, we consider another
setting with 20 regions and 10 targets. Fig. 6 again demonstrates that
directly updating the configuration after each cycle reduces the number of
cycles required to reach the local optimum. Finally, we note that the cost
iterations, even for variant (i), are not monotonic. This is due to the fact
that the step size may be too large and result in a cost increase in the
current iteration. This is also reflected in the oscillation of the plotted
monitoring duration. Fig. 7 shows that both variants converge to the same
cycle.

Figure 4: Depicting the optimization of the monitoring durations $\tau$ for
the small experiment utilizing both variant (i) (yellow) and (ii) (blue). The
first plot shows the cost per cycle, and the second plot shows the evolution
of $\tau_{1}$. The vertical lines indicate the cycles where (i) reaches steady
state. Figure 5: Depicting the optimized steady-state control law together
with the mean estimation errors. Figure 6: Same as Fig. 4 but for the larger
experiment. Vertical lines indicate times at which variant (i) updates the
$\tau_{k}$. Figure 7: The optimized configurations lead to the same cycle for
both methods (i) (yellow) and (ii) (blue).

##  VI Conclusion and Future Work

In this paper we presented a novel PM formulation by assuming that the agent
has piecewise smooth dynamics. With the goal of optimizing an average steady-
state mean estimation error, we designed a periodic agent trajectory in two
stages. First, in an offline approach, we created a global path planner and
computed a target visiting sequence. This yielded an initial cycle for the
agent to follow; this cycle was then optimized online. Simulations support
rapid convergence to the optimal configuration for a given visiting sequence,
at least under the settings considered.

While the approach we developed is a complete solution of the PM problem
defined, improvements on each level are possible. For instance, we currently
do not optimize the switching segments with respect to the global cost metric
– they are computed once and fixed afterwards. Their optimization in a global
sense could benefit the overall performance. This is not a straightforward
task, since even small alterations can change the traversed regions along the
switching path. It is conceivable that suddenly another target region is
visited along the altered switching path, which modifies the visiting sequence
and requires the addition of a new monitoring segment.

## References

  * [1] R. N. Smith, M. Schwager, S. L. Smith, B. H. Jones, D. Rus, and G. S. Sukhatme, “Persistent ocean monitoring with underwater gliders: Adapting sampling resolution,” Journal of Field Robotics, vol. 28, no. 5, pp. 714–741, 2011\.
  * [2] D. W. Casbeer, D. B. Kingston, R. W. Beard, and T. W. McLain, “Cooperative forest fire surveillance using a team of small unmanned air vehicles,” International Journal of Systems Science, vol. 37, no. 6, pp. 351–360, 2006\.
  * [3] S. C. Pinto, N. A. Vickers, F. Sharifi, and S. B. Andersson, “Tracking multiple diffusing particles using information optimal control,” in 2021 American Control Conference (ACC), pp. 4033–4038, IEEE, 2021.
  * [4] Y. Zhu and S. B. Andersson, “Control policy optimization for data harvesting in a wireless sensor network,” in 2022 IEEE 61st Conference on Decision and Control (CDC), pp. 7437–7442, IEEE, 2022.
  * [5] S. K. K. Hari, S. Rathinam, S. Darbha, K. Kalyanam, S. G. Manyam, and D. Casbeer, “Optimal UAV route planning for persistent monitoring missions,” IEEE Transactions on Robotics, vol. 37, no. 2, pp. 550–566, 2020.
  * [6] A. Nurkanović and M. Diehl, “NOSNOC: A software package for numerical optimal control of nonsmooth systems,” IEEE Control Systems Letters, 2022\.
  * [7] J. Hall, A. Nurkanovic, F. Messerer, and M. Diehl, “LCQPow: A solver for linear complementarity quadratic programs,” arXiv preprint arXiv:2211.16341, 2022.
  * [8] S. Karaman and E. Frazzoli, “Incremental sampling-based algorithms for optimal motion planning,” Robotics Science and Systems VI, vol. 104, no. 2, pp. 267–274, 2010.
  * [9] X. Lan and M. Schwager, “Planning periodic persistent monitoring trajectories for sensing robots in gaussian random fields,” in 2013 IEEE International Conference on Robotics and Automation, pp. 2415–2420, IEEE, 2013\.
  * [10] J. Hall, L. E. Beaver, C. G. Cassandras, and S. B. Andersson, “A bilevel optimization scheme for persistent monitoring,” in 2023 62nd IEEE Conference on Decision and Control (CDC), pp. 247–252, 2023.
  * [11] A. F. Filippov, Differential equations with discontinuous righthand sides: control systems, vol. 18.  Springer Science & Business Media, 2013.
  * [12] J. Cortes, “Discontinuous dynamical systems,” IEEE Control systems magazine, vol. 28, no. 3, pp. 36–73, 2008.
  * [13] X. Lan and M. Schwager, “A variational approach to trajectory planning for persistent monitoring of spatiotemporal fields,” in 2014 American Control Conference, pp. 5627–5632, IEEE, 2014.
  * [14] M. Ostertag, N. Atanasov, and T. Rosing, “Trajectory planning and optimization for minimizing uncertainty in persistent monitoring applications,” Journal of Intelligent & Robotic Systems, vol. 106, no. 1, pp. 1–19, 2022.
  * [15] S. C. Pinto, S. B. Andersson, J. M. Hendrickx, and C. G. Cassandras, “Optimal periodic multi-agent persistent monitoring of a finite set of targets with uncertain states,” in 2020 American Control Conference (ACC), pp. 5207–5212, IEEE, 2020.
  * [16] S. Welikala and C. G. Cassandras, “Event-driven receding horizon control for distributed persistent monitoring in network systems,” Automatica, vol. 127, p. 109519, 2021.
  * [17] S. C. Pinto, S. B. Andersson, J. M. Hendrickx, and C. G. Cassandras, “Multi-agent persistent monitoring of targets with uncertain states,” IEEE Transactions on Automatic Control, 2022.
  * [18] L. Zhao, W. Zhang, J. Hu, A. Abate, and C. J. Tomlin, “On the optimal solutions of the infinite-horizon linear sensor scheduling problem,” IEEE Transactions on Automatic Control, vol. 59, no. 10, pp. 2825–2830, 2014\.
  * [19] A. V. Rao, “A survey of numerical methods for optimal control,” Advances in the Astronautical Sciences, vol. 135, no. 1, pp. 497–528, 2009.
  * [20] S. P. Boyd and L. Vandenberghe, Convex optimization.  Cambridge university press, 2004.
  * [21] J. A. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “Casadi: a software framework for nonlinear optimization and optimal control,” Mathematical Programming Computation, vol. 11, no. 1, pp. 1–36, 2019.
  * [22] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical programming, vol. 106, pp. 25–57, 2006.