text:
\section{Introduction}

A {\bf complex sequence} of {\bf length} $N$ is an $N$-tuple
$$
\mathcal{X} = (\xi_0, \xi_1, \ldots, \xi_{N-1})
$$
where each $\xi_j$ is a complex number of modulus 1.
If furthermore each $\xi_j$ is a $D$th root of unity,
then the sequence is called a {\bf $D$-phase sequence},
or {\bf polyphase sequence} if we don't specify the value of $D$.
Throughout this article, lowercase Greek letters will denote complex numbers
and nonbold Latin letters will denote nonnegative integers, unless otherwise stated.

If $0 \le k \le N-1$, we define the {\bf acyclic autocorrelations} of $\mathcal{X}$ by
$$
\alf_k = \sum_{j=0}^{N-k-1} \ol{\xi_j} \xi_{j+k}
$$
and we define the {\bf cyclic autocorrelations} of $\mathcal{X}$ by
$$
\gam_k = \sum_{j=0}^{N-1} \ol{\xi_j} \xi_{j+k}
$$
where the bar denotes complex conjugation, and in the definition
of $\gam_k$, the addition in the subscript is mod $N$, or equivalently,
we assume $\xi_j$ is defined for $j \ge N$ by $\xi_{j+N}=\xi_j$.

Note that $\alf_k$ is a sum of $N-k$ terms, each of modulus 1.
Hence $\abs{\alf_k} \le N-k$.
Also note that if $1 \le k \le N-1$, we have $\gam_k = \alf_k + \ol{\alf_{N-k}}$.
It follows that if $\gam_k=0$, then
$\abs{\alf_k} = \abs{\alf_{N-k}} \le k$.

We can regard $\alf_k$ and $\gam_k$ as measuring resemblance between
the sequence $\mathcal{X}$ and a version of $\mathcal{X}$ that has been shifted by
$k$ positions (acyclically or cyclically respectively).
We have $\alf_0=\gam_0=N$, which we may call the {\bf trivial} autocorrelations.
Informally, we consider a sequence to be ``good'' if its nontrivial
autocorrelations are close to 0
(so it is ``uncorrelated'' with shifted versions of itself.)

How close to 0 can we make the cyclic autocorrelations,
and how close to 0 can we make the acyclic autocorrelations?
The first of those questions has an easier answer than the second.

We call $\mathcal{X}$ a {\bf perfect} sequence if $\gam_k=0$ for all $k\ne0$.
Many families of perfect sequences have been studied,
including Frank sequences \cite{FZ},
Chu sequences \cite{Chu},
and Milewski sequences \cite{Mil}.
There are perfect sequences of every length $N\ge2$.
For a good recent survey of perfect sequences,
see \cite[Section~4.1]{Sch16}.

We define the {\bf peak sidelobe level} (abbreviated ``PSL'') of $\mathcal{X}$ by
$$
\mathbf{P}(\mathcal{X}) = \max_{1\le k\le N-1} \abs{\alf_k},
$$
and we define the {\bf energy} of $\mathcal{X}$ by
$$
\mathbf{E}(\mathcal{X}) = \sum_{k=1}^{N-1} \abs{\alf_k}^2,
$$
which are two natural measures of the ``size'' of the acyclic autocorrelations.
A complex sequence with PSL at most 1
(i.e., with $\abs{\alf_k}\le1$ for all $k\ne0$)
is called a {\bf generalized Barker sequence}.
There exist generalized Barker sequences of all lengths $N\le70$
(see \cite{NC}), and it has been conjectured that they exist for all lengths.
See \cite[Question~4.10]{Sch16}
or \cite[p.~119]{Bor}.

A generalized Barker sequence of length $N$ has energy at most $N-1$.
But if we want an infinite family of complex sequences of increasing lengths $N$
that have small energy, the best known infinite families
(which are the Chu sequences and Frank sequences)
have energy growing like $O(N^{3/2})$. See \cite{Sch13}.

The {\bf merit factor} of a length $N$ complex sequence $\mathcal{X}$
is defined by
$$
\mathbf{F}(\mathcal{X}) = \frac{N^2}{2\mathbf{E}(\mathcal{X})}
$$
so asking for small energy is equivalent to asking for
large merit factor. The energy and merit factor are related to
the $L^4$ norm on the unit circle of the polynomial whose
coefficients are the $\xi_j$ (specifically, minimizing the energy
of the sequence is equivalent to minimizing the $L^4$ norm of the
polynomial). See \cite[Chapter~15]{Bor}.

We will construct a family of perfect polyphase sequences
that contains the Frank sequences, Chu sequences, and Milewski sequences
as special cases.
This is not the most general construction of this type.
In 1995, Mow provided a construction that
included all known infinite families of perfect polyphase sequences \cite{Mow}.
For more discussion of families of perfect polyphase sequences, see \cite{FD}.

Our family of sequences is not as general as Mow's,
but it includes several well-known families in one
surprisingly simple form.
We will refer to our sequences as ``LM sequences'', where LM could stand
for ``like Mow'' but also names two integer parameters we use.

Some facts are already known about acyclic autocorrelation
of families of perfect sequences. Turyn \cite{Tur}
showed that the PSL of the Frank sequence of length $N$
is asymptotically equal to $(1/\pi)\sqrt{N}$.
Mow and Li \cite{ML} showed that the PSL of the Chu
sequence of length $N$ is asymptotically equal to $C\sqrt{N}$
for a slightly larger constant.
The current author \cite{Mer} showed that the energy
of the Chu sequence of length $N$ is bounded above by
$(8/3\pi^{3/2})N^{3/2}$; this was improved by Schmidt \cite{Sch13}
who showed that the energy of the Chu sequence of length $N$
is asymptotically equal to $(1/\pi)N^{3/2}$,
and the energy of the Frank sequence of length $N$
is asymptotically equal to $(2/\pi^2)N^{3/2}$.
It appears that less is known about the PSL or energy
of Milewski sequences. For a good summary of acyclic autocorrelation
of polyphase sequences, see \cite[Section~4.2]{Sch16}.

In this article, we do not say much about the acyclic autocorrelations
of our sequences. As consequences of the proofs that our sequences are
perfect sequences, it will follow that the Frank sequence of length $N$
has PSL at most $\sqrt{N/2}$, and the Chu sequence of length $N$ has
PSL at most $\sqrt{N}$, which are slightly weaker than known results.
The current author conjectures that there exist polyphase sequences of
all lengths $N$ whose energy grows like $o(N^{3/2})$.
Perhaps further study of acyclic autocorrelations of general families
of perfect sequences will prove this conjecture.

We note that if $\mathcal{X}$ is a $D$-phase sequence, then each $\xi_j$ can be written as
$\zet_D^{p(j)}$, where $\zet_D = e^{2\pi i/D}$ and $p(j)$ belongs to the integers mod $D$.
We can specify the sequence by specifying the
values of $p(j)$.

\section{LM sequences}

Throughout the rest of this article, $L$ and $M$ are positive integers
such that $L$ divides $M$, and $N$ denotes $LM$.
Since $N=1$ is uninteresting, we assume $M\ge2$.
We will construct a $2M$-phase sequence of length $N$,
informally consisting of $M$ ``blocks'' each of length $L$.
As before, $\zet_D$ means $e^{2\pi i/D}$.

For all nonnegative integers $j$, we define
$$
s_j = \left\lfloor \frac{j}{L} \right\rfloor \qquad \mbox{and} \qquad
t_j = j - L\left\lfloor \frac{j}{L} \right\rfloor.
$$
Then $j=s_jL+t_j$.
One can verify that $s_{j+N} = s_j + M$ and $t_{j+N} = t_j$.

{\bf Proposition 2.1.} Let $L,M$ be as above, and let $A$ be an
integer with the same parity as $LM$.
Define a function on the nonnegative integers as follows:
$$
p(j) = 2s_jt_j + Ls_j^2 + As_j.
$$
Then $p$ is an integer-valued function that satisfies $p(j+N) \equiv p(j)$ mod $2M$.

{\it Proof.} We have
\begin{align*}
p(j+N) & = 2s_{j+N}t_{j+N} + Ls_{j+N}^2 + As_{j+N} \\
& = 2(s_j+M)t_j + L(s_j+M)^2 + A(s_j+M) \\
& = 2s_jt_j + 2Mt_j + Ls_j^2 + 2LMs_j + LM^2 + As_j + AM \\
& = p(j) + 2M(t_j+Ls_j) + M(LM+A)
\end{align*}
and since $LM+A$ is even, this proves the proposition.

It follows that if we define, for all $j$,
$$
\xi_j = \zet_{2M}^{p(j)},
$$
then we have $\xi_{j+N} = \xi_j$ for all $j$.

{\bf Definition 2.2.} We define the {\bf LM sequence} of length $N=LM$ to be
$$
\mathcal{X} = (\xi_0,\xi_1,\ldots,\xi_{N-1})
$$
where $\xi_j=\zet_{2M}^{p(j)}$ and $p(j)$ is as defined previously.

The sequence depends on our choice of $L$, $M$, and $A$.
(Recall we must have $L|M$, and $A \equiv LM$ mod 2.)
We single out some important special cases.

{\bf Special case (i).} If $M$ is even, choose $A=M$,
and if $M$ is odd, choose $A=L$. We then have
$$
p(j) = \begin{cases}
2s_jt_j + Ls_j^2 + Ms_j & \mbox{if $M$ is even}, \\
2s_jt_j + Ls_j^2 + Ls_j & \mbox{if $M$ is odd}.
\end{cases}
$$
{\bf Special case (ii).} If $M$ is even, choose $A=0$,
and if $M$ is odd, choose $A=L$. We then have
$$
p(j) = \begin{cases}
2s_jt_j + Ls_j^2 & \mbox{if $M$ is even}, \\
2s_jt_j + Ls_j^2 + Ls_j & \mbox{if $M$ is odd}.
\end{cases}
$$
{\bf Special case (i)a.} Consider the subcase of special case (i) where
$M=L$. Then $N=M^2$, and we have
$$
p(j) = 2s_jt_j + Ms_j^2 + Ms_j = 2s_jt_j + Ms_j(s_j+1).
$$
Since $s_j(s_j+1)$ is even, we have $p(j) \equiv 2s_jt_j$ mod $2M$.
We then have
$$
\xi_j = \zet_{2M}^{p(j)} = \zet_{2M}^{2s_jt_j} = \zet_M^{s_jt_j}
$$
which means $\mathcal{X}=(\xi_0,\ldots,\xi_{N-1})$ is the sequence of
length $M^2$ defined by
$$
\xi_j = \xi_{s_jM+t_j} = \zet_M^{s_jt_j}
$$
which is equivalent to the Frank sequence of length $M^2$ as defined in \cite{FZ} or \cite[Section~4.1]{Sch16}.

{\bf Special case (ii)a.} Consider the subcase of special case (ii) where
$L=1$. Then $N=M$, and we have
$$
p(j) = \begin{cases}
2s_jt_j + s_j^2 & \mbox{if $M$ is even}, \\
2s_jt_j + s_j^2 + s_j & \mbox{if $M$ is odd}.
\end{cases}
$$
If $L=1$, then $s_j=\lfloor j/1 \rfloor=j$ and $t_j=j-1j=0$, so the above becomes
$$
p(j) = \begin{cases}
j^2 & \mbox{if $M$ is even}, \\
j^2 + j & \mbox{if $M$ is odd}.
\end{cases}
$$
This means $\mathcal{X}=(\xi_0,\ldots,\xi_{N-1})$ is the sequence of length $M$
defined by
$$
\xi_j = \zet_{2M}^{p(j)} = \begin{cases}
\zet_{2M}^{j^2} & \mbox{if $M$ is even}, \\[1ex]
\zet_{2M}^{j^2+j} = \zet_M^{(j^2+j)/2} & \mbox{if $M$ is odd},
\end{cases}
$$
which is equivalent to the Chu sequence of length $M$ as defined in \cite{Chu} or \cite[Section~4.1]{Sch16}.

{\bf Special case (ii)b.} Consider the subcase of special case (ii) where
$L=G^H$ and $M=G^{H+1}$, where $H>0$. Then $N=G^{2H+1}$, and we have
$$
p(j) = \begin{cases}
2s_jt_j + G^H s_j^2 & \mbox{if $G$ is even}, \\
2s_jt_j + G^H (s_j^2+s_j) & \mbox{if $G$ is odd}.
\end{cases}
$$
This means $\mathcal{X}=(\xi_0,\ldots,\xi_{N-1})$ is the sequence of length $G^{2H+1}$
defined by
$$
\xi_j = \xi_{s_jG^H+t_j} = \begin{cases}
\zet_{2M}^{2s_jt_j + G^H s_j^2}
= \zet_{M}^{s_jt_j + \frac{G^H}{2} s_j^2} & \mbox{if $G$ is even}, \\[1ex]
\zet_{2M}^{2s_jt_j + G^H (s_j^2+s_j)}
= \zet_{M}^{s_jt_j + G^H\frac{s_j^2+s_j}{2}} &\mbox{if $G$ is odd},
\end{cases}
$$
or equivalently,
$$
\xi_j = \xi_{s_jG^H+t_j} =
\begin{cases}
\exp\big(\frac{2\pi i}{G^{H+1}}(s_jt_j+\frac{G^H}{2}s_j^2)\big)
& \mbox{if $G$ is even}, \\[1ex]
\exp\big(\frac{2\pi i}{G^{H+1}}(s_jt_j+G^H\frac{s_j^2+s_j}{2})\big)
& \mbox{if $G$ is odd},
\end{cases}
$$
which is equivalent to the Milewski sequence of length $G^{2H+1}$
as defined in \cite{Mil} or \cite[Section~4.1]{Sch16}.

\section{Useful lemmas}

As always, $\zet_D$ denotes $e^{2\pi i/D}$.
In this section, $x$ denotes a real number.

{\bf Lemma 3.1.} If $a$ and $k$ are any integers and $k$ is not a multiple of $D$, then
$$
\sum_{j=a}^{a+D-1} \zet_D^{kj} = 0.
$$
Proof: This sum of $D$ terms is invariant under multiplication by $\zet_D^k\ne1$.

Lemma. We have $\abs{1-e^{ix}}=2\abs{\sin\frac{x}{2}}$.

Proof: $\abs{1-e^{ix}}^2 = (1-\cos x)^2+\sin^2x = 2-2\cos x = 4\sin^2\frac{x}{2}$.

Lemma. If $0 < x \le\frac\pi2$ then
$\csc x \le \frac{\pi}{2x}$.

Proof: If $0 < x \le\frac\pi2$, we have
$\sin x \ge \frac{2}{\pi}x > 0$,
and taking reciprocals proves the lemma.

Lemma. If $\om=e^{ix}\ne1$, and $a<b$ are integers, then
$$
\abs{\om^a+\om^{a+1}+\cdots+\om^{b-1}} \le \abs{\csc\frac{x}2}.
$$
Proof:
If $S = \om^a + \om^{a+1} + \cdots + \om^{b-1}$, we have $S-\om S = \om^a-\om^b$
and so
$$
\abs{S} = \frac{\abs{\om^a-\om^b}}{\abs{1-\om}}
= \frac{\abs{e^{iax}-e^{ibx}}}{\abs{1-e^{ix}}}
\le \frac{2}{\abs{1-e^{ix}}} = \frac{2}{2\abs{\sin\frac{x}{2}}}
= \abs{\csc\frac{x}{2}}.
$$
{\bf Corollary 3.2.} If $k$ is not a multiple of $D$, then taking
$\om=\zet_D^k=e^{i(2\pi k/D)}\ne1$, we get
$$
\abs{\sum_{j=a}^{b-1} \zet_D^{kj}} \le \abs{\csc\frac{\pi k}{D}}.
$$
If furthermore we have $0 < k < D$, this becomes
$$
\abs{\sum_{j=a}^{b-1} \zet_D^{kj}} \le \csc\frac{\pi k}{D}.
$$
Definition. For positive integers $k$ and $D$,
let $\delta(k,D)$ be the distance from $k$
to the nearest multiple of $D$.
So if $k$ belongs to an interval of the
form $[jD,(j+\frac12)D]$, then
$k = jD+\delta(k,D)$, and if $k$ belongs
to an interval of the form
$[(j-\frac12)D,jD]$, then
$k = jD-\delta(k,D)$.
For example, $\delta(19,12)=\delta(29,12)=5$.

Fact. The function
$f(x)=\abs{\csc\frac{\pi x}{D}}$
is periodic with period $D$ and is symmetric
on the interval $[0,D]$, i.e.\ $f(D-x)=f(x)$.

{\bf Corollary 3.3.}
If $k$ is not a multiple of $D$,
and $k' = \delta(k,D)$, then
$$
\abs{\csc\frac{\pi k}{D}} = \csc\frac{\pi k'}{D}
$$
and furthermore, since
$0<k'\le\frac{D}{2}$, we have
$$
\abs{\csc\frac{\pi k}{D}} \le \frac{D}{2k'}.
$$

\section{Autocorrelation of LM sequences}

Throughout this section, $\mathcal{X} = (\xi_0, \ldots, \xi_{N-1})$
is the LM sequence of length $N$ defined previously.
So $L$, $M$, $A$, and $p(j)$ are as defined previously, and the autocorrelations
$\gam_k$ and $\alf_k$ are sums of terms of the form
$$
\ol{\xi_j} \xi_{j+k} = \zet_{2M}^{p(j+k)-p(j)}.
$$
Suppose $1 \le k \le N-1$. We want to show $\gam_k=0$,
and we want bounds on the size of $\alf_k$.

Note that either $k$ and $N-k$ are both multiples of $L$,
or $k$ and $N-k$ are both nonmultiples of $L$.

{\bf Proposition 4.1.} Let $\mathcal{X}=(\xi_0,\ldots,\xi_{N-1})$
be the LM sequence of length $N=LM$ defined previously.
Suppose $1 \le k \le N-1$, and suppose $k$ and $N-k$ are both
multiples of $L$. Then the autocorrelations $\gam_k$ and $\alf_k$
satisfy $\gam_k=0$ and $\abs{\alf_k} \le \sqrt{N/2}$.

{\it Proof.}
If $k$ and $N-k$ are both multiples of $L$, let $k = qL$
and let $N-k = rL$. So $q+r=M$, and $1 \le q \le M-1$.
We break the sums $\gam_k$ and $\alf_k$ into sums of $L$ terms:
\begin{align*}
\gam_k & = \sum_{j=0}^{ML-1} \ol{\xi_j} \xi_{j+k} = \sum_{i=0}^{M-1} \sum_{j=iL}^{(i+1)L-1} \ol{\xi_j} \xi_{j+k}, \\
\alf_k & = \sum_{j=0}^{rL-1} \ol{\xi_j} \xi_{j+k} = \sum_{i=0}^{r-1} \sum_{j=iL}^{(i+1)L-1} \ol{\xi_j} \xi_{j+k}.
\end{align*}
CLAIM 1: If $iL \le j \le (i+1)L-1$ and $k=qL$ as above, then
$$
p(j+k)-p(j) = 2qj + Lq^2 + Aq
$$
which is of the form $2qj+c_1$ where $c_1$ is independent of $i$ and $j$.

Claim 1 is proved in the appendix.

We then have
\begin{align*}
\gam_k & = \sum_{i=0}^{M-1} \sum_{j=iL}^{(i+1)L-1} \zet_{2M}^{p(j+k)-p(j)}
= \zet_{2M}^{c_1} \sum_{i=0}^{M-1} \sum_{j=iL}^{(i+1)L-1} \zet_{2M}^{2qj} \\
& = \zet_{2M}^{c_1} \sum_{j=0}^{ML-1} \zet_{2M}^{2qj}
= \zet_{2M}^{c_1} \sum_{i=0}^{L-1} \sum_{j=iM}^{(i+1)M-1} \zet_{2M}^{2qj} \\
& = \zet_{2M}^{c_1} \sum_{i=0}^{L-1} \sum_{j=iM}^{(i+1)M-1} \zet_{M}^{qj}
= \zet_{2M}^{c_1} \sum_{i=0}^{L-1} 0 = 0
\end{align*}
where we have used Lemma 3.1 and the fact that $q$ is not a multiple of $M$.
We also have
$$
\alf_k = \sum_{i=0}^{r-1} \sum_{j=iL}^{(i+1)L-1} \zet_{2M}^{p(j+k)-p(j)}
= \zet_{2M}^{c_1} \sum_{i=0}^{r-1} \sum_{j=iL}^{(i+1)L-1} \zet_{2M}^{2qj}
= \zet_{2M}^{c_1} \sum_{j=0}^{rL-1} \zet_{2M}^{2qj}
$$
which implies
$$
\abs{\alf_k} = \abs{\sum_{j=0}^{rL-1}\zet_{2M}^{2qj}}
= \abs{\sum_{j=0}^{rL-1}\zet_{M}^{qj}} \le \csc\frac{\pi q}{M}
= \csc\frac{\pi(k/L)}{M} = \csc\frac{\pi k}{N}
$$
where we have used Corollary 3.2 and the fact that $0<q<M$.

Since $\gam_k=0$, we know $\abs{\alf_k}=\abs{\alf_{N-k}}$.
Therefore when bounding $\abs{\alf_k}$, it suffices to consider
$k \le N/2$.
If $1 \le k \le \sqrt{N/2}$, then $\abs{\alf_k} \le k \le \sqrt{N/2}$.
If $\sqrt{N/2} \le k \le N/2$, then
$$
\abs{\alf_k} \le \csc\frac{\pi k}{N} \le \frac{N}{2k} \le \frac{N}{2\sqrt{N/2}} = \sqrt{\frac{N}{2}}
$$
which completes the proof of Proposition 4.1.

{\bf Proposition 4.2.} Let $\mathcal{X}=(\xi_0,\ldots,\xi_{N-1})$
be the LM sequence of length $N=LM$ defined previously.
Suppose $1 \le k \le N-1$, and suppose $k$ and $N-k$ are both
nonmultiples of $L$. Then the autocorrelation $\gam_k$
satisfies $\gam_k=0$.
Furthermore, in the special case $L=M$,
the autocorrelation $\alf_k$
satisfies $\abs{\alf_k}\le M = \sqrt{N}$.

{\it Proof.}
If $k$ and $N-k$ are both nonmultiples of $L$, let $k = qL+k_1$
and let $N-k = rL+k_2$, where $1 \le k_1,k_2 \le L-1$ and $k_1+k_2=L$.
So $q+r=M-1$.
We break the sums $\gam_k$ and $\alf_k$ into sums of $k_1$ terms
and sums of $k_2$ terms:
\begin{align*}
\gam_k & = \sum_{i=0}^{M-1} \sum_{j=iL}^{iL+k_2-1} \ol{\xi_j} \xi_{j+k}
+ \sum_{i=0}^{M-1} \sum_{j=iL+k_2}^{(i+1)L-1} \ol{\xi_j} \xi_{j+k}, \\
\alf_k & = \sum_{i=0}^{r} \sum_{j=iL}^{iL+k_2-1} \ol{\xi_j} \xi_{j+k}
+ \sum_{i=0}^{r-1} \sum_{j=iL+k_2}^{(i+1)L-1} \ol{\xi_j} \xi_{j+k}.
\end{align*}
CLAIM 2: If $iL \le j \le iL+k_2-1$ and $k=qL+k_1$ where $k_1,k_2$ are as above, then
$$
p(j+k)-p(j) = 2ik_1 + 2qj + 2qk_1 + Lq^2 + Aq
$$
which is of the form $2ik_1 + 2qj + c_2$ where $c_2$ is independent of $i$ and $j$.

CLAIM 3: If $iL+k_2 \le j \le (i+1)L-1$ and $k=qL+k_1$ where $k_1,k_2$ are as above, then
$$
p(j+k)-p(j) = -2ik_2 + 2\tq j - 2\tq k_2 + L\tq^2 + A\tq
$$
where $\tq=q+1$. This is of the form $-2ik_2+2\tq j+c_3$ where
$c_3$ is independent of $i$ and $j$.

Claims 2 and 3 are proved in the appendix.

We then have
\begin{align*}
\gam_k & = \sum_{i=0}^{M-1} \sum_{j=iL}^{iL+k_2-1} \zet_{2M}^{p(j+k)-p(j)}
+ \sum_{i=0}^{M-1} \sum_{j=iL+k_2}^{(i+1)L-1} \zet_{2M}^{p(j+k)-p(j)} \\
& = \zet_{2M}^{c_2} \sum_{i=0}^{M-1} \zet_{2M}^{2ik_1} \sum_{j=iL}^{iL+k_2-1} \zet_{2M}^{2qj}
+ \zet_{2M}^{c_3} \sum_{i=0}^{M-1} \zet_{2M}^{-2ik_2} \sum_{j=iL+k_2}^{(i+1)L-1} \zet_{2M}^{2\tq j} \\
& = \zet_{2M}^{c_2} \sum_{i=0}^{M-1} \zet_{2M}^{2ik_1} \sum_{j'=0}^{k_2-1} \zet_{2M}^{2q(iL+j')}
+ \zet_{2M}^{c_3} \sum_{i=0}^{M-1} \zet_{2M}^{-2ik_2} \sum_{j'=k_2}^{L-1} \zet_{2M}^{2\tq(iL+j')} \\
& = \zet_{2M}^{c_2} \sum_{i=0}^{M-1} \zet_{2M}^{2i(k_1+qL)} \sum_{j'=0}^{k_2-1} \zet_{2M}^{2qj'}
+ \zet_{2M}^{c_3} \sum_{i=0}^{M-1} \zet_{2M}^{2i(\tq L-k_2)} \sum_{j'=k_2}^{L-1} \zet_{2M}^{2\tq j'} \\
& = \zet_{2M}^{c_2} \cdot \sum_{i=0}^{M-1} \zet_{2M}^{2ik} \cdot \sum_{j'=0}^{k_2-1} \zet_{2M}^{2qj'}
+ \zet_{2M}^{c_3} \cdot \sum_{i=0}^{M-1} \zet_{2M}^{2ik} \cdot \sum_{j'=k_2}^{L-1} \zet_{2M}^{2\tq j'} \\
& = \zet_{2M}^{c_2} \cdot \sum_{i=0}^{M-1} \zet_{M}^{ik} \cdot \sum_{j'=0}^{k_2-1} \zet_{M}^{qj'}
+ \zet_{2M}^{c_3} \cdot \sum_{i=0}^{M-1} \zet_{M}^{ik} \cdot \sum_{j'=k_2}^{L-1} \zet_{M}^{\tq j'} \\
& = \zet_{2M}^{c_2} \cdot 0 \cdot \sum_{j'=0}^{k_2-1} \zet_{M}^{qj'}
+ \zet_{2M}^{c_3} \cdot 0 \cdot \sum_{j'=k_2}^{L-1} \zet_{M}^{\tq j'} = 0
\end{align*}
where we have used Lemma 3.1. (The hypotheses of Proposition 4.2 say that $k$
is not a multiple of $L$ and hence not a multiple of $M$.)
We also have, by similar manipulations,
\begin{align*}
\alf_k & = \sum_{i=0}^{r} \sum_{j=iL}^{iL+k_2-1} \zet_{2M}^{p(j+k)-p(j)}
+ \sum_{i=0}^{r-1} \sum_{j=iL+k_2}^{(i+1)L-1} \zet_{2M}^{p(j+k)-p(j)} \\
& = \cdots \\
& = \zet_{2M}^{c_2} \cdot \sum_{i=0}^{r} \zet_{M}^{ik} \cdot \sum_{j'=0}^{k_2-1} \zet_{M}^{qj'}
+ \zet_{2M}^{c_3} \cdot \sum_{i=0}^{r-1} \zet_{M}^{ik} \cdot \sum_{j'=k_2}^{L-1} \zet_{M}^{\tq j'}
\end{align*}
which implies
\begin{equation} \label{alfk}
\abs{\alf_k} = \abs{\sum_{i=0}^{r}\zet_{M}^{ik}} \abs{\sum_{j'=0}^{k_2-1}\zet_{M}^{qj'}}
+ \abs{\sum_{i=0}^{r-1}\zet_{M}^{ik}} \abs{\sum_{j'=k_2}^{L-1} \zet_{M}^{\tq j'}}.
\end{equation}
Now suppose we are in the special case $L=M$.
Note that $0 \le q \le M-1$. If $q=0$, then $k < L$, so $\abs{\alf_k} < k < L = M$.
If $q=M-1$, then $r=0$, so $N-k < L$, so $\abs{\alf_k} < N-k < L = M$.
So assume $1 \le q \le M-2$, implying $2 \le \tq \le M-1$.
Then neither $q$ nor $\tq$ is a multiple of $M$.
Lemma 3.1 then gives us
\begin{align*}
\sum_{j'=0}^{L-1} \zet_M^{qj'} & = \sum_{j'=0}^{M-1} \zet_M^{qj'} = 0, \\
\sum_{j'=0}^{L-1} \zet_M^{\tq j'} & = \sum_{j'=0}^{M-1} \zet_M^{\tq j'} = 0,
\end{align*}
which then implies
\begin{align*}
\abs{\sum_{j'=0}^{k_2-1} \zet_M^{qj'}}
& = \abs{\sum_{j'=k_2}^{L-1} \zet_M^{qj'}}
\le \min\{k_1,k_2\}, \\
\abs{\sum_{j'=0}^{k_2-1} \zet_M^{\tq j'}}
& = \abs{\sum_{j'=k_2}^{L-1} \zet_M^{\tq j'}}
\le \min\{k_1,k_2\}.
\end{align*}
Next, observe that we have $\delta(k,M) = \delta(k,L) = \min\{k_1,k_2\}$.
If we let $k'=\min\{k_1,k_2\}$, then Corollaries 3.2 and 3.3 give us
\begin{align*}
\abs{\sum_{i=0}^r \zet_M^{ik}} & \le \abs{\csc\frac{\pi k}{M}} \le \frac{M}{2k'}, \\
\abs{\sum_{i=0}^{r-1} \zet_M^{ik}} & \le \abs{\csc\frac{\pi k}{M}} \le \frac{M}{2k'},
\end{align*}
Then inequality (\ref{alfk}) implies
$$
\abs{\alf_k} \le \frac{M}{2k'} \cdot k' + \frac{M}{2k'} \cdot k' = M,
$$
which completes the proof of Proposition 4.2.

In the case $L \ne M$, it is less clear how to bound $\abs{\sum\zet_M^{qj'}}$
and $\abs{\sum\zet_M^{\tq j'}}$. We can bound them by $L$, but that does not
make it obvious whether $\abs{\alf_k}$ can be bounded by a multiple of $\sqrt{N}$.
A more careful analysis may be needed.

In summary, Propositions 4.1 and 4.2 together imply that every LM sequence
satisfies $\gam_k=0$ for all $k\ne0$, i.e.\ every LM sequence is a perfect sequence.
In the special case $L=1$ (which includes the Chu sequences),
$k$ is always a multiple of $L$, so Proposition 4.1
always applies, and the LM sequence has PSL at most $\sqrt{N/2}$.
In the special case $L=M$ (which includes the Frank sequences),
the LM sequence has PSL at most $\sqrt{N}$.

\section{Appendix}

As before, $L$ and $M$ are positive integers satisfying $L|M$,
and $A \equiv LM$ mod 2.
For all nonnegative integers $j$,
we define
\begin{align*}
s_j & = \left\lfloor \frac{j}{L} \right\rfloor, \\
t_j & = j - L \left\lfloor \frac{j}{L} \right\rfloor, \\
p(j) & = 2s_jt_j + Ls_j^2 + As_j.
\end{align*}
Proof of Claim 1.
If $iL \le j \le (i+1)L-1$ and $k=qL$, then
$$
(i+q)L \le j+k \le (i+q+1)L-1,
$$
so we have
\begin{align*}
s_j & = \left\lfloor \frac{j}{L} \right\rfloor = i, \\
s_{j+k} & = \left\lfloor \frac{j+k}{L} \right\rfloor = i+q, \\
t_j & = j - Li, \\
t_{j+k} & = (j+k) - L(i+q) = j - Li = t_j, \\
s_{j+k} - s_j & = q, \\
s_{j+k}^2 - s_j^2 & = (i+q)^2 - i^2 = 2iq + q^2, \\
s_{j+k}t_{j+k} - s_jt_j & = (s_{j+k}-s_j)t_j = q(j-Li).
\end{align*}
So then
\begin{align*}
p(j+k)-p(j) & = 2(s_{j+k}t_{j+k}-s_jt_j) + L(s_{j+k}^2-s_j^2) + A(s_{j+k}-s_j) \\
& = 2q(j-Li) + L(2iq+q^2) + Aq \\
& = 2qj + Lq^2 + Aq.
\end{align*}

Proof of Claim 2.
If $iL \le j \le iL+k_2-1$ and $k=qL+k_1$, and $k_1,k_2$ are positive integers
satisfying $k_1+k_2=L$, then
$$
(i+q)L + k_1 \le j+k \le (i+q+1)L - 1,
$$
so we have
\begin{align*}
s_j & = \left\lfloor \frac{j}{L} \right\rfloor = i, \\
s_{j+k} & = \left\lfloor \frac{j+k}{L} \right\rfloor = i+q, \\
t_j & = j - Li, \\
t_{j+k} & = (j+k) - L(i+q) = j - Li + k_1, \\
s_{j+k} - s_j & = q, \\
s_{j+k}^2 - s_j^2 & = (i+q)^2 - i^2 = 2iq + q^2, \\
s_{j+k}t_{j+k} - s_jt_j & = (i+q)(j-Li+k_1)-i(j-Li) \\
& = ik_1 + qj - qLi + qk_1.
\end{align*}
So then
\begin{align*}
p(j+k)-p(j) & = 2(s_{j+k}t_{j+k}-s_jt_j) + L(s_{j+k}^2-s_j^2) + A(s_{j+k}-s_j) \\
& = 2(ik_1+qj-qLi+qk_1) + L(2iq+q^2) + Aq \\
& = 2ik_1 + 2qj + 2qk_1 + Lq^2 + Aq.
\end{align*}

Proof of Claim 3.
If $iL+k_2 \le j \le (i+1)L-1$ and $k=qL+k_1$,
and $k_1,k_2$ are positive integers satisfying $k_1+k_2=L$, then
$$
(i+q+1)L \le j+k \le (i+q+1)L + k_1 - 1.
$$
Let $\tq = q+1$.
Then we have
\begin{align*}
s_j & = \left\lfloor \frac{j}{L} \right\rfloor = i, \\
s_{j+k} & = \left\lfloor \frac{j+k}{L} \right\rfloor = i+\tq, \\
t_j & = j - Li, \\
t_{j+k} & = (j+k) - L(i+q+1) = j - Li - k_2, \\
s_{j+k} - s_j & = \tq, \\
s_{j+k}^2 - s_j^2 & = (i+\tq)^2 - i^2 = 2i\tq + \tq^2, \\
s_{j+k}t_{j+k} - s_jt_j & = (i+\tq)(j-Li-k_2) - i(j-Li) \\
& = -ik_2 + \tq j - \tq Li - \tq k_2.
\end{align*}
So then
\begin{align*}
p(j+k)-p(j) & = 2(s_{j+k}t_{j+k}-s_jt_j) + L(s_{j+k}^2-s_j^2) + A(s_{j+k}-s_j) \\
& = 2(-ik_2 + \tq j - \tq Li - \tq k_2) + L(2i\tq + \tq^2) + A\tq \\
& = -2ik_2 + 2\tq j - 2\tq k_2 + L\tq^2 + A\tq.
\end{align*}

\bibliographystyle{amsplain}


text:
\section{Introduction}
\IEEEPARstart{W}{ith} merits of fast charging, high energy density, low cost, and long lifespan \cite{Ref1},\cite{Ref2}, Lithium-ion batteries (LiBs) provide a promising alternative to fossil fuels and have been applied to a wide range of applications, including electric vehicles (EV), personal mobility devices, portable instruments, smart-grid energy storage devices, etc. Accurate state of charge (SoC) estimation, which measures remaining capacity for LiBs [3], not only has a significant influence on the performance of its powered system, e.g., driving ranges of EVs, but also is important to assess the health status of LiBs and ensure a safe environment for LiBs powered systems.

Over the past decades, SoC estimation has attracted considerable attention both in industrial and academic community. In the early days, model-based methods \cite{Ref4},\cite{Ref5} have been widely used, which, however, require plentiful process knowledge that is time-consuming, expensive, and difficult to derive. Additionally, model-based methods are specific to details of the actual battery materials combinations and the fabrication process, resulting in poor generalization from one battery type to another.

With the rapid development in sensor and data acquisition technologies, a wealth of LiBs cycling data become available. Reflecting the inherent status using easily obtained data resources, data-driven SoC estimation methods have shown their superiorities in the era of big data and artificial intelligence. Taking current, voltage, and temperature as inputs, He et al. \cite{Ref6} proposed a method using artificial neural network to estimate SoC values. However, measurements were assumed to be time independent during model training, which might not hold true due to strong correlations between past and future samplings. To consider the time dependence of measurements, Li et al. \cite{Ref7} proposed an estimation method based on recurrent neural network (RNN), wherein hidden nodes were designed to memorize previous information by introducing feedback linkage. Although RNN outperforms traditional neural networks, it faces the vanishing descent problem due to an increasing burden of candidate states with the increase of time. Alternatively, long short-term memory (LSTM) network introduces a forget mechanism into the hidden nodes of RNN to overcome the vanishing descent problem by freeing memory from useless information. Chemali et al. \cite{Ref8} put forward a  LSTM-based method and achieved obvious improvements. However, a crucial implication of current methods [6]-[8] lies in that they are developed at a fixed ambient temperature.

Although researchers are fully motivated to develop more accurate SoC estimation models for LiBs, it is facing challenges due to various internal and external factors. Especially, LiBs are highly sensitive to ambient temperature [9], implying that the pre-trained model developed at a specific ambient temperature cannot work reliably at other ambient temperatures. However, it is common that the battery powered Internet of Thing devices or EVs have to operate over a wide range of ambient temperatures. For example, when an EV is driven in a long journey or rapidly (dis)charged, major ambient temperature changes may be experienced from day to night. The discharging efficiency will also drastically be reduced under extreme ambient temperatures, mainly low temperatures [10]. A SoC prediction system in the car that does not take this into account will show a poor prediction of SoC and driving range. Achievements dealing with varying temperatures with physical models have been reported [11], [12], however, they may suffer from the same challenges of model-based methods. To the best of the authors' knowledge, feasible solutions for data-driven SoC estimation have not been proposed yet.

Recently, transfer learning (TL) has been developed as a novel way to learn properties from source data and transfer them to the target data, which has achieved great success in image, audio, and text processing \cite{Ref13}. Inspired by this, it provides an opportunity for rapid modeling of SoC estimation from a comparatively small amount of data if estimation ability from known similar states (such as in this case other known temperatures) could be transferred to a new one. TL sequentially addresses three core questions: when to transfer, what to transfer, and how to transfer. Naturally, it raises three critical issues with respect to SoC estimation for LiBs taking temperature into account: 1) how to judge whether the previously developed model, i.e., the reference SoC estimation model, is still suitable or not; 2) once the estimation relationship has been changed, what kind of information in the reference model can be transferred; 3) how to achieve rapid SoC estimation modeling with less data at new temperatures.

In this paper, based on the in-depth analysis of LiBs data characteristics, the following recognition are gained: 1) temporal dynamics is similar within a small temperature region; 2) despite temporal dynamics may be dissimilar when temperature varies, consistent temporal dynamics, which remains the same at different temperatures, may exist due to the nature of the electrochemical process. Herein, the temporal dynamics refers to the correlations between past fluctuation and the future motion. Based on these recognitions, a TL-based SoC estimation method considering temporal dynamics is proposed. First, data processing is achieved through performing wavelet analysis on measurements to deal with the non-stationarity and multi-scale characteristics. Next, temporal correlations, which reflect the SoC estimation ability, are extracted using canonical variate analysis (CVA) as feature engineering. After that, canonical variates (CVs) are fed into LSTM network to develop the reference SoC estimation model, achieving higher prediction accuracy. The efficacy of estimation model in face of varying temperature is quantitatively evaluated by monitoring normal variation of CVs in a linear and explainable way. Furthermore, if the reference SoC estimation model becomes invalid, alarms will be issued promptly and a rapid model updating strategy is triggered by transferring consistent estimation ability. Estimation monitoring and rapid model updating minimize the required amount of new data and ensure a high adaptivity to temperature changes. The contributions of the proposed method are summarized below:
\begin{enumerate}[1)]
	\item A data-driven strategy is proposed for SoC estimation under varying temperatures;
	\item Estimation ability monitoring is designed to assess the efficacy of the reference model;
	\item  A rapid model updating strategy with less data is proposed to adapt to new temperatures via transferring consistent estimation ability.
\end{enumerate}

The rest of this article is structured as follows: Section II describes the discharge procedure and compares two machine learning (ML) methods used in this article. Next, details of the proposed methods are shown in Section III. In Section IV, efficacy of the proposed method is verified through a LiB benchmark. Finally, the work is concluded in the last section.
\section{Preliminary}
This section begins with the description of LiBs discharge procedure. After that, comparisons of two popular ML methods used in this article, i.e., CVA and LSTM, are given for better understanding.

\subsection{Discharging Procedure and Data Illustration}
Discharge procedure provides energy for LiBs powered devices. The available capacity is high at a full charge state, corresponding to a large value of SoC. Once discharging, SoC values will decrease until voltage signal reached a pre-set threshold. Fig. 1 depicts two widely adopted variables, i.e., current and voltage, which are obtained in an experimental discharge cycle at 25$^\circ$C under US06 dynamic testing of a Panasonic 18650PF cell \cite{Ref14}. Current signal varies between -20A to 10A, and voltage signal shows a descending trend. The true value of SoC is shown at the bottom of Fig. 1, which continuously decreases during the discharging.
\begin{figure}
	\centering
	\includegraphics[scale=0.55]{figure/FIG1_TII-20-3344.pdf}\\
	\caption{Measurements of a real discharge cycle with US06 testing at 25$^\circ$C of a Panasonic 18650PF cell \cite{Ref14}.}
	\vspace{-0.5cm}
	\label{MyFig1}
\end{figure}

\subsection{Comparisons between CVA and LSTM}
\begin{table*}[ht]
	\centering
	\tiny
	\renewcommand{\arraystretch}{1.2}
	\caption{Comprehensive comparisons between CVA and LSTM with respect to temporal dynamics.}
	\label{Table_1}
	\begin{center}
		\begin{threeparttable}
			\begin{tabular}{p{2.2cm} p{1.2 cm} p{1 cm} p{1 cm} p{1cm} p{2 cm} p{2.2 cm}}
				\hline
				\toprule
				\multicolumn{1}{c} {\diagbox {\textbf{Approach}}{\textbf{Property}}} & \multicolumn{1}{c} {\textbf{Model type}} & \multicolumn{1}{c} {\textbf{Data cost}} & \multicolumn{1}{c}{\textbf{Computation complexity}} & \multicolumn{1}{c}{\textbf{Solution type}} & {\textbf{What kind of temporal correlations can be captured?}} & {\textbf{Does the temporal correlation can be clearly described?}}\\
				\hline
				\multicolumn{1}{c}{CVA} & \multicolumn{1}{c}{White-box model} & \multicolumn{1}{c}{Low} & \multicolumn{1}{c}{SVD decomposition} & \multicolumn{1}{c}{Linear} & \multicolumn{1}{c}{Short-term} & \multicolumn{1}{c}{Yes} \\

				\multicolumn{1}{c}{LSTM} & \multicolumn{1}{c}{Black-box model} & \multicolumn{1}{c}{High} & \multicolumn{1}{c}{A multi-layer RNN} & \multicolumn{1}{c}{Nonlinear} & \multicolumn{1}{c}{Long-term}  & \multicolumn{1}{c}{No} \\
				\bottomrule
			\end{tabular}
		\end{threeparttable}
	\end{center}
\end{table*}

\begin{figure*}
	\centering
	\includegraphics[scale=0.35]{figure/FIG2_TII-20-3344.pdf}\\
	\caption{The structure and organization of the proposed method.}
	\label{MyFig2}
\end{figure*}

CVA is a dimensionality reduction method that maximizes the correlations between two datasets by finding linear combinations of high-dimensional process data \cite{Ref15}. A LSTM network is a sequence to sequence learning method, which learns how to relate an entire sequence with target state outputs in a state-space manner. To deal with the vanishing gradient problem of  traditional RNN, hidden nodes are designed as LSTM cells, where a forget gate is newly introduced and memory states are involved to remove redundant information. Due to limitation of pages, more information about these two algorithms can be found in \cite{Ref16}, \cite{Ref17}.

Although both methods are data-driven to analyze temporal dynamics, they present different working mechanisms and data requirements. Table \ref{Table_1} systematically presents the dissimilarities regarding both methods. First, CVA is a statistical ML method to find temporal correlations between variables in a linear and explainable way using latent variables. Differently, LSTM is a non-statistical ML model for long-period prediction. As a linear method, data requirements and computation complexity of CVA are much lower than that of LSTM. That is determined by the relatively simple structure of CVA. In contrast, LSTM owns a complex network, which may have multiple layers and thousands of nodes. As a result, the demands of hardware resource and amount of data to train LSTM are high.

In summary, CVA is good at revealing temporal correlations in a analytical expression with a small amount of data. And LSTM shows superiorities in finding long-period temporal prediction ability with amounts of data. To take advantages of the best of both methods for SoC estimation, an elaborate scheme is designed in this paper. The behind idea is as follows. To learn long-term temporal correlations of LiBs data, LSTM is used in its typical structure to serves as estimation model. For better estimation ability and overcoming the influence of varying temperatures, CVA is introduced as a feature extraction method and the purposes of these analytical features are threefold: 1) feeding features rather than original measurements benefits the accuracy improvement of LSTM based estimation ability; 2) extracted features allow for quantitative evaluation of the changes of estimation ability influenced by varying temperatures; 3) a prompt transfer learning works by classifying these features as consistent ones and inconsistent ones between the reference temperature and new temperatures.

\section{Methodology}
The basic structure and organization of the proposed method are shown in Fig. \ref{MyFig2} and each part is specified in the following subsections. First, a reference SoC estimation model is developed using a multi-layer LSTM network fed by latent variables extracted using CVA. Next, this reference model is evaluated by checking the variances of obtained latent variables. Finally, a TL-based rapid updating strategy is proposed if the reference model fails due to the influence of a change in temperature.

\subsection{The Reference SoC Estimation Model}
The reference SoC estimation model is developed at a reference temperature, where large amounts of data are available. LSTM is chosen for modeling by taking advantage of its ability in capturing long-period prediction. To further enhance LSTM, temporal features extracted by CVA are fed into LSTM instead of original measurements. The specific steps are:

\subsubsection{Data processing}
Considering voltage signal shows serious non-stationarity, it is improper to be directly applied in CVA and LSTM that are more suitable for stationary signals. As a multi-resolution decomposition method, wavelet analysis (WA) has been successfully applied in signal processing, which decomposes a non-stationary signal into a series of stationary components. The core idea behind WA is that components of signal at corresponding frequencies will be amplified and easy to recognize through a series of wavelet functions that are finite energy and smooth. Morlet et al. \cite{Ref18} proposed the discrete wavelet analysis (DWA) method, in which a signal was decomposed into an average value and a series details. Decomposition of current signal $\mathbf c_n$ and voltage signal $\mathbf v_n$ in discharge cycle $n$ using DWA is given below,
\begin{equation}
\begin{array}{l}
{\mathbf{c}_{n}=\mathbf{c}_{n,a}+\sum_{j=1}^{J_{c}} \mathbf{c}_{n,j}} \\ [1mm]
{\mathbf{v}_{n}=\mathbf{v}_{n,a}+\sum_{j=1}^{J_{v}} \mathbf{v}_{n,j}}
\end{array}
\end{equation}
where $\mathbf c_{n,a}$ and $\mathbf v_{n,a}$ are the average part reflecting the low-frequency information of current and voltage in the $n^{th}$ cycle, respectively; $\mathbf c_{n,j}$ and $\mathbf v_{n,j}$ are the $j^{th}$ detail reflecting high-frequency information of the $n^{th}$ current and voltage, respectively; $J_c$ and $J_v$ are the number of details of current and voltage, respectively.

The original two-dimensional data matrix $\mathbf X_n$ is extended into a new matrix with dimensions $K_n \times J_x$, in which $J_x=J_c+J_v+2$ and $K_n$ is the number of samplings of the $n^{th}$ cycle. For brevity, the new matrix is still denoted as $\mathbf X_n$.

\begin{figure}[tb]
\centering
\includegraphics[scale=0.55]{figure/FIG3_TII-20-3344.pdf}\\
\caption{Schematic strategy of the past data and the future data matrixes.}
\label{MyFig3}
\end{figure}

\subsubsection{CVA-based feature engineering} Although advanced deep learning networks are capable of developing estimation models with raw data,  feature engineering plays a vital role in further boosting the performance of data-driven models [19]-[21]. For instance, voltage falloff, a health indicator to average the discharging voltage with an equivalent interval of the whole discharging cycle, has successfully reduced the battery RUL estimation error [19], [20]. In this article, CVA-based feature engineering is preferred due to its ability to capture temporal correlation from time-series.

To relate the temporal property, past data vector $\mathbf x_{p,n}(k)=[\mathbf x_n(k-1), \mathbf x_n(k-2),..., \mathbf x_n(k-l)]$ and future data vector $\mathbf x_{f,n}(k)=[\mathbf x_n(k), \mathbf x_n (k+1),..., \mathbf x_n (k+h-1)]$ are formed at each sampling time, where $\mathbf x_n (k)$ is the $k^{th}$ row of $\mathbf X_n$; $p$ and $f$ indicate past information and future information, respectively. The parameters $l$ and $h$ relate to how many past data and future data will be used for modeling, which meets well with the purpose of the autocorrelation function. Thus, the specific values of these two parameters are calculated using the root summed squares of all variables from the normal data against certain confidence bound ($\pm 5\%$ here). Usually, these two parameters are given with the same value.

Then past data matrix $\mathbf X_{p,n}$ and future data matrix $\mathbf X_{f,n}$ are gained as ${\mathbf{X}_{p,n}=[\mathbf{x}_{p,n}(l), \mathbf{x}_{p,n}(l+1), ..., \mathbf{x}_{p,n}(K_n-h)]}$ and ${\mathbf{X}_{f,n}=[\mathbf{x}_{f,n}(l+1), \mathbf{x}_{f,n}(l+2), ..., \mathbf{x}_{f,n}(K_n-h+1)]}$. Fig. 3 illustrates the data arrangement procedure. Further, comprehensive past and future matrixes $\mathbf X_p$ and $\mathbf X_f$ are formed by putting $\mathbf X_{p,n}$ and $\mathbf X_{f,n}$ in each cycle together along variable-wise direction, respectively. With $\mathbf X_p$ and $\mathbf X_f$, an eigenvalue decomposition is performed below,
\begin{equation}
\mathbf \Sigma_{p p}^{-1/2} \mathbf \Sigma_{p f} \mathbf \Sigma_{f f}^{-1/2}=\mathbf{U}_{x} \boldsymbol{\Sigma}_{x} \mathbf{V}_{x}
\end{equation}
where $\mathbf \Sigma_{p p}$ and $\mathbf \Sigma_{f f}$ are the covariance matrixes of $\mathbf X_p$ and $\mathbf X_f$, respectively; $\mathbf \Sigma_{p f}$ is the cross-covariance matrix between $\mathbf X_p$ and $\mathbf X_f$; $\mathbf \Sigma_{x}$ is a diagonal matrix of nonnegative singular values; $\mathbf U_x$ and $\mathbf V_x$ are right-singular matrix and left-singular matrix, respectively.

In this way, CVs $\mathbf Z_x$ in $\mathbf X_p$ are calculated as linear combinations of $\mathbf X_p$ by retaining all features as follows,
\begin{equation}
\mathbf{Z}_{x}=\mathbf{J}_{x} \mathbf{X}_{p}=\mathbf{V}_{x}^{\mathrm{T}} \boldsymbol{\Sigma}_{p p}^{-1/2} \mathbf{X}_{p}
\end{equation}

Similarly, CVs $\mathbf Z_y$ in $\mathbf X_f$ are calculated as $\mathbf{Z}_{y}=\mathbf{J}_{y} \mathbf{X}_{f}=\mathbf{U}_{x}^{\mathrm{T}} \boldsymbol{\Sigma}_{f f}^{-1/2} \mathbf{X}_{f}$. It should be noted that both $\mathbf X_p$ and $\mathbf X_f$ should be normalized to zero mean and unit variance before being applied in Eq. (2).

\subsubsection{LSTM-based reference SoC estimation model}
A long-period prediction ability is achieved through a multi-layer LSTM network, which consists of one input layer, several LSTM hidden layers, and an output layer. Taking $\mathbf Z_x$ as inputs and the real values of SoC $\mathbf Y_x$ as ground truth, the reference SoC estimation model $\Gamma (\cdot)$ achieves sequential to sequential learning through minimizing the loss function $L$ between real values and estimations below,
\begin{equation}
\min \limits_{\Gamma(\cdot)} L=\sum_{k=1}^{K_{z}}\left(y(k)-\Gamma\left(\mathbf{z}_{x}(k)\right)\right)^{2}
\end{equation}
where $\mathbf z_x(k)$ is the $k^{th}$ row of $\mathbf Z_x$, $y(k)$ is the real values in $\mathbf Y_x$, and $K_z$ is the total number of samples in $\mathbf Z_x$.

In summary, the inputs of LSTM in reference SoC model are temporal feature $\mathbf Z_x$ and corresponding real SoC values as ground truth, and the well-trained model outputs predictions of SoC for new cycling data. The benefits of developing the hybrid architecture of CVA-based feature engineering and LSTM are twofold. On the one hand, the temporal correlations that contribute to estimation ability are specified. On the other hand, the advantages of the long-term learning ability of LSTM can be retained. It is worth noticing that interpretable insights on the correlation of individual features into the target output could be achieved by adopting interpretable LSTM \cite{Ref32}. During training, Adam optimization \cite{Ref22} is employed to update the weights and biases.

\subsection{Prediction Performance Evaluation}
Having a close relationship with estimation ability, temporal correlation is used as an indicator to evaluate SoC estimation performance by checking its variability. Therefore, $\mathbf Z_x$ is used to develop monitoring model based on multivariate statistical process analysis.
\subsubsection{Monitoring modelling of SoC estimation}
Diagonal elements of $\mathbf \Sigma_{x} = \operatorname{diag}(\alpha_1,\alpha_2,...,\alpha_{J_x l})$ in Eq. (2) reflect the correlation of pair-wise CVs in $\mathbf Z_x$ and $\mathbf Z_y$. For instance, $\alpha_1=\mathbf Z_{x,1}^T \mathbf Z_{y,1}$, where $\mathbf Z_{x,1}$ and $\mathbf Z_{y,1}$ are the first column of $\mathbf Z_x$ and $\mathbf Z_y$, respectively. According to the amplitude of correlation, $\mathbf Z_x$ are classified into two categories: system canonical variates (SCVs) and the residual variates (RVs). SCVs reflect main CVs with large correlations. And CVs with small correlations are grouped into RVs. Given there are $R$ SCVs, $\textbf Z_x$ can be decomposed as follows,
\begin{equation}
\begin{array}{l}
{\mathbf{Z}_{s}=\mathbf{J}_{s} \mathbf{X}_{p}=\mathbf{V}_{s}^{\mathrm{T}} \mathbf{\Sigma}_{p p}^{-1/2} \mathbf{X}_{p}} \\ [1mm]
{\mathbf{Z}_{r}=\mathbf{J}_{r} \mathbf{X}_{p}=(\mathbf{I}-\mathbf{V}_{s} \mathbf{V}_{s}^{\mathrm{T}}) \mathbf \Sigma_{p p}^{-1/2} \mathbf{X}_{p}}
\end{array}
\end{equation}
where $\mathbf Z_s$ contains SCVs, $\mathbf Z_r$ consists of RVs, $\mathbf V_s$ is the first $R$ columns of $\mathbf V_x$ in Eq. (2), and $\mathbf I$ is an identity matrix with dimension $J_xl$. Besides, $s$ and $r$ refer to the system subspace spanned by SCVs and the residual subspace spanned by RVs, respectively. The parameter $R$ stands for the number of dominant singular values, which is determined by finding a ``knee" point \cite{Ref23} from the singular matrix $\mathbf \Sigma_{x}$ calculated from Eq. (2).

To quantitatively measure the variance of SCVs and RVs, two monitoring statistics, i.e., indicators, are constructed using Hotelling's $T^2$ and squared prediction error ($SPE$) \cite{Ref24}. At sampling time $k$, they are calculated as follows,
\begin{equation}
\begin{array}{l}
{T^{2}(k)=\sum_{i=1}^{R} z_{s,i}(k)^{\mathrm{T}} z_{s,i}(k)} \\ [1mm]
{SPE(k)=\sum_{i=1}^{J_{x} l} z_{r,i}(k)^{\mathrm{T}} z_{r,i}(k)}
\end{array}
\end{equation}
where $z_{s,i}(k)$ is the $k^{th}$ row and the $i^{th}$ column in $\mathbf Z_s$, and $z_{r,i}(k)$ is the $k^{th}$ row and the $i^{th}$ column in $\mathbf Z_r$.

After checking the normality of $\mathbf Z_s$ and $\mathbf Z_r$ \cite{Ref25}, corresponding statistical control limits can be derived. If $\mathbf Z_s$ and $\mathbf Z_r$ follow Gaussian distribution, $\chi^2$ and $F$ distribution can be used to derive the statistical control limits $CL_s$ and $CL_r$ \cite{Ref26} with a given significance level, which is specified as 0.95 here. That is, the control limits will cover 95$\%$ normal samples. Otherwise, kernel density estimation \cite{Ref27} is suggested to estimate the probability distribution of $T^2$ and $SPE$ when $\mathbf Z_s$ and $\mathbf Z_r$ are non-Gaussian distribution, and then corresponding control limits can be calculated.

\subsubsection{Evaluation of the reference SoC estimation model}
For online monitoring, a data vector $\mathbf x_{new}(J_xl \times 1)$ is collected, which meets with the required window length of DWA. First, perform DWA on $\mathbf x_{new}$ according to Eq. (1) to obtain multi-solution information and then obtain the past data vector. Next, the normalized past data vector, which is still denoted as $\mathbf x_{new}$, is projected on the system subspace $\mathbf J_s$ and the residual subspace ${\mathbf J}_{r}$ given in Eq. (5) to get SCVs $\mathbf z_{new,s}=\mathbf J_s \mathbf x_{new,s}$ and the residuals $\mathbf z_{new,r}=\mathbf J_r \mathbf x_{new,r}$, respectively, in which $s$ and $r$ have the same meaning as that in Eq. (5). After that, the online monitoring statistics of the system subspace and the residual subspace are computed below,
\begin{equation}
\begin{aligned}
t^{2} &=\sum_{i=1}^{R} z_{new,s,i}^{\mathrm{T}} z_{new,s,i} \\
spe &=\sum_{i=1}^{J_xl} z_{new,r,i}^{\mathrm{T}} z_{new,r,i}
\end{aligned}
\end{equation}
where $z_{new,s,i}$ and $z_{new,r,i}$ are the $i^{th}$ element of $\mathbf z_{new,s}$ and $\mathbf z_{new,r}$, respectively.

\begin{figure}[htb]
	\centering
	\includegraphics[scale=0.75]{figure/FIG4_TII-20-3344.pdf}\\
	\caption{Flow chart of the transfer estimation learning at new temperatures.}
	\label{MyFig3}
\end{figure}

By sliding a moving window with fixed window length of DWA, $\mathbf x_{new}$ is updated with time evolution and monitoring of reference SoC model will proceed. By continuous comparing the online monitoring statistics with their control limits, performance of the reference SoC estimation model is judged according to the following rules:

Case I: If both online monitoring statistics are below their control limits, it indicates that SCVs and RVs both stay normal and estimation ability is not influenced by changes in temperature. Thus, the reference SoC estimation model trained in Eq. (4) can be directly adopted for estimation,
\begin{equation}
\hat{y}_{new}=\Gamma (\mathbf{z}_{new})
\end{equation}
where $\mathbf{z}_{new}$ is the collection of $\mathbf{z}_{new,s}$ and $\mathbf{z}_{new,r}$.

Case II: If at least one of online monitoring statistics is above its control limit in continuous three samplings, it means that partial CVs have been changed due to temperature variations. In this case, the current estimation relationship is improper and the reference SoC estimation model needs to be updated through TL, which will be shown in Subsection III.C.

\subsection{Transfer Estimation Learning at New Temperatures}
Once monitoring statistics are abnormal, it indicates that the reference SoC estimation model is no longer feasible to be used. Indeed, though temporal correlations may be dissimilar when temperature varies, consistent temporal correlations, which keep the same at different temperatures, may exist due to the fundamental nature of the electrochemical process in LiBs. Inspired by this, the basic idea of the proposed transfer estimation model is to inherit the common prediction ability from the source task meanwhile exploiting the ability of specific prediction in the target task.

The modeling procedures at the reference and the new temperatures are regarded as the source and the target task. Naturally, this raises two key issues when transfer process is triggered: 1) what to transfer, and 2) how to transfer. Fig. \ref{MyFig3} sketches the essence of the proposed transfer strategy, which includes two corresponding parts. The selection of consistent features from the reference SoC estimation model solves the first problem. The proposed transfer estimation model addresses the second one by combining consistent features from the reference temperature and the specific features at new temperatures.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.4]{figure/FIG5_TII-20-3344.pdf}\\
\begin{center}
\caption{Detailed illustration of the proposed method.}
\end{center}
\label{MyFig5}
\vspace{-0.5cm}
\end{figure*}

Let $\mathbf X_t$ a complete discharge cycle with current and voltage collected at new temperatures. $\mathbf X_t$ is decomposed using DWA and the past data matrix is arranged. The past data matrix is normalized with its own information and still denoted as $\mathbf X_t$. Correspondingly, the real values of SoC at new temperatures are denoted as $\mathbf Y_t$. To investigate the temporal dynamics of $\mathbf X_t$, CVA is performed on $\mathbf X_t$ in an analogous way to Eq. (2). CVs $\mathbf Z_t$ in target temperature are calculated by retaining all features below,
\begin{equation}
\mathbf{Z}_{t}=\mathbf{J}_{t} \mathbf{X}_{t}=\mathbf{V}_{t}^{\mathrm{T}} \mathbf {\Sigma}_{t t}^{-1/2} \mathbf{X}_{t}
\end{equation}
where $\mathbf J_t$ is projection matrix, $\mathbf V_t$ is left-singular matrix, and $\mathbf {\Sigma}_{t t}$ is covariance matrix of $\mathbf X_t$.

\begin{algorithm}[!ht]
	\scriptsize
	\caption{Selection of consistent features}
	\label{alg1}
	\begin{algorithmic}[1]
		\STATE Input: $\textbf Z_x$ and $\textbf Z_t$
		\STATE Initialize $q$ = 1
		\FOR{$q = 1, 2, ..., J_xl$}
		\STATE Retain the first $q$ CVs in $\textbf Z_x$ as $\textbf Z_{x,q}$
		\STATE Calculate $\textbf J_q$ and control limit $CL_q$ according to Eqs. (5) and (6) from $\textbf Z_{x,q}$
		\STATE Denote the first $q$ features in $\textbf Z_t$ as $\textbf Z_{t,q}$
		\STATE Calculate statistics $t_{t,q}^2$ from $\textbf Z_{t,q}$ at each sampling time according to Eq. (7)
		\IF    {$t_{t,q}^2 \leqslant CL_q$ holds for all sampling times}
		\STATE $q=q+1$
		\ELSE  [$t_{t,q}^2> CL_q$ holds at continuous three sampling times]
		\STATE $q=q-1$
		\ENDIF
		\ENDFOR
	\end{algorithmic}
\end{algorithm}

\subsubsection{Selection of consistent features}
Assuming the first $q$ features in $\mathbf Z_t$ are consistent with that of $\mathbf Z_x$, the monitoring statistics will be under the corresponding control limits calculated from the reference temperature. As a result, the proper value of $q$ can be determined by checking the variations of consistent features in $\mathbf Z_t$ with respect to the control limits derived from $\mathbf Z_x$. To select consistent features, a sequential selection strategy is given in {\textbf {Algorithm 1}}.

With the above procedures, the number of consistent features $q$ between the reference temperatures and target temperatures is determined. After that, consistent features $\mathbf Z_{x,q}$ in reference temperature are derived by separating the first $q$ rows in $\mathbf Z_x$. The consistent features in new temperatures is obtained as $\mathbf Z_{t,q}$, which is the collection of the first $q$ feature in $\textbf Z_t$. Correspondingly, the specific features $\mathbf Z_{t, \widetilde q}$ in new temperatures is computed by removing $\mathbf Z_{t,q}$ from $\textbf Z_t$.

\subsubsection{Development of transfer estimation model}
The proposed transfer estimation model consists of two independent LSTM networks. The first one is a multi-layer LSTM network that aims to learn shared prediction ability by training data pair $\{\mathbf Z_{x,q}, \mathbf Y_{sour}\}$, in which $\mathbf Y_{sour}$ is the real values of SoC in source task, i.e., the reference temperature. The other one is a single-layer LSTM network, which will be trained using a small amount of data pair $\{\mathbf Z_{t,\widetilde q}, \mathbf Y_t \}$, in which $\mathbf Y_t$ is the real values of SoC in target task (the new temperature). These two networks are independently designed for rapid TL purpose. The shared prediction ability $F_q(\cdot)$ is derived by minimizing estimation errors below,
\begin{equation}
\begin{array}{l}
{\min \limits_{F_{q}} \sum {({\hat{\mathbf Y}_{sour}}(k)-\mathbf{Y}_{sour}(k))^{2}}} \\
{{s.t.} \hat{\mathbf{Y}}_{sour}(k)=F_{q}(\mathbf{Z}_{x,q}(k))}
\end{array}
\end{equation}
where $\hat {\mathbf Y}_{sour}(k)$ is the estimation of $\mathbf Y_{sour}(k)$.

With $F_q(\cdot)$, the target task is represented as follows,
\begin{equation}
\begin{array}{l}
{\min \limits_{F_{t}} \sum (\hat{\mathbf{Y}}_{t}(k)-\mathbf{Y}_{t}(k))^{2}} \\
{{s.t.} \hat{\mathbf{Y}}_{t}(k)=\alpha_{1,k} F_{q}(\mathbf{Z}_{x,q}(k))+\alpha_{2,k} F_{t}(\mathbf{Z}_{t, \widetilde q}(k))}
\end{array}
\end{equation}
where $\hat {\mathbf Y}_t(k)$ is the estimation of $\mathbf Y_t(k)$; $\mathbf Z_{x,q}(k)$ and $\textbf Z_{t,\widetilde q}(k)$ are the $k^{th}$ row of $\mathbf Z_{x,q}$ and $\mathbf Z_{t,\widetilde q}$, respectively; $F_t(\cdot)$ is the specific prediction ability with respect to target temperature; $\alpha_{1,k}$ and $\alpha_{2,k}$ are real-time adjusting factors that both are initialized to be 0.5.

To balance the shared prediction ability and the specific prediction ability in Eq. (11), a strategy proposed in \cite{Ref28} is employed to update adjusting factors below,
\begin{equation}
\begin{array}{l}
{\alpha_{1,k}=\frac{{\alpha_{1,k} \Psi_{k}}(F_{q}(\mathbf{Z}_{x,q}(k)))}{\alpha_{1,k} \Psi_{k}(F_{q}(\mathbf{Z}_{x,q}(k)))+\alpha_{2, k} \Psi_{k}(F_{t}(\mathbf{Z}_{t,\widetilde q}(k)))}} \\ [2mm]
{\alpha_{2,k}=\frac{\alpha_{2,k} \Psi_{k}(F_{t}(\mathbf{Z}_{t,\widetilde q}(k)))}{\alpha_{1,k} \Psi_{k}(F_{q}(\mathbf{Z}_{x,q}(k)))+\alpha_{2,k} \Psi_{k}(F_{t}(\mathbf{Z}_{t,\widetilde q}(k)))}}
\end{array}
\end{equation}
where $\Psi_{k}(f(\cdot))=\exp (-\eta(f(\cdot)-y(k))^{2})$, $y(k)$ is the corresponding real value, and $\eta \in [0,1]$ is a coefficient.

From Eqs. (11) and (12), input data of the transfer estimation model include two parts. One is the consistent features $\mathbf Z_{x,q}$ with the corresponding real values of SoC $\mathbf Y_{sour}$ in the reference temperature as the label, and the other one is the specific features $\mathbf Z_{t,\widetilde q}$ with the corresponding real values $\mathbf Y_t$ in the new temperature as labels. Once the SoC transfer estimation model is well-trained, online estimations can be achieved for SoC values at the target temperature.
Here, the new sample at the $k^{th}$ sampling time is denoted as $\mathbf x_{t,new}(k)$, and then past data vector is arranged and normalized. After that, $\mathbf x_{t,new}(k)$ is projected on the subspace $\mathbf J_q$ expanded by consistent features and the remaining subspace. The consistent features and the specific features are obtained as $\mathbf z_{t,new,q}(k)=\mathbf J_{q} \mathbf x_{t,new}(k)$ and $\mathbf z_{t,new,\widetilde q}(k)=(\mathbf I - \mathbf J_{q}) \mathbf x_{t,new}(k)$. With $F_q(\cdot)$ and $F_t(\cdot)$, online estimation $\hat{y}_{t}(k)$ in new temperature is achieved below,
\begin{equation}
\hat{y}_{t}(k)=\alpha_{1}F_{q}(\mathbf z_{t,new,q}(k))+\alpha_{2}F_{t}(\mathbf z_{t,new,\widetilde q}(k))
\end{equation}
where $\alpha_1$ and $\alpha_2$ are the final values of adjusting factors.

For better understanding, Fig. 5 summarizes the main steps of the proposed method from off-line modelling and online estimation.

\begin{figure}[!ht]
\centering
\subfigure[]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=1.8in]{figure/FIG7a_TII-20-3344.pdf}
\end{minipage}%
}%
\subfigure[]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=1.8in]{figure/FIG7b_TII-20-3344.pdf}
\end{minipage}
}
\centering
\caption{Selection of propoer values for tuanble parameters at the reference temperature 10$^\circ$C with respect to (a) time lags $l$ and $h$ for reference estimation model and (b) the number of retained CVs $R$ for prediction performance evaluation model.}
\vspace{-0.5cm}
\end{figure}

\section{Results and discussions}
In this section, the proposed method is verified through a widely used LiB benchmark, which is provided by Dr. Kollmeyer \cite{Ref14}. It provides a series of drive cycles simulating an electric Ford F150 truck powered by a 35-kWh stack of 2.9 Ah Panasonic 18650PF cells (containing NCA oxide cathodes) at five discrete ambient temperatures (25$^\circ$C, 10$^\circ$C, 0$^\circ$C, -10$^\circ$C, and -20$^\circ$C). Ten drive cycles are available at each temperature, and each cycling data includes four variables, including current, voltage, battery surface temperature, and amp-hours. The first two variables are employed in the proposed method as input variables. The last variable amp-hours is measured to infer the real value of SoC. And the third variable is not used in our model to leave temperature as an external condition. All measurements are recorded at a frequency of 1Hz.

For fair comparisons, all simulations and data analysis are conducted on a desktop equipped with 16 processors of Intel Xeon E5-2620 v4 (20 MB cache, up to 2.10 GHz) and a multi-graphics processor unit of NVIDIA GeForce GTX 1080Ti at 11 GB. Programming software is Matlab with the version of 2019a. To evaluate the prediction accuracy, indices root mean square error ($RMSE$) and mean absolute error ($MAE$) are comprehensively used, which are given as follows,
\begin{equation}
\begin{array}{l}
{RMSE=\sqrt{\frac{1}{K} \sum_{k=1}^{K}(y_{t}(k)-\hat{y}_{t}(k))^{2}}} \\ [1mm]
{MAE=\frac{1}{K} \sum_{k=1}^{K}|y_{t}(k)-\hat{y}_{t}(k)|}
\end{array}
\end{equation}
where $\hat y_t(k)$ and $y_t(k)$ have the same meaning as that in Eq. (13), and $K$ is the number of sampling times in a testing cycle.

For all metrics, a smaller value indicates a better prediction accuracy [30], [31]. It should be noted that these indices are calculated with a complete cycling data.

\subsection{SoC Estimation Comparisons at Fixed Temperatures}
This part focuses on illustrating estimation performance of the proposed model at fixed temperatures. Five reference SoC models are developed corresponding to each temperature, where ten discharging cycles are available. Here, the first nine drive cycles are used for training, including one cycling data for validation. And the last cycling data is employed for testing. Considering the experiments are mainly conducted on 10$^\circ$C in \cite{Ref8}, comparisons are first completed at this temperature. Performing DWA on the current signal and the voltage signal of training data, each signal is decomposed into six components through retaining five details. Further, by conducting CVA on the decomposed data, a series of CVs are obtained according to Eqs. (2) and (3). In CVA, two crucial parameters $l$ and $h$ need to be determined. According to rules given in Section III.A, values of these two parameters are calculated using the root summed squares of all variables from the normal data against a certain confidence bound ($\pm 5 \%$ here). Therefore, both $l$ and $h$ are set to be 36 as shown in Fig. 6(a), where it is observed that the value steps into the defined confidence region \cite{Ref31}.

\begin{table*}[ht]
	\tiny
	\renewcommand{\arraystretch}{1.2}
	\caption{Comparisons between the proposed method and LSTM at five discrete temperatures concerning accuracy and network configuration.}
	\label{Table_2}
	\begin{center}
		\begin{threeparttable}
			\begin{tabular}{c c c c| c c c| c}
				\toprule
				\multicolumn{1}{c}{ \multirow{2}*{\textbf{Temp. ($^\circ$C)}}} & \multicolumn{1}{c}{\multirow{2}*{\textbf{Method}}} & \multicolumn{2}{c|}{\textbf{Indices}} & \multicolumn{3}{c|}{\textbf{Network configurations}} & \multicolumn{1}{c}{ \multirow{2}*{\textbf{Values of tunable parameters}}}\\
				\cline{3-7}
				\multicolumn{1}{c}{} & \multicolumn{1}{c}{} & $RMSE$ & $MAE$ & Iteration epochs & Training time & Network\\
				\hline
				\multirow{3}*{25} & The proposed method & 2.74 & 2.27 & 2000 & 105min34sec & L(50, 100)N(100) & $l=34,h=34,R=160$ \\
				& LSTM-based method & 5.46 & 4.29 &	15000 & 127min38sec & L(400)N(100) & - \\
				& Accuracy improvement & 49.82$\%$ & 47.09$\%$ & - & - & - & - \\
				\cline{1-8}
				\multirow{3}*{10} & The proposed method & 0.58& 0.46 & 2000 & 55min30sec &  L(50, 100)N(100) & $l=36,h=36,R=163$ \\
				& LSTM-based method & 1.41 & 1.12 & 15000 & 155min10sec & L(500)N(100) & -  \\
				& Accuracy improvement & 58.87$\%$ & 58.93$\%$ & - & - & - & - \\
				\cline{1-8}
				\multirow{3}*{0} & The proposed method & 1.91 &	1.70 & 2000 & 63min0sec &  L(50, 100)N(100)   & $l=40,h=40,R=165$  \\
				& LSTM-based method & 2.15 & 1.84 & 15000 & 135min3sec & L(500)N(100) & - \\
				& Accuracy improvement & 11.16$\%$ & 7.61$\%$ & - &	- & - & - \\
				\cline{1-8}
				\multirow{3}*{-10} & The proposed method &	2.89 & 2.17 & 3500 & 88min15sec &  L(150, 250 )N(100) & $l=35,h=35,R=155$ \\
				& LSTM-based method & 4.00 & 3.15 & 15000 & 116min45sec & L(400)N(100) & -   \\
				& Accuracy improvement & 27.75$\%$ & 31.11$\%$ & - & - & - & - \\
				\cline{1-8}
				\multirow{3}*{-20} & The proposed method & 4.94 & 3.95 & 2000 & 35min9sec &  L(50, 100)N(100) & $l=36,h=36,R=150$ \\
				& LSTM-based method & 6.53 & 5.14 & 15000	& 116min20sec & L(400)N(100) & - \\
				& Accuracy improvement & 24.35$\%$ & 23.15$\%$ & - & - & - & - \\
				\bottomrule
			\end{tabular}
			\begin{tablenotes}
				\footnotesize
				\item[-] denotes null.
			\end{tablenotes}
		\end{threeparttable}
	\end{center}
	\vspace{-0.8cm}
\end{table*}

\begin{figure}
	\centering
	\includegraphics[scale=0.55]{figure/FIG6_TII-20-3344.pdf}
	\caption{Prediction result comparisons between the proposed method and the LSTM-based method at 10$^\circ$C using the testing data.}
	\label{MyFig7}
	\vspace{-0.7cm}
\end{figure}

Inputting all CVs into the LSTM network with configurations given in Table II, Figs. 7(a) and 7(b) plot the predictions of SoC at 10$^\circ$C based on the proposed method and traditional LSTM-based method using testing discharge cycle, respectively. Fig. 7(c) further compares the absolute value of estimation error for both methods to give a better demonstration. Moreover, Table II summarizes the results at the other four temperatures with defined indices and specific values of tunable parameters. Values of each index at all temperatures have been greatly reduced using the proposed method. For 10$^\circ$C, index $RMSE$ is reduced by 58.87$\%$ from 1.41 (LSTM-based method) to 0.58 (the proposed method). The value is calculated through the ratio between reduced error and the original error, i.e., (1.41-0.58)/1.41. Accuracy improvements with respect to $RMSE$ are gained as 49.82$\%$, 11.16$\%$, 27.75$\%$, and 24.35$\%$ at 25$^\circ$C, 0$^\circ$C, -10$^\circ$C, and -20$^\circ$C, respectively. Therefore, prediction accuracy has been greatly improved and it proves that utilization of temporal dynamics enhances the precision of the results. Similar conclusions can be drawn with respect to index $MAE$.

\begin{table*}[ht]
	\tiny
	\renewcommand{\arraystretch}{1.2}
    	\centering
	\caption{Comparison of estimation errors using the proposed method with and without transfer learning.}
	\label{Table_3}
	\begin{center}
		\begin{threeparttable}
			\begin{tabular}{c c c c c| c c c| c}
				\toprule
				\multicolumn{1}{c}{\multirow{2}*{\textbf{Transferring case}}} & \multicolumn{1}{c}{\multirow{2}*{\textbf{Method}}} & \multicolumn{3}{c|}{\textbf{Indices}} & \multicolumn{3}{c|}{\multirow{1}*{\textbf{Transfer network configuration}}} & \multicolumn{1}{c}{\multirow{2}*{\textbf{Values for parameters}}} \\
				\cline{3-8}
				{} & {} & {$RMSE$} & {$MAE$} & {$H(\%)$} & {Iteration epochs} & {Training time} &{Network}\\
				\hline
				\multirow{3}*{From 10$^\circ$C to 25$^\circ$C} & {Model with TL\tnote{1}} & {6.93} & {5.73} & {86.44} & {100} & {1min45sec} & L(50)N(100) & {$q=200, \alpha_1=0.54, \alpha_2=0.46$}\\
				& {Model without TL\tnote{2}} & {33.77} & {29.37} & {-}  & {-} & {-}  & {-}  & {-} \\
				& {Accuracy improvement} & {79.48$\%$} & {80.49$\%$} & {-} & {-} & {-} & {-}  & {-} \\
				\cline{1-9}
				\multirow{3}*{From 10$^\circ$C to 0$^\circ$C} & {Model with TL\tnote{1}} & {6.00} & {4.50} & {56.84} & {100} & {1min10sec} & {L(50)N(100)} & {$q=108, \alpha_1=0.25, \alpha_2=0.75$}\\
				& {Model without TL\tnote{2}} & {37.33} & {32.05} & {-} & {-} & {-} & {-} & {-} \\
				& {Accuracy improvement} & {82.00$\%$} & {85.96$\%$} & {-} & {-} & {-} & {-}  & {-} \\
				\cline{1-9}
				\multirow{3}*{From 10$^\circ$C to -10$^\circ$C} & {Model with TL\tnote{1}} & {15.80} & 12.41 & {24.79} & {100}  & {1min23sec} & {L(50)N(100)} & {$q=45, \alpha_1=0.11, \alpha_2=0.89$}\\
				& {Model without TL\tnote{2}} & {36.08} & {31.88} & {-} & {-} & {-} & {-}  & {-} \\
				& {Accuracy improvement} & {56.21$\%$} & {61.73$\%$} & {-} & {-} & {-} & {-}  & {-} \\
				\cline{1-9}
				\multirow{3}*{From 10$^\circ$C to -20$^\circ$C} & {Model with TL\tnote{1}} & {21.59} & {16.96} & {8.28} & {100} & {1min9sec} & {L(50)N(100)} & {$q=15, \alpha_1=0.035, \alpha_2=0.965$}\\
				& {Model without TL\tnote{2}} & {36.03} & {41.53} & {-} & {-} & {-} & {-}  & {-} \\
				& {Accuracy improvement} & {40.08$\%$} & {59.16$\%$} & {-}& {-} & {-} & {-}  & {-} \\
				\bottomrule
			\end{tabular}
			\begin{tablenotes}
				\footnotesize
				\item[1] denotes the developed reference SoC model at 10$^\circ$C using transfer learning at new temperatures;
				\item[2] denotes the developed reference SoC model at 10$^\circ$C directly applied at new temperatures;
				\item[-] denotes null.
			\end{tablenotes}
		\end{threeparttable}
	\end{center}
\end{table*}

\begin{figure}[!ht]
	\centering
	\includegraphics[scale=0.55]{figure/FIG8_TII-20-3344.pdf}
	\caption{Monitoring results using US06 testing cycle at (a) 25$^\circ$C (b) 0$^\circ$C (c) -10$^\circ$C, and (d) -20$^\circ$C. (The red dashed line: control limits; the blue dotted line: monitoring statistics; the red dashed ellipse indicates abnormal samples)}
	\label{MyFig8}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.55]{figure/FIG9_TII-20-3344.pdf}\\
	\caption{Estimation results using the proposed method with TL at (a) 25$^\circ$C and (b) 0$^\circ$C when the reference temperature is 10$^\circ$C.}
	\label{MyFig9}
\end{figure}

The network configurations are available in Table II. For each reference temperature, the initial learning rate in the Adam optimization algorithm is given as 0.01, and batch size is selected as 32. To explain the way for configuring networks, L(50,100)N(100) at the reference temperature 10$^\circ$C is taken as an example. It means that the network that has three hidden layers with 50 cells in the first LSTM layer, 100 cells in the second LSTM layer, and the third fully-connected layer with 100 nodes. Finally, there is an output layer with one node to output the SoC estimation. We use dropout and L2 regularization to control model overfitting. The employing of dropout in the fully-connected layer will prevent overfitting by ignoring randomly selected neurons during training, hence reducing the sensitivity to individual neurons' specific weights. A 20$\%$ dropout rate is recommended as a good compromise between retaining model accuracy and preventing overfitting. Besides, we stop in the training process earlier when there is no improvement in the validation data set. Besides the provided guidance for tuning the network, automatic hyperparameter tuning methods, such as Bayesian optimization [33], particle swarm optimization [34], are helpful to alleviate the tedious tuning procedure. For data at different reference temperatures, the network architecture would be slightly different.

\subsection{Predication Performance Evaluation}
In this part, the influence of temperature on prediction ability is investigated using the proposed method. Taking 10$^\circ$C as the reference temperature, a monitoring model is established according to the procedures given in Subsection III.B. Same as the last section, the first eight drive cycles are used for training and the following one cycle is used for validation to determine the crucial parameters, i.e., the number of retained CVs $R$. The parameter $R$ is determined as 163 as shown in Fig. 6(b), where a knee point is observed. Besides, the specific values of $R$ for the other temperatures to develop monitoring model can be found in Table II.

After the monitoring model is completed, the last discharge cycle at 10$^\circ$C is used as a testing case to validate the developed model, which is tested as normal but not shown here for brevity. Moreover, Fig. 8(a) to Fig. 8(d) visualize the monitoring results at the other temperatures. For 25$^\circ$C (Fig. 8(a)), it is observed that $T^2$ statistics remains below the control limit while $SPE$ statistic reaches an abnormal value at the end of the discharge cycle. For 0$^\circ$C (Fig. 8(b)), the $SPE$ monitoring statistics becomes obviously abnormal already at the beginning of the discharge cycle. For -10$^\circ$C (Fig. 8(c)) and -20$^\circ$C (Fig. 8(d)), more serious abnormalities are observed for both $T^2$ and $SPE$. Abnormal samples are circled by red dashed ellipses. Therefore, changes of temperature exert unignorable influences on estimation ability. The larger the difference between reference temperature and new temperatures is the more dissimilar temporal dynamics will be.

\begin{table*}[ht]
    \tiny
    \centering
	\renewcommand{\arraystretch}{1.2}
	\caption{Estimation results and similarities using the proposed method from the reference temperature -10$^\circ$C to temperature -20$^\circ$C.}
	\label{table_V}
	\begin{center}
		\begin{threeparttable}
			\begin{tabular}{c c c c| c c c| c}
				\toprule
				\multicolumn{1}{c}{\multirow{2}*{\textbf{Transferring case}}} & \multicolumn{3}{c|}{\textbf{Indices}} & \multicolumn{3}{c|}{\multirow{1}*{\textbf{Transfer network configuration}}} & \multicolumn{1}{c}{\multirow{2}*{\textbf{Specific values for parameters}}}\\
				\cline{2-7}
				& {$RMSE$} & {$MAE$} & {$H(\%$)} & {Iteration epochs} & {Training time} & {Network} \\
				\hline
				{From -10$^\circ$C to -20$^\circ$C} & {9.01} & {7.65} & {55.80} & {100} & {2min6sec} & {L(50)N(100)} & {$q=110, \alpha_1=0.25, \alpha_2=0.75$}\\
				{From 10$^\circ$C to -20$^\circ$C}  & {21.59} & {16.96} &	{8.28} & {100} & 1min9sec & {L(50)N(100)} & {$q=15, \alpha_1=0.035, \alpha_2=0.965$}\\
				{Accuracy improvement} & {58.27$\%$} & {54.89$\%$} & {-} & {-} & {-}  & {-} & {-}\\
				\bottomrule
			\end{tabular}
			\begin{tablenotes}
				\footnotesize
				\item[-] denotes null.
			\end{tablenotes}
		\end{threeparttable}
	\end{center}
\end{table*}

\subsection{Transfer SoC Estimation at Varying Temperatures}
Choosing 10$^\circ$C as the reference temperature, the transferability of the proposed method is verified for the other four temperatures concerning SoC estimation. Here, ten discharging cycling data with real SoC values at 10$^\circ$C and the first discharge cycle with real SoC values at the other temperatures are used as the training dataset for four transfer models, respectively. And the second discharge cycle at each temperature is used as corresponding testing data. After the data processing and feature extraction using CVA as shown in the last two sections, we will make progress to select the consistent temporal features between 10$^\circ$C and the other four temperatures. According to the steps given in Subsection III.C, the number of consistent temporal features at 25$^\circ$C, 0$^\circ$C, -10$^\circ$C, and -20$^\circ$C are 200, 108, 45, and 15, respectively. For each temperature, a shared prediction model is trained using consistent temporal features selected at the reference temperature, i.e., 10$^\circ$C. Thereafter, specific temporal features at new temperatures are separated from off-line discharge cycle and then specific prediction models are trained. Table III summarizes the estimation results of the proposed method with TL. Besides, estimation results by directly applying the reference model into the other temperatures are given for comparison, which as to be expected, lead to very large prediction errors. In contrast, using the proposed method with TL estimation errors have been greatly reduced by 79.48$\%$, 82.00$\%$, 56.21$\%$, and 40.08$\%$ at 25$^\circ$C, 0$^\circ$C, -10$^\circ$C, and -20$^\circ$C, respectively. Especially, for 25$^\circ$C and 0$^\circ$C, the estimation results are shown in Fig. 9 for an intuitive presentation. Although errors exist due to the limitation of data amount, the estimation trends meet well with the real case. Besides, more accurate results can be derived if filter algorithms are applied as post-processing. It should be noticed that the network configurations in Table III are tuned for specific prediction ability $F_t(\cdot)$ in Eq. (11).

In addition to the reduction of estimation error, the proposed method contributes to further process understanding. The similarity between reference temperature and new temperatures is quantitatively accounted by taking advantage of consistent CVs. The index $H={\operatorname{var}\left(\mathbf{Z}_{x,q}^{T} \mathbf{Z}_{x,q}\right)}/{\operatorname{var}\left(\mathbf{Z}_{x}^{T} \mathbf{Z}_{x}\right)}$ is defined to reveal the similarity from the perspective of process variances, where $\operatorname{var}(\cdot)$ calculates the variance of the matrix. A larger value of $H$ indicates a higher similarity. As shown in Table III, high similarities are observed for 25$^\circ$C and 0$^\circ$C, which are more than 50$\%$ with respect to $H$. As a result, a better transfer estimation ability is achieved at these two cases. In contrast, similarities are low for -10$^\circ$C and -20$^\circ$C, leading to relatively poor estimation results.

Next, the efficacy of the proposed method is further verified through a transfer procedure from -10$^\circ$C to -20$^\circ$C. Choosing -10$^\circ$C as the reference temperature, the number of consistent temporal correlations between -20$^\circ$C and -10$^\circ$C is 110 using the first discharge cycle at -20$^\circ$C for training. The similarities between these two temperatures are calculated as shown in Table \ref{table_V}, which is 55.80$\%$. For better comparison, the results obtained from 10$^\circ$C to -20$^\circ$C are summarized in Table IV. Using the second discharge cycle at -20$^\circ$C as the testing data, a higher similarity is observed at the new reference temperature, i.e., -10$^\circ$C, and this contributes to estimation error reduced by as much as 58.27$\%$ for $RMSE$ in comparison with the result when the reference temperature is 10$^\circ$C. Similar results can be obtained with respect to index $MAE$.

\section{Conclusion}
In this article, a data-driven SoC estimation framework of Lithium-ion batteries is designed to consider the influences of varying ambient temperatures. The accuracy of reference SoC estimation model based on long short-term memory network has been greatly improved by taking advantage of temporal dynamics. Furthermore, it is more practical and effective in comparison with other algorithms through two achievements: 1) evaluating the influences of temperatures on SoC estimation ability by checking the variability of temporal features; 2) rapid updating the reference SoC estimation using advanced transfer learning when the estimation relationship has been changed. In this way, quantitative indicators contribute to automatically indicate the efficacy of the reference SoC estimation model. With a small number of cycles at new temperatures, estimation results are achieved by inheriting the consistent estimation ability. The proposed approach is exposed to extensive testing provided by a public benchmark and it has proven to be a powerful tool for SoC estimation without the limitation of battery type and material.

Exploiting transfer learning ability, the proposed method enables accurate SoC estimation from the reference temperature to new temperatures. However, efficacy of the proposed method may be challenged if the temperature experiences a large-scale variation, resulting in that new temperatures are far away from all reference temperatures. It may be promising to overcome the aforementioned challenge by introducing a hybrid estimation model with accessible domain knowledge.

\bibliographystyle{Bibliography/IEEEtranTIE}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction} \label{sec:Introduction}
The Black-Scholes equation has been used to price many financial instruments, with different boundary conditions defining different types of
instruments. These boundary conditions are not always conducive to an analytic solution of the Black-Scholes equation, and in practice numerical
methods are often used.

Symmetry methods can be used to solve difficult PDEs such as the Black-Scholes equation. However, boundary conditions can often cause problems.
Usually it is assumed that for a symmetry to admit a given boundary condition, the boundary must be invariant under the symmetry, and also the
function describing the boundary condition should be invariant. Ordinarily, since Lie symmetries are smooth, any non-smooth boundary
condition will not be admitted by a Lie symmetry, and so other methods would have to be used to solve the equation.

It has been shown by Goard \cite{Goard}  that the assumption just mentioned is overly restrictive. We will summarize this in
section~\ref{sec:The Invariant Surface Condition}. Unfortunately, the work of Goard does not overcome the problem arising from
the maturity condition of an option. As is shown in \cite{BlackScholes}, if $p$ represents time, $S$ represents the price of the
underlying asset and $V=V(S,p)$ is the price of a European call option, the boundary condition under consideration is that
when $p=T$ the option price $V(S,p)$, with strike price $K$, should be given by
	\begin{equation}
	V(S,T) = \max\{ S - K, 0\}=\begin{cases} S-K & \mbox{ if } S\geq K\\ 0 & \mbox{ if } S< K, \label{maturity}
		\end{cases}
	\end{equation}
which has a point at which the derivative with respect to $S$ does not exist. If one considers the solution surface to
the Black-Scholes equation, the condition (\ref{maturity}) forces the surface to have a fold or a crease at the point
$(K,T)$. Symmetry techniques do not accommodate such solutions, as one of the assumptions in using symmetry techniques
is that solutions are smooth.

If, however, we introduce a second boundary condition of the form
	\[
	V(S,p) = 0 \mbox{  when  } S=g(p)
	\]
where $g(p)$ is a function such that $g(T)=K$, we avoid the problem because we need not worry about values of $S$ less than $K$ for the first
boundary condition. The understanding here is that $S=g(\t)$ represents a lower barrier; if the option price falls below this value, the option
becomes worthless. In other words, we solve the Black-Scholes equation in the region given by $0\leq p \leq T$ and $S\geq g(p)$, and outside
this region the option price is taken to be zero.

\section{Notation and Preliminaries}
Let ${\bf x}=(x^1, x^2,\ldots , x^n)\in \mathbb{R}^n$ be the independent variable, with co-ordinates $x^i$. In the case of the Black-Scholes
equation, the independent variables are $p$ and $S$, where $S$ is the price of the underlying asset, and $p$ is time, which will be transformed to
${\bf x}=(x,t)$. Let
${\bf u}=(u^1, u^2,\ldots, u^m)\in \mathbb{R}^m$ be the dependent variable, with co-ordinates $u^{\alpha}$. For the Black-Scholes equation,
there is only one dependent variable, the value of the financial instrument in question, $V$, which will be transformed to ${\bf u}=u$. We use subscripts to indicate partial differentiation:
	\[
	u^{\alpha}_x = \frac{\p u^{\alpha}}{\p x}, \quad u^{\alpha}_{xx} = \frac{\p^2 u^{\alpha}}{\p x^2},
	\quad u^{\alpha}_{xt} = \frac{\p^2 u^{\alpha}}{\p x \p t},\mbox{ etc.}
	\]
We also use the notation $\u_{(1)}$ to represent the collection of all first derivatives of $\u$, and similarly
$\u_{(2)}$, $\u_{(3)}$, etc., represent collections of higher order derivatives.

An $n$-th order partial differential equation (PDE) may then be represented as
	\begin{equation}
	F\left(\x, \u, \u_{(1)}, \u_{(2)},\ldots, \u_{(n)} \right) = 0.\label{PDE}
	\end{equation}
A symmetry $X$ of such a PDE is a vector field $X=\xi^i\frac{\p}{\p x^i} + \phi^{\alpha}\frac{\p}{\p u^{\alpha}}$ (where repeated indices indicate
summation) that leaves solutions of the PDE invariant. In practice, this means that $X^{[n]}F|_{F=0} = 0$, where $X^{[n]}$ is the $n$-th
prolongation of $X$ (this is the vector field $X$ with extra terms added to show the action of $X$ on derivatives of $\u$). There are formulae
for calculating the prolongation coefficients but the details will not be relevant in this paper. The coefficients $\xi^i$ and $\phi^{\alpha}$
are assumed to be functions of $\x$ and $\u$ only (although in principle one may consider coefficients that are also functions of derivatives
of $\u$).

It is possible (for example, see \cite{OlverRosenau}) to use a symmetry to reduce either the order or the number of variables in a PDE, by solving the system
	\[
	\frac{d\,x^1}{\xi^1} = \cdots = \frac{d\, x^n}{\xi^n} = \frac{d\, u^1}{\phi^1} = \cdots = \frac{d\, u^m}{\phi^m}
	\]
and substituting solutions into the original equation $F=0$. For a PDE with boundary conditions, however, it is also necessary that the symmetry
satisfies the invariant surface condition.

\section{The Invariant Surface Condition} \label{sec:The Invariant Surface Condition}
We briefly outline the invariant surface condition and some of its implications in this section. For a more detailed explanation, we refer the reader
to \cite{Goard}. Here we assume that $\x=(x,t)$, $\u = (u)$, and
$X=\xi\frac{\p}{\p x} + \tau \frac{\p}{\p t} + \phi\frac{\p}{\p u}$. A solution $u=u(x,t)$ of
(\ref{PDE}) is invariant under the symmetry $X$ if and only if the invariant surface condition holds, i.e.
	\begin{equation}
	\xi(x,t,u)\frac{\p u}{\p x} + \tau(x,t,u)\frac{\p u}{\p t} = \phi(x,t,u).\label{ISC}
	\end{equation}

\subsection{Initial (or Terminal) Conditions}
If we impose the condition $u(x,T) = f(x)$ then substituting this condition into (\ref{ISC}), we get
	\begin{equation}
	\xi(x,T,f(x))f'(x) + \tau(x,T,f(x))\frac{\p u(x,T)}{\p t} = \phi(x,T,f(x)).\label{ISCi}
	\end{equation}
If $X$ leaves the boundary condition invariant, then this condition is automatically satisfied. On the other hand, if $X$ does not leave the
boundary condition invariant, then we solve (\ref{PDE}) for $u_t$, substitute it into (\ref{ISCi}), and solve (\ref{ISCi}) in one of two
ways: if $X$ is given then we can solve to find the most general $f$ allowed by the symmetry; if $f$ is given, then we can find the most
general symmetry of (\ref{PDE}) that admits solutions satisfying the boundary condition.

\subsection{Boundary Conditions}
Alternatively, we may wish to impose a boundary condition of the form
	\[
	u(x,t) = G(t) \mbox{  when  } x=g(t).
	\]
In this case the invariant surface condition (\ref{ISC}) becomes
	\[
	\xi(g(t),t,G(t))u_x + \tau(g(t),t,G(t))u_t = \phi(g(t),t,G(t)),
	\]
which, after some manipulation, amounts to
	\begin{equation}
	\xi_x u_x + \xi_u u^2_x + \xi u_{xx} + \tau_x u_t + \tau_u u_x u_t + \tau u_{xt} = \phi_x + \phi_u u_x. \label{barrier}
	\end{equation}
Again, this can be solved for the boundary condition, given $X$, or solved for $X$ given the boundary condition.

\section{Barrier Options}
Any option must satisfy the Black-Scholes equation, which we write here as
	\begin{equation}
	V_{p} + \frac{1}2 \sigma^2 S^2 V_{SS} + rS V_S - r V = 0. \label{BS}
	\end{equation}
where $V$ is the price of the option, $S$ is the price of the underlying commodity, and $p$ is time.

For a barrier option there will be both a terminal condition (the option matures when $p=T$) and a boundary condition, which we call a barrier. We
assume the barrier is as stated in section~\ref{sec:Introduction}, i.e. $V(S,p) = 0$ when $S=g(p)$ (so $G(p)=0$) and $g(T)=K$.

We could find symmetries and invariant solutions of this equation directly; however there are many parameters, and to simplify calculations,
we transform (\ref{BS}) to the Heat equation. The method below is not the only method; see Gazizov and Ibragimov \cite{GazizovIbragimov}.

\subsection{Conversion to the Heat Equation}
We break the transformation into steps:
\begin{enumerate}
	\item Let $t = \frac{\s^2}{2} (T - p)$. This reverses time, and (\ref{BS}) becomes
		\be
		-\frac{\s^2}{2} V_t + \frac{\s^2}{2} S^2 V_{SS} + rS V_S - rV = 0. \label{step1}
		\ee
	\item Let $S = K e^x$ (we assume that $S>0$; the factor $K$ introduced here simplifies computations when an initial condition is imposed),
	then (\ref{step1}) becomes
		\be
		-\frac{\s^2}{2}V_t + \frac{\s^2}{2}V_{xx} + \left(r-\frac{\s^2}{2} \right)V_x - rV = 0. \label{step2}
		\ee
	\item Let $w = e^{\a x}V$. By choosing $\a = \frac{1}2\left(\frac{2r}{\s^2} - 1\right)$, (\ref{step2}) transforms to
		\be
		-\frac{\s^2}{2}w_t + \frac{\s^2}{2}w_{xx} -\left( \frac{\s^2\a^2}{2} + r \right)w = 0. \label{step3}
		\ee
	\item Let $y=e^{\b t} w$. Choosing $\b = \frac{1}4\left(\frac{2r}{\s^2} + 1\right)^2$ gives us
	$-\frac{\s^2}{2}y_t + \frac{\s^2}{2}y_{xx} = 0$, i.e.
		\be
		y_t = y_{xx}.
		\ee
	\item Finally, let $u = \frac{1}K y$, to make the boundary conditions simpler.
\end{enumerate}

The boundary condition $V(S,T) = \max\{S-K,0 \}$ now becomes
	\[
	u(x, 0 ) = \max\left\{e^{(\a+1)x} - e^{\a x}, 0 \right\}=\begin{cases} e^{(\a+1)x} - e^{\a x} & \mbox{ if } x\geq 0\\ 0 & \mbox{ if } x< 0.
		\end{cases}
	\]

\subsection{Finite Dimensional Symmetries}
The Lie symmetries of the Heat equation $u_t = u_{xx}$ are well-known (for example, see \cite{Olver}):
	\begin{align*}
	X_1 &= 4xt \frac{\p}{\p x}  + 4t^2\frac{\p}{\p t} + u(-2t-x^2)\frac{\p}{\p u}\\
	X_2 &= x\frac{\p}{\p x} + 2t\frac{\p}{\p t} \\
	X_3 &= 2t\frac{\p}{\p x} - ux\frac{\p}{\p u} \\
	X_4 &= \frac{\p}{\p x} \\
	X_5 &= u\frac{\p}{\p u} \\
	X_6 &= \frac{\p}{\p t} \\
	X_{\infty} &= \psi(x,t) \frac{\p}{\p u} \mbox{ where } \psi_t = \psi_{xx}.
	\end{align*}
The most general finite-dimensional symmetry of the Heat equation is a linear combination of the first six of these:
	\begin{dmath*}
	X = [4c_1xt + c_2 x + 2c_3t + c_4]\frac{\p}{\p x} + [4c_1t^2 + 2c_2t + c_6]\frac{\p}{\p t} + [c_1u(-2t-x^2) - c_3ux - c_5 u]\frac{\p}{\p u}.
	\end{dmath*}
We now wish to apply the initial condition with this symmetry to (\ref{ISCi}) to find the most general symmetry that
allows the initial condition. However, due to the point of non-smoothness of the initial condition, we only consider
values $x>0$, bearing in mind that we will also be imposing a barrier condition. The condition
(\ref{ISCi}) becomes
	\begin{dmath*}
	\left( -c_1 x^2 - c_3 x - c_5 \right)\left( e^{(\a +1)x} - e^{\a x} \right) - (c_2 x + c_4) \left( (\a + 1)e^{(\a +1)x} - \a e^{\a x} \right)
			= c_6\left( (\a + 1)^2e^{(\a+1)x} - \a^2e^{\a x} \right)
	\end{dmath*}
which we solve for $c_1,\ldots, c_6$ by comparing coefficients of $e^{(\a+1)x}$, $e^{\a x}$, etc., to get
	\begin{align*}
	c_1 & = c_2 = c_3 = 0, \\
	c_4 & = -(2\a + 1) c_6, \\
	c_5 & = (\a^2 + \a)c_6.
	\end{align*}
Without loss of generality, we set $c_6 = 1$ so that the symmetry operator becomes
	\begin{dmath}
	X = -(2\a + 1)\frac{\p}{\p x} + \frac{\p}{\p t} - (\a^2 + \a) u \frac{\p}{\p u},
	\end{dmath}
and the system
	\[
	\frac{d\, x}{\xi} = \frac{d\, t}{\tau} = \frac{d\, u}{\phi}
	\]
becomes
	\[
	\frac{d\, x}{-2\a - 1} = \frac{d\, t}{1} = \frac{d\, u}{-\a^2 - \a}.
	\]
The left-hand pair can be solved to find
	\[
	I_1 = x + (2\a + 1)t,
	\]
and the right-hand pair can be solved to find
	\[
	u = I_2 e^{(-\a^2 - \a)t},
	\]
where $I_1$ and $I_2$ are constants of integration; the theory of differential equations tells us that these invariants are functionally dependent;
we express this as follows:
	\[
	u = h(I_1)e^{(-\a^2 - \a)t} = h(x + (2\a + 1)t)e^{(-\a^2 - \a)t}.
	\]
The heat equation $u_t = u_{xx}$ now becomes
	\[
	h''(I_1) = (2\a + 1)h'(I_1) - (\a^2 + \a)h(I_1),
	\]
which can be solved to give
	\[
	h(I_1) = A e^{(\a + 1)I_1} + B e^{\a I_1},
	\]
so that
	\[
	u = A e^{(\a + 1) x + (\a + 1)^2 t} + B e^{\a x + \a^2 t}.
	\]
Now $u(x,0) = A e^{(\a + 1)x} + B e^{\a x}$, and comparing this with the initial condition for $x\geq 0$, we see that $A=1$ and $B=-1$.

Next, we transform this back into the original variables $S$, $\t$ and $V$, to get:
	\[
	V = S - Ke^{-r(T-p)}.
	\]
$V$ does indeed solve the Black-Scholes equation, and satisfies the terminal condition for $S\geq K$; by inspection we see that the only
barrier allowed by this solution is
	\[
	S = K e^{-r(T-p)}.
	\]

\subsection{Infinite Dimensional Symmetries}
We now consider symmetries of the form
	\begin{dmath*}
	X = [4c_1xt + c_2 x + 2c_3t + c_4]\frac{\p}{\p x} + [4c_1t^2 + 2c_2t + c_6]\frac{\p}{\p t}
			+ [c_1u(-2t-x^2) - c_3ux - c_5 u + \psi(x,t)]\frac{\p}{\p u},
	\end{dmath*}
where $\psi_t = \psi_{xx}$, with the same initial condition as before, i.e. $u(x, 0 ) = \max\left\{e^{(\a+1)x} - e^{\a x}, 0 \right\}$;
the invariant surface condition (\ref{ISCi}) can now be written as
	\begin{dmath}
	\left( -c_1 x^2 - c_3 x - c_5 \right)\left( e^{(\a +1)x} - e^{\a x} \right) + \psi - (c_2 x + c_4) \left( (\a + 1)e^{(\a +1)x} - \a e^{\a x} \right)
			= c_6\left( (\a + 1)^2e^{(\a+1)x} - \a^2e^{\a x} \right) \label{ISCii}
	\end{dmath}

We start by considering functions $\psi$ of the form $\psi = A(t)e^{\a x} + B(t) e^{(\a + 1)x}$,
and since $\psi$ has to satisfy the Heat equation, we get
	\[
	\psi = k_1 e^{\a x + \a^2 t} + k_2 e^{(\a+1)x + (\a + 1)^2 t}.
	\]
Applying this to (\ref{ISCii}) and solving for arbitrary constants in a similar way to the previous subsection, we find that
	\begin{align*}
	c_1 & = c_2 = c_3 = 0, \\
	k_1 &= -\a c_4 -c_5 - \a^2 c_6, \\
	k_2 &= (\a + 1)c_4 + c_5 + (\a + 1)^2 c_6.
	\end{align*}
This gives rise to a three-dimensional sub-algebra generated by
	\begin{align*}
	X_a & = \frac{\p}{\p x} + \left(-\a e^{\a x + \a^2 t} + (\a + 1)e^{(\a+1)x + (\a + 1)^2 t}  \right)\frac{\p}{\p u},\\
	X_b & = \left(u - e^{\a x + \a^2 t} + e^{(\a+1)x + (\a + 1)^2 t}  \right)\frac{\p}{\p u},\\
	X_c & = \frac{\p}{\p t} + \left(-\a^2 e^{\a x + \a^2 t} + (\a + 1)^2e^{(\a+1)x + (\a + 1)^2 t}  \right)\frac{\p}{\p u}.
	\end{align*}

\subsection{Further Solutions}
We now use a linear combination of $X_a$, $X_b$ and $X_c$. The characteristic equation that we must solve is
	\[
	\frac{dx}{c_4} = \frac{dt}{c_6} = \frac{du}{-c_5 u + \psi(x,t)},
	\]
where
	\begin{dmath*}
	\psi = \left( -\alpha  c_4-c_5-\alpha ^2 c_6 \right) e^{\alpha x + \alpha ^2 t}
				 + \left((\alpha +1)c_4 + c_5 + (\alpha + 1)^2 c_6 \right) e^{(\alpha +1) x + (\alpha +1)^2 t}.
	\end{dmath*}
We solve $\frac{dx}{c_4} = \frac{dt}{c_6}$ to get $x=\frac{c_4}{c_6} t + I_1$, i.e. $I_1 = x- \frac{c_4}{c_6} t$. The equation
$\frac{dx}{c_4}= \frac{du}{-c_5 u + \psi(x,t)}$ can also be solved in one of three possible ways, depending on the arbitrary constant
$c_4$, $c_5$ and $c_6$. We assume that $\frac{c_5}{c_4}$ does not equal $\a$ or $\a + 1$ (otherwise, after a bit of work, we arrive at a
contradiction), and find
	\begin{dmath*}
	u = I_2 e^{-\frac{c_5}{c_4}x} - \frac{c_5 + \a^2 c_6 + \a c_4}{c_4 \a + c_5} e^{\a x + \a^2 t}
						+	\frac{(\a + 1)c_4 + c_5 + (\a+1)^2 c_6}{c_4(\a + 1) + c_5} e^{(\a + 1) x + (\a + 1)^2 t}.
	\end{dmath*}
We substitute this into the heat equation, and solve for $u$ using the fact that $I_2$ must be a function of $I_1$ to find that
	\begin{dmath*}
	u = A_1 e^{\frac{\left(c_4 - \sqrt{c_4^2-4 c_5 c_6}\right) \left(c_4 t-c_6 x\right)}{2 c_4^2 c_6^2}}
		+A_2 e^{\frac{\left(c_4^2-2 c_5 c_6 -\sqrt{c_4^4-4 c_4^2 c_5 c_6}\right) \left(c_4 t-c_6 x\right)}{2 c_4 c_6^2}}
		+\frac{\alpha ^2 c_6-\alpha  c_4+c_5}{\alpha  c_4+c_5} e^{\alpha x + \alpha^2 t}
		+ \frac{ (\alpha +1)^2 c_6+(\alpha +1) c_4-c_5 }{(\alpha +1) c_4+c_5} e^{ (\a + 1)x + (\a+1)^2 t}.
	\end{dmath*}
We require that $u(x,0) = e^{(\a+1)x} - e^{\a x}$ ($x>0$), and looking at the exponentials in our expression for $u$, we see that either $A_1=A_2=0$
(which leads to a contradiction), or
$c_4 = -(2 \alpha  +1)c_6$ and $c_5 = \alpha  (\alpha +1) c_6$.  We arrive at the solution
	\[
	u = e^{(\a+1)x + (\a+1)^2 t} - e^{\a x + \a^2 t},
	\]
which is the same solution that we found before.

\section{Conclusion}
We have found a solution to the Black-Scholes equation given by
	\[
	V = \begin{cases} S - Ke^{-r(T-p)} &\mbox{ if }\, 0\leq p \leq T,\, S\geq Ke^{-r(T-p)} \\ 0 &\mbox{ otherwise.}\end{cases}
	\]

At first glance, it may seem that the use of the infinite-dimensional symmetry $\psi\frac{\p}{\p u}$ does not gain any extra solutions;
however, we assumed a very specific form of $\psi$; other forms may generate other solutions. We note that although use of the
invariant surface condition was crucial to finding invariant solutions that satisfied the given boundary conditions, the invariant surface
condition is necessary, but not sufficient; the invariant surface condition restricted the choices of symmetries, but in some cases,
the boundary condition(s) restricted the choices further.

In conclusion, we have shown specifically how to use symmetry techniques to find prices of barrier options for a non-standard barrier. More generally,
we have demonstrated that PDEs with non-smooth boundary conditions can be solved using symmetry techniques if further boundary conditions are
imposed to remove any non-smooth points.