# Möbius Disjointness Conjecture on the product of
a circle and the Heisenberg nilmanifold

Jing Ma & Ronghui Wu  Jing Ma
School of Mathematics
Jilin University
Changchun 130012
P. R. China<EMAIL_ADDRESS>Ronghui Wu
School of Mathematics
Jilin University
Changchun 130012
P. R. China<EMAIL_ADDRESS>
###### Abstract.

Let $\mathbb{T}$ be the unit circle and $\Gamma\backslash G$ the
$3$-dimensional Heisenberg nilmanifold. We prove that the Möbius function is
linearly disjoint from a class of distal skew products on
$\mathbb{T}\times\Gamma\backslash G$. These results generalize a recent work
of Huang-Liu-Wang.

###### Key words and phrases:

Möbius Disjointness Conjecture, distal flow, skew product, Heisenberg
nilmanifold, measure complexity.

##  1\. Introduction

Let $(X,T)$ be a flow, where $X$ is a compact metric space and $T:X\rightarrow
X$ a continuous map. Let $\mu(n)$ be the Möbius function, that is $\mu(n)=0$
if $n$ is not square-free, and $\mu(n)=(-1)^{q}$ if $n$ is a product of $q$
distinct primes. $\mu$ is said to be linearly disjoint from $(X,T)$ if

$\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{n\leqslant N}\mu(n)f(T^{n}x)=0$

for any $f\in C(X)$ and any $x\in X$. The Möbius Disjointness Conjecture of
Sarnak [16, 17] states that the function $\mu$ is linearly disjoint from every
$(X,T)$ whose entropy is $0$. This conjecture has been established for many
cases, and we refer to the survey paper [4] for recent progresses. An
incomplete list for works related to the present paper is: Bourgain [1],
Bourgain-Sarnak-Ziegler [2], Green-Tao [7], Liu-Sarnak [12, 13], Wang [18],
Peckner [15], Huang-Wang-Ye [10], Litman-Wang [11], and Huang-Liu-Wang [9].

Distal flows are typical examples of zero-entropy flows, see Parry [14]. A
flow $(X,T)$ with a compatible metric $d$ is called distal if

$\inf_{n\geqslant 0}d(T^{n}x,T^{n}y)>0$

whenever $x\neq y.$ According to Furstenberg’s structure theorem of minimal
distal flows [6], skew products are building blocks of distal flows.

An example of distal flow is the skew product $T$ on the $2$-torus
$\mathbb{T}^{2}=(\mathbb{R}/\mathbb{Z})^{2}$ given by

(1.1) $T:(x,y)\mapsto(x+\alpha,y+h(x)),$

where $\alpha\in[0,1)$ and $h:\mathbb{T}\rightarrow\mathbb{R}$ is a continuous
function. For dynamical properties of this skew product, see for example
Furstenberg [5]. The Möbius disjointness for the skew product (1.1) was first
studied by Liu and Sarnak in [12, 13]. A result in [12] states that, if $h$ is
analytic with an additional assumption on its Fourier coefficients, then the
Möbius Disjointness Conjecture is true for the skew product
$(\mathbb{T}^{2},T).$ This result holds for all $\alpha$, as is not common in
the KAM theory. The aforementioned additional assumption was removed in Wang
[18]. It has been further generalized by Huang-Wang-Ye [10] to the case that
$h(x)$ is $C^{\infty}$-smooth.

Another example of distal flow is nilsystem. Let $G$ be a nilpotent Lie group
with a discrete cocomapct subgroup $\Gamma$. The group $G$ acts in a natural
way on the homogeneous space $\Gamma\backslash G.$ Fix $h\in G.$ Then the
transformation $T$ given by $T(\Gamma g)=\Gamma gh$ makes $(\Gamma\backslash
G,T)$ a nilsystem. The Möbius Disjointness Conjecture for these nilsystems was
proved by Green and Tao in [7].

Now let $G$ be the $3$-dimensional Heisenberg group with the cocompact
discrete subgroup $\Gamma,$ namely

$G=\tiny{\begin{pmatrix}1&\mathbb{R}&\mathbb{R}\\\ 0&1&\mathbb{R}\\\
0&0&1\end{pmatrix},\qquad\Gamma=\begin{pmatrix}1&\mathbb{Z}&\mathbb{Z}\\\
0&1&\mathbb{Z}\\\ 0&0&1\end{pmatrix}}.$

Then $\Gamma\backslash G$ is the $3$-dimensional Heisenberg nilmanifold. Let
$\mathbb{T}$ be the unit circle. Huang-Liu-Wang [9] define a skew product $S$
on $\mathbb{T}\times\Gamma\backslash G$

(1.2) $S:(t,\Gamma g)\mapsto\left(t+\alpha,\Gamma
g\footnotesize{\begin{pmatrix}1&\varphi_{2}(t)&\psi(t)\\\
0&1&\varphi_{1}(t)\\\ 0&0&1\end{pmatrix}}\right),$

where $\varphi_{1}$, $\varphi_{2}$ and $\psi$ are $C^{\infty}$-smooth periodic
functions with period $1$. They proved that the flow
$(\mathbb{T}\times\Gamma\backslash G,S)$ is distal, and in the case that
$\varphi_{1}=\varphi_{2}=\phi$, the skew product

$S_{0}:(t,\Gamma g)\mapsto\left(t+\alpha,\Gamma
g\footnotesize{\begin{pmatrix}1&\phi(t)&\psi(t)\\\ 0&1&\phi(t)\\\
0&0&1\end{pmatrix}}\right)$

is linearly disjoint with the Möbius function $\mu$.

The entropy of $(\mathbb{T}\times\Gamma\backslash G,S)$ is a zero, so the
Möbius Disjointness Conjecture is expected to hold on
$(\mathbb{T}\times\Gamma\backslash G,S)$. In this manuscript, we generalize
the result of Huang-Liu-Wang [9] by showing that the Möbius Disjointness
Conjecture holds on $(\mathbb{T}\times\Gamma\backslash G,S)$ if $\alpha$ is
rational. Let the skew product $T$ on $\mathbb{T}\times\Gamma\backslash G$ be
given by

(1.3) $T:(t,\Gamma g)\mapsto\left(t+\alpha,\Gamma
g\footnotesize{\begin{pmatrix}1&k\phi(t)&\psi(t)\\\ 0&1&\phi(t)\\\
0&0&1\end{pmatrix}}\right),$

where $\alpha\in[0,1)$, $k\in\mathbb{R}$ and $\phi,\psi$ be
$C^{\infty}$-smooth periodic functions from $\mathbb{R}$ to $\mathbb{R}$ with
period $1$. We also prove that the Möbius Disjointness Conjecture holds on
$(\mathbb{T}\times\Gamma\backslash G,T)$. In particular, the main results of
this manuscript are as follows.

######  Theorem 1.1.

Let $\mathbb{T}$ be the unit circle and $\Gamma\backslash G$ the
$3$-dimensional Heisenberg nilmanifold. Let $\alpha\in\mathbb{Q}\cap[0,1)$,
$\varphi_{1},\varphi_{2},\psi$ be $C^{\infty}$-smooth periodic functions from
$\mathbb{R}$ to $\mathbb{R}$ with period $1$. Let $S$ be the skew product on
$\mathbb{T}\times\Gamma\backslash G$ defined in (1.2). Then the Möbius
function $\mu$ is linearly disjoint from $(\mathbb{T}\times\Gamma\backslash
G,S)$.

######  Theorem 1.2.

Let $\mathbb{T}$ be the unit circle and $\Gamma\backslash G$ the
$3$-dimensional Heisenberg nilmanifold. Let $\alpha\in[0,1)$, $k\in\mathbb{R}$
and $\phi,\psi$ be $C^{\infty}$-smooth periodic functions from $\mathbb{R}$ to
$\mathbb{R}$ with period $1$ such that

(1.4) $\int_{0}^{1}\phi(t)dt=0.$

Let the skew product $T$ on $\mathbb{T}\times\Gamma\backslash G$ be given by
(1.3). Then the Möbius function $\mu$ is linearly disjoint from
$(\mathbb{T}\times\Gamma\backslash G,T)$.

Notations. We list some notations that we use in this paper. We write $e(x)$
for $e^{2\pi ix}$, and write $\|x\|$ for the distance between $x$ and the
nearest integer, that is

$\|x\|=\min_{n\in\mathbb{Z}}|x-n|.$

For positive $A,$ the notations $B=O(A)$ or $B\ll A$ mean that there exists a
positive constant $c$ such that $|B|\leqslant cA$. If the constant $c$ depends
on a parameter $b$, we write $B=O_{b}(A)$ or $B\ll_{b}A$. The notation
$A\asymp B$ means that $A\ll B$ and $B\ll A$. For a topological space $X,$ we
use $C(X)$ to denote the set of all continuous complex-valued functions on
$X$. If $X$ is a smooth manifold and $r\geqslant 1$ is an integer, then we use
$C^{r}(x)$ to denote the set of all $f\in C(X)$ that have continuous $r$-th
derivatives.

##  2\. Proof of Theorem 1.1

In this section, we prove Theorem1.1, where $\alpha$ is a rational number.

The Proposition 2.3 in [9] shows that the $\mathbb{C}$-linear subspace spanned
by two classes of functions $\mathcal{A}$ and $\mathcal{B}$ is dense in
$C(\mathbb{T}\times\Gamma\backslash G)$.

For integers $m,j$ with $0\leqslant j\leqslant m-1$, the functions $\psi_{mj}$
and $\psi^{*}_{mj}$ on $G$ are defined by

$\displaystyle\psi_{mj}\left(\footnotesize{\begin{pmatrix}1&y&z\\\ 0&1&x\\\
0&0&1\end{pmatrix}}\right)=e(mz+jx)\sum_{b\in\mathbb{Z}}e^{-\pi(y+b+\frac{j}{m})^{2}}e(mbx),$
$\displaystyle\psi^{*}_{mj}\left(\footnotesize{\begin{pmatrix}1&y&z\\\
0&1&x\\\
0&0&1\end{pmatrix}}\right)=ie(mz+jx)\sum_{b\in\mathbb{Z}}e^{-\pi(y+b+\frac{j}{m}+\frac{1}{2})^{2}}e\begin{pmatrix}\frac{1}{2}\begin{pmatrix}y+b+\frac{j}{m}\end{pmatrix}+\begin{pmatrix}mbx\end{pmatrix}\end{pmatrix}.$

Then $\psi_{mj}$ and $\psi^{*}_{mj}$ are $\Gamma$-invariant so that they can
be regarded as functions on $\Gamma\backslash G$.

Recall that there is a unique Borel probability measure on $\Gamma\backslash
G$ that is invariant under the right translations, and therefore
$L^{2}(\Gamma\backslash G)$ can be defined. Let $V_{0}$ be the subspace of
$L^{2}(\Gamma\backslash G)$ consisting of all functions $f\in
L^{2}(\Gamma\backslash G)$ satisfying

$f\left(\Gamma g\footnotesize{\begin{pmatrix}1&0&z\\\ 0&1&0\\\
0&0&1\end{pmatrix}}\right)=f(\Gamma g)$

for any $g\in G$ and $z\in\mathbb{R}$. Set $C_{0}=V_{0}\cap C(\Gamma\backslash
G)$.

######  Lemma 2.1.

([9, Proposition 2.3]) Let $\mathcal{A}$ be the subset of
$f\in\mathrm{C}(\mathbb{T}\times\Gamma\backslash G)$ satisfying

$f:\left(t,\Gamma\footnotesize{\begin{pmatrix}1&y&z\\\ 0&1&x\\\
0&0&1\end{pmatrix}}\right)\mapsto
e(\xi_{1}t+\xi_{2}x+\xi_{3}y)\psi\left(\Gamma\footnotesize{\begin{pmatrix}1&y&z\\\
0&1&x\\\ 0&0&1\end{pmatrix}}\right),$

where $\xi_{1},\xi_{2},\xi_{3}\in\mathbb{Z}$ and
$\psi=\psi_{mj},\bar{\psi}_{mj},\psi^{*}_{mj}$ or $\bar{\psi}^{*}_{mj}$ for
some $0\leqslant j\leqslant m-1$. Here $\bar{\psi}_{mj}$ and
$\bar{\psi}^{*}_{mj}$ stand for the complex conjugates of $\psi_{mj}$ and
$\psi^{*}_{mj}$, respectively. Let $\mathcal{B}$ be subset of
$f\in\mathrm{C}(\mathbb{T}\times\Gamma\backslash G)$ satisfying $f((t,\Gamma
g))=f_{1}(t)f_{2}(\Gamma g)$ with $f_{1}\in\mathrm{C}(\mathbb{T})$ and
$f_{2}\in\mathrm{C_{0}}(\Gamma\backslash G)$. Then the $\mathbb{C}$-linear
subspace spanned by $\mathcal{A}\cup\mathcal{B}$ is dense in
$C(\mathbb{T}\times\Gamma\backslash G)$.

Following [9], we should separately consider the two cases, namely
$f\in\mathcal{A}$ and $f\in\mathcal{B}$. The case $f\in\mathcal{A}$ will be
handled by Fourier analysis and a classical result of Hua [8], which is a
generalization of Davenport [3].

######  Lemma 2.2.

([8]) Let
$f(x)=\alpha_{d}x^{d}+\alpha_{d-1}x^{d-1}+\cdots+\alpha_{1}x+\alpha_{0}\in\mathbb{R}[x]$
and $0\leqslant a<q$. Then, for arbitrary $A>0$,

$\sum_{\begin{subarray}{c}n\leqslant N\\\ n\equiv a\hskip-4.83694pt\mod
q\end{subarray}}\mu(n)e(f(n))\ll_{A}\frac{N}{\log^{A}N},$

where the implied constant may depend on $A$, $q$ and $d$, but is independent
of $\alpha_{d},\dots,\alpha_{0}$.

For the case that $f\in\mathcal{A}$ we have the following proposition.

######  Proposition 2.3.

Let $(\mathbb{T}\times\Gamma\backslash G,S)$ be as in Theorem 1.1 with
$\alpha\in\mathbb{Q}\cap[0,1).$ Let $\mathcal{A}$ be as above. Then, for any
$(t_{0},\Gamma g_{0})\in\mathbb{T}\times\Gamma\backslash G,$ any
$f\in\mathcal{A}$ and any $A>0$,

$\sum_{n\leqslant N}\mu(n)f(S^{n}(t_{0},\Gamma
g_{0}))\ll_{A}\frac{N}{\log^{A}N},$

where the implied constant depends on $A$ and $\alpha$ only.

###### Proof.

For simplicity, we only consider a typical $f\in\mathcal{A}$ defined by

(2.1) $f\left(t,\Gamma\footnotesize{\begin{pmatrix}1&y&z\\\ 0&1&x\\\
0&0&1\end{pmatrix}}\right)=e(t+x+y+z)\sum_{l\in\mathbb{Z}}e^{-\pi(y+l)^{2}}e(lx).$

A general $f$ can be treated in the same way.

To compute $f(S^{n}(t_{0},\Gamma g_{0}))$ via $\eqref{2}$, we define

$\displaystyle S_{1}(n;t)$ $\displaystyle=\sum_{l=0}^{n-1}\varphi_{1}(\alpha
l+t),\qquad S_{2}(n;t)=\sum_{l=0}^{n-1}\varphi_{2}(\alpha l+t),$
$\displaystyle S_{3}(n;t)$ $\displaystyle=\sum_{l=0}^{n-1}\psi(\alpha l+t),\
\qquad S_{4}(n;t)=\sum_{l=0}^{n-2}\varphi_{2}(\alpha
l+t)\sum_{j=l+1}^{n-1}\varphi_{1}(\alpha j+t)$

for $n\geqslant 1$, $t\in\mathbb{T}$ and set
$S_{1}(0;t)=S_{2}(0;t)=S_{3}(0;t)=S_{4}(0;t)=S_{4}(1;t)=0$ for simplicity. A
straightforward calculation gives

(2.2) $S^{n}:(t_{0},\Gamma g_{0})\mapsto(t_{0}+n\alpha,\Gamma g_{n}),$

where

(2.3)
$g_{n}:=g_{0}{\footnotesize{\begin{pmatrix}1&S_{2}(n;t_{0})&S_{3}(n;t_{0})+S_{4}(n;t_{0})\\\
0&1&S_{1}(n;t_{0})\\\ 0&0&1\end{pmatrix}}}.$

Now write

$g_{0}=\footnotesize{\begin{pmatrix}1&y_{0}&z_{0}\\\ 0&1&x_{0}\\\
0&0&1\end{pmatrix}},\qquad g_{n}=\footnotesize{\begin{pmatrix}1&y_{n}&z_{n}\\\
0&1&x_{n}\\\ 0&0&1\end{pmatrix},}$

where we may assume $x_{0},y_{0},z_{0}\in[0,1)$ without loss of generality, so
that $\eqref{4'}$ gives

$\displaystyle x_{n}$ $\displaystyle=x_{0}+S_{1}(n;t_{0}),$ $\displaystyle
y_{n}$ $\displaystyle=y_{0}+S_{2}(n;t_{0}),$ $\displaystyle z_{n}$
$\displaystyle=z_{0}+y_{0}S_{1}(n;t_{0})+S_{3}(n;t_{0})+S_{4}(n;t_{0}).$

Substituting $\eqref{4'}$ into $\eqref{3'}$ and combining with (2.1), we
obtain that

(2.4) $\begin{split}f(S^{n}(t_{0},\Gamma
g_{0}))&=f\left(t_{0}+n\alpha,\Gamma\footnotesize{\begin{pmatrix}1&y_{n}&z_{n}\\\
0&1&x_{n}\\\ 0&0&1\end{pmatrix}}\right)\\\
&=e(t_{0}+n\alpha+x_{n}+y_{n}+z_{n})\sum_{l\in\mathbb{Z}}e^{-\pi(y_{n}+l)^{2}}e(lx_{n}).\end{split}$

To analyze $\eqref{5'}$, we define
$\omega:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ by

$\omega(u,v)=\sum_{l\in\mathbb{Z}}e^{-\pi(v+y_{0}+l)^{2}}e\big{(}l\big{(}u+x_{0}\big{)}\big{)}.$

Then $\omega$ is an double periodic analytic function which is bounded on
$\mathbb{R}\times\mathbb{R}$. With this function $\omega$ we can rewrite (2.4)
as

(2.5) $\begin{split}f(S^{n}(t_{0},\Gamma
g_{0}))=\rho&\omega\big{(}S_{1}(n;t_{0}),S_{2}(n;t_{0})\big{)}\\\ &\times
e\Big{(}n\alpha+(y_{0}+1)S_{1}(n;t_{0})+S_{2}(n;t_{0})+S_{3}(n;t_{0})+S_{4}(n;t_{0})\Big{)},\end{split}$

where $\rho:=e(t_{0}+x_{0}+y_{0}+z_{0})$.

Recall that $\alpha\in\mathbb{Q}\cap[0,1)$ in the present situation, so that
we can write $\alpha=a/q$ with $0\leqslant a<q$ and $(a,q)=1.$ Thus for any
periodic function $h$ with period $1,$ we have
$h(l_{1}\alpha+t_{0})=h(l_{2}\alpha+t_{0})$ whenever $l_{1}\equiv
l_{2}\hskip-4.2679pt\mod q$. For any $0\leqslant b<q$ and any periodic
function $h$ with period $1$, define

$\gamma(h,b)=\sum_{l=0}^{b-1}h(l\alpha+t_{0})$

and set $\gamma(h)=\gamma(h,q)/q$. For any $n\equiv b\hskip-4.2679pt\mod q$,
write $n-1\equiv b^{\prime}\hskip-4.2679pt\mod q$, $0\leqslant
b,b^{\prime}<q$. Then

$\displaystyle S_{1}(n;t_{0})$
$\displaystyle=(n-b)\gamma(\varphi_{1})+\gamma(\varphi_{1},b),$ $\displaystyle
S_{2}(n;t_{0})$
$\displaystyle=(n-b)\gamma(\varphi_{2})+\gamma(\varphi_{2},b),$ $\displaystyle
S_{3}(n;t_{0})$ $\displaystyle=(n-b)\gamma(\psi)+\gamma(\psi,b)$

and

$\displaystyle S_{4}(n;t_{0})$
$\displaystyle=\sum_{l=0}^{n-2}\varphi_{2}(\alpha
l+t_{0})\sum_{j=l+1}^{n-1}\varphi_{1}(\alpha j+t_{0})$
$\displaystyle=\Big{(}\frac{n-1-q-b^{\prime}}{q}+\cdots+1\Big{)}S_{1}(q;t_{0})S_{2}(q;t_{0})+\sum_{l=0}^{b^{\prime}-1}\varphi_{2}(\alpha
l+t_{0})\sum_{j=l+1}^{b^{\prime}}\varphi_{1}(\alpha j+t_{0})$
$\displaystyle\hskip
17.07164pt+(n-1-b^{\prime})q^{-1}\sum_{l=0}^{q-1}\varphi_{2}(\alpha
l+t_{0})\sum_{j=l+1}^{q+b^{\prime}}\varphi_{1}(\alpha j+t_{0})$

are all real-valued polynomials in $n$ of degree lease than $3$ with
coefficients depending on $\alpha$, $b$, $b^{\prime}$ and $m$. It follows from
this and (2.5) that

$f(S^{n}(t_{0},\Gamma
g_{0}))=\rho\omega\big{(}S_{1}(n;t_{0}),S_{2}(n;t_{0})\big{)}e\big{(}P_{b,b^{\prime}}(n)\big{)},$

where $P_{b,b^{\prime}}(n)$ is a real-valued polynomial in $n$ of degree lease
than $3$ with coefficients depending on $\alpha,b,b^{\prime}$ and $m$, and
then

(2.6) $\sum_{n\leqslant N}\mu(n)f(S^{n}(t_{0},\Gamma
g_{0}))=\rho\sum_{b=0}^{q-1}\sum_{\begin{subarray}{c}n\leqslant N\\\ n\equiv
b\hskip-4.2679pt\mod q\end{subarray}}\mu(n)\tilde{\omega}(n)e(P(n,b)),$

where $\tilde{\omega}(n)=\omega\big{(}S_{1}(n;t_{0}),S_{2}(n;t_{0})\big{)}$ is
a bounded smooth weight. Hence, by Lemma 2.2 we have

(2.7) $\sum_{\begin{subarray}{c}n\leqslant N\\\ n\equiv b\hskip-5.12149pt\mod
q\end{subarray}}\mu(n)e(P(n,b))\ll_{A}\frac{N}{\log^{A}N}$

for arbitrary $A>0$, where the implied constant depending on $q$ (hence on
$\alpha$) and $A$ only. Substituting this to (2.6) we obtain the desired
estimate. ∎

The other case $f\in\mathcal{B}$ can be reduced to the case of skew products
on $\mathbb{T}^{3}$ which is already known as [9, Corollary 3.2].

######  Lemma 2.4.

([9, Corollary 3.2]) Let $\alpha\in\mathbb{R}$, and let
$h_{1},h_{2}:\mathbb{T}\rightarrow\mathbb{R}$ be $C^{\infty}$-smooth
functions, $T:\mathbb{T}^{3}\rightarrow\mathbb{T}^{3}$ be given by
$T(x,y,z)=(x+\alpha,y+h_{1}(x),z+h_{2}(x))$. Then the Möbius Disjointness
Conjecture holds for $(\mathbb{T}^{3},T)$.

For the case that $f\in\mathcal{B}$ we have the following proposition.

######  Proposition 2.5.

Let $\mathcal{B}\subset C(\mathbb{T}\times\Gamma\backslash G)$ as above. Let
$S$ be as in Theorem 1.1 and let $f\in\mathcal{B}$. Then, for any
$(t_{0},\Gamma g_{0})\in\mathbb{T}\times\Gamma\backslash G,$

$\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{n=1}^{\mathrm{N}}\mu(n)f(S^{n}(t_{0},\Gamma
g_{0}))=0.$

###### Proof.

Let $\tilde{S}:\mathbb{T}^{3}\rightarrow\mathbb{T}^{3}$ be given by

$\tilde{S}:(t,x,y)\mapsto(t+\alpha,x+\varphi_{1}(t),y+\varphi_{2}(t)),$

and $\pi$ be the projection of $\mathbb{T}\times\Gamma\backslash G$ onto
$\mathbb{T}^{3}$ given by

$\pi:\left(t,\Gamma\footnotesize{\begin{pmatrix}1&y&z\\\ 0&1&x\\\
0&0&1\end{pmatrix}}\right)\mapsto(t,x,y).$

Then we have $\pi\circ S=\tilde{S}\circ\pi,$ and hence
$(\mathbb{T}^{3},\tilde{S})$ is a topological factor of
$(\mathbb{T}\times\Gamma\backslash G,S)$.

Since $f\in\mathcal{B}$, there are some $f_{1}\in C(\mathbb{T})$ and $f_{2}\in
V_{0}\cap C(\Gamma\backslash G)$ such that $f(t,\Gamma g)=f_{1}(t)f_{2}(\Gamma
g)$. Then for any $z^{\prime}\in\mathbb{R}$,

$f_{2}\left(\Gamma\footnotesize{\begin{pmatrix}1&y&z\\\ 0&1&x\\\
0&0&1\end{pmatrix}}\right)=f_{2}\left(\Gamma\footnotesize{\begin{pmatrix}1&y&z+z^{\prime}\\\
0&1&x\\\ 0&0&1\end{pmatrix}}\right)$

which shows that $f_{2}$ is independent of the $z$-component and induces a
well-defined continuous function $\tilde{f_{2}}\in C(\mathbb{T}^{2})$ given by

$\tilde{f_{2}}(x,y)=f_{2}\left(\Gamma\footnotesize{\begin{pmatrix}1&y&z\\\
0&1&x\\\ 0&0&1\end{pmatrix}}\right)$

for any $z\in\mathbb{R}.$ Define $\tilde{f}(t,x,y)\in C(\mathbb{T}^{3})$ by

$\tilde{f}(t,x,y)=f_{1}(t)\tilde{f_{2}}(x,y).$

Then $f(t,\Gamma g)=\tilde{f}\circ\pi(t,\Gamma g)$ for any $(t,\Gamma
g)\in\mathbb{T}\times\Gamma\backslash G.$ Hence

$f(S^{n}(t_{0},\Gamma g_{0}))=\tilde{f}\circ\pi\circ S^{n}(t_{0},\Gamma
g_{0})=\tilde{f}\circ\tilde{S}^{n}\circ\pi(t_{0},\Gamma g_{0})$

for any $n\geqslant 1$, and the sequence $\\{f(S^{n}(t_{0},\Gamma
g_{0}))\\}_{n\geqslant 1}$ is observed in $(\mathbb{T}^{3},\tilde{S})$. The
desired result follows from Lemma 2.4. ∎

_Proof of Theorem 1.1_. Combining Propositions 2.3 and 2.5 we get that Theorem
1.1 holds.∎

##  3\. Measure complexity

To prove Theorem 1.2 for irrational $\alpha$, we will use the criterion given
by Huang-Wang-Ye in [10]. To do this, we need the concept of measure
complexity. In this section, we will collect some concepts and facts from [10]
without proof.

Let $(X,T)$ be a flow, and $M(X,T)$ the set of all $T$-invariant Borel
probability measures on $X$. A metric $d$ on $X$ is said to be compatible if
the topology induced by $d$ is the same as the given topology on $X$. For a
compatible metric $d$ and an $n\in\mathbb{N},$ define

(3.1) $\bar{d}_{n}(x,y)=\frac{1}{n}\sum_{j=0}^{n-1}d(T^{j}x,T^{j}y)$

for $x,y\in X.$ Then, for $\varepsilon>0$, let

$B_{\bar{d}_{n}}(x,\varepsilon)=\\{y\in X:\bar{d}_{n}(x,y)<\varepsilon\\},$

with which we can further define, for $\rho\in M(X,T),$

$s_{n}(X,T,d,\rho,\varepsilon)=\min\\{m\in\mathbb{N}:\exists
x_{1},...,x_{m}\in X,\ \textrm{s. t.}\
\rho(\cup_{j=1}^{m}B_{\bar{d}_{n}}(x_{j},\varepsilon))>1-\varepsilon\\}.$

Let $(X,d,T,\rho)$ be as above, and let $\\{u(n)\\}_{n\geqslant 1}$ be an
increasing sequence satisfying $1\leqslant u(n)\rightarrow\infty$ as
$n\rightarrow\infty.$ We say that the measure complexity of $(X,d,T,\rho)$ is
weaker than $u(n)$ if

$\liminf_{n\rightarrow\infty}\frac{s_{n}(X,T,d,\rho,\varepsilon)}{u(n)}=0$

for any $\varepsilon>0$. In view of Lemma 3.2 this property is independent of
the choice of compatible metrics. Hence we can say instead that the measure
complexity of $(X,T,\rho)$ is weaker than $u(n).$ We say the measure
complexity of $(X,T,\rho)$ is sub-polynomial if the measure complexity of
$(X,T,\rho)$ is weaker than $n^{\tau}$ for any $\tau>0$.

We will use the equivalent condition provided by Huang-Wang-Ye [10] to derive
that the Möbius Disjointness Conjecture holds for our flow
$(\mathbb{T}\times\Gamma\backslash G,T)$ as in Theorem 1.2 for $\alpha$
irrational. The following lemma is the main theorem of [10] which gives a
criterion of the Möbius Disjointness Conjecture.

######  Lemma 3.1.

([10, Theorem 1.1]) If the measure complexity of $(X,T,\rho)$ is sub-
polynomial for any $\rho\in M(X;T)$, then the Möbius Disjointness Conjecture
holds for $(X,T)$.

Let $(X,T)$ and $(Y,S)$ be two flows, $d$ and $d^{\prime}$ the metrics on $X$
and $Y$, respectively. Let $\rho\in M(X,T)$, $\nu\in M(Y,S)$,
$\mathcal{B}_{X}$ and $\mathcal{B}_{Y}$ be the Borel $\sigma$-algebras of $X$
and $Y$, respectively. We say $(X,\mathcal{B}_{X},T,\rho)$ is measurably
isomorphic to $(Y,\mathcal{B}_{Y},S,\nu)$, if there exist $X^{\prime}\subset
X$, $Y^{\prime}\subset Y$ with $\rho(X^{\prime})=\nu(Y^{\prime})=1$,
$TX^{\prime}\subset X^{\prime}$, $SY^{\prime}\subset Y^{\prime}$, and an
invertible measure-preserving map $\theta:X^{\prime}\rightarrow Y^{\prime}$
such that $\theta\circ T(x)=S\circ\theta(x)$ for any $x\in X^{\prime}$. The
following lemma is important when calculating the measure complexity.

######  Lemma 3.2.

([10, Proposition 2.2]) Let $\\{u(n)\\}_{n\geqslant 1}$ be an increasing
sequence satisfying $1\leqslant u(n)\rightarrow\infty$ as $n\rightarrow 1$.
Assume that $(X,\mathcal{B}_{X},T,\rho)$ is measurably isomorphic to
$(Y,\mathcal{B}_{Y},S,\nu)$. Then the measure complexity of $(X,d,T,\rho)$ is
weaker than $u(n)$ if and only if the measure complexity of
$(Y,d^{\prime},S,\nu)$ is weaker than $u(n)$.

Now we need to choose a proper metric on $\mathbb{T}\times\Gamma\backslash G$.
The following facts can be found in Sections 2 and 5 in Green-Tao [7], which
we will directly state without proof.

The lower central series filtration $G_{\bullet}$ on $G$ is the sequence of
closed connected subgroups

$G=G_{1}\supseteq G_{2}\supseteq G_{3}=\\{\textup{id}_{G}\\},$

where

$G_{2}=[G_{1},G]=\begin{pmatrix}\begin{smallmatrix}1&0&\mathbb{R}\\\ 0&1&0\\\
0&0&1\end{smallmatrix}\end{pmatrix}$

and idG is the identity element of $G$. Let $\mathfrak{g}$ be the Lie algebra
of $G$. Let

$X_{1}=\begin{pmatrix}\begin{smallmatrix}0&0&0\\\ 0&0&1\\\
0&0&0\end{smallmatrix}\end{pmatrix},\qquad
X_{2}=\begin{pmatrix}\begin{smallmatrix}0&1&0\\\ 0&0&0\\\
0&0&0\end{smallmatrix}\end{pmatrix},\qquad
X_{3}=\begin{pmatrix}\begin{smallmatrix}0&0&1\\\ 0&0&0\\\
0&0&0\end{smallmatrix}\end{pmatrix}.$

Then $\mathcal{X}=\\{X_{1},X_{2},X_{3}\\}$ is a Mal’cev basis adapted to
$G_{\bullet}$. The corresponding Mal’cev coordinate map
$\kappa:G\rightarrow\mathbb{R}^{3}$ is given by

(3.2) $\kappa\left(\footnotesize{\begin{pmatrix}\begin{smallmatrix}1&y&z\\\
0&1&x\\\ 0&0&1\end{smallmatrix}\end{pmatrix}}\right)=(x,y,z-xy)$

and the metric $d_{G}$ on $G$ is defined by

$d_{G}(g_{1},g_{2}):=\inf\Big{\\{}\sum\limits_{i=0}^{n-1}\min(|\kappa(h^{-1}_{i-1}h_{i})|,|\kappa(h^{-1}_{i}h_{i-1})|)\
:\ h_{0},\dots,h_{n}\in G,h_{0}=g_{1},h_{n}=g_{2}\Big{\\}}$

from which we see that $d_{G}$ is left-invariant. By (3.2), we can get

(3.3)
$\left|\kappa\left(\footnotesize{\begin{pmatrix}\begin{smallmatrix}1&y&z\\\
0&1&x\\\
0&0&1\end{smallmatrix}\end{pmatrix}}\right)\right|\leqslant|x|+|y|+|z|$

with $x,y\in[0,1)$. The above metric on $G$ descends to a metric on
$\Gamma\backslash G$ given by

$d_{\Gamma\backslash G}(\Gamma g_{1},\Gamma
g_{2}):=\inf\big{\\{}d_{G}(g_{1}^{\prime},g_{2}^{\prime})\ :\
g_{1}^{\prime},g_{2}^{\prime}\in G,\Gamma g_{1}=\Gamma g_{1}^{\prime},\Gamma
g_{2}=\Gamma g_{2}^{\prime}\big{\\}}$.

It can be proved that $d_{\Gamma\backslash G}$ is indeed a metric on
$\Gamma\backslash G$. Since $d_{G}$ is left-invariant,

(3.4) $d_{\Gamma\backslash G}(\Gamma g_{1},\Gamma
g_{2})=\inf_{\gamma\in\Gamma}d_{G}(g_{1},\gamma g_{2}).$

Finally, take $d_{\mathbb{T}}$ to be the canonical Euclidean metric on
$\mathbb{T}$, and $d=d_{\mathbb{T}\times\Gamma\backslash G}$ the
$l^{\infty}$-product metric of $d_{\mathbb{T}}$ and $d_{\Gamma\backslash G}$
given by

(3.5) $d\big{(}(t_{1},\Gamma g_{1}),(t_{2},\Gamma
g_{2})\big{)}=\max\big{(}d_{\mathbb{T}}(t_{1},t_{2}),d_{\Gamma\backslash
G}(\Gamma g_{1},\Gamma g_{2})\big{)}.$

In view of Lemma 3.2, the choice of compatible metrics does not affect the
measure complexity. Thus the above choice of $d$ is admissible.

##  4\. Proof of Theorem 1.2

In this section, we assume that $\alpha$ is irrational and prove Theorem 1.2
for irrational $\alpha$. Then combine with Theorem 1.1 we get Theorem 1.2.

Let

$\alpha=[0;a_{1},a_{2},\dots,a_{i},\dots]=\frac{1}{a_{1}+\frac{1}{a_{2}+\frac{1}{a_{3}+\cdots}}}$

be the continued fraction expansion of $\alpha$. Then the expansion is
infinite since $\alpha$ is irrational. Let
$l_{i}/q_{i}=[0;a_{1},a_{2},\dots,a_{i}]$ be the $i$-th convergent of
$\alpha$. Let $\mathcal{Q}=\\{q_{i}:i\geqslant 1\\}.$ For $B>2,$ define

$\mathcal{Q}^{\flat}(B)=\\{q_{i}\in\mathcal{Q}:q_{i+1}\leqslant
q_{i}^{B}\\}\cup\\{1\\},\qquad\mathcal{Q}^{\sharp}(B)=\\{q_{i}\in\mathcal{Q}:q_{i+1}>q_{i}^{B}>1\\}.$

Furthermore, we define

$M_{1}(B)=\bigcup_{q_{i}\in\mathcal{Q}^{\sharp}(B)}\\{m\in\mathbb{Z}:q_{i}\leqslant|m|<q_{i+1},q_{i}|m\\}\cup\\{0\\}$

and define $M_{2}(B)=\mathbb{Z}\backslash M_{1}(B)$. Now expand $\phi$ into
Fourier series

$\phi(t)=\sum_{m\in\mathbb{Z}}\hat{\phi}(m)e(mt),$

and further decompose $\phi$ as

$\begin{split}\phi(t)&=\phi_{1}(t)+\phi_{2}(t)\\\ &:=\sum_{m\in
M_{1}(B)}\hat{\phi}(m)e(mt)+\sum_{m\in
M_{2}(B)}\hat{\phi}(m)e(mt).\end{split}$

We call $\phi_{1}$ and $\phi_{2}$ the resonant and non-resonant part of
$\phi,$ respectively. Let $\eta(t)=\phi^{2}(t).$ We can do the same
decompositions for $\eta$ and $\psi,$ getting

(4.1) $\eta(t)=\eta_{1}(t)+\eta_{2}(t),\qquad\psi(t)=\psi_{1}(t)+\psi_{2}(t).$

Note that the above decompositions of $\phi$, $\eta$ and $\psi$ depend on the
parameter $B$, though we do not make it explicit.

In our proof of Theorem 1.2 for irrational $\alpha$ we still need a lemma from
[9].

######  Lemma 4.1.

([9, Lemma 4.2]) Let $B>2$ and let $\\{a(m)\\}_{m\in\mathbb{Z}}$ be a sequence
such that $|a(m)|\ll m^{-2B}$. Then the series

$\sum_{m\in M_{2}(B)}\frac{a(m)}{e(m\alpha)-1}$

is absolutely convergent.

Define

(4.2) $g_{\phi}(t)=\sum_{m\in
M_{2}(B)}\hat{\phi}(m)\frac{e(mt)}{e(m\alpha)-1}.$

Since $\phi$ is assumed to be $C^{\infty}$-smooth, we have
$\hat{\phi}(m)\ll|m|^{-2B}$ for any $B>0$. Therefore, by Lemma 4.1, $g_{\phi}$
is a continuous periodic function with period $1$ and satisfies

(4.3) $g_{\phi}(t+\alpha)-g_{\phi}(t)=\phi_{2}(t).$

Similarly, there exist continuous periodic functions $g_{\eta}$ and $g_{\psi}$
with period $1$ such that

(4.4)
$\eta_{2}(t)=g_{\eta}(t+\alpha)-g_{\eta}(t),\qquad\psi_{2}(t)=g_{\psi}(t+\alpha)-g_{\psi}(t).$

To investigate the resonant part, define

(4.5) $\Phi_{n}(t)=\sum_{l=0}^{n-1}\phi_{1}(l\alpha+t),\quad
H_{n}(t)=\sum_{l=0}^{n-1}\eta_{1}(l\alpha+t),\quad\Psi_{n}(t)=\sum_{l=0}^{n-1}\psi_{1}(l\alpha+t)$

for $n\in\mathbb{N}$ and $t\in\mathbb{T}$. For $n=q_{i}$, Huang-Liu-Wang [9]
proved the following proposition.

######  Lemma 4.2.

([9, Lemma 4.3]) Let $B>2$. Then there exists a positive constant
$C_{1}=C_{1}(B)$ depending on $B$ only, such that the three inequalities

$|\Phi_{q_{i}}(t)-q_{i}\hat{\phi}(0)|\leqslant C_{1}q_{i}^{-B+1};\ \
|H_{q_{i}}(t)-q_{i}\hat{\eta}(0)|\leqslant C_{1}q_{i}^{-B+1};\ \
|\Psi_{q_{i}}(t)-q_{i}\hat{\psi}(0)|\leqslant C_{1}q_{i}^{-B+1}$

hold simultaneously for all $t\in\mathbb{T}$ and all
$q_{i}\in\mathbb{Q}^{\sharp}(B)$.

Now we are ready to prove Theorem 1.2 for irrational $\alpha$ and the skew
product $T$ defined in (1.3). We will prove the following proposition.

######  Proposition 4.3.

Let $(\mathbb{T}\times\Gamma\backslash G,T)$ be as in Theorem 1.2 with
$\alpha$ irrational. Then the measure complexity of
$(\mathbb{T}\times\Gamma\backslash G,T)$ is sub-polynomial for any $\rho\in
M(\mathbb{T}\times\Gamma\backslash G,T)$.

###### Proof.

Fix $\tau>0.$ We want to show that, for any $\varepsilon>0,$

$\liminf_{n\rightarrow\infty}\frac{s_{n}(\mathbb{T}\times\Gamma\backslash
G,T,d,\rho,\varepsilon)}{n^{\tau}}=0.$

Without loss of generality, we assume that both $\tau$, $\varepsilon<10^{-2}$,
and $\tau^{-1}$, $\varepsilon^{-1}$ are integers. Set $B=8\tau^{-1}+1$.

Firstly, assume that $\mathcal{Q}^{\sharp}(B)$ is infinite. Construct a
transformation $R:\mathbb{T}\times\Gamma\backslash
G\rightarrow\mathbb{T}\times\Gamma\backslash G$ as

(4.6) $R:(t,\Gamma g)\mapsto\left(t,\Gamma
g\footnotesize{\begin{pmatrix}1&kg_{\phi}(t)&\frac{1}{2}kg_{\phi}^{2}(t)-\frac{1}{2}kg_{\eta}(t)+g_{\psi}(t)\\\
0&1&g_{\phi}(t)\\\ 0&0&1\end{pmatrix}}\right).$

Write $T_{1}=R^{-1}\circ T\circ R.$ Then a straightforward calculation gives

$\begin{split}T_{1}(t,\Gamma g)&=R^{-1}\circ T\circ R(t,\Gamma g)\\\
&=R^{-1}\circ T\left(t,\Gamma
g\footnotesize{\begin{pmatrix}1&kg_{\phi}(t)&\frac{1}{2}kg_{\phi}^{2}(t)-\frac{1}{2}kg_{\eta}(t)+g_{\psi}(t)\\\
0&1&g_{\phi}(t)\\\ 0&0&1\end{pmatrix}}\right)\\\
&=R^{-1}\hskip-2.84526pt\left(\hskip-2.84526ptt+\alpha,\Gamma
g\footnotesize{\begin{pmatrix}1&k\phi(t)+kg_{\phi}(t)&\psi(t)+kg_{\phi}(t)\phi(t)+\frac{1}{2}kg_{\phi}^{2}(t)-\frac{1}{2}kg_{\eta}(t)+g_{\psi}(t)\\\
0&1&\phi(t)+g_{\phi}(t)\\\ 0&0&1\end{pmatrix}}\hskip-2.84526pt\right)\\\
&=\left(t+\alpha,\Gamma
g\footnotesize{\begin{pmatrix}1&k\phi(t)+kg_{\phi}(t)-kg_{\phi}(t+\alpha)&\omega\\\
0&1&\phi(t)+g_{\phi}(t)-g_{\phi}(t+\alpha)\\\
0&0&1\end{pmatrix}}\right),\end{split}$

where

$\begin{split}\omega&=\frac{1}{2}kg_{\phi}^{2}(t+\alpha)+\frac{1}{2}kg_{\phi}^{2}(t)-kg_{\phi}(t)g_{\phi}(t+\alpha)+\frac{1}{2}kg_{\eta}(t+\alpha)\\\
&\hskip 22.76219pt-
g_{\psi}(t+\alpha)-k\phi(t)g_{\phi}(t+\alpha)+\psi(t)+k\phi(t)g_{\phi}(t)-\frac{1}{2}kg_{\eta}(t)+g_{\psi}(t).\end{split}$

So $\omega$ can be simplified as

(4.7)
$\begin{split}\omega&=\frac{1}{2}k(g_{\phi}(t+\alpha)-g_{\phi}(t))^{2}+\frac{1}{2}k(g_{\eta}(t+\alpha)-g_{\eta}(t))\\\
&\hskip 130.88268pt-k\phi(t)(g_{\phi}(t+\alpha)-g_{\phi}(t))+\psi_{1}(t)\\\
&=\frac{1}{2}k\phi_{2}^{2}(t)+\frac{1}{2}k\eta_{2}(t)-k\phi(t)\phi_{2}(t)+\psi_{1}(t)\\\
&=\frac{1}{2}k(\phi(t)-\phi_{2}(t))^{2}-\frac{1}{2}k\phi^{2}(t)+\frac{1}{2}k\eta_{2}(t)+\psi_{1}(t)\\\
&=\frac{1}{2}k\phi_{1}^{2}(t)-\frac{1}{2}k\eta_{1}(t)+\psi_{1}(t),\end{split}$

where we have used (4.2), (4.4) and (4.3). It follows that

(4.8) $T_{1}:(t,\Gamma g)\mapsto\left(t+\alpha,\Gamma
g\footnotesize{\begin{pmatrix}1&k\phi_{1}(t)&\frac{1}{2}k\phi_{1}^{2}(t)-\frac{1}{2}k\eta_{1}(t)+\psi_{1}(t)\\\
0&1&\phi_{1}(t)\\\ 0&0&1\end{pmatrix}}\right)$

and then

$T_{1}^{n}:(t,\Gamma g)\mapsto\left(t+n\alpha,\Gamma
g\footnotesize{\begin{pmatrix}1&k\Phi_{n}(t)&\frac{1}{2}k\Phi_{n}^{2}(t)-\frac{1}{2}kH_{n}(t)+\Psi_{n}(t)\\\
0&1&\Phi_{n}(t)\\\ 0&0&1\end{pmatrix}}\right).$

Clearly, $R$ is a homeomorphism on $\mathbb{T}\times\Gamma\backslash G$. Hence
by Lemma 3.2, we need only to show that the measure complexity of
$(\mathbb{T}\times\Gamma\backslash G,T_{1},v)$ is weak than $n^{\tau}$, where
$v=\rho\circ R$.

Let $C_{1}=C_{1}(B)>0$ be the constant in Lemma 4.2. The functions
$\phi_{1}(t),$ $\eta_{1}(t)$ and $\psi_{1}(t)$ are Lipschitz continuous, and
therefore there exists $L>0$ such that

(4.9) $\begin{split}&|\phi_{1}(t_{1})-\phi_{1}(t_{2})|\leqslant
L\|t_{1}-t_{2}\|,\\\ &|\eta_{1}(t_{1})-\eta_{1}(t_{2})|\leqslant
L\|t_{1}-t_{2}\|,\\\ &|\psi_{1}(t_{1})-\psi_{1}(t_{2})|\leqslant
L\|t_{1}-t_{2}\|\end{split}$

for any $t_{1},t_{2}\in\mathbb{T}.$ We also assume that $L$ is large enough
such that $L>\varepsilon^{-1}.$ Moreover, since $\phi_{1}(t),$ $\eta_{1}(t)$
and $\psi_{1}(t)$ are continuous, there exists a constant $C_{2}>0$ such that

$|\phi_{1}(t)|\leqslant C_{2},\qquad|\eta_{1}(t)|\leqslant
C_{2},\qquad|\psi_{1}(t)|\leqslant C_{2}$

for all $t\in\mathbb{T}.$ Since $q_{i}\rightarrow\infty$ as
$i\rightarrow\infty$, there exists a constant $I_{0}>0$ such that
$(C_{1}+C_{2})/q_{i}<\varepsilon$ for all $i\geqslant I_{0}$. For $i\geqslant
I_{0}$, define

$F_{1}(i)=\Big{\\{}t=\frac{j\varepsilon}{Lq_{i}}\in\mathbb{T}\ :\
j=0,1,\dots,\frac{Lq_{i}}{\varepsilon}-1\Big{\\}}$

and

$F_{2}(i)=\left\\{\Gamma
g=\footnotesize{\begin{pmatrix}1&j_{2}(q_{i}^{2}L)^{-1}&j_{3}(q_{i}^{2}L)^{-1}\\\
0&1&j_{1}(q_{i}^{2}L)^{-1}\\\ 0&0&1\end{pmatrix}}\in\Gamma\backslash G\ :\
j_{1},j_{2},j_{3}=0,1,\dots,q_{i}^{2}L-1\right\\}.$

Let

$F(i)=\\{(t,\Gamma g)\in\mathbb{T}\times\Gamma\backslash G\ :\ t\in
F_{1}(i),\;\Gamma g\in F_{2}(i)\\}.$

Then $\sharp F(i)=\varepsilon^{-1}L^{4}q_{i}^{7}.$

Now assume that $q_{i}\in\mathcal{Q}^{\sharp}(B)$ with $i\geqslant I_{0}$ and
set

(4.10) $n_{i}=q_{i}^{B-1}.$

Then any positive integer $m\leqslant n_{i}$ can be uniquely written as

(4.11) $m=a_{m}q_{i}+b_{m}$

with $0\leqslant b_{m}<q_{i}$ and $a_{m}\leqslant q_{i}^{B-2}.$ By the
definition of $F(i),$ for any

$(t,\Gamma g)=\left(t,\Gamma\footnotesize{\begin{pmatrix}1&y&z\\\ 0&1&x\\\
0&0&1\end{pmatrix}}\right)\in\mathbb{T}\times\Gamma\backslash G$

with $x,y,z\in[0,1),$ there exists

$(t^{*},\Gamma
g^{*})=\left(t^{*},\Gamma\footnotesize{\begin{pmatrix}1&y^{*}&z^{*}\\\
0&1&x^{*}\\\ 0&0&1\end{pmatrix}}\right)\in F(i)$

such that $\|t-t^{*}\|\leqslant\varepsilon/(Lq_{i})$ and

(4.12) $\max\\{|x-x^{*}|,\ |y-y^{*}|,\
|z-z^{*}|\\}\leqslant\frac{1}{q_{i}^{2}L}.$

We want to show that $d(T_{1}^{m}(t,\Gamma g),T_{1}^{m}(t^{*},\Gamma g^{*}))$
is small for any $m\leqslant n_{i}$, where $n_{i}$ is as in $\eqref{14}$. Let

$Y(m)=\footnotesize{\begin{pmatrix}1&k\Phi_{m}(t)&\frac{1}{2}k\Phi_{m}^{2}(t)-\frac{1}{2}kH_{m}(t)+\Psi_{m}(t)\\\
0&1&\Phi_{m}(t)\\\ 0&0&1\end{pmatrix}}$

and

$Y^{*}(m)=\footnotesize{\begin{pmatrix}1&k\Phi_{m}(t^{*})&\frac{1}{2}k\Phi_{m}^{2}(t^{*})-\frac{1}{2}kH_{m}(t^{*})+\Psi_{m}(t^{*})\\\
0&1&\Phi_{m}(t^{*})\\\ 0&0&1\end{pmatrix}}.$

Then we have

$T_{1}^{m}(t,\Gamma g)=(t+\alpha m,\Gamma gY(m)),\quad T_{1}^{m}(t^{*},\Gamma
g^{*})=(t^{*}+\alpha m,\Gamma g^{*}Y^{*}(m)).$

Therefore, by our choice of the metric on $\mathbb{T}\times\Gamma\backslash G$
in (3.5), we have

(4.13) $d\big{(}T_{1}^{m}(t,\Gamma g),T_{1}^{m}(t^{*},\Gamma
g^{*})\big{)}\leqslant\max\big{(}\|t-t^{*}\|,d_{\Gamma\backslash G}(\Gamma
gY(m),\Gamma g^{*}Y^{*}(m))\big{)}.$

Since the term $\|t-t^{*}\|$ can be arbitrarily small as
$q_{i}\rightarrow\infty$, it remains only to bound the last term in (4.13). By
the triangle inequality and (3.4) we get

$\displaystyle d_{\Gamma\backslash G}(\Gamma gY(m),\Gamma g^{*}Y^{*}(m))$
$\displaystyle\leqslant d_{\Gamma\backslash G}(\Gamma g^{*}Y(m),\Gamma
gY(m))+d_{\Gamma\backslash G}(\Gamma g^{*}Y^{*}(m),\Gamma g^{*}Y(m))$ (4.14)
$\displaystyle\leqslant d_{G}(g^{*}Y(m),gY(m))+d_{G}(g^{*}Y^{*}(m),g^{*}Y(m))$
$\displaystyle=d_{G}(g^{*}Y(m),gY(m))+d_{G}(Y^{*}(m),Y(m)),$

where the last equality follows from the left invariance of $d_{G}$.
Furthermore, from the definition of $d_{G}$ we get

(4.15)
$\begin{split}d_{G}(g^{*}Y(m),gY(m))&\leqslant|\kappa(Y(m)^{-1}g^{-1}g^{*}Y(m))|,\\\
d_{G}(Y^{*}(m),Y(m))&\leqslant|\kappa(Y(m)^{-1}Y^{*}(m))|,\end{split}$

where $\kappa$ is the Mal’cev coordinate map (3.2), and $|\cdot|$ is the
$l^{\infty}$-norm on $\mathbb{R}^{3}.$

A straightforward calculation gives

$\displaystyle Y(m)^{-1}$ $\displaystyle g^{-1}g^{*}Y(m)$
$\displaystyle=\footnotesize{\begin{pmatrix}1&-k\Phi_{m}(t)&\frac{1}{2}k\Phi_{m}^{2}(t)+\frac{1}{2}kH_{m}(t)-\Psi_{m}(t)\\\
0&1&-\Phi_{m}(t)\\\
0&0&1\end{pmatrix}\begin{pmatrix}1&y^{*}-y&z^{*}-z+xy-x^{*}y\\\ 0&1&x^{*}-x\\\
0&0&1\end{pmatrix}}$ $\displaystyle\hskip
159.33542pt\times\footnotesize{\begin{pmatrix}1&k\Phi_{m}(t)&\frac{1}{2}k\Phi_{m}^{2}(t)-\frac{1}{2}kH_{m}(t)+\Psi_{m}(t)\\\
0&1&\Phi_{m}(t)\\\ 0&0&1\end{pmatrix}}$
$\displaystyle=\footnotesize{\begin{pmatrix}1&y^{*}-y&z^{*}-z+y(x-x^{*})+\Phi_{m}(t)(kx-
kx^{*}-y+y^{*})\\\ 0&1&x^{*}-x\\\ 0&0&1\end{pmatrix}}.$

Since $x,y\in[0,1)$, from (3.3) we get

$|\kappa(Y(m)^{-1}g^{-1}g^{*}Y(m))|\leqslant\big{(}\max(|k|,1)|\Phi_{m}(t)|+2\big{)}(|x-x^{*}|+|y-y^{*}|+|z-z^{*}|).$

By Lemma 4.2, and noticing that $\hat{\phi}(0)=0$ as in (1.4), we get

$|\Phi_{q_{i}}|\leqslant C_{1}q_{i}^{-B+1}.$

Hence by the definition of $\Phi_{n}(t)$ and $\eqref{15}$, we obtain

$\begin{split}|\Phi_{m}(t)|&\leqslant\sum_{r=0}^{a_{m}-1}|\Phi_{q_{i}}(t+rq_{i}\alpha)|+\sum_{l=0}^{b_{m}}|\phi_{1}(t+(a_{m}q_{i}+l)\alpha)|\leqslant\frac{C_{1}a_{m}}{q_{i}^{B-1}}+C_{2}q_{i}\leqslant\frac{C_{1}}{q_{i}}+C_{2}q_{i}.\end{split}$

Thus by $\eqref{16},$ $\eqref{18}$ and $\eqref{19}$, we obtain

(4.16)
$\begin{split}d_{G}(g^{*}Y(m),gY(m))&\leqslant\big{(}\max(|k|,1)|\Phi_{m}(t)|+2\big{)}(|x-x^{*}|+|y-y^{*}|+|z-z^{*}|)\\\
&\leqslant\big{(}\max(|k|,1)(1+C_{2}q_{i})+2\big{)}\frac{3}{q_{i}^{2}L}\leqslant
6\big{(}\max(|k|,1)+1\big{)}\varepsilon.\end{split}$

The treatment of $d_{G}(Y^{*}(m),Y(m))$ is similar. We calculate that

$\displaystyle Y^{-1}(m)Y^{*}(m)$
$\displaystyle=\footnotesize{\begin{pmatrix}1&-k\Phi_{m}(t)&\frac{1}{2}k\Phi_{m}^{2}(t)+\frac{1}{2}kH_{m}(t)-\Psi_{m}(t)\\\
0&1&-\Phi_{m}(t)\\\ 0&0&1\end{pmatrix}}$ $\displaystyle\hskip
17.07164pt\times\footnotesize{\begin{pmatrix}1&k\Phi_{m}(t^{*})&\frac{1}{2}k\Phi_{m}^{2}(t^{*})-\frac{1}{2}kH_{m}(t^{*})+\Psi_{m}(t^{*})\\\
0&1&\Phi_{m}(t^{*})\\\ 0&0&1\end{pmatrix}}$
$\displaystyle=\footnotesize{\begin{pmatrix}1&k(-\Phi_{m}(t)+\Phi_{m}(t^{*}))&\tilde{\omega}\\\
0&1&-\Phi_{m}(t)+\Phi_{m}(t^{*})\\\ 0&0&1\end{pmatrix}},$

where

$\tilde{\omega}=\frac{1}{2}k(\Phi_{m}(t^{*})-\Phi_{m}(t))^{2}+\frac{1}{2}k(H_{m}(t)-H_{m}(t^{*}))+(\Psi_{m}(t^{*})-\Psi_{m}(t)).$

By Lemma 4.2, (4.9) and (4.11), we get

$\displaystyle|\Phi_{m}(t^{*})-\Phi_{m}(t)|$
$\displaystyle\leqslant\sum_{r=0}^{a_{m}-1}\Big{(}|\Phi_{q_{i}}(t^{*}+rq_{i}\alpha)|+|\Phi_{q_{i}}(t+rq_{i}\alpha)|\Big{)}$
$\displaystyle\hskip
22.76219pt+\sum_{l=0}^{b_{m}}|\phi_{1}(t^{*}+(a_{m}q_{i}+l)\alpha)-\phi_{1}(t+(a_{m}q_{i}+l)\alpha)|$
$\displaystyle\leqslant\frac{C_{1}}{q_{i}}+q_{i}L\|t^{*}-t\|<2\varepsilon.$

Similarly, we can also obtain $|H_{m}(t^{*})-H_{m}(t)|\leqslant 2\varepsilon$
and $|\Psi_{m}(t^{*})-\Psi_{m}(t)|\leqslant 2\varepsilon$. Therefore, by the
definition of $d_{G}$ and (3.3), we get

(4.17) $d_{G}(Y^{*}(m),Y(m))\leqslant|\kappa
Y(m)^{-1}Y^{*}(m)|<8\max(|k|,1)\varepsilon.$

From $\eqref{17},$ $\eqref{18},$ $\eqref{20}$ and $\eqref{21},$ we conclude
that

$d\big{(}T_{1}^{m}(t,\Gamma g),T_{1}^{m}(t^{*},\Gamma
g^{*})\big{)}\leqslant\big{(}14\max(|k|,1)+6\big{)}\varepsilon$

for all $m\leqslant n_{i}$. Here, and in what follows, $n_{i}$ is as in
(4.10). So, by the definition of $\bar{d}_{n}$ in (3.1), we have

$\displaystyle\bar{d}_{n_{i}}((t,\Gamma g),(t^{*},\Gamma g^{*}))$
$\displaystyle=\frac{1}{n_{i}}\sum_{m=0}^{n_{i}-1}d(T_{1}^{m}(t,\Gamma
g),T_{1}^{m}(t^{*},\Gamma
g^{*}))\leqslant\big{(}14\max(|k|,1)+6\big{)}\varepsilon.$

This means that $\mathbb{T}\times\Gamma\backslash G$ can be covered by $\sharp
F(k)=\varepsilon^{-1}L^{4}q_{i}^{7}$ balls of radius
$\delta:=(14\max(|k|,1)+6)\varepsilon$ under the metric $\bar{d}_{n_{i}}$,
since $(t,\Gamma g)$ can be chosen arbitrarily. It follows that

$s_{n_{i}}(\mathbb{T}\times\Gamma\backslash
G,T_{1},d,v,\delta)\leqslant\varepsilon^{-1}L^{4}q_{i}^{7}.$

Since $\mathcal{Q}^{\sharp}(B)$ is infinite, we can let $q_{i}$ tend to
infinity along $\mathcal{Q}^{\sharp}(B),$ getting

$\begin{split}\liminf_{n\rightarrow\infty}\frac{s_{n}(\mathbb{T}\times\Gamma\backslash
G,T_{1},d,v,\delta)}{n^{\tau}}&\leqslant\liminf_{\begin{subarray}{c}i\rightarrow\infty\\\
q_{i}\in\mathcal{Q}^{\sharp}(B),\ i\geqslant
I_{0}\end{subarray}}\frac{s_{n_{i}}(\mathbb{T}\times\Gamma\backslash
G,T_{1},d,v,\delta)}{n_{i}^{\tau}}\\\
&\leqslant\liminf_{\begin{subarray}{c}i\rightarrow\infty\\\
q_{i}\in\mathcal{Q}^{\sharp}(B),\ i\geqslant
I_{0}\end{subarray}}\frac{\varepsilon^{-1}L^{4}q_{i}^{7}}{q_{i}^{8+\tau}}\leqslant\liminf_{\begin{subarray}{c}i\rightarrow\infty\\\
q_{i}\in\mathcal{Q}^{\sharp}(B),\ i\geqslant
I_{0}\end{subarray}}\frac{\varepsilon^{-1}L^{4}}{q_{i}}=0.\end{split}$

Since $\varepsilon$ can be arbitrary small, this means that the measure
complexity of $(\mathbb{T}\times\Gamma\backslash G,T,\rho)$ is weaker that
$n^{\tau}$ when $\mathcal{Q}^{\sharp}(B)$ is infinite.

Finally, we deal with the case that $\mathcal{Q}^{\sharp}(B)$ is finite. Now
$M_{1}(B)$ is also finite. Hence the conclusion of Lemma 4.1 still holds if we
replace $M_{2}(B)$ by $\mathbb{Z}\backslash\\{0\\}.$ Hence the functions
$\tilde{g}_{\phi}(t),$ $\tilde{g}_{\eta}(t),$ and $\tilde{g}_{\psi}(t)$
defined by

$\tilde{g}_{\phi}(t)=\sum_{m\neq
0}\frac{\hat{\phi}(m)e(mt)}{e(m\alpha)-1},\quad\tilde{g}_{\eta}(t)=\sum_{m\neq
0}\frac{\hat{\eta}(m)e(mt)}{e(m\alpha)-1},\quad\tilde{g}_{\psi}(t)=\sum_{m\neq
0}\frac{\hat{\psi}(m)e(mt)}{e(m\alpha)-1}$

are continuous and periodic with period $1$. Thus we can write

$\displaystyle\phi(t)=\tilde{g}_{\phi}(t+\alpha)-\tilde{g}_{\phi}(t),$
$\displaystyle\eta(t)=\hat{\eta}(0)+\tilde{g}_{\eta}(t+\alpha)-\tilde{g}_{\eta}(t),$
$\displaystyle\psi(t)=\hat{\psi}(0)+\tilde{g}_{\psi}(t+\alpha)-\tilde{g}_{\psi}(t).$

Notice that $\hat{\phi}(0)=0,$ and so there is no constant term in the first
equation. Similarly to (4.6), we define
$\tilde{R}:\mathbb{T}\times\Gamma\backslash
G\rightarrow\mathbb{T}\times\Gamma\backslash G$ by

$\tilde{R}:(t,\Gamma g)\mapsto\left(t,\Gamma
g\begin{pmatrix}1&k\tilde{g}_{\phi}(t)&\frac{1}{2}k\tilde{g}_{\phi}^{2}(t)-\frac{1}{2}k\tilde{g}_{\eta}(t)+\tilde{g}_{\psi}(t)\\\
0&1&\tilde{g}_{\phi}(t)\\\ 0&0&1\end{pmatrix}\right).$

Then $\tilde{T}_{1}:=\tilde{R}^{-1}\circ T\circ\tilde{R}$ is given by

$\tilde{T}_{1}:(t,\Gamma g)\mapsto\begin{pmatrix}t+\alpha,\Gamma
g\begin{pmatrix}1&0&-\frac{1}{2}k\hat{\eta}(0)+\hat{\psi}(0)\\\ 0&1&0\\\
0&0&1\end{pmatrix}\end{pmatrix}$

as in $\eqref{12}$. Again by Lemma 3.2, the measure complexity of
$(\mathbb{T}\times\Gamma\backslash G,T,\rho)$ is weaker than $n^{\tau}$ if and
only if the measure complexity of $(\mathbb{T}\times\Gamma\backslash
G,\tilde{T}_{1},v)$ is weaker than $n^{\tau},$ where $v=\rho\circ\tilde{R}$.
However, $d$ is invariant under $\tilde{T}_{1}.$ So we have for any
$n\geqslant 1$ and $\varepsilon>0$ that

$s_{n}(\mathbb{T}\times\Gamma\backslash
G,\tilde{T}_{1},d,v,\varepsilon)=s_{1}(\mathbb{T}\times\Gamma\backslash
G,\tilde{T}_{1},d,v,\varepsilon).$

Since $\mathbb{T}\times\Gamma\backslash G$ is compact, we have
$s_{1}(\mathbb{T}\times\Gamma\backslash
G,\tilde{T}_{1},d,v,\varepsilon)<\infty$ and consequently

$\lim_{n\rightarrow\infty}\frac{s_{n}(\mathbb{T}\times\Gamma\backslash
G,\tilde{T}_{1},d,v,\varepsilon)}{n^{\tau}}=\lim_{n\rightarrow\infty}\frac{s_{1}(\mathbb{T}\times\Gamma\backslash
G,\tilde{T}_{1},d,v,\varepsilon)}{n^{\tau}}=0.$

Hence the measure complexity of $(\mathbb{T}\times\Gamma\backslash G,T,\rho)$
is also weaker than $n^{\tau}$ if $\mathcal{Q}^{\sharp}(B)$ is finite. The
proof is complete. ∎

_Proof of Theorem 1.2_. Note that the the skew product $T$ in Theorem 1.2 is a
special case of the skew product $S$ in Theorem 1.1. Theorem 1.2 follows from
Theorem 1.1, Proposition 4.3 and Lemma 3.1. ∎

Acknowledgments This work is supported by the National Natural Science
Foundation of China (Grant No. 11771252).

## References

  * [1] J. Bourgain, _On the correlation of the Möbius function with rank-one systems_ , J. Anal. Math. 120 (2013), 105-130.
  * [2] J. Bourgain, P. Sarnak and T. Ziegler, _Disjointness of Möbius from horocycle flows_ , in: From Fourier Analysis and Number Theory to Radon Transforms and Geometry, in: Dev. Math., vol. 28, Springer, New York, 2013, 67-83.
  * [3] H. Davenport, _On some infinite series involving arithmetical functions_ , II, Quart. J. Math. 8 (1937), 313-350.
  * [4] S. Ferenczi, J. Kulaga-Przymus and M. Lemanczyk, _Sarnak’s conjecture: What’s new_ , in: Ergodic Theory and Dynamical Systems in Their Interactions with Arithmetics and Combinatorics, Lecture Notes in Math., vol. 2213, Springer, Cham, 2018, 163-235.
  * [5] H. Furstenberg, _Strict ergodicity and transformation of the torus_ , Amer. J. Math., 83 (1961), 573-601.
  * [6] H. Furstenberg, _The structure of distal flows_ , Amer. J. Math., 85 (1963), 477-515.
  * [7] B. Green and T. Tao, _The Möbius function is strongly orthogonal to nilsequences_ , Ann. Math., 175 (2012), 541-566.
  * [8] L. K. Hua, _Additive theory of prime numbers_ , Transl. Math. Monogr. 13, Amer. Math. Soc., Providence, 1965.
  * [9] W. Huang, J. Y. Liu and K. Wang, _Möbius disjointness for skew products on a circle and a nilmanifold_ , Discrete Contin. Dyn. Syst. 41 (2021), 3531–3553.
  * [10] W. Huang, Z. Wang and X. Ye, _Measure complexity and M $\ddot{o}$obius disjointness_, Adv. Math., 347 (2019), 827-858.
  * [11] M. Litman and Z. Wang, _Möbius disjointness for skew products on the Heisenberg nilmanifold_ , Proc. Amer. Math. Soc., 147 (2019), 2033-2043.
  * [12] J. Liu and P. Sarnak, _The Möbius function and distal flows_ , Duke Math. J., 164 (2015), 1353-1399.
  * [13] J. Liu and P. Sarnak, _The Möbius disjointness conjecture for distal flows_ , in Proceedings of the Sixth International Congress of Chinese Mathematician, Vol. I, Adv. Lect. Math. (ALM) 36, Int. Press, Somerville, MA, 2017, 327-335.
  * [14] W. Parry, _Zero entropy of distal and related transformations_ , Topological Dynamics, (Symposium, Colorado State Univ., Ft. Collins, Colo., 1967), 383-389.
  * [15] R. Peckner, _Möbius disjointness for homogeneous dynamics_ , Duke Math. J., 167 (2018), 2745-2792.
  * [16] P. Sarnak, _Three lectures on the Möbius function_ , randomness and dynamics, IAS Lecture Notes, 2009; http://publications.ias.edu.
  * [17] P. Sarnak, _Möbius randomness and dynamics_ , Not. S. Afr. Math. Soc., 43 (2012), 89-97.
  * [18] Z. Wang, _Möbius disjointness for analytic skew products_ , Invent. Math., 209 (2017), 175-196.


# Ubiquitous Robot Control Through Multimodal Motion Capture Using Smartwatch
and Smartphone Data

Fabian C Weigend1, Neelesh Kumar2, Oya Aran2, and Heni Ben Amor1 1Weigend, F &
Ben Amor, H are affiliated with SCAI, Arizona State University $\\{$fweigend,
<EMAIL_ADDRESS>2Kumar, N & Aran, O. are affiliated with Corporate
Functions-R&D, Procter and Gamble<EMAIL_ADDRESS>
###### Abstract

We present an open-source library for seamless robot control through motion
capture using smartphones and smartwatches. Our library features three modes:
Watch Only Mode, enabling control with a single smartwatch; Upper Arm Mode,
offering heightened accuracy by incorporating the smartphone attached to the
upper arm; and Pocket Mode, determining body orientation via the smartphone
placed in any pocket. These modes are applied in two real-robot tasks,
showcasing placement accuracy within 2 cm compared to a gold-standard motion
capture system. WearMoCap stands as a suitable alternative to conventional
motion capture systems, particularly in environments where ubiquity is
essential. The library is available at: www.github.com/wearable-motion-
capture.

##  I INTRODUCTION

As consumer wearables become increasingly ubiquitous, IMU-based motion capture
utilizing smartwatch and smartphone data emerges as a widely accessible
solution [1, 2]. While motion capture from consumer wearables may offer lower
accuracy compared to optical and specialized IMU-based systems, their key
advantage lies in their ubiquity and familiarity. Most users carry these
devices with them consistently, enabling pose tracking anytime and anywhere
[1].

Recently, we showcased the potential of utilizing smartwatch data for robot
control [3, 4]. In this work, we compile these approaches into a library named
Wearable Motion Capture (WearMoCap). It offers three operational modes
tailored to different requirements of precision and portability. We evaluate
WearMoCap in two real-robot tasks.

##  II WearMoCap MODES

We begin by outlining the pose estimation methodology for the three motion
capture modes as depicted in Figure 1.

####  II-1 Watch Only

Utilizing the optimal neural network architecture proposed by [3], an LSTM is
employed to estimate the orientations of the lower arm, ${\bf
q}_{\mathrm{la}}$, and the upper arm, ${\bf q}_{\mathrm{ua}}$, derived from a
sequence of watch sensor data, including calibrated IMU, orientation, and
pressure readings. Subsequently, positional values are computed through
forward kinematics. The Watch Only mode requires the user to maintain a
constant forward-facing direction post-calibration.

####  II-2 Upper Arm

In this mode, the user wears a watch and straps a phone to their upper arm
(see Figure 1). Utilizing an LSTM with three layers of 128 neurons, WearMoCap
estimates ${\bf q}_{\mathrm{la}}$ and ${\bf q}_{\mathrm{ua}}$ from the
combined sensor data from both devices. Joint positions are then estimated
through forward kinematics. In contrast to Watch Only mode, users are free to
turn and change their body orientation.

####  II-3 Pocket

This mode streams the sensor data from the watch and a phone in the user’s
pocket. WearMoCap employs a Differentiable Ensemble Kalman Filter [5] to
update an ensemble of states from previous estimates and the sensor data from
both devices. For the detailed methodology, we refer to [4]. Each ensemble
member characterizes the orientation of the lower arm, ${\bf
q}_{\mathrm{la}}$, the upper arm, ${\bf q}_{\mathrm{ua}}$, and additionally
captures the rotation around the up-axis of the hip, ${\bf q}_{\mathrm{hi}}$.
Joint positions are determined through forward kinematics.

Figure 1: Our three WearMoCap modes track the arm pose with varying accuracy
and flexibility. We evaluated all modes in real-robot tasks.

##  III REAL-ROBOT TASKS

We assessed the practicality of WearMoCap in two human-robot tasks approved by
the ASU IRB as STUDY00018521.

###  III-A Handover

In the handover task, a Universal Robot 5 (UR5) picked up an object from a
table and handed it over to a human subject. Successful completion of this
task relied on WearMoCap accurately tracking the position of the human hand.
Human subjects were seated on a rotating chair facing the UR5. As illustrated
in Step 1 of Figure 2, the tabletop area between the robot and the subject was
divided into three distinct zones. Within these zones, subjects performed two
handovers: one with the hand positioned at a lower and another with the hand
at a higher elevation. Subjects repeated these handovers in all other zones in
a random order. While the subject’s orientation remained fixed in the Watch
Only mode, they had the flexibility to adjust their orientation by rotating
the chair in the other two modes.

To mitigate potential learning effects, we randomized the order of tracking
modes. We recorded the handover distance as the disparity between the hand
position and the cube when the participant completed the handover with a voice
command (Step 5). Additionally, we recorded handover time, representing the
duration taken from Step 2 to Step 5.

Figure 2: Handover Task: Human subjects used WearMoCap to perform handovers
with a UR5. A voice command completed the handover and the robot let go of the
cube. Figure 3: Intervention Task: The UR5 sorted the green and red cubes into
the bin. Human subjects use WearMoCap to interrupt the robot and to place the
green cube at a target location.

###  III-B Intervention

For this task, the UR5 robot autonomously picked up colored cubes (green or
red) and dropped them at target locations of the corresponding color. As
depicted in Figure 3 Step 1, the human subject stood facing the robot, and
there were three possible target locations for the green cube. However, the
robot was not trained to correctly place green cubes. Whenever the robot
selected a green cube, the subject intervened with a voice command and
directed the robot with WearMoCap to place the cube correctly.

Five subjects performed this task with each WearMoCap mode randomly assigned
and for each of the three green target locations. We used OptiTrack to record
the placement distance as the distance between the placed green cube and the
center of its target location. Additionally, we recorded the task completion
time as the duration from issuing the ”Follow me” command in Step 4 until Step
5.

##  IV Results And Conclusions

As our baseline, subjects completed all tasks using OptiTrack. Table I
compares baseline objective task metrics to the performance achieved with
WearMoCap modes. Notably, the Watch Only mode exhibits higher error rates, as
indicated by its greater intervention placement distance (+5.2 cm) and
handover distance (+4.5 cm). Conversely, the Upper Arm mode demonstrates
superior accuracy, with deviations below +2 cm in both tasks. Also, the Pocket
mode surpasses the Watch Only mode in distance metrics, attributable to its
additional degree of freedom for fine-tuning positioning. However, this extra
flexibility resulted in longer task completion times, as subjects had to
navigate adjustments in arm motion alongside changes in body orientation.

TABLE I: Summarized Differences to OptiTrack Baseline Task | Method | Dist. (cm) | Time (s) | Trials
---|---|---|---|---
Handover | OptiTrack |  $6.8$${}\pm{}$ | $1.6$ |  $9.2$${}\pm{}$ | $3.2$ | 24
Watch Only |  +$4.5$${}\pm{}$ | $9.7$ |  +$3.3$${}\pm{}$ | $8.1$ | 24
Pocket |  +$2.2$${}\pm{}$ | $3.6$ |  +$3.4$${}\pm{}$ | $6.3$ | 24
Upper Arm |  +$1.9$${}\pm{}$ | $3.7$ |  +$0.5$${}\pm{}$ | $5.2$ | 24
Intervent. | OptiTrack |  $2.4$${}\pm{}$ | $1.5$ |  $17.5$${}\pm{}$ | $4.5$ | 15
Watch Only |  +$5.2$${}\pm{}$ | $6.0$ |  +$10.5$${}\pm{}$ | $7.8$ | 15
Pocket |  +$2.9$${}\pm{}$ | $4.3$ |  +$11.4$${}\pm{}$ | $10.7$ | 15
Upper Arm |  +$1.7$${}\pm{}$ | $5.2$ |  +$4.0$${}\pm{}$ | $5.6$ | 15

Overall, WearMoCap proved suitable for the evaluated real-robot tasks.
However, variations in task performance among the three modes revealed
application trade-offs.

##### Watch Only

Solely utilizing a smartwatch offers convenience in terms of availability and
setup. However, our findings exhibit considerable increases in placement
deviations and completion times. The suitability of this mode depends on the
task at hand; while the human can mitigate final centimeter deviations during
handovers, the prediction inaccuracies are impractical for high-fidelity
teleoperation control, such as pick-and-place tasks.

##### Upper Arm

Wearing an upper-arm fitness strap is an additional step compared to the other
modes. Yet, our results confirm that an additional IMU on the upper arm
enables accurate pose tracking, as it has been previously established by [6,
7]. With the widespread availability of smart devices, this mode emerges as a
ubiquitous and promising alternative to OptiTrack for robot teleoperation.

##### Pocket

The Pocket mode simply requires users to stow the phone in a pocket and allows
to move around freely. This stands in contrast to the Watch Only mode, where
users must maintain a fixed forward-facing direction. Our results underscore
that the additional tracking of body orientation facilitated more precise
control than in Watch Only mode. Still, results are less precise than with the
Upper Arm mode. The Pocket mode, therefore, balances the precision and
convenience of the other two modes.

In conclusion, WearMoCap stands as a promising alternative to state-of-the-art
motion capture systems. When weighing its trade-offs for suitable tasks, it
effectively transforms smartwatches and smartphones into versatile interfaces,
facilitating intuitive robot control anytime and anywhere.

## References

  * [1] V. Mollyn, R. Arakawa, M. Goel, C. Harrison, and K. Ahuja, “Imuposer: Full-body pose estimation using imus in phones, watches, and earbuds,” in _Proceedings of the 2023 CHI Conference on Human Factors in Computing Systems_ , ser. CHI ’23. Association for Computing Machinery, 2023.
  * [2] M. Liu, S. Yang, W. Chomsin, and W. Du, “Real-time tracking of smartwatch orientation and location by multitask learning,” in _Proceedings of the 20th ACM Conference on Embedded Networked Sensor Systems_ , ser. SenSys ’22. Association for Computing Machinery, 2023, p. 120–133.
  * [3] F. C. Weigend, S. Sonawani, M. Drolet, and H. B. Amor, “Anytime, anywhere: Human arm pose from smartwatch data for ubiquitous robot control and teleoperation,” in _2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2023, pp. 3811–3818.
  * [4] F. C. Weigend, X. Liu, S. Sonawani, N. Kumar, V. Vasudevan, and H. B. Amor, “iRoCo: Intuitive robot control from anywhere using a smartwatch,” 2024.
  * [5] X. Liu, G. Clark, J. Campbell, Y. Zhou, and H. B. Amor, “Enhancing state estimation in robots: A data-driven approach with differentiable ensemble kalman filters,” in _2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2023, pp. 1947–1954.
  * [6] C. Yang, J. Chen, and F. Chen, “Neural learning enhanced teleoperation control of baxter robot using imu based motion capture,” in _2016 22nd International Conference on Automation and Computing (ICAC)_. IEEE, 2016, pp. 389–394.
  * [7] V. Joukov, J. Ćesić, K. Westermann, I. Marković, D. Kulić, and I. Petrović, “Human motion estimation on lie groups using imu measurements,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 1965–1972.


# Computing an Optimal Control Policy for an Energy Storage

Pierre Haessig∗†, Thibaut Kovaltchouk†, Bernard Multon†, Hamid Ben Ahmed†,
Stéphane Lascaud‡ * Corresponding author<EMAIL_ADDRESS>SATIE
CNRS laboratory - ENS Rennes, Bruz, France‡ LME department - EDF R&D,
Écuelles, FranceCopyright © 2014 Pierre Haessig et al. This is an open-access
article distributed under the terms of the Creative Commons Attribution
License, which permits unrestricted use, distribution, and reproduction in any
medium, provided the original author and source are credited.
http://creativecommons.org/licenses/by/3.0/

###### Abstract

We introduce StoDynProg, a small library created to solve Optimal Control
problems arising in the management of Renewable Power Sources, in particular
when coupled with an Energy Storage System. The library implements generic
Stochastic Dynamic Programming (SDP) numerical methods which can solve a large
class of Dynamic Optimization problems.

We demonstrate the library capabilities with a prototype problem: smoothing
the power of an Ocean Wave Energy Converter. First we use time series analysis
to derive a stochastic Markovian model of this system since it is required by
Dynamic Programming. Then, we briefly describe the “policy iteration”
algorithm we have implemented and the numerical tools being used. We show how
the API design of the library is generic enough to address Dynamic
Optimization problems outside the field of Energy Management. Finally, we
solve the power smoothing problem and compare the optimal control with a
simpler heuristic control.

###### Index Terms:

Stochastic Dynamic Programming, Policy Iteration Algorithm, Autoregressive
Models, Ocean Wave Energy, Power Smoothing.

##  1 Introduction to Power Production Smoothing

Electric power generated by renewable sources like wind, sun or ocean waves
can exhibit a strong _variability_ along time. Because on an electricity grid
the energy production must match the consumption, this variability can be an
issue for the grid stability. Yet most of the time, fluctuations of renewable
power sources are absorbed without trouble thanks to regulation mechanisms
which make flexible generation units adjust their production in real-time.
Therefore, the production-consumption equilibrium can be maintained.

However, there are cases where fluctuations may be considered too strong to be
fed directly to the grid so that an _energy storage system_ , acting as a
_buffer_ , may be required to smooth out the production. The schematic of the
system considered in this article is given on figure 1.

Figure 1: Power smoothing with an Energy Storage: an example of an Optimal
Control problem.

###  1.1 Smoothing with an Energy Storage

Electricity generation from _ocean waves_ (with machines called Wave Energy
Converters) is an example where the output power can be _strongly
fluctuating_. This is illustrated on figure 2 where the output power
$P_{prod}(t)$ from a particular wave energy converter called SEAREV is
represented during 100 seconds.

We just mention that this production time series does not come from
measurements but from an hydro-mechanical simulation from colleagues since the
SEAREV is a big 1 MW - 30 meters long machine which is yet to be built
[Ruellan-2010].

The oscillations of $P_{prod}(t)$ at a period of about 1.5 s comes from the
construction of the SEAREV: in short, it is a floating _double-pendulum_ that
oscillates with the waves. Also, because _ocean waves have a stochastic
behavior_ , the amplitude of these oscillations is irregular.

Figure 2: Smoothing the Ocean Power injected to the grid using an Energy
Storage controlled by the simple linear law (3). The storage buffers the
difference between the two powers.

Therefore an energy storage absorbing a power $P_{sto}$ can be used to smooth
out the power $P_{grid}$ injected to the electricity network:

$P_{grid}(t)=P_{prod}(t)-P_{sto}(t)$ (1)

The energy of the storage then evolves as:

$E_{sto}(k+1)=E_{sto}(k)+P_{sto}(k)\Delta t$ (2)

expressed here in discrete time ($\Delta t=0.1\text{ s}$ throughout this
article), without accounting for losses. The storage energy is bounded: $0\leq
E_{sto}\leq E_{rated}$, where $E_{rated}$ denotes the storage capacity which
is set to 10 MJ in this article (i.e. about 10 seconds of reserve at full
power)

It is a control problem to choose a power smoothing law. We present the
example of a linear feedback control:

$P_{grid}(t)=\frac{P_{max}}{E_{rated}}E_{sto}(t)$ (3)

where $P_{max}$ is the rated power of the SEAREV (1.1 MW). This law gives
“good enough” smoothing results as it can be seen on figure 2.

The performance of the smoothing is greatly influenced by the _storage sizing_
(i.e. the choice of the capacity $E_{rated}$). This question is not addressed
in this article but was discussed by colleagues [Aubry-2010]. We also don’t
discuss the choice of the storage _technology_ , but it is believed that
super-capacitors would be the most suitable choice. Because energy storage is
very expensive (~20 k€/kWh or ~5 k€/MJ for supercaps), there is an interest in
studying how to make the best use of a given capacity to avoid a costly over-
sizing.

###  1.2 Finding an Optimal Smoothing Policy

Control law (3) is an example of heuristic choice of policy and we now try to
go further by finding an _optimal_ policy.

Optimality will be measured against a _cost function_ $J$ that penalizes the
average variability of the power injected to the grid:

$J=\frac{1}{N}\mathbb{E}\left\\{\sum_{k=0}^{N-1}c(P_{grid}(k))\right\\}\quad\text{with
$N\rightarrow\infty$}$ (4)

where $c$ is the _instantaneous cost_ (or penalty) function which can be
$c(P_{grid})=P_{grid}^{2}$ for example. Expectation $\mathbb{E}$ is needed
because the production $P_{prod}$ is a stochastic input, so that the output
power $P_{grid}$ is also a random variable.

This minimization problem falls in the class of _stochastic dynamic
optimization_. It is _dynamic_ because decisions at each time-step cannot be
taken independently due to the coupling along time introduced by the evolution
of the stored energy (2). To describe the dynamics of the system, we use the
generic notation

$x_{k+1}=f(x_{k},u_{k},\varepsilon_{k})$ (5)

where $x,u,\varepsilon$ are respectively _state_ variables, _control_
variables and _perturbations_. State variables are the “memory” of the system.
The stored energy $E_{sto}$ is here the only state variable, but more will
appear in section 2.2. Control variables are the ones which values must be
chosen at each instant to optimize the cost $J$. The injected power $P_{grid}$
is here the single control variable.

Dynamic optimization (also called _optimal control_) is addressed by the
Dynamic Programming method [Bertsekas-2005] which yields a theoretical
analysis of the solution structure. Indeed, once all state variables (i.e.
“memories”) of the system are identified, the optimum of the cost $J$ is
attained by a “state feedback” policy, that is a policy where the control is
chosen as _a function of the state:_

$P_{grid}(t)=\mu(x(t))$ (6)

The goal is then to find the _optimal_ feedback function $\mu$. Since
$E_{sto}$ is a state variable, policy (3) is in fact a special case of (6).
Since $\mu$ has no special structure in the general case1, it will be
_numerically computed on a grid_ over the state space. We cover the algorithm
for this computation in section 3.

####  1.2.1 Prerequisite

Dynamic Programming does require that stochastic perturbations are
_independent_ random variables (i.e. the overall dynamical model must be
Markovian) and this is not true for the $P_{prod}(k)$ time series. Therefore
we devote section 2 to the problem of expressing $P_{prod}$ as a discrete-time
Markov process, using _time series analysis_. This will yield new state
variables accounting for the dynamics of $P_{prod}$.

##  2 Stochastic Model of a Wave Energy Production

We now take a closer look at the $P_{prod}$ time series. A 1000 s long
simulation is presented on figure 4, along with a zoom to better see the
structure at short time scales. An histogram is also provided which shows that
$P_{prod}$ is clearly _non-gaussian_. This precludes the direct use of
“standard” time series models based on Autoregressive Moving Average (ARMA)
models [Brockwell-1991].

However, we can leverage the knowledge of the inner working of the SEAREV.
Indeed, by calling $\Omega$ the rotational speed of the inner pendulum with
respect to the hull, we know that the output power is:

$P_{prod}=T_{PTO}(\Omega).\Omega$ (7)

where $T_{PTO}$ is the torque applied to the pendulum by the electric machine
which harvests the energy (PTO stands for “Power Take Off”). Finding the best
$T_{PTO}$ command is actually another optimal control problem which is still
an active area of research in the Wave Energy Conversion community
[Kovaltchouk-2013]. We use here a “viscous damping law, with power leveling”,
that is $T_{PTO}(\Omega)=\beta.\Omega$. This law is applied as long as it
yields a power below $P_{max}$. Otherwise the torque is reduced to level the
power at 1.1 MW as can be seen on figure 4 whenever the speed is more than 0.5
rad/s.

Thanks to equation (7), we can thus model the speed $\Omega$ and then deduce
$P_{prod}$. Modeling the speed is much easier because it is quite Gaussian
(see fig. 4) and has a much more regular behavior which can be captured by an
ARMA process.

###  2.1 Autoregressive Model of the Speed

Within the ARMA family, we restrict ourselves to the autoregressive (AR)
processes because we need a Markovian model. The equation of an AR(p) model
for the speed is:

$\Omega(k)=\phi_{1}\Omega(k-1)+\dots+\phi_{p}\Omega(k-p)+\varepsilon(k)$ (8)

where $p$ is the order of the model and $\varepsilon(k)$ is a series
independent random variables. Equation (8) indeed yields a Markovian process,
using the lagged observations of the speed $\Omega(k-1),\dots,\Omega(k-p)$ as
state variables.

AR(p) model fitting consists in _selecting_ the order $p$ and _estimating_ the
unknown coefficients $\phi_{1},\dots,\phi_{p}$ as well as the unknown variance
of $\varepsilon$ which we denote $\sigma_{\varepsilon}^{2}$.

Figure 3: Autocorrelation function (acf) of the speed data, compared with the
acf from two AR(2) models, fitted with two different methods.

Figure 4: Speed & Power time series from a 1000 seconds SEAREV simulation
(sample Em_1.txt). The gray rectangle time interval is enlarged in the middle
panel. Distribution histogram on the right.

####  2.1.1 Order selection

is generally done using _information criterions_ such as AIC or BIC
[Brockwell-1991], but for this modeling problem, we restrict ourselves to the
smallest order which can reproduce the _decaying oscillations_ of the
autocorrelation function. Autocorrelation (acf) of the speed is plotted on
figure 3 where we can see that a model of order $p=2$ can indeed reproduce the
autocorrelation up to about 15 s of time lags (a 1st order model would only
yield an exponential decay without oscillations). 15 s is thought to be the
time horizon of interest when using a 10 MJ/1.1 MW energy storage.

Keeping the model order low is required to maintain the dimension of the
overall state vector under 3 or 4. The underlying issue of an exponentially
growing complexity will appear in section 3 when solving the Dynamic
Programming equation.

####  2.1.2 Parameter estimation

Once the order is selected, we have to estimate coefficients $\phi_{1}$,
$\phi_{2}$ and $\sigma_{\varepsilon}^{2}$. “Classical” fitting methodology
[Brockwell-1991] is based on Conditional Maximum Likelihood Estimators (CMLE).
This method is readily available in GNU R with the arima routine or in Python
with statsmodel.tsa.ar_model.

However, we have plotted the autocorrelation of the estimated AR(2) model on
figure 3 to show that CMLE is _not appropriate:_ oscillations of the acf
clearly decay too slowly compared to the data acf.

The poor adequacy of this fit is actually a consequence of our choice of a low
order model which implies that the AR(2) process can only be an _approximation
of the true process_. Statistically speaking, our model is _misspecified_ ,
whereas CMLE is efficient for correctly specified models only. This problem
has been discussed in the literature [McElroy-2013] and has yielded the
“Multi-step ahead fitting procedure”.

Being unfamiliar with the latter approach, we compute instead
$\phi_{1},\phi_{2}$ estimates which _minimize the difference_ between the
theoretical AR(2) acf and the data acf. The minimization criterion is the sum
of the squared acf differences over a range of lag times which can be chosen.
We name this approach the “multi-lags acf fitting” method. Minimization is
conducted with fmin from scipy.optimize.

method |  $\hat{\phi}_{1}$ |  $\hat{\phi}_{2}$ |  $\hat{\sigma}_{\varepsilon}$
---|---|---|---
CMLE |  1.9883 (.0007) |  -0.9975 (.0007) |  0.00172
fit on 15 s |  1.9799 |  -0.9879 |  0.00347
TABLE II: AR(2) fitting results from the two methods (along with standard
error when available).

The result of this acf fitting over time lags up to 15 s (i.e. 150 lags) is
shown on figure 3 while numerical estimation results are given in table I.

With the model obtained from this multi-lags method, we can simulate speed and
power trajectories and check that they have a “realistic behavior”. We can
thus infer that the dynamic optimization algorithm should make appropriate
control decisions out of it. This will be discussed in section 3.3. Going
further, it would be interesting to study the influence of the AR parameters
(including order $p$) on the dynamic optimization to see if the “multi-lags
acf fitting” indeed brings an improvement of the final cost function $J$.

###  2.2 Reformulation as a state-space model

The AR(2) model is a state-space model with state variables being the lagged
observations of the speed $\Omega(k-1)$ and $\Omega(k-2)$. In order to get a
model with a better “physical interpretation” we introduce the variable
$A_{k}=(\Omega_{k}-\Omega_{k-1})/\Delta t$ which is the backward discrete
derivative of $\Omega$. As the timestep gets smaller $A_{k}$ comes close to
the acceleration (in rad/s2) of the pendulum. Using $(\Omega,A)$ as the state
vector, we obtain the following state-space model:

$\begin{split}\begin{pmatrix}\Omega_{k}\\\
A_{k}\end{pmatrix}=&\begin{bmatrix}\phi_{1}+\phi_{2}&-\phi_{2}\Delta t\\\
(\phi_{1}+\phi_{2}-1)/\Delta
t&-\phi_{2}\end{bmatrix}\begin{pmatrix}\Omega_{k-1}\\\ A_{k-1}\end{pmatrix}\\\
+&\begin{bmatrix}1\\\ 1/\Delta t\end{bmatrix}\varepsilon_{k}\end{split}$ (9)

We now have a stochastic Markovian model for the power production of the
SEAREV. Taken together with state equation of the storage (2) and algebraic
relations (1) and (7), we have a Markovian model of the overall system. The
state vector $x=(E_{sto},\Omega,A)$ is of dimension 3 which is just small
enough to apply the Stochastic Dynamic Programming method.

##  3 Optimal storage control with Dynamic Programming

###  3.1 The Policy Iteration Algorithm

We now give an overview of the _policy iteration_ algorithm that we
implemented to solve the power smoothing problem described in the
introduction. Among the different types of dynamic optimization problems, it
is an “infinite horizon, average cost per stage problem” (as seen in (4)).
While at first this cost equation involves a summation over an infinite number
of instants, the Dynamic Programming approach cuts this into two terms: the
present and the whole future. In the end, the optimization falls back to
solving this equilibrium equation:

$\begin{split}J+\tilde{J}(x)=\min_{u\in
U(x)}\underset{w}{\mathbb{E}}\Big{\\{}\underbrace{c(x,u,w)}_{\text{instant
cost}}\\\ +\underbrace{\tilde{J}(f(x,u,w))}_{\text{cost of the
future}}\Big{\\}}\end{split}$ (10)

where $J$ is the minimized average cost and $\tilde{J}$ is the transient (or
differential) cost function, also called _value function_.

Note that eq. (10) is a functional equation for $\tilde{J}$ which should be
solved for _any value_ of the state $x$ in the state space. In practice, it is
solved in a _discrete grid_ that must be chosen so that the variations of
$\tilde{J}$ are represented with enough accuracy. Also, the optimal policy
$\mu$ appears implicitly as the _argmin_ of this equation, that is the optimal
control $u$ for each $x$ value of the state grid.

####  3.1.1 Equation solving

The simplest way to solve eq. (10) is to iterate the right-hand side, starting
with a zero value function. This is called _value iteration_.

A more efficient approach is _policy iteration_. It starts with an initial
policy (like the heuristic linear (3)) and gradually improves it with a two
steps procedure:

  * 1.

policy evaluation: the current policy is evaluated, which includes computing
the average cost (4) and the so-called _value function_

  * 2.

policy improvement: a single step of optimization with policy iteration yields
a improved policy. Then this policy must be again evaluated (step 1).

The policy evaluation involves solving the equilibrium equation without the
minimization step but with a fixed policy $\mu$ instead:

$\begin{split}J_{\mu}+\tilde{J}_{\mu}(x)=\underset{w}{\mathbb{E}}\Big{\\{}\underbrace{c(x,\mu(x),w)}_{\text{instant
cost}}\\\ +\underbrace{\tilde{J}_{\mu}(f(x,\mu(x),w))}_{\text{cost of the
future}}\Big{\\}}\end{split}$

It can be solved by iterating the right-hand side like for policy iteration
but much faster due to the absence of minimization. In the end, a few policy
improvement iterations are needed to reach convergence. More details about the
value and policy iteration algorithms can be found in [Bertsekas-2005]
textbook for example. The conditions for the convergence, omitted here, are
also discussed.

###  3.2 StoDynProg library description

We have created a small library to _describe_ and _solve_ optimal control
problems (in discrete time) using the Stochastic Dynamic Programming (SDP)
method. It implements the value iteration and policy iteration algorithms
introduced above. Source code is available on GitHub
https://github.com/pierre-haessig/stodynprog under a BSD 2-Clause license.

####  3.2.1 Rationale for a library, benefits of Python

Because the SDP algorithms are in fact quite simple (they can be written with
one set of nested for loops) we were once told that they should be written
from scratch for each new problem. However we face in our research in energy
management several optimization problems with slight structural differences so
that code duplication would be unacceptably high. Thus the motivation to write
a unified code that can handle all our use cases, and hopefully some others’.

StoDynProg is pure Python code built with numpy for multi-dimensional array
computations. We also notably use an external multidimensional interpolation
routine by Pablo Winant (see 3.2.5 below).

The key aspect of the flexibility of the code is its ability to handle
problems of _arbitrary dimensions_ (in particular the state space and the
control space). This impacts particularly the way to iterate over those
variables. Our code makes thus a heavy use of Python tuple packing/unpacking
machinery and itertools.product to iterate on rectangular grids of arbitrary
dimension.

####  3.2.2 API description

StoDynProg provides two main classes: SysDescription and DPSolver.

####  3.2.3 SysDescription

holds the description of the discrete-time dynamic optimization problem.
Typically, a user writes its dynamics function (the Python implementation of
$f$ in (5)) and handles it to a SysDescription instance:

from stodynprog import SysDescription# SysDescription object with proper
dimensions# of state (2), control (1) and perturbation (1)mysys =
SysDescription((2, 1, 1))def my_dyn(x1, x2, u, w): ’dummy dynamics’ x1_next =
x1 + u + w x2_next = x2 + x1 return (x1_next, x2_next)# assign the dynamics
function:mysys.dyn = my_dynWe use here a setter/getter approach for the dyn
property. The same is used to describe the cost function ($c$ in (4)). We
believe the property approach enables simplified user code compared to a class
inheritance mechanism. With some inspiration of Enthought traits, the setter
has a basic validation mechanism that checks the signature of the function
being assigned (with getargspec from the inspect module).

####  3.2.4 DPSolver

holds parameters that tunes the optimization process, in particular the
discretized grid of the state. Also, it holds the code of the optimization
algorithm in its methods. We illustrate here the creation of the solver
instance attached to the previous system:

from stodynprog import DPSolver# Create the solver for ‘mysys‘ system:dpsolv =
DPSolver(mysys)# state discretizationx1_min, x1_max, N1 = (0, 2.5, 100)x2_min,
x2_max, N2 = (-15, 15, 100)x_grid = dpsolv.discretize_state(x1_min, x1_max,
N1, x2_min, x1_max, N2)Once the problem is fully described, the optimization
can be launched by calling dpsolv.policy_iteration with proper arguments about
the number of iterations.

For more details on StoDynProg API usage, an example problem of _Inventory
Control_ is treated step-by-step in the documentation (created with Sphinx).

####  3.2.5 Multidimensional Interpolation Routine

StoDynProg makes an intensive use of a multidimensional interpolation routine
that is not available in the “standard scientific Python stack”. Interpolation
is needed because the algorithm manipulates two scalar functions which are
discretized on a grid over the state space: the value function $\tilde{J}$ and
feedback policy $\mu$. Thus, functions are stored as $n$-d arrays, where $n$
is the dimension of the state vector ($n=3$ for ocean power smoothing
example). In the course of the algorithm, the value function needs to be
evaluated between grid points, thus the need for interpolation.

####  3.2.6 Requirements and Algorithm Selection

No “fancy” interpolation method is required so linear interpolation is a good
candidate. Speed is very important because it is called many times. Also, it
should accept vectorized inputs, so that interpolation of multiple points can
be done efficiently in one call. We assert that the functions will be stored
on a _rectangular grid_ which should simplify interpolation computations. The
most stringent requirement is _multidimensionality_ (for $0\leq n\leq 4$)
which rules out most available tools.

We have evaluated 4 routines (details available in an IPython Notebook within
StoDynProg source tree):

  * •

LinearNDInterpolator class from scipy.interpolate

  * •

RectBivariateSpline class from scipy.interpolate

  * •

map_coordinates routine from scipy.ndimage

  * •

and MultilinearInterpolator class written by Pablo Winant within its Dolo
project [Winant-2010] for Economic modelling (available on
https://github.com/albop/dolo).

The most interesting in terms of performance and off-the-shelf availability is
RectBivariateSpline which exactly meets our needs expect for
multidimensionality because it’s limited to $n=2$. LinearNDInterpolator has no
dimensionality limitations but works with unstructured data and so does not
leverage the rectangular structure. Interpolation time was found 4 times
longer in 2D, and unacceptably long in 3D. Then map_coordinates and
MultilinearInterpolator were found to both satisfy all our criterions but the
latter being consistently 4 times faster (both 2D and 3D case). Finally we
also selected MultilinearInterpolator because it can be instantiated to retain
the data once and then called several time. We find the usage of this object-
oriented interface more convenient than functional interface of
map_coordinates.

###  3.3 Results for Searev power smoothing

Figure 5: Storage control policy: Power injected to the grid as a function of
speed and acceleration, for 7 levels of stored energy between empty and full.

Figure 6: Comparison of the power smoothing behavior between the _heuristic_
(dark blue) and _optimized_ (light blue) storage management policies (storage
capacity of 10 MJ). Stored energy on the bottom panel.

We have applied the policy iteration algorithm to the SEAREV power smoothing
problem introduced in section 1. The algorithm is initialized with the linear
storage control policy (3). This heuristic choice is then gradually improved
by each policy improvement step.

####  3.3.1 Algorithm parameters

About 5 policy iterations only are needed to converge to an optimal strategy.
In each policy iteration, there is a policy evaluation step which requires
1000 iterations to converge. This latter number is dictated by the time
constant of the system (1000 steps $\leftrightarrow$ 100 seconds) and 100
seconds is the time it takes for the system to “decorrelate”, that is loose
memory of its state (both speed and stored energy).

We also need to decide how to discretize and bound the state space of the
{SEAREV + storage} system:

  * •

for the stored energy $E_{sto}$, bounds are the natural limits of the storage:
$E_{sto}\in[0,E_{rated}]$. A grid of 30 points yields precise enough results.

  * •

for the speed $\Omega$ and the acceleration $A$, there are no natural bounds
so we have chosen to limit the values to $\pm 4$ _standard deviations_. This
seems wide enough to include most observations but not too wide to keep a good
enough resolution. We use grids of 60 points to keep the grid step small
enough.

This makes a state space grid of $30\times 60\times 60\approx 110\text{k}$
points. Although this number of points can be handled well by a present
desktop computer, this simple grid size computation illustrates the commonly
known weakness of Dynamic Programming which is the “Curse of Dimensionality”.
Indeed, this size grows exponentially with the number of dimensions of the
state so that for practical purpose state dimension is limited to 3 or 4. This
explains the motivation to search a low order model for the power production
time series in section 2.

####  3.3.2 Algorithm execution time

With the aforementioned discretization parameters, policy evaluation takes
about 10 s (for the 1000 iterations) while policy improvement takes 20 s (for
one single value iteration step). This makes 30 s in total for one policy
iteration step, which is repeated 5 times. Therefore, the optimization
converges in about 3 minutes. This duration would grow steeply should the grid
be refined.

As a comparison of algorithm efficiency, the use of _value iteration_ would
takes much longer than _policy iteration_. Indeed, it needs 1000 iterations,
just like policy evaluation (since it is dictated by the system’s
“decorrelation time”) but each iteration involves a costly optimization of the
policy so that it takes 20 s. This makes altogether 5 hours of execution time,
i.e. 100 times more than policy iteration!

As possible paths to improve the execution time, we see, at the _code level_ ,
the use of more/different vectorization patterns although vectorized
computation is already used a lot. Maybe the use of Cython may speed up
unavoidable loops but this may not be worth the loss of flexibility and the
decrease in coding speed. Optimization at the _algorithm level_ , just as
demonstrated with “policy vs. value iteration”, is also worth investigating
further. In the end, more use of Robert Kern’s line_profiler will be needed to
decide the next step.

####  3.3.3 Output of the computation

The policy iteration algorithm solves equation (10) and outputs the minimized
cost $J$ and two arrays: function $\tilde{J}$ (transient cost) and function
$\mu$ (optimal policy (6)), both expressed on the discrete state grid (3d
grid).

We focus on $\mu$ which yields the power $P_{grid}$ that should be injected to
the grid for any state of the system. Figure 5 is a Mayavi surface plot which
shows $P_{grid}(\Omega,A)$ for various levels of $E_{sto}$. Observations of
the result are in agreement with what can be expected from a reasonable
storage control:

  * •

the more energy there is in the storage, the more power should be injected to
the grid (similar to the heuristic control (3)).

  * •

the speed and acceleration of the SEAREV also modulates the injected power,
but to a lesser extent. We may view speed and acceleration as approximate
measurements of the _mechanical energy_ of the SEAREV. This energy could be a
hidden influential state variable, in parallel with the stored energy.

  * •

the injected power is often set between 0.2 and 0.3 MW, that is _close to the
average_ power production.

Such observations show that the algorithm has _learned_ from the SEAREV
behavior to take sharper decisions compared to the heuristic policy it was
initialized with.

####  3.3.4 Qualitative analysis of the trajectory

To evaluate the storage control policy, we simulate its effect on the sample
SEAREV data we have (instead of using the state space model used for the
optimization). The only adaptation required for this trajectory simulation is
to transform the _policy array_ ($\mu$ known on the state grid only) into a
_policy function_ ($\mu$ evaluable on the whole state space). This is achieved
using the same n-dimensional interpolation routine used in the algorithm.

A simulated trajectory is provided on figure 6 to compare the effect of the
optimized policy with the heuristic linear policy (3). As previously said, the
storage capacity is fixed at $E_{rated}=10\text{ MJ}$ or about 9 seconds of
charge/discharge at the rated power.

Positive aspect, the optimized policy yields an output power that is generally
closer to the average (thin gray line) than the linear policy. This better
smoothing of the “peaks and valleys” of the production is achieved by a better
usage of the available storage capacity. Indeed, the linear policy generally
under-uses the higher levels of energy.

As a slight negative aspect, the optimized policy yields a “spiky” output
power in the situations of high production (200–220 s). In this situation, the
output seems worse that the linear policy. We connect this underperformance to
the linear model (9) used to represent the SEAREV dynamics. The linearity
holds well for small movements but not when the speed is high and the pendulum
motion gets very abrupt (acceleration high above 4 standard deviations which
contradicts the Gaussian distribution assumption). Since the control
optimization is based on the linear model, the resulting control law cannot
appropriately manage these non-linear situations. Only an upgraded model would
genuinely solve this problem but we don’t have yet an appropriate low-order
non-linear model of the SEAREV. One quick workaround to reduce the power peaks
is to shave the acceleration measurements (not demonstrated here).

####  3.3.5 Quantitative assessment

We now numerically check that the optimized policy brings a true enhancement
over the linear policy. We simulate the storage with the three 1000 s long
samples we have and compute the power variability criterion2 for each.

Figure 7 shows the standard deviation for each sample in three situations:
without storage (which yields the natural standard deviation of the SEAREV
production), with a storage controlled by the linear policy and finally the
same storage controlled by the optimized policy. Sample Em_1.txt was used to
fit the state space model (9) but we don’t think this should introduce a too
big bias because of the low model order.

Beyond the intersample variability, we can see the consistent improvement
brought by the optimized law. Compared to the linear policy, the standard
deviation of the injected power is reduced by about 20 % (27 %, 16 %, 22 % for
each sample respectively). We can conclude that the variability of the
injected power is indeed reduced by using the Dynamic Programming.

Figure 7: Effect of optimizing the storage control on three SEAREV production
time series. Standard deviation compared to the heuristic linear control case
is reduced by about 20 %.

Because the RMS deviation criterion used in this article is not directly
limited or penalized in current grid codes, there is no financial criterion to
decide whether the observed deviations are acceptable or not. Therefore we
cannot conclude if the ~20% reduction of the variability brought by optimal
control is valuable. Nevertheless, there exists criterions like the "flicker"
which are used in grid codes to set standards of power quality. Flicker, which
is way more complicated to than an additive criterion like (4) could be used
to put an economic value on a control strategy. This is the subject of ongoing
research.

##  4 Conclusion

With the use of standard Python modules for scientific computing, we have
created StoDynProg, a small library to solve Dynamic Optimization problems
using Stochastic Dynamic Programming.

We have described the mathematical and coding steps required to apply the SDP
method on an example problem of realistic complexity: smoothing the output
power of the SEAREV Wave Energy Converter. With its generic interface,
StoDynProg should be applicable to other Optimal Control problems arising in
Electrical Engineering, Mechanical Engineering or even Life Sciences. The only
requirement is an appropriate mathematical structure (Markovian model), with
the “Curse of Dimensionality” requiring a state space of low dimension.

Further improvements on this library should include a better source tree
organization (make a proper package) and an improved test coverage.
\raisebox{10.00002pt}{\hypertarget{id13}{}}\hyperlink{id4}{1}\raisebox{10.00002pt}{\hypertarget{id13}{}}\hyperlink{id4}{1}footnotetext:
In the special case of a linear dynamics and a quadratic cost (“LQ control” ),
the optimal feedback is actually a _linear_ function. Because of the state
constraint $0\leq E_{sto}\leq E_{rated}$, the storage control problem falls
outside this classical case.
\raisebox{10.00002pt}{\hypertarget{id14}{}}\hyperlink{id12}{2}\raisebox{10.00002pt}{\hypertarget{id14}{}}\hyperlink{id12}{2}footnotetext:
Instead of using the exact optimization cost (4) (average quadratic power in
MW2), we actually compute the standard deviation (in MW). It is mathematically
related to the quadratic power and we find it more readable.

## References

  * [Aubry-2010] J. Aubry, P. Bydlowski, B. Multon, H. Ben Ahmed, and B. Borgarino. _Energy Storage System Sizing for Smoothing Power Generation of Direct Wave Energy Converters_ , 3rd International Conference on Ocean Energy, 2010.
  * [Bertsekas-2005] D. P. Bertsekas, _Dynamic Programming and Optimal Control_ , Athena Scientific, 2005.
  * [Brockwell-1991] P. J. Brockwell, and R. A. Davis. _Time Series: Theory and Methods_ , Springer Series in Statistics, Springer, 1991.
  * [Kovaltchouk-2013] T. Kovaltchouk, B. Multon, H. Ben Ahmed, F. Rongère, J. Aubry, and A. Glumineau. _Influence of control strategy on the global efficiency of a Direct Wave Energy Converter with electric Power Take-Off_ , EVER 2013 conference, 2013.
  * [McElroy-2013] T. McElroy, and M. Wildi. _Multi-step-ahead estimation of time series models_ , International Journal of Forecasting, 29: 378–394, 2013.
  * [Ruellan-2010] M. Ruellan, H. Ben Ahmed, B. Multon, C. Josset, A. Babarit, and A. Clément. _Design Methodology for a SEAREV Wave Energy Converter_ , IEEE Trans. Energy Convers, 25: 760–767, 2010.
  * [Winant-2010] P. Winant. _Dolo, a python library to solve global economic models_ , http://albop.github.io/dolo, 2010.