#  Large deviations for Markov jump processeswith uniformly diminishing rates

Andrea Agazzi  Department of Mathematics, Duke University, 120 Science Dr,
Durham, NC 27710, USA<EMAIL_ADDRESS>,  Luisa Andreis  ,  Robert I. A.
Patterson  and  D. R. Michiel Renger  Weierstraß-Institut für Angewandte
Analysis und Stochastik, Mohrenstrasse 39, 10117 Berlin, Germany  {andreis,
patterson<EMAIL_ADDRESS>
###### Abstract.

We prove a large-deviation principle (ldp) for the sample paths of jump Markov
processes in the small noise limit when, possibly, _all_ the jump rates vanish
_uniformly_ , but slowly enough, in a region of the state space. We further
discuss the optimality of our assumptions on the decay of the jump rates. As a
direct application of this work we relax the assumptions needed for the
application of ldps to, _e.g._ , Chemical Reaction Network dynamics, where
vanishing reaction rates arise naturally particularly the context of mass
action kinetics.

##  1\. Introduction

###  1.1. Large deviations of Markov jump processes

We study a family of $d$-dimensional Markov jump processes
$\\{X^{v}\\}_{v\in\mathbb{N}}$ with state space $(v^{-1}\mathbb{Z})^{d}$,
deterministic initial condition $X^{v}(0)=x_{0}^{v}\in(v^{-1}\mathbb{Z})^{d}$
and generator:

$(\mathcal{L}^{v}f)(x):=v\sum_{r\in\mathcal{R}}\Lambda^{v}_{r}(x)\big{[}f(x+v^{-1}{\gamma^{r}})-f(x)\big{]}.$
(1)

Here $\mathcal{R}$ is the finite set of possible jumps,
$\gamma^{r}\in\mathbb{Z}^{d}$ are the fixed jump vectors, and
$v\Lambda^{v}_{r}\leavevmode\nobreak\ :\leavevmode\nobreak\
(v^{-1}\mathbb{Z})^{d}\to[0,\infty)$ the associated jump rates. The parameter
$v$ controls the noise in the system, and the scaling is chosen so that
$\Lambda^{v}_{r}(x)$ converge as $v\to\infty$. Under this scaling it is known
that the paths $X^{v}$ concentrate on solutions of the _fluid limit_ ODE
[Kurtz1970]:

$\frac{\mathrm{d}}{\mathrm{d}t}x(t)=\sum_{r\in\mathcal{R}}\lambda_{r}(x(t))\gamma^{r}\,,$
(2)

where the continuous rates $\lambda_{r}\leavevmode\nobreak\
:\leavevmode\nobreak\ \mathbb{R}^{d}\to[0,\infty)$ are the limits of
$\Lambda^{v}_{r}(x)$ as $v\to\infty$.

The process (1) and the ODE (2) are used as microscopic and macroscopic models
for a wide range of applications. For example, in the context of chemical
reactions $X^{v}$ denotes the concentrations (number of molecules per unit
volume) of $d$ species being transformed by a set of reactions $\mathcal{R}$.
Here, for each reaction $r\in\mathcal{R}$, the vectors $\gamma^{r}$ encode
which species are removed and created when a reaction $r$ occurs, and
$\lambda_{r}(x)$ are the reaction rates [Kurtz1972]. In that setting $v$ can
be interpreted as the volume size over which the concentrations are averaged.
Other typical applications using similar models include biological systems
involving predator-prey interaction, birth/cell division and death, biological
fitness models, as well as epidemiological models. In these settings, large-
deviation techniques are often applied to simplify the dynamical landscape of
the complex, high dimensional microscopic model (1) while retrieving
quantitative information about the random fluctuations around the mean (2).
This information can be used for example to study non-equilibrium
thermodynamics [esposito19, MielkePattersonPeletierRenger2015], to speed up
simulations of rare events [delmoral04, vde04, vde19], or to study spontaneous
transitions between metastable states [freidlin12] also in the multi-scale
setting [popovic19].

The classical proof of the Large-Deviation Principle (ldp) uses a tilting,
also called a change-of-measure technique. The main challenge there is that
the tilting can only be performed around sufficiently regular paths, whereas
the large-deviation principle needs to be proven for any non-typical path.
Therefore, the large-deviation lower bound requires an approximation argument,
either for the random process or for the rate functional. A particular
challenge in either case is to approximate a path without changing its
starting point, which is required when proving the large-deviation principle
under a deterministic initial condition. This becomes more difficult if the
jump rates vanish in some regions of the state space, which is however an
inherent property of the models used in many application domains. For example,
in the context of chemical kinetics, it is natural to assume that the rate of
a chemical reaction vanishes when the concentration of one or more of the
reactants approaches $0$. Similarly in the context of infectious disease
models, the rate of spread of a virus is usually modelled as a linear function
of the infected population.

######  Example 1.1.

The problem is nicely illustrated by the simplest model for autocatalysis or
cell division, in chemical notation: $\mathsf{A}\to 2\mathsf{A}$. In this case
there is only one species and one jump, so we may write the linear jump rate
as $\Lambda^{v}(x)=\lambda(x)=x$, starting from a concentration with one
particle, i.e. $X^{v}(0)=x_{0}^{v}=1/v$. Clearly, the process converges to the
solution of $\dot{x}(t)=x(t)$ with initial condition $x(0)=x_{0}=0$, that is
$x(t)\equiv 0$. In other words, the process is expected to stick to the
degenerate set $\partial\mathcal{S}=\\{x=0\\}$, which corresponds to the
boundary of the state space $\mathbb{R}_{+}$. However, it can be calculated
(as an application of [Kordecki97, Lemma 2.1]) that
$\lim_{v\to\infty}v^{-1}\log\mathbb{P}\big{(}X^{v}(t)\geq\delta\big{)}=\delta\log(1-e^{-t})$
for any fixed time $t>0$ and $\delta>0$. Although this is a large-deviation
result about the marginal $X^{v}(t)$, it suggests that the paths can also
“escape” from the boundary with finite large-deviation cost. On the other
hand, we note that choosing $x_{0}^{v}\equiv 0$ implies that $X^{v}(t)\equiv
0$ almost surely. Hence, whether or not escape is possible depends strongly on
the initial condition. This is a similar principle as what is sometimes called
a “well-prepared initial condition” in $\Gamma$-convergence theory
[Mielke2016].

For more general models of the type (1), one expects that as long as the jump
rates $\lambda_{r}(x)$ do not vanish too rapidly approaching the degenerate
points and when starting at a well-chosen initial condition near such points,
then the process will be able to escape with finite large-deviation cost, and
the large-deviation principle should still hold. In this paper, we show that
this is indeed the case. More specifically, denoting by $\partial\mathcal{S}$
the set where some – or possibly _all_ – of the reaction rates $\lambda$
vanish, we prove that when the rates near $\partial{\mathcal{S}}$ decrease
slower than
$\exp\left(-1/\mathrm{dist}(x,\partial{\mathcal{S}})^{\alpha}\right)$, i.e.,

$\mathrm{dist}(x,\partial{\mathcal{S}})^{\alpha}\log\lambda(x)\to
0\qquad\text{as }\quad\mathrm{dist}(x,\partial{\mathcal{S}})\to 0,$ (3)

the large-deviation principle holds under the assumption $\alpha\in[0,1)$.
Furthermore, we show that these conditions on the decay of the rates close to
the degenerate set are, in some sense, necessary and not just sufficient.

###  1.2. Literature and approach

Early papers [Feng1994dynamic, Leonard1995] establishing a sample-path large-
deviation principle for jump Markov processes mimicked the Dawson–Gärtner
approach [Dawson1987], where one first derives an abstract large-deviation
result for the empirical measure on paths, and then contracts it to obtain a
large-deviation principle for the path of the empirical measure. Another
approach, now considered the classic tilting technique, was first used in
[ShwartzWeiss95] assuming that all the jump rates are uniformly bounded away
from zero in the domain of interest. In [ShwartzWeiss05], the authors relaxed
this condition by assuming the existence _a subset_ of jumps with rates that
are uniformly bounded away from zero and that “push” the process away from the
degeneracy region. This assumption is further relaxed in [dupuis16, agazzi18,
agazzi182, anderson18], requiring only the existence of a _sequence_ of jumps
that sequentially transport the process away from the problematic region.

Recent works have taken steps to generalise these assumptions to the uniformly
vanishing case (when all rates can vanish in some region of state space), in
the context of chemical kinetics [pattersonrenger19] and in the context of
infectious disease models [pardoux16, pardoux17]. These papers give sufficient
conditions on the models at hand to bypass the technical difficulties
encountered in the proof of the ldp lower bound when some of the jump rates
are not bounded away from zero. We mention that the work [pattersonrenger19]
assumes “sufficiently random” initial conditions to bypass the problem, but we
shall focus on a deterministic initial condition.

The problem of vanishing rates is addressed more completely in [pardoux16,
pardoux17], where the authors obtain large-deviation estimates for vanishing
rates when the microscopic initial condition allows escaping the degenerate
set with positive – although vanishing in $v$ – probability. Their approach is
based on a careful adaptation of the standard tilting argument to obtain the
ldp lower bound for processes. In particular, to bypass the problem of jump
rates vanishing in some regions of state space, the change of measure
performed by the authors depends on the large-deviation scaling parameter $v$,
which is inversely proportional to the jump size. This replaces the problem of
escaping to an $\mathcal{O}(1)$ distance from these degeneracy regions
uniformly in $v$ to the one of escaping to $\mathcal{O}(1/\sqrt{v})$. Their
result allows for jump rates that behave as in (3) for $\alpha\in[0,1/2)$.

In this paper we bypass the change-of-measure technique altogether
establishing more direct and concise bounds on the process level. The main
challenge in realizing this strategy is to identify a set of paths occurring
with sufficient probability to recover the ldp lower bound while allowing for
simple estimation of such probability. The core of the approach consists of
showing separate estimates on the total number of jumps and on the types of
jumps for paths in such set. This allows us to extend the assumption (3) to
any $\alpha\in[0,1)$ while covering a larger family of processes than the
existing literature. In addition, we provide a counterexample showing that our
upper bound for the exponent $\alpha$ is optimal: If the rates of the process
decay as (3) for $\alpha\geq 1$ with sufficient uniformity, as we make precise
below, the process will no longer be able to escape the degenerate set with
finite large-deviation cost.

###  1.3. Outline

The paper is structured as follows. In Section 2 we introduce our notation, we
list our assumptions and we state our main result, Theorem LABEL:th:main,
namely the ldp. We also illustrate the generality of our result with some
examples. The proof of Theorem LABEL:th:main is split into two sections: In
Section LABEL:Sec:lower we prove the ldp lower bound, while Section
LABEL:Sec:upper deals with the ldp upper bound. Finally, in Section
LABEL:Sec:counterex, we discuss the optimality of our assumptions on the decay
of the rates and in Section LABEL:Sec:constants we summarize the key
quantities of the proof.

##  2\. Notation and results

We start by giving a concrete example of the properties of systems we aim to
generalize in this paper, introducing some important quantities in an
intuitive way:

######  Example 2.1 (Mass action kinetics).

In the context of chemical kinetics, one indexes the dimension of state space
with a set of _species_ $\\{\mathsf{S}_{i}\\}$ representing the chemical
compounds in the system of interest. To describe interactions between
different compounds one defines reactions $r\in\mathcal{R}$ via
$\gamma^{r,\mathrm{in}},\gamma^{r,\mathrm{out}}\in\mathbb{N}_{0}^{d}$ and

$\sum_{i=1}^{d}\gamma^{r,\mathrm{in}}_{i}\,\mathsf{S}_{i}\longrightarrow\sum_{i=1}^{d}\gamma^{r,\mathrm{out}}_{i}\,\mathsf{S}_{i}.$

This is to be understood as saying $\gamma^{r,\mathrm{in}}_{i}$ copies of
species $\mathsf{S}_{i}$ are consumed in each $r$-reaction while
$\gamma^{r,\mathrm{out}}_{i}$ copies are produced. The reaction rate
$\Lambda_{r}^{v}$ is specified via a rate constant $k_{r}\geq 0$ as

$\Lambda_{r}^{v}(x):=k_{r}\frac{1}{v^{\sum_{i}|\gamma_{i}^{r,\mathrm{in}}|}}\prod_{i=1}^{d}\binom{vx_{i}}{\gamma^{r,\mathrm{in}}_{i}}\gamma^{r,\mathrm{in}}_{i}!\qquad\forall
x\in\left(v^{-1}\mathbb{N}_{0}\right)^{d}\,,$

where $\binom{\cdot}{\cdot}$ denotes the binomial coefficient. The jump vector
is $\gamma^{r}:=\gamma^{r,\mathrm{out}}-\gamma^{r,\mathrm{in}}$. These rates
are bounded from above on compact sets, and they converge to
$\lambda_{r}(x)=k_{r}\prod_{i=1}^{d}x_{i}^{\gamma_{i}^{r,\mathrm{in}}}$ as
$v\to\infty$. It is easy to see that $\Lambda_{r}^{v}(x)=0$ whenever
$vx_{i}<\gamma^{r,\mathrm{in}}$, so that $X^{v}\in(v^{-1}\mathbb{N}_{0})^{d}$
almost surely. Consequently, the degenerate set of the limiting dynamics is a
subset of the boundary $\partial\mathbb{R}_{\geq 0}^{d}$. As we shall
illustrate in examples below, different choices of reactions result in $X^{v}$
being confined on subsets of $(v^{-1}\mathbb{N}_{0})^{d}$.

We start by defining the set of reachable points of the process. Throughout,
we fix a sequence of deterministic initial conditions $\\{x_{0}^{v}\\}$. By
the potentially degenerate character of the stochastic dynamics at hand, we
reduce the state space $(v^{-1}\mathbb{Z})^{d}$ to the set of reachable points
of the process with that initial condition
$x_{0}^{v}\in(v^{-1}\mathbb{Z})^{d}$:

$\mathcal{S}_{v}:=\\{x\in\left(v^{-1}\mathbb{Z}\right)^{d}\colon\mathbb{P}[\exists
t\geq 0,\,X^{v}(t)=x\mid X^{v}(0)=x_{0}^{v}]>0\\}\,.$

Note that, by definition, $\Lambda^{v}_{r}(x)=0$ whenever
$x+v^{-1}\gamma^{r}\notin{\mathcal{S}}_{v}$ for any $x\in\mathcal{S}_{v}$.
Assuming that in the limit $v\to\infty$, the initial values
$x_{0}^{v}\in\left(v^{-1}\mathbb{Z}\right)^{d}$ converge to
$x_{0}\in\mathbb{R}^{d}$, we write
$\displaystyle{\mathcal{S}=\bigcap_{n\in\mathbb{N}}\overline{\bigcup_{v\geq
n}\mathcal{S}_{v}}}$ where the raised line indicates topological closure. We
assume throughout that $\mathcal{S}$ is compact, but discuss how to relax this
assumption in Remark LABEL:r:exptight. We associate to $\mathcal{S}$ the set
of jumps

$\mathcal{R}_{\geq 0}:=\left\\{r\in{\mathcal{R}}\colon\exists
x\in\mathcal{S},\leavevmode\nobreak\ \lambda_{r}(x)>0\,\right\\}$

Notice that, depending on the sequence of initial conditions, the same Markov
process may have a different state space $\SS_{v}$ and different set of jumps
${\mathcal{R}}_{\geq 0}$. We refer to Example 2.3 for a situation where
$\mathcal{R}_{\geq 0}\neq\mathcal{R}$. However, by abuse of notation, we will
drop the index $\geq 0$ and refer to this set simply as ${\mathcal{R}}$.

Finally, we define the _degenerate_ set – also referred to as “boundary” from
its topological characterization in many application domains – as
$\partial\mathcal{S}:=\\{x\in\mathcal{S}\leavevmode\nobreak\
:\leavevmode\nobreak\ \exists r\in{\mathcal{R}},\leavevmode\nobreak\
\lambda_{r}(x)=0\\}$. This represents the set of points where the limiting
process is degenerate, i.e., where the classical proof of the large-deviation
principle will not immediately apply. Observe that this is a slight abuse of
both notation and terminology, since this degenerate set $\partial\mathcal{S}$
may be different from the actual topological boundary of the set
$\mathcal{S}$.

The following example clarifies the role of the sequence of initial conditions
on the resulting state space.

######  Example 2.2.

The mass action kinetics model $\mathsf{A}\leftrightarrow\mathsf{B}$ (see
Example 2.1 for definition of the rates and jump vectors) with initial
conditions $x_{0}^{v}=(0,1+1/v)$ results in
$\mathcal{S}_{v}=\\{x\in(v^{-1}\mathbb{Z}_{\geq 0})^{2}\leavevmode\nobreak\
:\leavevmode\nobreak\ x_{1}+x_{2}=1+1/v\\}$ and
$\mathcal{S}=\\{x\in\mathbb{R}_{\geq 0}^{2}\leavevmode\nobreak\
:\leavevmode\nobreak\ x_{1}+x_{2}=1\\}$.

######  Example 2.3.

The nontrivial effect of different sequences of initial conditions is captured
by the system

$\displaystyle\mathsf{B}\leftrightarrow
2\mathsf{B}\qquad\mathsf{A}\leftrightarrow\mathsf{2A+B},$

with mass action kinetics (see Example 2.1). For this model, the sequence
$x_{0}^{v}=(1/v,0)$ results in $\mathcal{S}_{v}=(v^{-1}\mathbb{Z}_{\geq
0})^{2}\setminus\\{0\\}$ and $\SS=\mathbb{R}_{\geq 0}^{2}$. However, if
$x_{0}^{v}=(0,1/v)$ we have $\mathcal{R}_{\geq 0}=\\{B\to 2B,2B\to B\\}$ and
the dynamics are restricted to
$\mathcal{S}_{v}=\\{0\\}\times(v^{-1}\mathbb{Z}_{\geq 0})$ resulting in
$\SS=\\{x\in\mathbb{R}_{\geq 0}^{2}\leavevmode\nobreak\ :\leavevmode\nobreak\
x_{1}=0\\}$.

###  2.1. Assumptions

To ensure existence of the limit, we require the reaction rates to satisfy
some conditions.

######  Assumption 1 (Convergence and regularity of rates).

We assume the following.

  1. a)

There exists a collection of non-negative functions
$\\{\lambda_{r}\\}_{r\in{\mathcal{R}}}$, Lipschitz continuous on a
neighborhood of ${\mathcal{S}}$ in $\mathbb{R}^{d}$, such that

$\lim_{v\to\infty}\sup_{x\in{\mathcal{S}}_{v}}\sum_{r\in{\mathcal{R}}}\left\lvert{{\Lambda_{r}^{v}(x)}-\lambda_{r}(x)}\right\rvert=0.$

  2. b)

There exists $\aleph>0$ such that for all $r\in{\mathcal{R}}$, $v>0$ and
$x\in\mathcal{S}_{v}$ with $\Lambda^{v}_{r}(x)>0$, we have

$\frac{\Lambda_{r}^{v}(x)}{\lambda_{r}(x)}\geq\aleph.$

As we outline in Section LABEL:s:path the proof of our main theorem is based
on the construction of short linear paths moving the process away from the
boundary $\partial{\mathcal{S}}$. We now introduce notation to decompose the
state space into subsets, in each of which the linear path will be fixed. More
precisely, following a standard approach first presented in [ShwartzWeiss05],
we cover the state space $\mathcal{S}$ with the relative interior of finitely
many convex, compact sets $\\{{\mathcal{A}}_{i}\\}_{i\in\mathcal{I}}$ with
${\mathcal{A}}_{i}\subseteq\mathcal{S}$ for all $i\in\mathcal{I}$. We then
define $\partial\mathcal{A}_{i}:=\partial{\mathcal{S}}\cap{\mathcal{A}}_{i}$
and let $\mathcal{I}^{\mathrm{bd}}\subseteq\mathcal{I}$ be the subset of
indices for which $\partial{\mathcal{A}}_{i}\neq\varnothing$.

We now present the assumptions for the lower bound. We assume that, whenever
the process starts from an initial condition close to $\partial{\mathcal{S}}$
(where possibly all the rates are zero), one can identify a finite sequence of
favorable jumps, which we call the _escape sequence_ , that push the process
away from the boundary. We further crucially assume that the rates of such
favorable jumps do not decay too fast as we approach the boundary. This is
captured by the following counterexample.

######  Example 2.4.

Consider the family of Markov jump processes $\\{X^{v}\\}$ with generator

$\mathcal{L}^{v}f(x):=ve^{-\frac{k}{x}}\left(f(x+v^{-1})-f(x)\right)\qquad\text{for
}f\leavevmode\nobreak\ :\leavevmode\nobreak\
v^{-1}\mathbb{N}_{0}\to\mathbb{R}\,,$ (4)

for any $k>0$. The above process, which for small $x$ is a time-changed
version of the autocatalytic process introduced in Example 1.1, has only one
possible jump in the positive direction with rate
$\Lambda^{v}(x)=e^{-\frac{k}{x}}$ s.t. $\lim_{\varrho\to
0}\varrho\left(\inf_{x\colon x\geq\varrho}\log\lambda(x)\right)=-k\neq 0\,.$
For the sequence of initial conditions $x_{0}^{v}=1/v$, we have
${\mathcal{S}}_{v}=v^{-1}\mathbb{N}$ and ${\mathcal{S}}=\mathbb{R}_{\geq 0}$.
Then, for any $w>0$ and $\varepsilon\in(0,w/2)$ the probability of observing a
realization of $X^{v}$ in an $\varepsilon$-neighborhood of the path $z(s)=sw$
on the interval $s\in[0,1]$ can be trivially estimated as

$\mathbb{P}\left[\sup_{t\in[0,1]}|X_{t}^{v}-z(t)|\leq\varepsilon\right]\leq\mathbb{P}[X_{1}^{v}\geq
w-\varepsilon]\leq\mathbb{P}[X_{1}^{v}\geq\widetilde{w}]$

for $\widetilde{w}=\min(k,w-\varepsilon)/2$. Denoting by $\tau_{i}$ the
waiting time between the $i-1$-th and $i$-th jump of the Poisson process
$X_{t}^{v}$ at $x\in\mathcal{S}_{v}$, we further have

$\mathbb{P}[X_{1}^{v}\geq\widetilde{w}]\leq\prod_{i=1}^{\lfloor
v\widetilde{w}\rfloor}\mathbb{P}[\tau_{i}\leq 1]=\prod_{i=1}^{\lfloor
v\widetilde{w}\rfloor}1-\exp[-(ve^{-kv/i})]\leq\exp\left(\sum_{i=1}^{\lfloor
v\widetilde{w}\rfloor}\left(\log v-kv/i\right)\right).$

The rough estimate above yields

$\frac{1}{v}\log\mathbb{P}[X_{1}^{v}\geq\widetilde{w}]\leq\frac{1}{v}\lfloor
v\widetilde{w}\rfloor\log v-k\sum_{i=1}^{\lfloor
v\widetilde{w}\rfloor}\frac{1}{i}<(\widetilde{w}-k)\log
v-k(1+\log\widetilde{w})$ (5)

which approaches $-\infty$ as $v\to\infty$.

As the example above shows, sufficiently fast decay of the rates of the
process $X^{v}$ implies the divergence of the large-deviation cost of any
nontrivial path starting on the boundary $\partial{\mathcal{S}}$. We now
proceed to give sufficient assumptions guaranteeing that this does not happen
in the general setting. In particular, to capture the idea of escaping a
boundary in the higher dimensional setting, for each $\mathcal{A}_{i}$ we
define directions $w_{i}$ with some structural properties (Assumption 2 a))
allowing to construct linear paths that leave such boundaries. We assume that
these paths can be realized as a sequence of jumps $\mathcal{E}_{i}$ whose
rates do not decay too fast (Assumption 2 b)-c)), as to avoid for the
realization of such path to have an infinite large-deviation cost. Denoting
throughout by $\mathcal{B}_{\varrho}(x)$ the Euclidean ball of radius
$\varrho$ in $\mathbb{R}^{d}$ and by $|A|$ the Lebesgue measure of the set
$A$, we summarize such assumptions below:

######  Assumption 2 (Escape).

There exist constants
$\varepsilon,\varepsilon^{\prime},\varepsilon^{\prime\prime}>0$ such that for
each $j\in\mathcal{I}$, the following holds:

  1. a)

If $j\in{\mathcal{I}}^{bd}$ there is a $w_{j}\in\mathbb{R}^{d}$ with
$\|w_{j}\|=1$ and $\kappa_{j}\in(0,1)$ such that whenever
$x\in\mathcal{A}_{j}$ and
$\inf_{y\in\partial{\mathcal{S}}}\|x-y\|<\varepsilon^{\prime}$ and
$t\in(0,\varepsilon)$

    1. i)

$t\mapsto\inf_{y\in\partial{\mathcal{S}}}\|x+t\,w_{j}-y\|$ is increasing, and

    2. ii)

$\mathcal{B}_{t\kappa_{j}}(x+t\,w_{j})\cap\partial{\mathcal{S}}=\varnothing$.

We write $\kappa_{-}=\min_{j\in\mathcal{I}^{bd}}\\{\kappa_{j}\\}$. If
$j\in{\mathcal{I}}\setminus{\mathcal{I}}^{bd}$ we choose $w_{j}=0$.

  2. b)

There exists a finite sequence ${\mathcal{E}}_{j}:=(r_{1},\dots,r_{n_{j}})$ of
jumps in ${\mathcal{R}}$ with

${\limsup_{v\to\infty}\frac{1}{v}\left|{\log\Lambda_{r_{k}}^{v}\left(x_{0}^{v}+v^{-1}\sum_{i=1}^{k-1}\gamma^{r_{i}}\right)}\right|=0,}\qquad
k=1,\dotsc n_{j}\,,$

and $\sum_{i=1}^{n}\gamma^{r_{i}}=\alpha_{j}w_{j}$ for some $\alpha_{j}>0$ .

  3. c)

Defining
$Z_{0}^{v}:=\\{x_{0}^{v}+v^{-1}\sum_{i=1}^{k}\gamma^{r_{i}}\leavevmode\nobreak\
:\leavevmode\nobreak\ k\in(0,\dots,n_{j}-1)\\}$, for all
$r\in{\mathcal{E}}_{j}$ and $T>0$

$\lim_{\varrho\to
0}\sup_{x\in{\mathcal{A}}_{j}\cup\bigcup_{v}Z_{0}^{v}}\int_{0}^{\varrho}\\!\left|\log\lambda_{r}(x+sw_{j})\right|\mathrm{d}s=0\,.$

  4. d)

Let $\mathcal{W}_{j,\kappa^{\prime\prime}}:=\\{w_{j}+y\leavevmode\nobreak\
:\leavevmode\nobreak\ \|y\|<\kappa^{\prime\prime}\\}$. There exists
$\kappa^{\prime\prime}<\kappa_{-}/3$ such that for any
$x\in\mathcal{A}_{j}\cup\bigcup_{v}Z_{0}^{v}$ we have that for all
$r\in\mathcal{R}$ with $\lambda_{r}(x)<\varepsilon^{\prime\prime}$ the rates
$\lambda_{r}(\,\cdot\,)$ are nondecreasing along paths $t\mapsto x+tw$ for any
$w\in\mathcal{W}_{j,\kappa^{\prime\prime}}$, for $t\in(0,\varepsilon)$.

It is readily verified Assumption 1 and 2 are satisfied by mass action
kinetics rates on a convex domain [agazzi18].

######  Remark 2.5.

While Assumption 2 c) is natural in terms of our proof, we note that it is
automatically satisfied whenever there exists $\alpha\in[0,1)$ such that

$\lim_{\varrho\to
0}\varrho^{\alpha}\left(\inf_{x\in{\mathcal{A}}_{j}{\cup\bigcup_{v}Z_{0}^{v}}\colon
d(x,\partial{\mathcal{A}}_{j})>\varrho}\log\lambda_{r}(x)\right)\to
0\quad\hbox{for all }r\in\mathcal{E}_{j}\,,{\hbox{\rm(6)}\kern
5.69046pt\kern-5.69046pt}:decay{=}e5{6}$


# Time-Series Regeneration with Convolutional Recurrent Generative Adversarial
Network for Remaining Useful Life Estimation

Xuewen Zhang, Yan Qin, Chau Yuen, Lahiru Jayasinghe, and Xiang Liu This work
was supported by the A*STAR-NTU-SUTD Joint Research Grant on Artificial
Intelligence Partnership under Grant RGANS1906 and in part supported by the
National Natural Science Foundation of China under Grant 61903327.
(Corresponding author: Yan Qin)X.W. Zhang and X. Liu are with the Department
of Software and Microelectronics, Peking University, Peking, 100871 China. (
e-mail<EMAIL_ADDRESS>xliu@ss.pku.edu.cn).Y. Qin, C. Yuen, and L.
Jayasinghe are with the Engineering Product Development Pillar of Singapore
University of Technology and Design, 8 Somapah Road, 487372 Singapore.
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

For health prognostic task, ever-increasing efforts have been focused on
machine learning-based methods, which are capable of yielding accurate
remaining useful life (RUL) estimation for industrial equipment or components
without exploring the degradation mechanism. A prerequisite ensuring the
success of these methods depends on a wealth of run-to-failure data, however,
run-to-failure data may be insufficient in practice. That is, conducting a
substantial amount of destructive experiments not only is high costs, but also
may cause catastrophic consequences. Out of this consideration, an enhanced
RUL framework focusing on data self-generation is put forward for both non-
cyclic and cyclic degradation patterns for the first time. It is designed to
enrich data from a data-driven way, generating realistic-like time-series to
enhance current RUL methods. First, high-quality data generation is ensured
through the proposed convolutional recurrent generative adversarial network
(CR-GAN), which adopts a two-channel fusion convolutional recurrent neural
network. Next, a hierarchical framework is proposed to combine generated data
into current RUL estimation methods. Finally, the efficacy of the proposed
method is verified through both non-cyclic and cyclic degradation systems.
With the enhanced RUL framework, an aero-engine system following non-cyclic
degradation has been tested using three typical RUL models. State-of-art RUL
estimation results are achieved by enhancing capsule network with generated
time-series. Specifically, estimation errors evaluated by the index score
function have been reduced by 21.77$\%$, and 32.67$\%$ for the two employed
operating conditions, respectively. Besides, the estimation error is reduced
to zero for the Lithium-ion battery system, which presents cyclic degradation.

###### Index Terms:

Multivariate time-series generation, degradation pattern learning, generative
adversarial network, convolutional recurrent network, remaining useful life
estimation.

##  I Introduction

Remaining useful life (RUL) estimation plays a vital role in condition-based
maintenance, which estimates the useful life left on diversely daily used
components and industrial systems in advance before they thoroughly break down
[1]. RUL estimation enables early warning for a maintainer to avoid unplanned
downtime of equipment that may cause catastrophic consequences [2]. For
instance, estimation of remaining hours for running aero-engines contributes
to the maintenance plan and avoids accidents; estimation of available
discharging cycles for Lithium-ion batteries (LiB) powered devices helps to
timely replace failed batteries. With the ever-increasing strict requirements
of operation safety, providing accurate RUL estimation becomes a significant
concern in both academic and industrial societies.

Nowadays, advanced machine learning methods, especially the convolutional
neural network (CNN) [3] and the recurrent neural network [4], have been the
mainstream for RUL estimation. As a challenging issue in many industrial
scenarios, like machinery, aerospace engineering, energy storage, etc., RUL
has already attracted substantial attention in each research field. According
to the data presentation of failure components, the current RUL methods could
be broadly classified into two categories regardless of associated
disciplines, in which one is the non-cyclic degradation pattern, and the other
one is cyclic degradation. Non-cyclic degradation refers to that a particular
object continuously degrades over time, which is one-time. Unlike the non-
cyclic degradation pattern, cyclic degradation follows an inconsistently
repetitive but similar degradation behavior for the same object, leading to
performance indicator deterioration from cycle to cycle. For RUL estimation of
non-cyclic degradation pattern, Zheng et al. [5] applied deep long short-term
memory network (LSTM) into RUL estimation to learn the long-term degradation
trend of faulty components. Experimental results supported to conclude that
LSTM outperformed other traditional methods, e.g., support vector machine,
with respect to RUL estimation. Wu et al. [6] proposed a similar RUL
estimation structure to [5], but differential features of measurements were
additionally fed to present extra degradation information. Jayasinghe et al.
[7] incorporated temporal convolutions into LSTM to comprehensively consider
both long-term and short-term time dependencies. Assigning weights for
different layers, Chen et al. [8] put forward an attention LSTM-based RUL
estimation model for machinery process by jointly considering features at each
layer rather than the final features. Recently, the capsule network (CapsNet)
has been applied for RUL estimation, which has the ability to present features
in high-dimension space [9]. For RUL estimation of cyclic-degradation pattern,
the mainstream machine learning algorithms have been reported. Taking the LiB
system as an example, informative measurements, like charging/discharging
voltage, charging/discharging current, and impedance, can be used for the RUL
estimation by learning the regression relationship between these features and
the degradation performance [10], [11]. In additional to the regression way to
perform RUL estimation, it is feasible to treate the RUL estimation in an
autoregression way due to the cyclic degradation nature through learning past
performance to forecast the future ones [12], [13]. GAN is capable of
forecasting by constructing the generator as an autoregressive model. For
instance, GAN is employed for stock price prediction [14], the typhoon
trajectory [15], and failure frame [16].

Despite the employment of various advanced deep learning algorithms, the
essence of these models [3]-[16] lies in regressing the concerned performance
indicator, which may be high-cost or time-consuming to measure, on available
sensor data. In other words, a regression relationship will be developed for
two groups of specific variables. The methods mentioned above inferred faulty
evolution patterns from measured time-series, assisting preventive maintenance
and quantitative indication. Therefore, the acquisition of time-series data is
the primary but fundamental step to obtain satisfied estimation. This step
involves a large amount of data to train a good neural network to avoid
overfitting and keep generalization ability. In the above literature, the data
are run-to-failure time-series that record long-term degradation procedure of
machines or components from health status to failure, indicating a whole time-
series as a sample. Although it is not difficult to collect plenty of data in
the era of Big Data, collecting enough run-to-failure time-series is still an
obstacle in practice. Taking the turbofan engine as an example, it may take
months even years for a turbofan engine to naturally experience initial
degradation until it completely steps into failure status. Practically,
failure status is undesirable and prohibited to avoid safety violation and
catastrophic damage. Hence parts are replaced much before failure due to
safety reasons. As a result, data shortage in terms of failure cycles is an
ever-increasing challenge for advanced machine learning-based RUL estimation
to successfully fulfill the task of health prognosis.

As an alternative to physical experiments, data generation technology is a
promising solution in enriching data, creating realistically synthetic data
through learning data distribution from real sequences. Among data generation
algorithms, the most important one is the generative adversarial network (GAN)
proposed in 2014 [17]. It succeeds in introducing a novel adversarial
modelling framework, which includes a generative network and a discriminative
network. These two networks construct a two-player minimax game, using CNNs to
sidestep the calculation of intractable likelihood functions. Since GAN has
been proposed, its superiorities in data generation have been noticed and
attracted extensive studies in image synthesis and natural language
processing. Recently, the success of GAN in synthesizing sequences raises
attention [18]. Li et al. [19] applied GAN to increase available tokens for
discrete sequences. Ramponi et al. [20] proposed a time-conditional GAN for
sequence augmentation, which introduced a deconvolutional neural network in
generative step and CNN in the discriminative network. Since the above-
mentioned methods [18]-[20] are not designed for time-series generation, the
characteristics of temporal correlations between past samples and future data,
which is the nature of industrial time-series, cannot be ensured if they are
directly applied. Besides, industrial time-series are always multivariate
rather than a single variable, presenting correlations between variables. To
solve these problems, Esteban et al. [21] attempted to generate medical time-
series by replacing CNN with LSTM in the original GAN framework to capture
long-term temporal correlation. However, this method may not be suitable for
run-to-failure data, in which inherent complexity and strong dynamics are
observed. That is, multivariate in degradation procedure simultaneously
present serious non-stationarity and long-term temporal correlation within
each cycle. In addition, some industrial time-series are cyclic in nature,
e.g., the battery system that is charged and discharged over repeating cycles.
These characteristics cause challenges to generate high-quality time-series
using current methods. Moreover, to the best of the authors’ knowledge, few
studies have been reported about the hybrid RUL model with time-series
regeneration, which is a non-destructive and safe way to overcome the data
shortage problem of RUL estimation.

To address the aforementioned problems, an enhanced RUL estimation framework
is proposed, which combines the proposed time-series generation method, named
convolutional recurrent GAN (CR-GAN), with current RUL estimation models.
First, both generator and discriminator are improved in the proposed CR-GAN.
In the generator, a two-channel fusion model is given to simultaneously learn
both spatial and temporal evolution nature of the concerned multivariate time-
series. To achieve this, CNN and LSTM are combined to ensure learning ability,
and initial multivariate time-series are generated according to the hidden
relationship with consideration of external information. Next, the
discriminator employs LSTM to verify the performance of the generated data,
acting as an evaluator. The interactions between the generator and the
discriminator ensure that the generated data are realistic (i.e., similar to
the real data). Finally, with augmented time-series, the proposed RUL
framework is tested and compared in actual industrial scenarios. The
contributions of the proposed method are summarized below:

1) An enhanced RUL estimation concept is first proposed with the idea of time-
series self-generation to overcome the shortage of run-to-failure data in a
low-cost, rapid, and safe way, improving the current RUL estimation models
with more data;

2) CR-GAN model is designed to ensure the quality of generated time-series,
serving both non-cyclic and cyclic degradation patterns in a unified
framework;

3) We evaluate our proposed method with an aero-engine system [22] and a LiB
battery system [23] both from NASA to demonstrate its effectiveness in
reducing RUL estimation error with generated time-series.

The remainder of this article is structured as follows: In Section II, the RUL
estimation problem is classified according to different degradation patterns,
following with the motivation. The proposed method is detailed in Section III,
including data preparation, specific steps of the enhanced RUL estimation
framework, and online implementation. Section IV illustrates the efficacy of
the proposed method in typical systems. Conclusions are drawn in the last
section.

Figure 1: The main components and flowchart of (a) the traditional methods and
(b) the proposed method.

##  II Problem Formulation and Motivation

This section summarizes current RUL estimation problems into two groups
according to different degradation patterns, including non-cyclic degradation
and cyclic degradation. After that, the motivation for the proposed method is
presented.

###  II-A Problem Formulation of RUL Estimation

####  II-A1 RUL estimation of non-cyclic degradation pattern

Non-cyclic degradation refers to that a certain object continuously degrades
over time, which is one-time. For example, when an aero-engine begins to fail,
its RUL decreases until it totally breaks down. Thereby, the non-cyclic
degradation of an object constructs a two-dimensional data matrix. Assuming
$J$ process variables are measured at discrete time instance
$k=1,2,...,K_{i}$, a matrix $\textbf{X}_{i}(J\times K_{i})$ is formed, where
$K_{i}$ is the length of the $i^{th}$ object. Correspondingly, a vector
recording the remaining life is $\textbf{Y}_{i}=[K_{i},K_{i}-1,...,1]$, which
decreases linearly with the increase of time. Given there are $I$ run-to-
failure objects, a series of pair-wise data matrixes
$\\{\textbf{X}_{1},\textbf{Y}_{1}\\}$,$\\{\textbf{X}_{2},\textbf{Y}_{2}\\}$,…,$\\{\textbf{X}_{I},\textbf{Y}_{I}\\}$
are prepared, each of which corresponds to a certain object. For RUL
estimation, the goal is to find a degradation pattern ${\phi_{nc}}$ from
objects 1, 2, $...,I$ by minimizing the estimation error as follows,

$\min\limits_{\phi_{nc}}\sum_{i=1}^{I}\sum_{k=1}^{K_{i}}\|\textbf{Y}_{i}(k)-{\phi_{nc}}(\textbf{X}_{i}(1:k))\|_{2}^{2}$
(1)

where $\textbf{X}_{i}(1:k)$ is time-series ranged from time 1 to $k$, and
$\textbf{Y}_{i}(k)$ is the $k^{th}$ sample in $\textbf{Y}_{i}$.

####  II-A2 RUL estimation of cyclic degradation pattern

Different from the non-cyclic degradation pattern, cyclic degradation follows
an inconsistently repetitive but similar degradation behavior for the same
object, leading to performance indicator deterioration from cycle to cycle.
For instance, after a discharging cycle, LiB’s real capacity may slightly
drop. After thousands of charging and discharging cycles, the capacity loss
will cumulate, and its real capacity will reach the end of life (EOL)
threshold, indicating the replacement of LiB. For such a kind of degradation
pattern, assuming $J$ process variables are measured at discrete time instance
$k=1,2,...,K_{c}$ in the $c^{th}$ cycle, a data matrix $\textbf{X}_{c}(J\times
K_{c})$ is formed, where $K_{c}$ is the cycle length. Correspondingly, the
performance indicator is recorded as $Y_{c}$, which presents a cumulative
influence of the whole cycle. By collecting $C$ run-to-failure cycles, a
series of pair-wise cycling data $\\{\textbf{X}_{1},Y_{1}\\}$,
$\\{\textbf{X}_{2},Y_{2}\\}$, …, $\\{\textbf{X}_{C},Y_{C}\\}$, each of which
corresponds to a certain cycle. For RUL estimation, the goal is to find a
relationship ${\phi_{c}}$ from cycles 1, 2, $...,C$ by minimizing the
estimation error below,

$\min\limits_{\phi}\sum_{c=1}^{C}\|Y_{c}-{\phi_{c}}(\textbf{X}_{c})\|_{2}^{2}$
(2)

Due to cyclic degradation nature, it is promising to consider RUL estimation
in an autoregression modeling way through directly learning past performance
to forecast the future ones. However, this way may have a real-time
requirement on performance indicators. Since this article mainly focuses on
the regression way through easy-measured sensor data to predict the concerning
performance, the autoregressive-based RUL estimation models are not included
in this classification.

###  II-B Motivation

From the above analysis, data matrixes of both degradation patterns can be
arranged as a three-dimensional matrix $\textbf{X}(N\times J\times K_{n})$,
where the symbol $N$ indicates the number of components $I$ in non-cyclic
degradation or cycles $C$ for the same object in cyclic degradation, and
$K_{n}$ is the corresponding sample length. For readability,
$\textbf{X}(N\times J\times K_{n})$ is used as a general form in the following
analysis.

Strong dynamics and non-stationary characteristics are observed for both types
of degradation patterns, resulting in varying data distribution of process
variables. Although the number of variables $J$ is fixed, it is more
reasonable to consider the variable dimension combining sampling time. That
is, the original data matrix with dimensions $J\times K$ is extended to
$1\times JK$ by treating a complete run-to-failure time-series as a basic data
unit. Usually, $JK$ is greatly large than $N$, resulting in a shortage of
training cycles. Instead of increasing the amount of cycling data physically,
it is more affordable to enrich data through a data-driven way. If we could
model the degradation nature, it is possible to construct a realistic
simulation to mimic it. The above analysis is the basic idea of the proposed
method, which meets well with the purpose of GAN.

Figure 2: The data presentation of (a) a single variable, (b) non-cyclic
degradation pattern, and (c) cyclic degradation pattern.

##  III RUL Estimation with Time-series Generation

In this section, the organization and highlights of the proposed method are
illustrated in Fig. 1, in which five sequential parts are given. First, data
preparations are conducted on both non-cyclic and cyclic degradation systems
to transform run-to-failure cycles into a three-dimensional data matrix
without changing their physical structure. Next, data segmentation is
performed to increase the data amount. After that, the inherent evolution
trend of degradation patterns is learned through the proposed CR-GAN, enabling
the generation of similar degradation trajectories with the same nature. Then
generated temporal sequences combined with the real run-to-failure data
sequences are mixed as new inputs to train RUL estimation models, improving
its estimation ability. At last, when a new time-series is available, the
well-trained RUL estimation model can be adopted for online prediction.

###  III-A Data Preparation

The degradation system always consists of dozens of slowly varying variables,
as shown in Fig. 2(a). As mentioned in Section II.A, the data structure in a
complete run-to-failure object following a non-cyclic degradation pattern is
denoted as $\textbf{X}_{i}(J\times K_{i})$ as shown in Fig. 2(b). A three-
dimensional form can be arranged by putting data matrixes together without
considering the order of objects. However, for cyclic degradation, cycling
data should be arranged according to the ascending order of cycle number.
Therefore, a series of sequential two-dimension data matrixes are extended to
a three-dimensional data matrix X as shown Fig. 2(c), which takes spatial
correlation over cycles into account inherently.

Denote the maximum and the minimum of each variable among the data matrix X as
$\textbf{M}_{max}=[M_{1},M_{2},...,M_{J}]$ and
$\textbf{V}_{min}=[V_{1},V_{2},...,V_{J}]$, respectively. Each sample is
scaled using minimax normalization into range [0, 1] as follows,

$\displaystyle{\overline{x}}(n,j,k)=\frac{x(n,j,k)-V_{j}}{M_{j}}$ (3)

where $x(n,j,k)$ is the $j^{th}$ variable in the $k^{th}$ time index of the
$n^{th}$ cycle, $V_{j}$ is the minimum of the $j^{th}$ variable, $M_{j}$ is
the maximum of the $j^{th}$ variable, and $\overline{x}(n,j,k)$ is the
normalized sample in data matrix $\overline{\textbf{X}}$.

###  III-B Data Segmentation

To enable the basic training of the CR-GAN model, data segmentation is
performed after normalization. The $n^{th}$ complete run-to-failure time-
series $\mathbf{X}_{n}$ in $\overline{\textbf{X}}$ is divided into a series of
overlapped local trajectories. Each local trajectory is obtained by sliding
the complete time-series along the time direction with a fixed window size
$L$. A complete time-series is transformed into many locals trajectories to
enrich the training data in a step-wise way. Taking engine data [22] as an
example, Fig. 3 illustrates the procedure of data segmentation with step size
1 using the first two variables. By performing the data segmentation for all
run-to-failure time-series, the number of $N$ complete time-series is enlarged
to $N(K-L+1)$, where $K$ is the total length of a run-to-failure time-series.
For brevity, data matrixes after segmentation are still denoted as
$\overline{\textbf{X}}$. As a tunable parameter, the value of $L$ is suggested
to choose 80$\%\sim$90$\%$ of $K$, enriching the amounts of time-series from
10 to 20 times without changing of the degradation behavior.

Figure 3: The illustration of data segmentation for a complete run-to-failure
time-series from engine data.

###  III-C The CR-GAN Model

Two kinds of characteristics need to be carefully analyzed in generating high-
quality run-to-failure time-series: the first one is temporal correlations
between the past samples and future data in each cycle, and the second one is
the similarities and dissimilarities between different cycles along cycle
direction. These characteristics are ensured by a novel Convolutional
Recurrent GAN (CR-GAN) designed with the combination of LSTM and CNN.

(a) The detailed structure of generative network

(b) The detailed structure of discriminative network

Figure 4: The detailed network structure for (a) Generator and (b)
Discriminator.

####  III-C1 The structure of CR-GAN Model

The generative model $G$ introduces both convolutional ability and long-term
prediction ability to generate run-to-failure time-series. As shown in Fig.
4(a), the generator integrates convolutional layers and LSTM layers together.
Specifically, CNN layers focus on describing spatial features of the sequence,
and LSTM layers aim at learning temporal correlations between the past samples
and the future samples. Moreover, the fusion of spatial and temporal
information is achieved by putting these two networks into an integrated
network. Here, the noise vector is assumed to follow Gaussian distribution,
which has been widely adopted in most literature [14]-[21], [24]-[27]. Through
the temporal learning ability, the generator outputs a series of generated
time-series for further evaluation of the discriminator. Fig. 4(b) describes
the network structure of the discriminator, where the quality of generated
time-series is evaluated by classifying generated time-series and real time-
series. The discriminator is implemented using the LSTM network and the fully-
connected layer. For concisely, only one convolutional layer or one LSTM layer
have been presented in both the generator and the discriminator, and more
layers can be easily stacked to achieve better performance during practical
training.

For non-cyclic data, data distributions of variables are varying along time
direction for the same object. However, for cyclic data, the data distribution
of a variable in a specific cycle is treated as a data unit, and it slowly
changes over cycles. Thereby, the time index of non-cyclic data and cycle
index in cyclic data could serve as an auxiliary indicator to achieve better
RUL estimation. That is, the time index can be served as the condition vector
$\mathbf{Y}$.

By introducing condition information $\mathbf{Y}$ into Generator $G$ and
Discriminator $D$, CR-GAN replaces $D(\textbf{x})$ and $G(\textbf{z})$ with
$D(\textbf{x}|\mathbf{Y})$ and $G(\textbf{z}|\mathbf{Y})$, respectively. At
this point, both the generator and the discriminator take extra information Y
as a condition, which is fed into the two models as additional input as shown
in Figs. 4(a) and 4(b). Through the continual interaction between
Discriminator and Generator, unqualified time-series will be filtered out via
iterative solution, which is specified in the following part.

####  III-C2 Iterative solution of CR-GAN Model

The loss of Discriminator consists of two parts. The first one is the
probability of generated samples belonging to the real data, which is as large
as possible. The second one is the probability of generated samples belonging
to the synthetic data, which is as small as possible. As a result, the
objective function of Discriminator is constructed as follows,

$\displaystyle maxV_{D}$ $\displaystyle=\mathbb{E}_{\textbf{x}\sim
p_{\text{data}}(\textbf{x})}[\log D(\textbf{x}|\mathbf{Y})]$ (4)
$\displaystyle+\mathbb{E}_{\textbf{z}\sim
p_{\textbf{z}}(\textbf{z})}[\log(1-D(G(\textbf{z}|\mathbf{Y})))$

where $p_{z}(\textbf{z})$ is the aprior distribution of noise variables z,
$G(\textbf{z})$ is the generated data using Generator $G$, $D(\textbf{x})$
presents the probability that x came from the given real data, and
$D(G(\textbf{z}))$ stands for the probability that generated data come from
the distribution of the given real data.

Input: Normalized real run-to-failure time-series $\overline{\mathbf{X}}$,
noise vector $\mathbf{z}$ follows Gaussian distribution

Output: Generator $\mathbf{\Theta}_{g}$ and Discriminator
$\mathbf{\Theta}_{d}$

1 Initialize $\mathbf{\Theta}_{g}$ of Generator and $\mathbf{\Theta}_{d}$ of
Discriminator

2 repeat

3 repeat//Update the discriminator network

4 Sample $m$ noise vector $\mathbf{z}_{1},\mathbf{z}_{2},...,\mathbf{z}_{m}$
from the prior distribution $P_{z}$ and pick up time index ${\mathbf{Y}}_{m}$
from known range;

5 Generate $m$ time-series with generative network
$\tilde{\mathbf{x}}_{m}=G(\mathbf{z}_{m}|{\mathbf{Y}}_{m})$;

6 Update discriminator parameter $\mathbf{\Theta}_{d}$ to maximize,

7
$V_{D}=\frac{1}{m}\sum_{i=1}^{m}logD(\mathbf{x}_{i}|\mathbf{Y}_{i})+\frac{1}{m}\sum_{i=1}^{m}log(1-D(\tilde{\mathbf{x}}_{i}|{\mathbf{Y}}_{i}))$

8
$\mathbf{\Theta}_{d}\leftarrow\mathbf{\Theta}_{d}+\nabla\tilde{V}_{D}(\mathbf{\Theta}_{d})$

9 until  _Given iterations (here is 3)_ ;

10 repeat//Update the generator network

11 Sample another $m$ noise vector
$\mathbf{z}_{1},\mathbf{z}_{2},...,\mathbf{z}_{m}$ from the prior distribution
$P_{z}$ and pick up time index ${\mathbf{Y}}_{m}$ from known range;

12 Update generator parameter $\mathbf{\Theta}_{g}$ to minimize,

13
$V_{G}=\frac{1}{m}\sum_{i=1}^{m}logD(\mathbf{x}_{i}|\mathbf{Y}_{i})+\frac{1}{m}\sum_{i=1}^{m}log(1-D(G(\mathbf{z}_{i}|\tilde{\mathbf{Y}}_{i})))$

14
$\mathbf{\Theta}_{g}\leftarrow\mathbf{\Theta}_{g}+\nabla\tilde{V}_{G}(\mathbf{\Theta}_{g})$

15 until  __;

16until  _$\mathbf{\Theta}_{g}$ and $\mathbf{\Theta}_{d}$ converged_;

Algorithm 1 Iterative solution of CR-GAN

The discriminator is trained according to the given number of iterations.
Next, the generator loss $V_{G}$ is minimized by keeping weights of the
discriminator constant as follows,

$minV_{G}=\mathbb{E}_{\textbf{z}\sim
p_{\textbf{z}}(\textbf{z})}[\log(1-D(G(\textbf{z})))]$ (5)

The training of CR-GANs involves optimizations of $V_{D}$ and $V_{G}$
iteratively. Conceptually, the training procedure corresponds to a minimax
two-player game as below,

$\displaystyle\min_{G}\max_{D}V(D,G)$
$\displaystyle=\mathbb{E}_{\textbf{x}\sim p_{\text{data}}(\textbf{x})}[\log
D(\textbf{x}|\mathbf{Y})]$ (6) $\displaystyle+\mathbb{E}_{\textbf{z}\sim
p_{\textbf{z}}(\textbf{z})}[\log(1-D(G(\textbf{z}|\mathbf{Y})))]$

where the ideal objective is the discriminator cannot distinguish the
generated samples from the training samples.

The solving procedure of CR-GAN is specified in Algorithm 1, in which
Discriminator and Generator are iteratively executed until converge. Feeding
$\overline{\textbf{X}}$ into the proposed CR-GAN with additional information
Y, a series of qualified generated run-to-failure cycles are generated as
${\textbf{X}_{g}}(N_{g}\times J\times K)$, where $N_{g}$ is the number of
generated cycles.

Remark: Due to the randomness of the input noise vector, irrational time-
series are possibly generated with CR-GAN. The RUL estimation ability may
deteriorate if a false degradation pattern is learned. To filter out these
irrational time-series, we employ a “Train on Real, Test on Synthetic (TSTR)”
concept [21] to recognize the unqualified time-series. In this approach, the
real run-to-failure time-series are trained as a supervised model using CNN.
Then we input the generated time-series using CR-GAN into the trained
classification model. The irrational time-series will be identified if it does
not belong to the class of real time-series. These irrational time-series will
be discarded and not used in subsequential modeling.

###  III-D Enhanced RUL Model with Data Generation

With the generated data, the enhanced RUL framework is constructed. By merging
${\textbf{X}_{g}}$ into $\overline{\textbf{X}}$ along cycle-direction, a new
data matrix ${\textbf{X}_{n}}$ is formed with dimensions $N_{n}\times J\times
K$, in which $N_{n}=N_{g}+N$. Selecting a proper RUL estimation model, and
training the model structure with the matrix ${\textbf{X}_{n}}$. Also
${\textbf{X}_{n}}$ can be directly fed into existing RUL estimation models
with slight adjustments of tunable parameters, such as the number of nodes in
LSTM or CNN layers. The RUL estimation relationship ${\phi_{nc}(\cdot)}$ or
${\phi_{c}(\cdot)}$ as shown in Eqs. 1 and 2 are achieved when the iteration
epoches are reached or the iteration accuracy are met.

###  III-E Online Estimation using the Enhanced RUL Model

If a new data vector $\textbf{x}_{new}(J\times 1)$ from a non-cyclic component
or a new cycling data matrix $\textbf{X}_{new}(J\times K_{new})$ from cyclic
component is upcoming, all variables are first normalized into the range
$[0,1]$ with the normalization information obtained in Eq. 3, in which
$K_{new}$ is the cycle length of the new cycle. For clarity, the normalized
data information is still denoted as $\textbf{x}_{new}$ and
$\textbf{X}_{new}$, respectively.

With well-trained model ${\phi_{nc}(\cdot)}$ and ${\phi_{c}(\cdot)}$, the
estimated RUL is calculated by inputting $\textbf{x}_{new}$ and
$\textbf{X}_{new}$ below,

$RUL_{new}=\left\\{\begin{array}[]{ll}{\phi_{nc}}(\textbf{x}_{{new}})&\text{Non-
cyclic degradation}\\\ {\phi_{c}}(\textbf{X}_{{new}})&\text{Cyclic
degradation}\end{array}\right.$ (7)

where $RUL_{new}$ is the remaining times of non-cyclic component or the
remaining cycles of a cyclic component.

To measure performance of the enhanced RUL estimation model, the indices
scoring function ($SF$) [5] and root mean square error ($RMSE$) [28] [29] are
used. $RMSE$ equally weights both early estimation and delay estimation, which
is defined as follows,

$RMSE=\sqrt{\frac{1}{N_{t}}\sum_{c=1}^{N_{t}}(d_{c})^{2}}$ (8)

where $d_{c}$ is the difference between $RUL_{new}$ and the real useful life
left, and $N_{t}$ is the number of components in non-cyclic degradation or
testing cycles in cyclic degradation.

Specifically, index $SF$ is a piece-wise function calculated as follows,

$SF=\left\\{\begin{array}[]{ll}\sum_{c=1}^{N_{t}}{(e^{-\frac{d_{c}}{13}}-1)}&{d_{c}<0}\\\
\sum_{c=1}^{N_{t}}{(e^{\frac{d_{c}}{10}}-1)}&{d_{c}\geq 0}\end{array}\right.$
(9)

Index $SF$ gives a large punishment for the case when lagged estimation is
observed. With these two metrics, performance of the enhanced RUL estimation
method can be quantitatively compared with the results without generated data.
Both values of $SF$ and $RMSE$ are the smaller the better. Considering $RMSE$
and $SF$ are overall metrics by averaging or summing all weighted estimation
errors, it is recommended to further analyze the distribution of estimation
errors. For two methods with a similar value of $RMSE$, the method who has
more estimation errors distributed in the low-error zones will be preferred.

##  IV Results and discussions

To demonstrate the efficacy of the proposed framework, three cases are
employed in this section, one of which is a numerical example and the other
two are practical ones. The numerical case focuses on verifying the quality of
the generated time-series from an intuitive view. The other two practical
cases comprehensively evaluate the proposed framework by simultaneously
considering the quality of generated time-series and RUL estimation accuracy
from non-cyclic degradation using aero-engine open dataset [22] and cyclic
degradation with LiB open dataset [23].

For fair comparisons, all simulations are conducted on a desktop equipped with
6 processors of Intel Xeon E5-2620 v4 and a multi-graphics processor unit of
NVIDIA GeForce GTX 1080Ti. Programming language is Python 3.6 with scientific
computing library “Pandas (0.25.3)”, “Numpy (1.16.4)”, “Scipy (1.2.1)”,
“Matplotlib (3.0.3)”, “Scikit-learn (0.20.3)”, “Tensorflow (1.15.0)”, and
“Keras (2.3.1)”.

###  IV-A Time-series Generation with A Numerical Case

Two kinds of simulation signals are constructed as real time-series, including
a sinusoidal wave and a smooth signal. The $i^{th}$ generated sinusoidal time-
series $S_{i}(k)$ is generated with varying amplitude and frequency, which is
given below,

${S_{i}(k)=A_{r}sin(\omega_{i}k+\varphi_{i})}$ (10)

where $\omega_{i}$ is frequency randomly selected from range [1, 5], $A_{i}$
is amplitude randomly picked up from range [0.1, 0.9], and $\varphi_{i}$ is
phase randomly given between [$-\pi,\pi$].

A more complex smooth signal is constructed with time-varying variations.
Repeating the process ten times with the same length of 16, and then the
corresponding cycling time-series are obtained with dimensions of $10\times
1\times 16$. Fig. 5 plots the real waves of the generated time-series using
the proposed CR-GAN for the sinusoidal and the smooth signal, respectively. In
this case, the quality of the generated sine waves is examined by visual
inspection, which is almost the same as given real waves.

Figure 5: Example of real curves and generated curves for (a) sinusoidal
signal and (b) smooth signal.

Figure 6: Performance comparisons between the proposed method and its counterparts for generated variables with (a) Dataset FD001 and (b) Dataset FD003. TABLE I: Description of variables collected in engine system. No. | Symbol | Variable description | Unit
---|---|---|---
1 | $S_{1}$ | Total temperature at fan inlet |  oC
2 | $S_{2}$ | Total temperature at low pressure compressor outlet |  oC
3 | $S_{3}$ | Total temperature at high pressure compressor outlet |  oC
4 | $S_{4}$ | Total temperature at low pressure turbine outlet |  oC
5 | $S_{5}$ | Pressure at fan inlet | Pa
6 | $S_{6}$ | Total pressure in bypass-duct | Pa
7 | $S_{7}$ | Total pressure at high pressure compressor outlet | Pa
8 | $S_{8}$ | Physical fan speed | r/min
9 | $S_{9}$ | Physical core speed | r/min
10 | $S_{10}$ | Engine pressure ratio | 1
11 | $S_{11}$ | Static pressure at high pressure compressor outlet | Pa
12 | $S_{12}$ | Ratio of fuel flow to Ps30 | 1
13 | $S_{13}$ | Corrected fan speed | r/min
14 | $S_{14}$ | Corrected core speed | r/min
15 | $S_{15}$ | Bypass ratio | 1
16 | $S_{16}$ | Burner fuel-air ratio | 1
17 | $S_{17}$ | Bleed enthalpy | 1
18 | $S_{18}$ | Demanded fan speed | r/min
19 | $S_{19}$ | Demanded corrected fan speed | r/s
20 | $S_{20}$ | High pressure turbine coolant bleed | lbm/s
21 | $S_{21}$ | Low pressure turbine coolant bleed | lbm/s

###  IV-B RUL Estimation for Non-cyclic Degradation Pattern

In this section, a more complex multivariate dataset, i.e., commercial modular
aviation propulsion system simulation (C-MAPSS) [22] [30], provided by NASA is
used for further evaluation. C-MAPSS provides complete degradation procedures
of the simulated engine system from normal status to faulty status, which is a
typical non-cyclic degradation pattern. According to the number of operating
conditions, four datasets in C-MAPSS are classified into two groups, one of
which (FD001 and FD003) follows single data distribution and the other one is
multi-distribution (FD002 and FD004). Since the proposed CR-GAN does not
consider the characteristics of multiple distributions, only the datasets.
FD001 and FD003 are chosen for analysis. For both FD001 and FD003, there are
100 training engines and 100 testing engines, each of which contains 21
variables listed in Table I [22]. Besides, real RUL values of training engines
are known as a priori information.

TABLE II: Specific Configuration of Generator in the proposed CR-GAN. Layer | Network | Parameters
---|---|---
1 | Convolution layer | Filter=18, Kernel=2, Strides=1, Activation=Relu
Max-pooling layer | Pool size=2, Strides=2
2 | Convolution layer | Filter=18, Kernel=2, Strides=1, Activation=Relu
Max-pooling layer | Pool size=2, Strides=2
3 | Convolution layer | Filter=18, Kernel=2, Strides=1, Activation=Relu
Max-pooling layer | Pool size=2, Strides = 2
4 | Fully connected layer | Unit =100, Activation =Relu, Dropout=0.2
5 | LSTM layer | Cell=100, Dropout=0.2
6 | Fully connected layer | Unit=1, Activation=Sigmoid
TABLE III: Specific Configuration of Discriminator in proposed CR-GAN Layer | Network | Parameters
---|---|---
1 | LSTM layer | Cell=100, Dropout=0.2
2 | Fully connected layer | Unit=1, Activation=Sigmoid
TABLE IV: Comparison of time complexity between the proposed method CR-GAN with its counterparts. Method | Time complexity
---|---
RCGAN | $O(2W)$
T-CGAN | $O(\sum_{l=1}^{4}(s_{l}^{2}n_{l}m_{l}^{2}))$
CR-GAN (Proposed) | $O(3\sum_{l=1}^{3}(s_{l}^{2}n_{l}m_{l}^{2})+4W)$

####  IV-B1 Multivariate time-series generation

To verify the proposed method quantitatively, a “Train on Synthetic data, Test
on Real data” provided in [21] is employed for quantitative comparison.
Specifically, we use the synthetic data generated by CR-GAN to train a
classifier model and then test it on a set of real examples. Besides, two
competitive algorithms, RCGAN [21] and T-CGAN [20], are employed for
comparison. Considering only seven variables are used for time-series
generation in RCGAN and T-CGAN, the same variables are employed as model input
for a fair comparison. The employed variables in FD001 and FD003 can be
observed from Fig. 6(a) and Fig. 6(b), respectively. The specific
configurations of the proposed CR-GAN are listed in Tables II and III.

Fig. 6(a) summarizes all methods for FD001 using a classification evaluation
index $AUROC$, which signifies how much a model could distinguish different
classes. By checking the generated time-series of each variable, it is
observed that the value of $AUROC$ for CR-GAN is very close to 1, indicating a
high-quality generation ability. For RCGAN, it achieves the best result
(95$\%$) for variable $S_{21}$, however, it is still worse than that of the
proposed method (99$\%$). Here, $S_{21}$ stands for the twenty-first variable,
and the subscript is used in the same way for other variables. Besides, the
result of Variable $S_{7}$ in RCGAN is less than $60\%$, yielding a poor
accuracy. Similar results are observed for T-CGAN, which is also much worse
than the proposed method. The corresponding results of FD003 are summarized in
Fig. 6(b), and almost all values concerning index $AUROC$ are higher than
90$\%$, illustrating that the proposed model performs well.

TABLE V: Comparisons between current RUL methods and their enhanced methods
for engine data.

Dataset | Approach | Indices
---|---|---
$SF$ | $RMSE$
FD001 | TCMN | 1220.00 | 23.57
TCMN with CR-GAN⋆ | 1027.00 | 18.39
Accuracy improvement ($\%$) | 15.82 | 21.98
LSTM | 338.00 | 16.14
LSTM with CR-GAN⋆ | 301.06 | 15.30
Accuracy improvement ($\%$) | 10.93 | 5.21
CapsNet | 276.34 | 12.58
CapsNet with CR-GAN⋆ | 216.17 | 12.49
Accuracy improvement ($\%$) | 21.77 | 0.72
FD003 | TCMN | 1300.00 | 21.70
TCMN with CR-GAN⋆ | 1257.00 | 17.92
Accuracy improvement ($\%$) | 3.31 | 17.42
LSTM | 852.00 | 16.18
LSTM with CR-GAN⋆ | 257.32 | 15.58
Accuracy improvement ($\%$) | 69.80 | 3.71
CapsNet | 283.81 | 11.71
CapsNet with CR-GAN⋆ | 191.10 | 12.06
Accuracy improvement ($\%$) | 32.67 | -2.98

  * $\star$

denotes the proposed methods.

Figure 7: Comparison on dataset FD001 for (a) Estimation error of TCMN and
TCMN with CR-GAN, (b) Error distribution of TCMN and TCMN with CR-GAN, (c)
Estimation error of LSTM and LSTM with CR-GAN, (d) Error distribution of LSTM
and LSTM with CR-GAN, (e) Estimation error of CapsNet and CapsNet with CR-GAN,
and (f) Error distribution of CapsNet and CapsNet with CR-GAN.

With the network architecture presented in Tables II and III, the overall time
complexity of the proposed method is calculated by a linear combination of the
generator and the discriminator. Specifically, when the input channel is one,
the time complexity of the $l^{th}$ convolutional layer is estimated as
$O(s_{l}^{2}n_{l}m_{l}^{2})$ [31], where $n_{l}$ is the number of filters in
the $l^{th}$ layer (18 here), $s_{l}$ is the spatial size of the filter (2
here), and $m_{l}$ is the spatial size of the output feature map. Since the
RNN/LSTM is local in space and time, the complexity of an LSTM/RNN layer per
time step is equal to $O(W)$, where $W$ is the number of weight. Finally, the
time complexity of combining the generator and the discriminator is calculated
as $O(3\sum_{l=1}^{3}(s_{l}^{2}n_{l}m_{l}^{2})+4W)$ through the linear
superposition of all convolutional layers and LSTM layers. The time complexity
of RC-GAN and T-CGAN can be computed in the same way. As shown in Table IV, it
is observed that the time complexity of the proposed method is a litter higher
than that of RC-GAN and T-CGAN since more CNN and LSTM layers are used. In
comparison with RCGAN and T-CGAN, the increased training time of the proposed
method is not significant. Moreover, once the offline model has been well
trained, it could be directly adopted in the online estimation stage with
rapid computing.

####  IV-B2 Enhanced RUL estimation with generated time-series

With the CR-GAN model, performances of the enhanced RUL estimation is verified
by feeding the generated time-series into three typical RUL estimation methods
[5], [7], [9] with model re-training. The same model structure listed in
Tables II and III are adopted.

Figure 8: Comparison on dataset FD003 for (a) Estimation error of TCMN and
TCMN with CR-GAN, (b) Error distribution of TCMN and TCMN with CR-GAN, (c)
Estimation error of LSTM and LSTM with CR-GAN, (d) Error distribution of LSTM
and LSTM with CR-GAN, (e) Estimation error of CapsNet and CapsNet with CR-GAN,
and (f) Error distribution of CapsNet and CapsNet with CR-GAN.

It should be noted that all variables are used for RUL estimation modeling.
This is different from the purpose of data generation in the last subsection,
which only considers fair comparison with existing methods. First, data
preparation is arranged to remove redundant variables from both training and
testing data. By checking mean and variance of each variable, variables with
constant value are Variables $S_{1},S_{5},S_{6},S_{10},S_{16},S_{18},S_{19}$
in FD001 and Variables $S_{1},S_{5},S_{10},S_{16},S_{18},S_{19}$ in FD003,
which are non-informative. Next, the degradation of the engine system
typically starts after a certain degree of usage, resulting in the less
importance of the normal data samplings. According to apriori knowledge, the
samplings in the last 125 hours in each cycle have remained for data
generation, and the window size $L$ is chosen as 100 here. In this way, the
number of original time-series is enlarged 26 times to enable basic training
of CR-GAN. With data normalization in the way given in Eq. 6, the inputs of
CR-GAN is normalized data matrixes with dimensions $2600\times 14\times\ 100$
and dimensions $2600\times 15\times\ 100$ for FD001 and FD003, respectively.

TABLE VI: RUL estimation error comparison between the enhanced LSTM and
original LSTM for cyclic degradation battery data.

Battery Name | Method | Number of Generated Time-series | Actual RUL | Predicted RUL | RUL Error | RMSE
---|---|---|---|---|---|---
B0005 | LSTM | – | 74 | 94 | 20 | 0.113
LSTM with CR-GAN⋆ | 50 | 78 | 4 | 0.0822
100 | 76 | 2 | 0.0814
150 | 73 | -1 | 0.0563
B0006 | LSTM | – | 58 | 106 | 48 | 0.1216
LSTM with CR-GAN⋆ | 50 | 77 | 19 | 0.0852
100 | 73 | 15 | 0.0699
150 | 57 | -1 | 0.0722
B0018 | LSTM | – | 46 | 66 | 20 | 0.0607
LSTM with CR-GAN⋆ | 50 | 46 | 0 | 0.0441
100 | 46 | 0 | 0.0343
150 | 46 | 0 | 0.0239

  * $\star$

denotes the proposed methods.

(a) Battery B0005

(b) Battery B0006

(c) Battery B0018

Figure 9: Accuracy comparison between LSTM and its enhanced framework for (a)
Battery B0005, (b) Battery B0006, and (c) Battery B0018.

After filtering out irrational time-series, a total number of 200 run-to-
failure time-series with a length of 100 is generated using CR-GAN. These
synthetic time-series are mixed with original time-series to form a new
training dataset for the RUL estimation model. Here, three baseline RUL
estimation methods, LSTM [5], temporal convolution neural network (TCMN) [7],
and CapsNet [9], are employed for comparison. With the given steps through
Sections III. A to III. E, Table V summarizes the estimation results of all
enhanced methods using two criteria, including $SF$ and $RMSE$. For the
enhanced TCMN, the value of $RMSE$ is reduced from 23.57 to 18.39 for FD001
and from 21.70 to 17.92 for FD003, respectively. With respect to $SF$, the
value is also significantly reduced using the proposed framework. For the
enhanced LSTM algorithm, the estimation error of $RMSE$ is reduced from 16.14
to 15.30 for FD001 and from 16.18 to 15.58 for FD003, respectively. CapsNet
[32] proposed by Hinton in 2017 aims at extracting vector-shaped features to
model the hierarchical relationship. The state-of-art results of $RMSE$ index
reported by CapsNet are 12.59 and 11.71 for FD001 and FD003, respectively [9].
The proposed framework helps CapsNet further reducing the value of the index
$SF$ from 276.34 to 216.17 for FD001 and from 283.81 to 191.10 for FD003,
respectively. More precisely, in Table V, improvements are calculated in the
ratio that the reduced errors are divided by the values of its reference
model.

As an auxiliary evaluation to $RMSE$ and $SF$, the estimation error of the
specific engine and the corresponding distribution are shown in Figs. 7 and 8
for FD001 and FD003, respectively. The engine index here is sorted according
to the ascending of real RUL. It provides an intuitive way to observe the
performance by analyzing each engine. It is observed that adding generated
time-series will drive the estimations closer to their ground truth. Further,
the distribution of estimation errors is shown in the histogram for the first
time to account for the percentage of error band with a specific interval,
which is four here. For the similar values of $RMSE$ and $SF$, a method is
regarded better if more estimation errors are distributed in the low-value
zones. In Figs. 7(b), 7(d), and 7(f), it is observed that error distributions
of the proposed method for a higher ratio in low-value zones in comparison
with their counterparts. A similar conclusion can be drawn for FD003 through
analyzing Figs. 8(b), 8(d), and 8(f).

###  IV-C RUL Estimation for Cyclic Degradation Pattern

As a reusable component, LiB suffers from performance degradation due to
inherent electrochemical reactions during repetitive usage. A LiB needs to be
replaced when its real capacity drops to a certain threshold, i.e., EOL, which
is usually determined as 70$\%$ of its rated capacity [33]. Correspondingly,
RUL is calculated as the difference between predicted EOL and the real value.
In practice, discharging current and discharging voltage could be measured for
capacity inference for RUL estimation. Having the data characteristics of
cycling run-to-failure time-series, LiB is an ideal benchmark to test the
performance of the proposed method concerning cyclic degradation pattern.
Here, LiB data provided by NASA Ames Prognostics Center of Excellence [23] are
employed for testing the performance of the enhanced RUL framework.

####  IV-C1 Multivariate time-series generation

Datasets B0005, B0006, and B0018 are selected for analysis, which have been
used in [34] that adopts LSTM for RUL estimation. For cycling data in each
dataset, current signal, voltage signal, and ground truth capacity are
employed as the generated variables, and the time index is used as the
condition information. In [34], two datasets and the first 50 cycles in the
third one are used for training in turn, and the remaining ones are used for
testing. According to steps given in Section III, 150 cycling data located in
the first 50 cycles are generated. Here, the same configurations listed in
Tables II and III are adopted with $L=200$.

####  IV-C2 Enhanced RUL estimation with generated time-series

Employing the same LSTM model structure in [34], the RUL estimation ability of
the proposed framework is further tested. Taking current and voltage signals
as input variables, the estimation model outputs the prediction of capacity.
Moreover, to reveal the influences of data amount on RUL estimation accuracy,
three different amounts of generated time-series are considered, including 50,
100, and 150. According to Table VI, two conclusions can be clearly drawn.
First, both $RMSE$ and RUL errors have been significantly reduced with the
introduction of generated time-series. Especially, RUL estimation error has
been reduced to zero for Battery B0018 when the number of generated time-
series is 150 cycles. Second, the more generated data used for model re-
training, the smaller RUL error will be, which can be observed from the value
of index $RMSE$. To intuitively show the estimation results, Fig. 9 compares
LSTM and its enhanced model when 150 generated cycles are used for model
training. The y-label of the top subplot of Fig. 9 presents the state of
health (SoH), which indicates the health status of a battery by informing the
remaining capacity or point that has been reached in the life cycle. Besides,
the absolute error distribution of each model with respect to SoH is shown in
Fig. 9 using the boxplot.

From the above analysis for both types of degradation systems, RUL estimation
errors have significantly been reduced with the aid of regeneration of run-to-
failure time-series, demonstrating the efficacy and success of the proposed
framework.

##  V Conclusion and Future work

To effectively address the challenge of data shortage of run-to-failure time-
series, this article puts forward the idea of time-series generation for the
first time in improving RUL estimation. It is a self-learning enhancement
framework aiming at integrating time-series generation into current RUL
estimation models. The specificity of this article is to synthesize realistic
multivariate time-series with a novel network named convolutional recurrent
conditional generative adversarial network (CR-GAN), yielding a wealth of data
to escalate the training procedure of existing RUL models. Moreover, the
degradation nature of various systems is well studied, producing the
classification of non-cyclic and cyclic degradation patterns. By comparing the
current RUL estimation models with their enhanced versions, practical
experiment results on both aero-engine system, and Lithium-ion aging clearly
demonstrate the strengths of the proposed framework in terms of diverse
evaluation indices. Therefore, it would be interesting to apply this idea for
more industrial scenarios where it is economically expensive or dangerous to
gain enough run-to-failure data.

This article verifies the feasibility of employing CR-GAN for generating time-
series for better RUL estimation. However, there are still some possible
directions deserving more efforts:

1) Run-to-failure time-series possess plenty of statistical properties, such
as slow-varying trends, etc. It is possible that the objective function in the
original GAN is not able to preserve all the statistics of the real time-
series. It may be promising to consider proper statistical constraints in the
objective function of the original GAN to obtain more reliable time-series.

2) Gaussian distribution has been widely selected to present the given input
noise vector in the image processing field. Although Gaussian distribution is
continually adopted for time-series generation, it is interesting to
investigate the performance when other noise statistics are concerned.

3) Although irrational generated time-series could be filtered out using a
post-processing way, it is recommended to consider a direct way during
adversarial learning. Especially, irrational time-series may be easily
recognized by experienced engineers equipped with plentiful domain knowledge.
Therefore, it may provide a feasible solution to filter out irrational time-
series with a discriminator considering domain knowledge.

## References

  * [1] Y.G. Lei, N.P. Li, L. Guo, N.B. Li, T. Yan, and J. Lin, “Machinery health prognostics: A systematic review from data acquisition to RUL prediction,” Mech. Syst. Signal Process., vol. 104, pp. 799-834, 2018.
  * [2] A.K.S. Jardine, D. Lin, and D. Banjevic, “A review on machinery diagnostics and prognostics implementing condition-based maintenance,” Mech. Syst. Signal Process., vol. 20, no. 7, pp. 1483-1510, 2006.
  * [3] R. Liu, B. Yang, and A. G. Hauptmann, “Simultaneous bearing fault recognition and remaining useful life prediction using joint loss convolutional neural network,” IEEE Trans. Ind. Informat., vol. 16, no. 1, pp. 87-96, 2020.
  * [4] H. Miao, B. Li, C. Sun, and J. Liu, “Joint learning of degradation assessment and RUL prediction for aeroengines via dual-task deep LSTM networks,” IEEE Trans. Ind. Informat., vol. 15, no. 9, pp. 5023-5032, 2019.
  * [5] S. Zheng, K. Ristovski, A. Farahat, and C. Gupta, “Long short-term memory network for remaining useful life estimation,” in Proc. IEEE Int. Conf. Progn. Health Manag., 2017, pp. 88-95.
  * [6] Y. Wu, M. Yuan, S. Dong, L. Lin, and Y. Liu, “Remaining useful life estimation of engineered systems using vanilla LSTM neural networks,” Neurocomputing, vol. 275, pp. 167-179, 2018.
  * [7] L. Jayasinghe, T. Samarasinghe, C. Yuen, J.C.N. Low, and S.S. Ge, “Temporal convolutional memory networks for remaining useful life estimation of industrial machinery,” in IEEE Proc. IEEE Int. Conf. Ind. Technol., 2019, pp. 915-920.
  * [8] Z. Chen, M. Wu, R. Zhao, F. Guretno, R. Yan, and X. Li, “Machine remaining useful life prediction via an attention based deep learning approach,” IEEE Trans. Ind. Electron., doi: 10.1109/TIE.2020.2972443.
  * [9] R.P. Andrés, L.D. Enrique, and P. Rodrigo, “A novel deep capsule neural network for remaining useful life estimation,” Proc. Inst. Mech. Eng. O J. Risk Reliab., vol. 234, no. 1, pp. 151-167, 2020.
  * [10] Y. Zhang, R. Xiong, H. He, and M.G. Pecht, “Long short-term memory recurrent neural network for remaining useful life prediction of Lithium-ion batteries,” IEEE Trans. Veh. Technol., vol. 67, no. 7, pp. 5695-5705, 2018.
  * [11] K.A. Severson, P.M. Attia, N. Jin, et al., “Data-driven prediction of battery cycle life before capacity degradation,” Nat. Energy, vol. 4, pp. 383-391, 2019.
  * [12] R. Razavi-Far, S. Chakrabarti, M. Saif, and E. Zio, “An integrated imputation-prediction scheme for prognostics of battery data with missing observations,” Expert Syst. Appl., vol. 115, pp. 709-723, 2019.
  * [13] R. Razavi-Far, S. Chakrabarti, M. Saif, E. Zio, and V. Palade, “Extreme learning machine based prognostics of battery life,” Int. J. Artif. Intell. Tools, vol. 27, no. 8, 2018.
  * [14] K. Zhang, G. Zhong, J. Dong, S. Wang, and Y. Wang, “Stock market prediction based on generative adversarial network,” Procedia Comput. Sci., vol. 147, pp. 400-406, 2019.
  * [15] M. Rüttgers, S. Lee, S. Jeon, and D. You, “Prediction of a typhoon track using a generative adversarial network and satellite images,” Sci. Rep., vol. 9, no. 6057, pp. 1-15, 2019.
  * [16] Y. Kwon and M. Park, “Predicting future frames using retrospective cycle GAN,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., 2019, pp. 1811-1820.
  * [17] I.J. Goodfellow, J.P. Abadie, M. Mirza, et al, “Generative adversarial networks,” in Proc. Int. Conf. Neural Inf. Process Syst., 2014, pp. 2672–2680.
  * [18] L.T. Yu, W.N. Zhang, J. Wang, and Y. Yu, “Seqgan: Sequence generative adversarial nets with policy gradient,” in Proc. AAAI Conf. Artif. Intell., 2017, pp. 2852-2858.
  * [19] J. Li, W. Monroe, T. Shi, S. Jean, A. Ritter, and D. Jurafsky, “Adversarial learning for neural dialogue generation,” in Proc. Conf. Empirical Methods in Natural Language Processing, 2017, pp. 2157-2169.
  * [20] G. Ramponi, P. Protopapas, M. Brambilla, and R. Janssen, “T-CGAN: conditional generative adversarial network for data augmentation in noisy time series with irregular sampling,” ArXiv e-preprint arXiv:1811.08295, 2018.
  * [21] C. Esteban, S.L. Hyland, and G. Rätsch, “Real valued (Medical) time series generation with recurrent conditional GANs,” ArXiv e-prints arXiv:1706.02633, 2017.
  * [22] E. Ramasso and A. Saxena, “Performance benchmarking and analysis of prognostic methods for CMAPSS datasets,” Int. J. Prog. Health Manage., vol. 5, pp. 1-15, 2014.
  * [23] B. Saha and K. Goebel, “Battery data set,” NASA Ames Prognostics Data Repository (http://ti.arc.nasa.gov/project/prognostic-data-repository), NASA Ames Research Center, 2007.
  * [24] M. Kahng, N. Thorat, D.H. Chau, F.B. Viégas, and M. Wattenberg, “GAN Lab: understanding complex deep generative models using interactive visual experimentation,” IEEE Trans. Vis. Comput. Graph., vol. 25, no. 1, pp. 310-320, 2019.
  * [25] C. You, G. Li, Y. Zhang, X. Zhang, and H. Shan, “CT Super-resolution GAN constrained by the identical, residual, and cycle learning ensemble (GAN-CIRCLE),” IEEE Trans. Medical Imaging, vol. 39, no. 1, pp. 188-203, 2020.
  * [26] H. Zhang, T. Xu, H. Li, S. Zhang, X. Wang, X. Huang, and D. Metaxas, “StackGAN: text to photo-realistic image synthesis with stacked generative adversarial networks,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., 2017, pp. 5908-5916.
  * [27] T. Karras, S. Laine, and T. Aila, “A style-based generator architecture for generative adversarial networks,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., 2019, pp. 4401-4410.
  * [28] W. Shao, Z. Ge, Z. Song, and J. Wang, “Semi-supervised robust modeling of multimode industrial processes for quality variable prediction based on student’s t mixture model,” IEEE Trans. Ind. Informat., vol. 16, no. 5, pp. 2965-2976, 2020.
  * [29] B. Song, X.G. Zhou, H.B. Shi, and Y. Tao, “Performance indicator oriented concurrent subspace process monitoring method,” IEEE Trans. Ind. Electron., vol. 66, pp. 5535-5545, 2019.
  * [30] W. Peng, Z. Ye, and N. Chen, “Joint online RUL prediction for multivariate deteriorating systems,” IEEE Trans. Ind. Informat., vol. 15, no. 5, pp. 2870-2878, 2019.
  * [31] E. Tsironi, P. Barros, C. Weber, and S. Wermter, “An analysis of convolutional long short-term memory recurrent neural networks for gesture recognition,” Neurocomputing, vol. 268, pp. 76-86, 2017.
  * [32] S. Sabour, N. Frosst, and G.E. Hinton,“ Dynamic routing between capsules,” in Prof. Int. Conf. Neural Inf. Process Syst., 2017, pp. 1-11.
  * [33] E. Chemali, P.J. Kollmeyer, M. Preindl, R. Ahmed, and A. Emadi, “Long short-term memory networks for accurate state-of-charge estimation of Li-ion batteries,” IEEE Trans. Ind. Electron., vol. 65, no. 8, pp. 6730-6739, 2018.
  * [34] J. Qu, F. Liu, Y. Ma, and J. Fan, “A neural-network-based method for RUL prediction and SOH monitoring of Lithium-Ion battery,” IEEE Access, vol. 7, pp. 87178-87191, 2019.