# A Generalization of Varnavides’s Theorem

Asaf Shapira School of Mathematics, Tel Aviv University, Tel Aviv 69978,
Israel. Email<EMAIL_ADDRESS>Supported in part by ERC Consolidator Grant
863438 and NSF-BSF Grant 20196.

###### Abstract

A linear equation $E$ is said to be sparse if there is $c>0$ so that every
subset of $[n]$ of size $n^{1-c}$ contains a solution of $E$ in distinct
integers. The problem of characterizing the sparse equations, first raised by
Ruzsa in the 90’s, is one of the most important open problems in additive
combinatorics. We say that $E$ in $k$ variables is abundant if every subset of
$[n]$ of size $\varepsilon n$ contains at least $\text{poly}(\varepsilon)\cdot
n^{k-1}$ solutions of $E$. It is clear that every abundant $E$ is sparse, and
Girão, Hurley, Illingworth and Michel asked if the converse implication also
holds. In this note we show that this is the case for every $E$ in $4$
variables. We further discuss a generalization of this problem which applies
to all linear equations.

##  1 Introduction

Turán-type questions are some of the most well studied problems in
combinatorics. They typically ask how “dense” should an object be in order to
guarantee that it contains a certain small substructure. In the setting of
graphs, this question asks how many edges an $n$-vertex graph should contain
in order to force the appearance of some fixed graph $H$. For example, a
central open problem in this area asks, given a bipartite graph $H$, to
determine the smallest $T=T_{H}(\varepsilon)$ so that for every $n\geq T$
every $n$-vertex graph with $\varepsilon{n\choose 2}$ edges contains a copy of
$H$ (see [3] for recent progress). A closely related question which also
attracted a lot of attention, is the supersaturation problem, introduced by
Erdős and Simonovits [5] in the 80’s. In the setting of Turán’s problem for
bipartite $H$, the supersaturation question asks to determine the largest
$T^{*}_{H}(\varepsilon)$ so that every $n$-vertex graph with
$\varepsilon{n\choose 2}$ edges contains at least
$(T^{*}_{H}(\varepsilon)-o_{n}(1))\cdot n^{h}$ labelled copies of $H$, where
$h=|V(H)|$ and $o_{n}(1)$ denotes a quantity tending to $0$ as
$n\rightarrow\infty$. One of the central conjectures in this area, due to
Sidorenko, suggests that $T^{*}_{H}(\varepsilon)=\varepsilon^{m}$, where
$m=|E(H)|$ (see [4] for recent progress).

We now describe two problems in additive number theory, which are analogous to
the graph problems described above. We say that a homogenous linear equation
$\sum^{k}_{i=1}a_{i}x_{i}=0$ is invariant if $\sum_{i}a_{i}=0$. All equations
we consider here will be invariant and homogenous. Given a fixed linear
equation $E$, the Turán problem for $E$ asks to determine the smallest
$R=R_{E}(\varepsilon)$ so that for every $n\geq R$, every
$S\subseteq[n]:=\\{1,\ldots,n\\}$ of size $\varepsilon n$ contains a solution
to $E$ in distinct integers. For example, when $E$ is the equation $a+b=2c$ we
get the Erdős–Turán–Roth problem on sets avoiding $3$-term arithmetic
progressions (see [7] for recent progress). Continuing the analogy with the
previous paragraph, we can now ask to determine the largest
$R^{*}_{E}(\varepsilon)$ so that every $S\subseteq[n]$ of size $\varepsilon n$
contains at least $(R^{*}_{E}(\varepsilon)-o_{n}(1))\cdot n^{k-1}$ solutions
to $E$, where $k$ is the number of variables in $E$. We now turn to discuss
two aspects which make the arithmetic problems more challenging than the graph
problems.

Let us say that $E$ is sparse if there is $C=C(H)$ so that111It is easy to see
that this definition is equivalent to the one we used in the abstract.
$R_{E}(\varepsilon)\leq\varepsilon^{-C}$. The first aspect which makes the
arithmetic landscape more varied is that while in the case of graphs it is
well known (and easy) that for every bipartite $H$ we have
$T_{H}(\varepsilon)=\text{poly}(1/\varepsilon)$, this is no longer the case in
the arithmetic setting. Indeed, while Sidon’s equation $a+b=c+d$ is sparse, a
well known construction of Behrend [1] shows that $a+b=2c$ is not
sparse222More precisely, it shows that in this case
$R_{E}(\varepsilon)\geq(1/\varepsilon)^{c\log 1/\varepsilon}$. Here and
throughout this note, all logarithms are base $2$.. The problem of determining
which equations $E$ are sparse is a wide open problem due to Ruzsa, see
Section $9$ in [9].

Our main goal in this paper is to study another aspect which differentiates
the arithmetic and graph theoretic problems. While it is easy to translate a
bound for $T_{H}(\varepsilon)$ into a bound for $T^{*}_{H}(\varepsilon)$ (in
particular, establishing that
$T^{*}_{H}(\varepsilon)\geq\text{poly}(\varepsilon)$ for all bipartite $H$),
it is not clear if one can analogously transform a bound for
$R_{E}(\varepsilon)$ into a bound for $R^{*}_{E}(\varepsilon)$. The first
reason is that while we can average over all subsets of vertices of graphs, we
can only average over “structured” subsets of $[n]$. This makes is hard to
establish a black-box reduction/transformation between $R_{E}(\varepsilon)$
and $R^{*}_{E}(\varepsilon)$. The second complication is that, as we mentioned
above, we do not know which equations are sparse. This makes it hard to
directly relate these two quantities. Following [6], we say that $E$ is
abundant if $R^{*}_{E}(\varepsilon)\geq\varepsilon^{C}$ for some $C=C(E)$.
Clearly, if $E$ is abundant then it is also sparse. Girão, Hurley, Illingworth
and Michel [6] asked if the converse also holds, that is, if one can transform
a polynomial bound for $R_{E}(\varepsilon)$ into a polynomial bound for
$R^{*}_{E}(\varepsilon)$. Our aim in this note is to prove the following.

######  Theorem 1.1.

If an invariant equation $E$ in $4$ variables is sparse, then it is also
abundant. More precisely, if $R_{E}(\varepsilon)\leq\varepsilon^{-C}$ then
$R^{*}_{E}(\varepsilon)\geq\frac{1}{2}\varepsilon^{8C}$ for all small enough
$\varepsilon$.

Given the above discussion it is natural to extend the problem raised in [6]
to all equations $E$.

######  Problem 1.2.

Is it true that for every invariant equation $E$ there is $c=c(E)$, so that
for all small enough $\varepsilon$

$R^{*}_{E}(\varepsilon)\geq 1/R_{E}(\varepsilon^{c})\;.$

It is interesting to note that Varnavides [11] (implicitly) gave a positive
answer to Problem 1.2 when $E$ is the equation $a+b=2c$. In fact, essentially
the same argument gives a positive answer to this problem for all $E$ in $3$
variables. Hence, Problem 1.2 can be considered as a generalization of
Varnavides’s Theorem. Problem 1.2 was also implicitly studied previously in
[2, 8]. In particular, Kosciuszko [8], extending earlier work of Schoen and
Sisask [10], gave direct lower bounds for $R^{*}_{E}$ which, thanks to [7],
are quasi-polynomially related to those of $R_{E}$.

The proof of Theorem 1.1 is given in the next section. For the sake of
completeness, and as a preparation for the proof of Theorem 1.1, we start the
next section with a proof that Problem 1.2 holds for equations in $3$
variables. We should point that a somewhat unusual aspect of the proof of
Theorem 1.1 is that it uses a Behrend-type [1] geometric argument in order to
find solutions, rather than avoid them.

##  2 Proofs

In the first subsection below we give a concise proof of Varnavides’s Theorem,
namely, of the fact that Problem 1.2 has a positive answer for equations with
$3$ variables. In the second subsection we prove Theorem 1.1.

###  2.1 Proof of Varnavides’s theorem

Note that for every equation $E$, there is a constant $C$ such that for every
prime $p\geq Cn$ every solution of $E$ with integers $x_{i}\in[n]$ over
$\mathbb{F}_{p}$ is also a solution over $\mathbb{R}$. Since we can always
find a prime $Cn\leq p\leq 2Cn$, this means that we can assume that $n$ itself
is prime333The factor $2C$ loss in the density of $S$ can be absorbed by the
factor $c$ in Problem 1.2. and count solutions over $\mathbb{F}_{n}$. So let
$S$ be a subset of $\mathbb{F}_{n}$ of size $\varepsilon n$ and let
$R=R_{E}(\varepsilon/2)$. For $b=(b_{0},b_{1})\in(\mathbb{F}_{n})^{2}$ and
$x\in[R]$ let $f_{b}(x)=b_{1}x+b_{0}$ and444Since $f([R])$ is a subset of
$[R]$ (rather than $S$), it might have been more accurate to denote $f([R])$
by $f^{-1}([R])$ but we drop the $-1$ to make the notation simpler.
$f_{b}([R])=\\{x\in[R]:f_{b}(x)\in S\\}$. Pick $b_{0}$ and $b_{1}$ uniformly
at random from $\mathbb{F}_{n}$ and note that for any $x\in[R]$ the integer
$f_{b}(x)$ is uniformly distributed in $\mathbb{F}_{n}$. Hence,

$\mathbb{E}|f_{b}([R])|=\varepsilon R\;.$

It is also easy to see that for every $x\neq y$ the random variables
$f_{b}(x)$ and $f_{b}(y)$ are pairwise independent. Hence

$\mathrm{Var}|f_{b}([R])|\leq\varepsilon R\;.$

Therefore, by Chebyshev’s Inequality we have

$\mathbb{P}\left[|f_{b}([R])|\leq\frac{\varepsilon}{2}R\right]\leq\frac{\varepsilon
R}{\varepsilon^{2}R^{2}/4}\leq 1/2\;.$

In other words, at least $n^{2}/2$ choices of $b$ are such that
$|f_{b}([R])|\geq\frac{\varepsilon}{2}R$. By our choice of $R$ this means that
$f_{b}([R])$ contains $3$ distinct integers $x_{1},x_{2},x_{3}$ which satisfy
$E$ and such that $f_{b}(x_{i})\in S$. Note that if $x_{1},x_{2},x_{3}$
satisfy $E$ then so do $f_{b}(x_{1}),f_{b}(x_{2}),f_{b}(x_{3})$. Let us denote
the triple $(f_{b}(x_{1}),f_{b}(x_{2}),f_{b}(x_{3}))$ by $s_{b}$. We have thus
obtained $n^{2}/2$ solutions $s_{b}$ of $E$ in $S$. To conclude the proof we
just need to estimate the number of times we have double counted each solution
$s_{b}$. Observe that for every choice of $s_{b}=\\{s_{1},s_{2},s_{3}\\}$ and
distinct $x_{1},x_{2},x_{3}\in[R]$, there is exactly one choice of
$b=(b_{0},b_{1})\in(\mathbb{F}_{n})^{2}$ for which $b_{1}x_{i}+b_{0}=s_{i}$
for every $1\leq i\leq 3$. Since $[R]$ contains at most $R^{2}$ solutions of
$E$ this means that for every solution $s_{1},s_{2},s_{3}\in S$ there are at
most $R^{2}$ choices of $b$ for which $s_{b}=\\{s_{1},s_{2},s_{3}\\}$. We
conclude that $S$ contains at least $n^{2}/2R^{2}$ distinct solutions, thus
completing the proof.

###  2.2 Proof of Theorem 1.1

As in the proof above, we assume that $n$ is a prime and count the number of
solutions of the equation $E:~{}\sum^{4}_{i=1}a_{i}x_{i}=0$ over
$\mathbb{F}_{n}$. Let $S$ be a subset of $\mathbb{F}_{n}$ of size $\varepsilon
n$, and let $d$ and $t$ be integers to be chosen later and let $X$ be some
subset of $[t]^{d}$ to be chosen later as well. For every
$b=(b_{0},\ldots,b_{d})\in(\mathbb{F}_{n})^{d+1}$ and
$x=(x_{1},\ldots,x_{d})\in X$ we use $f_{b}(x)$ to denote
$b_{0}+\sum^{d}_{i=1}b_{i}x_{i}$ and $f_{b}(X)=\\{x\in X:f_{b}(x)\in S\\}$. We
call $b$ good if $|f_{b}(X)|\geq\varepsilon|X|/2$. We claim that at least half
of all possible choices of $b$ are good. To see this, pick
$b=(b_{0},\ldots,b_{d})$ uniformly at random from $(\mathbb{F}_{n})^{d+1}$,
and note that for any $x\in X$ the integer $f_{b}(x)$ is uniformly distributed
in $\mathbb{F}_{n}$. Hence,

$\mathbb{E}|f_{b}(X)|=\varepsilon|X|\;.$

It is also easy to see that for every $x\neq y\in X$ the random variables
$f_{b}(x)$ and $f_{b}(y)$ are pairwise independent. Hence

$\mathrm{Var}|f_{b}(X)|\leq\varepsilon|X|\;.$

Therefore, by Chebyshev’s Inequality we have555We will make sure $|X|\geq
8/\varepsilon$.

$\mathbb{P}\left[\left|f_{b}(X)\right|\leq\frac{\varepsilon}{2}\left|X\right|\right]\leq
4/\varepsilon|X|\leq 1/2\;,$ (2.1)

implying that at least half of the $b$’s are good. To finish the proof we need
to make sure that every such choice of a good $b$ will “define” a solution
$s_{b}$ in $S$ in a way that $s_{b}$ will not be identical to too many other
$s_{b^{\prime}}$. This will be achieved by a careful choice of $d$, $t$ and
$X$.

We first choose $X$ to be the largest subset of $[t]^{d}$ containing no three
points on one line. We claim that

$|X|\geq t^{d-2}/d\;.$ (2.2)

Indeed, for an integer $r$ let $B_{r}$ be the points $x\in[t]^{d}$ satisfying
$\sum^{d}_{i=1}x^{2}_{i}=r$. Then every point of $[t]^{d}$ lies on one such
$B_{r}$, where $1\leq r\leq dt^{2}$. Hence, at least one such $B_{r}$ contains
at least $t^{d-2}/d$ of the points of $[t]^{d}$. Furthermore, since each set
$B_{r}$ is a subset of a sphere, it does not contain three points on one line.

We now turn to choose $t$ and $d$. Let $C$ be such that
$R_{E}(\varepsilon)\leq(1/\varepsilon)^{C}$. Set $a=\sum^{4}_{i=1}|a_{i}|$ and
pick $t$ and $d$ satisfying

$(1/\varepsilon)^{2C}\geq
t^{d}\geq\left(\frac{2dt^{2}a^{d}}{\varepsilon}\right)^{C}\;.$ (2.3)

Taking $t=2^{\sqrt{\log 1/\varepsilon}}$ and $d=2C\sqrt{\log 1/\varepsilon}$
satisfies666Recalling (2.2), we see that since $C\geq 1$ (indeed, a standard
probabilistic deletion method argument shows that if an equation has $k$
variables, then $C(E)\geq 1+\frac{1}{k-2}$), we indeed have $|X|\geq
8/\varepsilon$ as we promised earlier. the above for all small enough
$\varepsilon$. Note that by (2.3) and our choice of $C$ we have
$R_{E}\left(\frac{\varepsilon}{2dt^{2}a^{d}}\right)\leq t^{d}$.

Let us call a collection of $4$ vectors $x^{1},x^{2},x^{3},x^{4}\in X$ helpful
if they are distinct, and they satisfy $E$ in each coordinate, that is, for
every $1\leq i\leq d$ we have $\sum^{4}_{j=1}a_{j}x^{j}_{i}=0$. We claim that
for every good $r$, there are useful $x^{1},x^{2},x^{3},x^{4}\in f_{r}(X)$. To
see this let $M$ denote the integers $1,\ldots,(at)^{d}$ and note that (2.2)
along with the fact that $r$ is good implies that

$|f_{r}(X)|\geq\varepsilon|X|/2\geq\frac{\varepsilon
t^{d}}{2t^{2}d}=\frac{\varepsilon}{2dt^{2}a^{d}}\cdot|M|$ (2.4)

Now think of every $d$-tuple $x\in X$ as representing an integer $p(x)\in[M]$
written in base $at$. So we can also think of $f_{r}(X)$ as a subset of $[M]$
of density at least $\varepsilon/2dt^{2}a^{d}$. By (2.3), we have

$M=(at)^{d}\geq t^{d}\geq
R_{E}\left(\frac{\varepsilon}{2dt^{2}a^{d}}\right)\;,$

implying that there are distinct $x^{1},x^{2},x^{3},x^{4}\in f_{r}(X)$ for
which $\sum^{4}_{j=1}a_{j}\cdot p(x^{j})=0$. But note that since the entries
of $x^{1},x^{2},x^{3},x^{4}$ are from $[t]$ there is no carry when evaluating
$\sum^{4}_{j=1}a_{j}\cdot p(x^{j})$ in base $at$, implying that
$x^{1},x^{2},x^{3},x^{4}$ satisfy $E$ in each coordinate. Finally, the fact
that $\sum_{j}a_{j}=0$ and that $x_{i}^{1},x_{i}^{2},x_{i}^{3},x_{i}^{4}$
satisfy $E$ for each $1\leq i\leq d$ allows us to deduce that

$\sum^{4}_{j=1}a_{j}\cdot
f_{b}(x^{j})=\sum^{4}_{j=1}a_{j}\cdot(b_{0}+\sum^{d}_{i=1}b_{i}x^{j}_{i})=\sum^{d}_{i=1}b_{i}\cdot(\sum^{4}_{j=1}a_{j}x^{j}_{i})=0\;,$

which means that $f_{b}(x^{1}),f_{b}(x^{2}),f_{b}(x^{3}),f_{b}(x^{4})$ forms a
solution of $E$. So for every good $b$, let $s_{b}$ be (some choice of)
$f_{b}(x^{1}),f_{b}(x^{2}),f_{b}(x^{3}),f_{b}(x^{4})\in S$ as defined above.
We know from (2.1) that at least $n^{d+1}/2$ of all choices of $b$ are good,
so we have thus obtained $n^{d+1}/2$ solutions $s_{b}$ of $E$ in $S$. To
finish the proof we need to bound the number of times we have counted the same
solution in $S$, that is, the number of $b$ for which $s_{b}$ can equal a
certain $4$-tuple in $S$ satisfying $E$.

Fix $s=\\{s_{1},s_{2},s_{3},s_{4}\\}$ and recall that $s_{b}=s$ only if there
is a helpful $4$-tuple $x^{1},x^{2},x^{3},x^{4}$ (as defined just before
equation (2.4)) such that $f_{b}(x^{i})=s_{i}$. We claim that for every
helpful $4$-tuple $x^{1},x^{2},x^{3},x^{4}$, there are at most $n^{d-2}$
choices of $b=(b_{0},\ldots,b_{d})$ for which $s_{b}=s$. Indeed recall that by
our choice of $X$ the vectors $x^{1},x^{2},x^{3}$ are distinct and do not lie
on one line. Hence they are affine independent777That is, if we turn these
three $d$-dimensional vectors into $(d+1)$-dimensional vectors, by adding a
new coordinate whose value is $1$, we get three linearly independent vectors.
over $\mathbb{R}$. But since the entries of $x^{i}$ belong to $[t]$ and $t\leq
1/\varepsilon$ we see that for large enough $n$ the vectors
$x^{1},x^{2},x^{3}$ are also affine independent over $\mathbb{F}_{n}$. This
means that the system of three linear equations

$b_{0}+b_{1}x^{1}_{1}+\ldots+b_{d}x^{1}_{d}=s_{1}$
$b_{0}+b_{1}x^{2}_{1}+\ldots+b_{d}x^{2}_{d}=s_{2}$
$b_{0}+b_{1}x^{3}_{1}+\ldots+b_{d}x^{3}_{d}=s_{3}$

(in $d+1$ unknowns $b_{0},\ldots,b_{d}$ over $\mathbb{F}_{n}$) has only
$n^{d-2}$ solutions, implying the desired bound on the number of choices of
$b$. Since $|X|\leq t^{d}\leq(1/\varepsilon)^{2C}$ by (2.3) we see that $X$
contains at most $(1/\varepsilon)^{8C}$ helpful $4$-tuples. Altogether this
means that for every $s_{1},s_{2},s_{3},s_{4}\in S$ satisfying $E$, there are
at most $(1/\varepsilon)^{8C}n^{d-2}$ choices of $b$ for which $s_{b}=s$.
Since we have previously deduced that $S$ contains at least
$\frac{1}{2}n^{d+1}$ solutions $s_{b}$, we get that $S$ contains at least
$\frac{1}{2}\varepsilon^{8C}n^{3}$ distinct solutions, as needed.

Acknowledgment: I would like to thank Yuval Wigderson and an anonymous referee
for helpful comments and suggestions.

## References

  * [1] F. A. Behrend, On sets of integers which contain no three terms in arithmetical progression, Proc. Nat. Acad. Sci. U. S. A. 32 (1946), 331–-332.
  * [2] T. F. Bloom, Translation invariant equations and the method of Sanders, Bull. Lond. Math. Soc. 44 (2012), 1050-–1067.
  * [3] B. Bukh, Extremal graphs without exponentially-small bicliques, Duke Math. J., to appear.
  * [4] D. Conlon, J. H. Kim, C. Lee and J. Lee, Some advances on Sidorenko’s conjecture, J. London Math. Soc. 98 (2018), 593–608.
  * [5] P. Erdős and M. Simonovits, Supersaturated graphs and hypergraphs, Combinatorica 3 (1983), 181-–192.
  * [6] A. Girão, E. Hurley, F. Illingworth and L. Michel, Abundance: asymmetric graph removal lemmas and integer solutions to linear equations, arXiv:2310.18202.
  * [7] Z. Kelley and R. Meka, Strong bounds for 3-progressions, arXiv:2302.05537.
  * [8] T. Kosciuszko, Counting solutions to invariant equations in dense sets, arXiv:2306.08567.
  * [9] I. Z. Ruzsa, Solving a linear equation in a set of integers I, Acta Arith. 65 (1993), 259–282.
  * [10] T. Schoen and O. Sisask, Roth’s theorem for four variables and additive structures in sums of sparse sets, Forum Math. Sigma 4 (2016) e5, 28.
  * [11] P. Varnavides, On certain sets of positive density, J. London Math. Soc. 34 (1959) 358–360.


# Renormalization Group as a Koopman Operator

###### Abstract

Koopman operator theory is shown to be directly related to the renormalization
group. This observation allows us, with no assumption of translational
invariance, to compute the critical exponents $\eta$ and $\delta$, as well as
ratios of critical exponents, of classical spin systems from single
observables alone. This broadens the types of problems that the
renormalization group framework can be applied to and establish universality
classes of. In addition, this connection may allow for a new, data-driven way
in which to find the renormalization group fixed point(s), and their relevant
and irrelevant directions.

## Renormalization Group as a Koopman Operator

William T. Redman

University of California, Santa Barbara

<EMAIL_ADDRESS>
##  1 Introduction

Understanding the behavior of multi-component complex systems is a central
problem in all areas of science. In physics, the study of phase transitions
and critical phenomena led to the development of the renormalization group
(RG), one of the most successful tools in theoretical physics [1, 2, 3, 4]. By
turning disparate physical problems into dynamical systems problems, it
enabled a coherent understanding of the universality of critical exponents and
allowed a controlled, albeit difficult, manner in which to compute them.

Given its success, there have been attempts in applying the RG framework to
problems beyond traditional physics. A recent example was the discovery of an
exact mapping between the RG and deep neural networks (DNNs) based on
Restricted Boltzmann machines (RBMs) [5]. However, a challenge in directly
importing ideas and techniques from RG theory to these systems is that they do
not have the same symmetries as the problems that RG theory was developed on
have (e.g. translational invariance). Therefore, the extension of the RG
framework to compute the critical exponents, and thereby define the
universality classes, of non-symmetric systems is an important open question
that likely requires a new approach.

A recent advance in the study of nonlinear dynamical systems is Koopman
operator theory (KOT). In addition to being based on a foundation of rigorous
mathematical theory, there are a number of effective, data-driven algorithms
inspired by KOT that extract the dynamics underlying a wide range of physical
systems. KOT has been indirectly linked with the RG in several different
fields. First, both the RG and KOT have been used to successfully discover the
dynamics of partial differential equations [6, 7, 8]. Second, the RG has been
used to extend principal component analysis (PCA) and perform dimensionality
reduction [9], while KOT has been shown to be related to PCA [10, 11]. Third,
recent work applying normal form theory to the RG equations gave a unified
approach to characterize the nonlinear generalizations of scaling functions
[12]. KOT has connections to normal form theory [13, 14, 15]. Finally, the RG
has been shown to have an exact mapping to RBM DNNs [5], and KOT has been used
to develop novel machine learning techniques [16, 17, 18].

These connections between the RG and KOT motivated us to look more closely at
the two. Does there exist a direct link between them? As we show below, the RG
is, in fact, a Koopman operator.

This paper is organized as follows. We begin by providing a brief overview of
basic KOT. This will allow us to show that, by definition, the RG is a Koopman
operator. We include references to a number of recent theoretical and applied
research, which provide more detail about KOT, for the interested reader. We
then use a KOT inspired algorithm to compute two critical exponents, and two
ratios of critical exponents, of the Ising and three-state Potts models in 2D.
We show that our algorithm requires only measuring a single observable, and
does not rely on explicit calculations of how applying the RG changes coupling
constants. We compare these results to the standard Monte Carlo RG (MCRG)
approach. We find that our method performs similarly, while requiring less
information about the exact RG flow, and is considerably faster. We end by
outlining possible uses of this Koopman RG. In particular, we suggest that it
can be used to compute the universality classes of more complex systems, such
as RBM DNNs and spin glasses. Finally, we note that, by using more advanced
KOT methods, we may be able to compute the RG fixed point(s), and their
relevant (unstable) and irrelevant (stable) directions, in a data-driven way.

##  2  KOT and its connection to the RG operator

KOT is a spectral dynamical systems theory that was first developed by Bernard
Koopman in 1931 in the context of classical mechanics [19], and then later
expanded upon by Koopman and John von Neumann in 1932 [20]. It has seen a
great increase in attention over the past two decades as a wave of new data-
driven methods [21, 22, 23, 24, 15] and underlying mathematical theory [15,
26, 25] has allowed it to be applied to the dynamics of fluids [27, 28, 29,
30], and to the study of power grids [31, 32], logistics [33], urban
insurgency [34], and building energy [35].

The key insight in KOT is that there exists an infinite dimensional linear
operator, the Koopman operator (which is related to the composition operator),
whose spectrum provides information on the dynamics of nonlinear systems. The
Koopman operator, $U^{t}$, is defined, to be the time evolution operator of a
given observable $\mathbf{g}$

$U^{t}\mathbf{g}\left(\mathbf{x}_{0}\right)=\mathbf{g}\left(\mathbf{S}^{t}\left(\mathbf{x}_{0}\right)\right)$
(1)

where $\mathbf{S}^{t}$ is the dynamics that act on $\mathbf{g}$ [15]. Time
here can either be discrete or continuous.

The block spin RG is defined to be a map in the infinite dimensional space of
possible Hamiltonians with coupling constants $\left(K_{1},K_{2},...\right)$,
which we refer to as K-space from now on [36]. While K-space is infinite
dimensional, because we are ultimately interested in real systems that are
assumed to have a finite interaction length, we will restrict ourselves to
assuming that K-space is finite dimensional, with $d_{K}$ dimensions, and with
a “vanilla” $L^{2}$ function space on it. This additionally allows us to avoid
problems in nonlinear functional analysis that come with considering an
infinite dimensional space.

The block spin RG transformation, $R_{b}$, acting on a Hamiltonian with
coupling constants $\mathbf{K}_{0}=(K_{1},K_{2},...,K_{d_{K}})$, is equivalent
to a coarse grained Hamiltonian with new coupling constants. This can be done
$n$ times

$R_{b}^{n}H(\mathbf{K}_{0})=H(T^{n}(\mathbf{K}_{0}))$ (2)

where $T$ is the transformation from one point in K-space to another,
following the chosen blocking procedure. Note that, when considering the
evolution of $H(\mathbf{K}_{0})$ in this way (with the RG iteration number
acting as a “time”), we are defining a dynamical system. It is from this that
can connect the RG to KOT.

Eqs. 1 and 2 immediately imply that, by definition, the RG is a Koopman
operator in K-space. Because the Hamiltonian defines the value of all of the
observables of interest (such as magnetization), we can apply KOT methods to
their measured quantities to gain information about $R_{b}$. Although we only
consider the block spin RG here, the Wilsonian (momentum space) RG, which
integrates over continuous degrees of freedom, is a continuous Koopman
operator.

The finite section method (also commonly known as the Galerkin projection) is
a simple algorithm that is frequently used to compute the Koopman operator
from time series data [15, 22]. It finds the approximate Koopman operator,
$\tilde{U}$, from the data matrix $F$, which is comprised of the first $n$
time points of $m$ observables $\mathbf{g}_{1},...,\mathbf{g}_{m}$. That is
$F=\Big{[}\mathbf{g}_{1}|...|\mathbf{g}_{m}\Big{]}$. In particular,
$\tilde{U}$ is given by

$\tilde{U}=F^{+}F^{\prime}$ (3)

where $F^{\prime}$ is the data matrix shifted one time step forward and
$F^{+}$ is the Moore-Penrose pseudoinverse of $F$ (i.e.
$F^{+}=(F^{\dagger}F)^{-1}F^{\dagger}$). In practice, a single observable is
often used to populate the data matrix, with each column being time delayed
from the others [37]. This allows for a rich, informative spectrum even from a
single observable.

##  3 Results

Having found that the RG is a Koopman operator in K-space, we explored whether
we could successfully apply tools from KOT to calculate the critical exponents
of the 2D Ising model and the 2D three-state Potts model.

These classical spin models have a Hamiltonian of the form

$H(\mathbf{K})=-\sum_{i}K_{1}s_{i}-\sum_{\langle
i,j\rangle}K_{2}s_{i}s_{j}-...$ (4)

where the $\\{s_{i}\\}$ are the spins of the system, and the $K_{i}$ are the
strengths of the different interaction types (e.g. $K_{2}$ is the strength of
the nearest neighbor coupling). In these models, $K_{1}$ is interpreted as an
external magnetic field, and is often referred to as $h$ instead.

We considered the magnetization, $m$, which scales as

$m\sim t^{\beta}$ (5)

near the critical manifold, where $t$ is the difference between the critical
temperature, $T_{c}$, and actual temperature (i.e. $t=(T_{c}-T)/T_{c}$) [36].
Because the magnetization is a function of the Hamiltonian, there exists a
Koopman operator, $U_{m}$, that evolves the magnetization in RG iteration
time. This is a proxy for the block spin RG, $R_{b}$, and its spectrum is
related to $\beta$ (see below).

To numerically approximate $U_{m}$, we first equilibrated spin systems using
standard Monte Carlo approaches [38]. We then performed block spin
renormalizations, with $b^{2}$ spins in each block, $n_{\text{R}}$ times,
measuring the magnetization after each blocking (see Supplemental Material for
more details about our procedure). This gave $n_{\text{R}}+1$ values of $m$.
We used the finite section method, Eq. 3, to get

$\tilde{U}_{m}=\left(m_{0},...,m_{n_{\text{R}}-1}\right)^{+}\left(m_{1},...,m_{n_{\text{R}}}\right)^{\dagger}$
(6)

where $m_{j}$ is the magnetization of the spin system after $j$ block spin
renormalizations.

To see how the spectrum of $\tilde{U}_{m}$ is related to $\beta$, note that
near the critical temperature

$m_{k}\sim t_{k}^{\beta}\sim b^{\beta k/\nu}t^{\beta}$ (7)

where $\nu$ is the critical exponent related to the change in $t$ from
applying the block spin RG (i.e. $t_{k}\sim b^{k/\nu}t$) and to the
correlation length (i.e. $\xi\sim t^{-\nu}$) [36]. Therefore,

$\begin{split}\left(m_{1},...,m_{n_{\text{R}}}\right)&\sim\left(b^{\beta/\nu}t^{\beta},...,b^{\beta
n_{R}/\nu}t^{\beta}\right)\\\
&=b^{\beta/\nu}\left(t^{\beta},...,b^{\beta(n_{R}-1)/\nu}t^{\beta}\right)\\\
&=b^{\beta/\nu}\left(m_{0},...,m_{n_{\text{R}}-1}\right)\end{split}$ (8)

and so, by Eq. 6, we have that

$\tilde{U}_{m}\sim b^{\beta/\nu}$ (9)

Taking $\log_{b}$ of Eq. 9 gives $\beta/\nu$. This ratio of critical exponents
is informative, as the larger it is, the more strongly the magnetization
changes as we apply $R_{b}$ and move away from the critical manifold.
Additionally, it tells us how the magnetization and correlation length compare
in their sensitivities to changes in $t$. Note that, even if we had enriched
our construction of the Koopman operator (e.g. by including time delayed
versions of the magnetization [37]), the additional eigenvalues of
$\tilde{U}_{m}$ wouldn’t be related to other critical exponents. This is
because critical exponents are determined by the flow along the relevant
direction(s) of a given RG fixed point, and the 2D Ising model (as well as the
2D three-state Potts model) has a single relevant dimension, in the $h=0$
case, for its non-trivial fixed point.

The error in estimating $\beta/\nu$ using this method, as a function of
initial condition in K-space, is shown in Fig. 1a. We compared these results
to those obtained by using the standard Monte Carlo RG (MCRG) method [39, 40,
41, 38]. This method approximates the RG transformation near the fixed point
Hamiltonian, $H^{*}$, by the linearization

$T_{\alpha\beta}=\left[\partial K_{\alpha}^{(n+1)}/\partial
K_{\beta}^{(n)}\right]_{H^{*}}$ (10)

$\partial K_{\alpha}^{(n+1)}/\partial K_{\beta}^{(n)}$ can be solved for by
using the chain rule and certain identities that require computing spin-spin
correlations [42, 40]. If the matrix $T$ is constructed by using only even
interactions in the fixed point Hamiltonian (i.e. $K_{2},K_{4},$ etc.), $\nu$
is related to the largest eigenvalue of $T$. If $T$ is instead constructed by
using only odd interactions (or in the case of [40], only $K_{2}$ and $h$),
the critical exponent $\eta$ is related to the largest eigenvalue of $T$. The
remaining critical exponents are found using standard critical exponent
relationships [36] (see Supplemental Material).

Figure 1: Error in approximating $\beta/\nu$ for the 2D Ising model. (a) Error
in approximating $\beta/\nu$ using the KOT finite section method. (b) Error in
approximating $\beta/\nu$ using the MCRG method, while fixing $\eta=0.25$
[36]. This limits our computation to be over even K-space. (c) RG iteration
time averaged magnetization in the $(K_{2},K_{4})$ plane. This serves as a
proxy for the RG flow. (d) same as (c), but for the flow in the $(K_{2},h)$
plane. (e) Error in approximating $\beta/\nu$ using the MCRG method, setting
$h=0$, as was done in (a). For all error subplots (a, b, and e), a maximum
value of 100% (yellow) was set.

We note that the MCRG requires knowledge of the RG flow in both the odd and
even K-spaces to compute $\beta/\nu$ (as both $\eta$ and $\nu$ are needed to
solve for $\beta$), whereas the finite section method does not. If we ignore
some information, for instance, if we fix the odd interaction to be $0$ (as we
did when we using the finite section method), there is large error over a
considerable range of K-space when computing $\beta/\nu$ (Fig. 1e). This is
because the flow in the odd K-space is different from that in the even K-space
(Fig. 1c,d). Fixing the value of $\eta$ to it’s exact value of $0.250$ [36],
and searching only over even K-space, the MCRG computed $\beta/\nu$ is
slightly better as compared to our Koopman RG method, but is accurate over a
larger range of K-space (Fig. 1b). Similar results are seen when evaluating
$\beta/\nu$ in the 2D three-state Potts model (see Supplemental Material Fig.
S1).

While we have plotted $\beta/\nu$ because of it’s direct relation to the
observable we recorded across RG iteration time, we can also use standard
critical exponent relationships to compute another critical exponent ratio, as
well as two “full” critical exponents (see Supplemental Material for more
details).

$\gamma/\nu=d-2\beta/\nu$ (11) $\eta=2-d+2\beta/\nu$ (12)
$\delta=d\left(\beta/\nu\right)^{-1}-1$ (13)

where $\gamma$ is critical exponent related to the scaling of the magnetic
susceptibility near $T_{C}$, $\eta$ is the critical exponent related to the
behavior of the two point correlation function right at $T_{C}$, and $\delta$
is the critical exponent related to the scaling of the energy of the system as
a function of the magnetization near $T_{c}$. Evaluating Eqs. 11-13 using the
$\beta/\gamma$ found at $K_{2}=0.3947,K_{4}=0$, and $h=0$ from 100 different
simulations, gives mean values of $\beta/\nu=0.1269$, $\gamma/\nu=1.74632$,
$\eta=0.2537$, and $\delta=14.7652$. These compare well to the true values
$\beta/\nu=1/8$, $\gamma/\nu=7/4$, $\eta=1/4$, and $\delta=15$ [36].

Finally, our KOT method does not require any explicit calculation of how
applying the block spin RG changes K. This property shows the power of
recognizing the RG as a Koopman operator: simply by recording an observable as
we apply the RG, we can recover properties of the system in the form of
(ratios of) critical exponents. This removes the constraint of only being able
to work on systems where solutions for the linearization of $T$, $\partial
K_{\alpha}^{(n+1)}/\partial K_{\beta}^{(n)}$, are available. Therefore, the
Koopman RG allows for the calculation of (ratios of) critical exponents even
in systems that are not translationally invariant, as long as a coarse
graining procedure exists. Because there have been some that have been
successfully used (e.g. [43]), this is not a particularly strong requirement.

##  4 Discussion

We started by connecting the renormalization group to a powerful nonlinear
dynamical systems theory, Koopman operator theory. This insight was suggested
by the fact that both KOT and the RG have been shown to be related to the
asymptotics of partial differential equations [6, 7, 8], principal component
analysis [10, 11, 9], normal form theory [13, 14, 15, 12], and machine
learning [5, 16, 17, 18]. We showed that, when viewed as a dynamical system in
“RG iteration time”, by definition, the RG is a Koopman operator in the space
of coupling constants. The fact that KOT has been successfully applied to
understanding the dynamics of a wide range of systems in a data-driven manner
[27, 28, 21, 29, 32, 33, 34, 35] led us to investigate whether we could import
methods from KOT to evaluate critical exponents, and whether those methods
would afford us benefits that the standard Monte Carlo RG [38, 39, 40, 41]
approach cannot.

We showed that, for both the 2D Ising model and the 2D three-state Potts
model, the Koopman RG approach, which made use of the finite section method,
allowed for an evaluation of $\beta/\nu$ that was close to as good as that of
the MCRG on the same part of K-space, yet required less information about the
exact RG flow (Fig. 1 and Fig. S1). This is because the MCRG requires
knowledge of the RG flow in both the odd and the even K-spaces, whereas the
KOT RG method does not. Additionally, our method was significantly faster, as
we only had to measure the magnetization at each RG iteration time point,
whereas the MCRG required computing spin-spin correlations of all spin pairs
over a lengthy time period (see Supplemental Material). We also showed that,
by using critical exponent relationships [36], we were also able to compute
$\gamma/\nu$, $\eta$, and $\delta$ to good agreement with their true values.

While the robust and fast calculation of these (ratios of) critical exponents
is encouraging, the real advantage in using the Koopman RG approach is that it
does not rely on any explicit formulation of the RG transformation in terms of
$\mathbf{K}$, like the MCRG does. This greatly widens the range of systems
that we can apply the RG framework to. Applying our method to problems that
are not translationally invariant and, therefore, have been largely unstudied
using numerical real space RG methods, will be a direction of future work. The
example of deep neural networks is particularly exciting [5]. By computing
their (ratios of) critical exponents, we can establish universality classes.
Such universality classes could offer a principled way in which to explore how
architectures, learning rules, and data sets affect the performance of DNNs.
More generally, the perspective taken here has already inspired work on
accelerating the training of DNNs via KOT [44, 45]. The Koopman RG could also
be applied to spin glasses. Because spin glasses are not translationally
invariant, MCRG based methods require a mapping to an effective, nonrandom
Hamiltonian [46]. This mapping is tricky, and it is not clear that it does not
affect the results (e.g. Table 1 of [47]). While there is a plethora of work
computing the critical exponents of spin glasses via other approaches
(including many that specifically compute $\eta$), there is much disagreement
(as reviewed in [47, 48]). Whether these differences arise from nonlinear
corrections to the RG equations is something that we hope to explore using the
Koopman RG.

Finally, while our results followed from applying tools from KOT, there is
significantly more to the theory than what we have used here. For example, KOT
allows us to represent the coarse grained Hamiltonian as

$R_{b}H(\textbf{K}_{0})=\sum_{k}\mu_{k}\mathbf{v}_{k}\phi_{k}\left(\textbf{K}_{0}\right)$
(14)

where $\textbf{v}_{k}$ is the kth Koopman mode, $\phi_{k}$ is the kth Koopman
eigenfunction, and $\mu_{k}$ is the kth eigenvalue [15, 22]. As we did
earlier, we can consider the evolution of $H(\textbf{K}_{0})$ in RG iteration
time. The same representation, after $n$ applications of the RG, then gives

$R_{b}^{n}H(\textbf{K}_{0})=\sum_{k}\mu_{k}^{n}\mathbf{v}_{k}\phi_{k}\left(\textbf{K}_{0}\right)$
(15)

From this, we see that the dynamics of the RG flow, in each of the directions
of K-space defined by the Koopman modes, is determined by the magnitude its
corresponding $\mu_{k}$ [22]. While the accurate computation of the Koopman
modes and eigenfunctions is difficult (although doable [25, 26, 22]), this
representation, in principle, allows us to find the fixed points, and their
relevant and irrelevant directions, of the RG flow all from measuring a single
observable. Because these are not known for many systems, and because recent
work has shown that even relatively simple models can have non-trivial fixed
points with exotic features [49], we believe this is an area the Koopman RG
can play an important role in.

Our paper highlights the power that applying techniques developed in the field
of nonlinear dynamical systems offers when working with the RG. We see our
work very much in the same spirit as the recent success in using normal form
theory to classify nonlinear generalizations of scaling functions [12]. Given
the connection of KOT with normal form theory, and nonlinear systems in
general [13, 14, 15], it would be interesting to see whether our method could
be used to provide numerical predictions as to what universality families
systems belong to.

KOT is an exciting, powerful, and growing theory, and we hope that this
illustrates its potential as a tool for physicists.

## Acknowledgements

We would like to thank Prof. Igor Mezic for the introduction to, and helpful
discussions on, KOT. We also thank Dean Huang for his advice on early attempts
at this problem and Akshunna S. Dogra for discussions on DNNs. We thank the
referees for their constructive feedback, which strengthened the paper. The
author is supported by a Chancellor Fellowship from UCSB.

## References

  * [1] L.P. Kadanoff, Scaling laws for Ising models near $T_{c}$, Physics Physique Fizika 2, 263 (1966).
  * [2] K.G. Wilson, Renormalization group and critical phenomena. I. Renormalization group and the Kadanoff scaling picture, Phys. Rev. B 4, 3174 (1971).
  * [3] K.G. Wilson, Renormalization group and critical phenomena. II. Phase-space cell analysis of critical behavior, Phys. Rev. B 4, 3184 (1971).
  * [4] K.G. Wilson, The renormalization group: Critical phenomena and the Kondo problem, Rev. Mod. Phys. 47, 773 (1975).
  * [5] P. Mehta and D.J. Schwab, An exact mapping between the variational renormalization group and deep learning, arXiv (2014).
  * [6] L.Y. Chen, N. Goldenfeld, and Y. Oono, Renormalization group theory for global asymptotic analysis, Phys. Rev. Lett. 73, 1311 (1994).
  * [7] L.Y. Chen, N. Goldenfeld, and Y. Oono, Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory, Phys. Rev. E 54, 376 (1996).
  * [8] J. N. Kutz, J. L. Proctor, and S. L. Brunton, Koopman theory for partial differential equations, arXiv (2016).
  * [9] S. Bradde and W. Bialek, PCA meets RG, Journal of Statistical Physics 167, 462 (2017).
  * [10] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition, Journal of Neuroscience Methods 258, 1 (2016).
  * [11] S. Klus, F. Nüske, P. Koltai, H. Wu, I. Kevrekidis, C.Schutte,and F.Noe, Data-driven model reduction and transfer operator approximation, Journal of Nonlinear Science 28, 985 (2018).
  * [12] A. Raju, C. B. Clement, L. X. Hayden, J. P. Kent-Dobias, D. B. Liarte, D. Z. Rocklin, and J. P. Sethna, Normal form for renormalization groups, Phys. Rev. X 9, 021014 (2019).
  * [13] A. Mauroy and I. Mezic, Global stability analysis using the eigenfunctions of the Koopman operator, IEEE Transactions on Automatic Control 61, 3356 (2016).
  * [14] R. Mohr and I. Mezic, Koopman principle eigenfunctions and linearization of diffeomorphisms, arXiv (2016).
  * [15] I. Mezic, Spectral Koopman Operator Methods in Dynamical Systems (Springer, 2019) (in press).
  * [16] Y. Kawahara, Dynamic mode decomposition with reproducing kernels for Koopman spectral analysis, in Advances in Neural Information Processing Systems 29, edited by D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett (Curran Associates, Inc., 2016) pp. 911–919.
  * [17] N. Takeishi, Y. Kawahara, and T. Yairi, Learning Koopman invariant subspaces for dynamic mode decomposition, in Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17 (Curran Associates Inc., USA, 2017) pp. 1130–1140.
  * [18] N. Takeishi, Y. Kawahara, Y. Tabei, and T. Yairi, Bayesian dynamic mode decomposition, in Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17 (2017) pp. 2814–2821.
  * [19] B. O. Koopman, Hamiltonian systems and transformation in Hilbert space, Proceedings of the National Academy of Sciences 17, 315 (1931)
  * [20] B. O. Koopman and J. v. Neumann, Dynamical systems of continuous spectra, Proceedings of the National Academy of Sciences 18, 255 (1932)
  * [21] I. Mezic, Analysis of fluid flows via spectral properties of the Koopman operator, Annual Review of Fluid Mechanics 45, 357 (2013)
  * [22] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science 25, 1307 (2015).
  * [23] K. Fujii and Y. Kawahara, Dynamic mode decomposition in vector-valued reproducing kernel Hilbert spaces for extracting dynamical structure among observables, Neural Networks 117, 94 (2019).
  * [24] N. Takeishi, Kernel learning for data-driven spectral analysis of Koopman operators, in Proceedings of The Eleventh Asian Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 101, edited by W. S. Lee and T. Suzuki (PMLR, Nagoya, Japan, 2019) pp. 956–971.
  * [25] I. Mezic, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics 41, 309 (2005).
  * [26] I. Mezic, Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry, Journal of Nonlinear Science 10.1007/s0032-019-09598-5 (2019).
  * [27] C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter, and S. S. Henningson, Spectral analysis of nonlinear flows, Journal of Fluid Mechanics 641, 115 (2009).
  * [28] M. Budisíc, R. Mohr, and I. Mezic, Applied Koopmanism, Chaos: An Interdisciplinary Journal of Nonlinear Science 22, 047510 (2012).
  * [29] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113, 3932 (2016),
  * [30] H. Arbabi and I. Mezic, Study of dynamics in post transient flows using Koopman mode decomposition, Phys. Rev. Fluids 2, 124402 (2017).
  * [31] Y. Susuki, I. Mezíc, F. Raak, and T. Hikihara, Applied Koopman operator theory for power systems technology, Nonlinear Theory and Its Applications, IEICE 7, 430 (2016).
  * [32] M. Korda, Y. Susuki,and I. Mezic, Power grid transient stabilization using koopman model predictive control, IFAC-PapersOnLine 51, 297 (2018).
  * [33] J. Hogg, M. Fonoberova, I. Mezic, and R. Mohr, Koopman mode analysis of agent-based models of logistics processes, PLOS ONE 14, 1 (2019).
  * [34] M. Fonoberova, I. Mezic, J. Mezic, and R. Mohr, An agent-based model of urban insurgence: Effect of gathering sites and Koopman mode analysis, PLOS ONE 13, 1 (2018).
  * [35] M. Georgescu and I. Mezíc, Building energy modeling: A systematic approach to zoning and model reduction using Koopman mode analysis, Energy and Buildings 86, 794 (2015).
  * [36] N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group (Levant Books, 2005).
  * [37] F. Takens, Detecting strange attractors in turbulence. In: Rand D., Young LS. (eds) Dynamical Systems and Turbulence, Warwick 1980. Lecture Notes in Mathematics, vol 898. Springer, Berlin, Heidelberg (1981)
  * [38] M. Newman and G. Barkema, Monte Carlo Methods in Statistical Physics, Vol. 1st (Clarendon Press, 1999)
  * [39] S. K. Ma, Renormalization Group by Monte Carlo Methods, Phys. Rev. Lett. 37, 461 (1976).
  * [40] R. H. Swendsen, Monte Carlo Renormalization Group, Phys. Rev. Lett. 42, 859 (1979).
  * [41] R. H. Swendsen, Monte Carlo Calculation of Renormalized Coupling Parameters, Phys. Rev. Lett. 52, 1165 (1984).
  * [42] T. Niemeyer and J. van Leeuwn, Wilson theory for 2-dimensional Ising spin systems, Physica 71, 17 (1974).
  * [43] M. Newman and D. Watts, Renormalization group analysis of the small-world network model, Physics Letters A 263, 341 (1999).
  * [44] A. S. Dogra, Dynamical Systems and Neural Networks, arXiv:2004.11826 (2020).
  * [45] A. S. Dogra and W. T. Redman, Optimizing Neural Networks via Koopman Operator Theory, arXiv:2006.02361 (2020).
  * [46] J.-S. Wang and R.H. Swendsen, Monte Carlo renormalization-group study of Ising spin glasses, Physical Review B 37, 7745 (1988)
  * [47] T. Jorg and H.G. Katzgraber, Universality and universal finite-size scaling functions in four-dimensional Ising spin glasses, Physical Review B 77, 214426 (2008).
  * [48] H.G. Katzgraber, M. Korner, and A.P. Young, Universality in three-dimensional Ising spin glasses: A Monte Carlo study, Physical Review B 73, 224432 (2006).
  * [49] J.T. Young, A.V. Gorshkov, M. Foss-Feig, and M.F. Maghrebi, Nonequilibrium Fixed Points of Coupled Ising Models, Physical Review X 10, 011039 (2020)


# A practical optimal control approach for two-speed actuators

Alexandre Girard and H. Harry Asada This work was supported in part by The
Boeing Company, the Fonds Québécois de la recherche sur la nature et les
technologies (FQRNT) and the Natural Sciences and Engineering Research Council
of Canada (NSERC). A. Girard and H. H. Asada are with Department of Mechanical
Engineering, Massachusetts Institute of Technology, Cambridge, MA, USA, email
<EMAIL_ADDRESS>
###### Abstract

This paper addresses the closed-loop control of an actuator with both a
continuous input variable (motor torque) and a discrete input variable (mode
selection). In many applications, robots have to bear large loads while moving
slowly and also have to move quickly through the air with almost no load,
leading to conflicting requirements for their actuators. An actuator with
multiple gear ratios, like in a powertrain, can address this issue by allowing
an effective use of power over a wide range of output speed. However, having
discrete modes of operation adds complexity to the high-level control and
planning. Here a controller for two-speed actuators that automatically select
both the best gear ratio and the motor torque is developed. The approach is
to: first derive a low-dimensional model, then use dynamic programming to find
the best actions for all possible situations, and last use regression analysis
to extract simplified global feedback laws. This approach produces simple
practical nearly-optimal feedback laws. A controller that globally minimizes a
quadratic cost function is derived for a two-speed actuator prototype, global
stability is proven and performance is demonstrated experimentally.

##  I INTRODUCTION

Robotic systems often operate in distinctively different torque-speed
conditions. For instance, a legged robot has to move its leg forward quickly
through the air and, once touching the ground, it has to bear a large load
[1]. Those challenging extremum requirements often lead to the use of
oversized and inefficient actuators, which inhibit the performance of robots
that carry their actuators and energy source [2] [3], like wearable and mobile
robots. Actuators using multiple gear ratios, like in a powertrain, can
address this issue by allowing a more effective use of power over a wide range
of speeds. However, gear-shifting is a very non-linear process and its use in
a robotic context raise many issues. One challenge is the high-level control
of a system with discrete modes of operation. The issue comes from the fact
that the system has input variables of different natures: continuous variables
(motor torques) and discrete variables (mode selections), hence most control
engineering tools are not suited to tackle this problem.

For simple scenarios like steady state operations, the mode selection can be
based on simple principles. For instance, for a two-speed electric motor at a
steady speed, the best operating mode for power output could be selected based
on torque-speed curves (see Fig. 1a). Alternatively, for maximizing
accelerations, the reduction ratio could be selected based on the best motor-
load inertia match. However, finding the best mode selection policy for a
general task where the whole trajectory is important cannot be done using such
point-by-point optimums. For instance, to globally minimize the energy use of
a car, it can be advantageous to coast with the transmission in the neutral
mode and use the engine only intermittently. To find such solutions the global
behavior and whole trajectories must be considered. Furthermore, the problem
gets harder for robotic systems because of the multi-DOF aspects. For
instance, a manipulator using two two-speed actuators would have $2^{2}$
operating modes, with their force properties illustrated by 4 possible
manipulability ellipsoids at Fig. 1b. Hence, except for some basic scenario it
is not trivial to find good feedback policies for the mode selection in a
robotics context and an approach generalizable to complex situations is thus
needed.

(a) Two-speed actuator

(b) 2 DOF manipulator

Figure 1: Force-speed characteristics of different operating modes

This paper address the design of global feedback laws, for both continuous and
discrete inputs, that lead to a desired global behavior, which will be
formulated as reaching a target while minimizing a cost function. A practical
computational approach that put simplicity before absolute optimality is
proposed. The technique could be used on any low-dimensional system for a wide
range of control objectives. Here, this paper focus on demonstrating the
approach for the task of driving a 11.4 kg (25 lbs) mass for point-to-point
motions using the two-speed actuator shown at Fig. 2.

Figure 2: Dual-speed dual-motor linear actuator prototype

This actuator use a special dual-motor architecture to enable fast and
seamless gear shifting, which is critical in robotic applications. Motor 1 and
motor 2 have reduction ratios of 4:1 and 72:1 respectively and are coupled to
a ball screw (20 mm lead) through a 3-ports gearbox. The motor in use can be
selected by closing or opening the brake. In previous work, low-level
controllers were developed to control each mode of operation independently and
transitions from one mode to another [2]. In this paper, a high-level
controller that selects the operating mode and the torque command is designed
(see Fig. 3). The organization of the paper is as follow: In section I-A,
related control techniques are briefly reviewed. In section II, global
feedback laws for the two-speed actuator are derived using dynamic
programming. In section III, experimentation results using the prototype are
presented and discussed.

Figure 3: Controller architecture

###  I-A Related works

Discrete mode of operations are present in many system, such as powertrains
with multiple gear ratios, pneumatic or hydraulic systems equipped with on/off
valves and special actuators [4][3][5]. Considering the global behavior of
such systems, including the mode selection, lead to a hybrid dynamical model.
Most optimal control techniques are based on either variational approaches or
some form of gradient descent to find a trajectory that minimize a cost
function [6]. Hence those techniques cannot be used directly to optimize
discrete variables. An interesting approach to get around this problem is to
use the switching instants as optimization parameters instead [7][8]. However,
to use this approach a sequence of operating modes must be predefined first.
Mixed-integer programming has also been used to generate optimal open-loop
trajectories of dynamical system with both continuous and discrete input
variables [9] [10]. However, open-loop trajectories can be unstable and new
trajectories must be computed for each initial condition. Hence, from a
practical point of view, techniques generating feedback policies are
preferable. Most of the analytical results regarding feedback control of
hybrid systems are for specific cases, for instance the optimal feedback law
of linear hybrid systems with linear constraints for a quadratic cost function
have been shown to have a particular form [11]. One computational technique
that generate feedback laws and that can be used for non-linear systems with
any kind of constraints is dynamic programming [12]. Two disadvantages of the
techniques are however that it only works for low-dimensional systems (so
called curse of dimensionality) and also that the resulting feedback laws are
in the form of a look-up table.

##  II Control laws synthesis

###  II-A Methodology overview

The proposed approach is to discretize the continuous control problem into a
graph search problem where the discrete input actions can be considered
naturally, and then use dynamic programming to find global optimal feedback
policies. First, a simplified low-dimensional model of the controlled system
is derived to keep the problem tractable. Second, the optimal control problem
is solved using dynamic programming. Third, the resulting optimal policies are
approximated with simple feedback laws, to conduct a stability analysis and
also for the ease of implementation.

###  II-B Modeling

Using the previously developed low-level controller, when the selected mode is
changed the system automatically transit from one mode to another quickly and
seamlessly [2]. Hence, for the purpose of designing the high-level controller,
a simple model that assumes transitions are instantaneous is used, see Fig. 4.
While this model does not capture the transition behavior, it has the
advantage of capturing most of the global behavior with only two state
variables. However, this modeling simplification can be unsuited for loads
exhibiting fast dynamics

Figure 4: Lumped parameter simplified model

Two input variables are used; $u_{2}\in\\{1,2\\}$ a discrete variable
representing the mode selection and $u_{1}\in\Re^{1}$ representing the
continuous torque command sent to the motor currently in use. Hence, the model
is formulated with a hybrid form:

$\displaystyle\ddot{x}=$ $\displaystyle
f(\dot{x},u_{1},u_{2})=\textstyle\frac{1}{m_{r,i}}\left[-b_{r,i}\dot{x}+\frac{R_{i}}{L_{o}}u_{1}\right]\displaystyle,\;i=u_{2}$
(1)

where $m_{r,i}$ and $b_{r,i}$ are the reflected mass and damping of each mode
of operation:

$\displaystyle\textstyle
m_{r,i}=m_{o}+\left(\frac{R_{i}}{L_{o}}\right)^{2}J_{i}\quad\quad
b_{r,i}=b_{o}+\left(\frac{R_{i}}{L_{o}}\right)^{2}b_{i}$ (2)

In state-space form:

$\displaystyle\textstyle\underline{\dot{x}}=F(\underline{x},\underline{u})=\left[\begin{array}[]{c}\dot{x}\\\
\text{eq.(1)}\end{array}\right]\quad\underline{x}=\left[\begin{array}[]{c}x\\\
\dot{x}\end{array}\right]\;\underline{u}=\left[\begin{array}[]{c}u_{1}\\\
u_{2}\end{array}\right]$ (9)

###  II-C Problem formulation

The control problem of obtaining the global desired behavior is formulated as
minimizing a scalar cost $J$ that is a function of the state trajectory
$\underline{x}(:)$ and the inputs trajectory $\underline{u}(:)$, while
constraining both states and inputs to be in their respective domains:

$\displaystyle\operatornamewithlimits{min}\limits_{\underline{u}(:)}\;$
$\displaystyle J(\underline{x}(:),\underline{u}(:))=\int
g(\underline{x}(t),\underline{u}(t))dt$ (10) $\displaystyle s.t.\quad$
$\displaystyle\underline{\dot{x}}=F(\underline{x},\underline{u})$ (11)
$\displaystyle u_{1}\in[-0.02,0.02]\quad(Nm)$ (12) $\displaystyle
u_{2}\in\\{1,2\\}$ (13) $\displaystyle x_{1}\in[-150,150]\quad(mm)$ (14)
$\displaystyle x_{2}\in[-500,500]\quad(mm/sec)$ (15)

An additional constraint is imposed, when the output speed is greater than 20
mm/sec only the operating mode $u_{1}=1$ can be used (to limit motor speeds to
8000 RPM). Then three different additive cost functions are investigated; a
quadratic cost function $g_{q}$ where error is weighted against control
effort, a function corresponding $g_{t}$ to the minimal time problem and a
function $g_{e}$ corresponding to thermal losses in the motors, leading to a
simplified minimum wasted energy problem (neglecting mechanical losses).

$\displaystyle g_{q}(\underline{x},\underline{u})$
$\displaystyle=w_{1}\;x^{2}+w_{2}\;\dot{x}^{2}+w_{3}\;u_{1}^{2}$ (16)
$\displaystyle g_{t}(\underline{x},\underline{u})$ $\displaystyle=1$ (17)
$\displaystyle g_{e}(\underline{x},\underline{u})$
$\displaystyle=u_{1}^{2}\quad\propto\quad\dot{Q}_{motor}=rI^{2}$ (18)

###  II-D Dynamic programming parameters

A dynamic programming algorithm, also known as value iteration, is used to
find both the optimal cost-to-go function and optimal policies for an infinite
horizon. The optimal control problem is numerically solved for the different
three cost functions. The discretization parameters are as follow: the time
step is 0.02 sec, the state space is discretized into an even 501 x 501 grid
and the continuous torque is discretized into 51 discrete control options, for
a total of 102 possible control actions including the mode selection.
Artificial out-of-bound and on-target termination states are included to
guarantee convergence [13].

###  II-E Numerical results

Fig. 5-7 illustrate the numerical results of the dynamic programming algorithm
for the three different cost functions. The absence of color indicates states
with no solution (a constraint will be violated for any possible control
actions). Fig. 8 shows the closed loop behavior of the system in the phase
plane when the optimal policy is applied.

(a) Minimum time

(b) Quadratic cost

(c) Minimum energy

Figure 5: Optimal cost-to-go $J^{*}$

(a) Minimum time

(b) Quadratic cost

(c) Minimum energy

Figure 6: Optimal policy for the continuous torque command $u_{1}$ [Nm]

(a) Minimum time

(b) Quadratic cost

(c) Minimum energy

Figure 7: Optimal policy for the mode selection $u_{2}$

(a) Minimum time

(b) Quadratic cost

(c) Minimum energy

Figure 8: Closed loop behavior with the optimal policy illustrated in the
phase plane

#### Minimum time

For the minimum time problem, the optimal policy is a bang-bang law for
$u_{1}$ and always using highly-geared mode when possible. Note that the bang-
bang switching curve accounts for the fact the large gear ratio will be used
during the final part of the trajectory.

#### Quadratic cost

For the quadratic cost, the mode selection optimal policy is almost as simple
as the minimum time problem except for small features in quadrant II and IV.
The more interesting result comes from the continuous torque control law, the
gains when using the large reduction ratio are larger than those when using
the small reduction ratio. This results in the controller taking action mainly
at low speed when its actions have the biggest impacts on the system, and lead
to a highly non-linear closed-loop behavior. This is the opposite of what
would have been obtained using feedback linearization: large gains with the
small gear and small gains with the large gear to linearize the global
behavior. Hence implementing feedback linearization would have led to a poor
performance considering this quadratic cost metric.

#### Minimum energy

For the minimum energy controller, interestingly the mode selection policy is
not trivial even for this simple linear model. This shows that it does not
take much complexity to have non-trivial optimal policies for hybrid systems.
Here, the large reduction ratio is used almost only for braking, and in
quadrant II and IV the small reduction ratio is used even at low speed to
coast with low viscous resistance. Also globally the gains are much lower than
the other controllers except for zones where it is necessary to use energy to
stay in the domain.

###  II-F Control laws approximation

In this section, simplified algebraic control laws are extracted from the
numerical map of optimal control actions for the quadratic cost controller.
The dynamic programming algorithm outputs a list of optimal actions that is
assumed are samples of a map of optimal actions
$\underline{u}^{*}=\pi^{*}(\underline{x})$. The goal is then to find a
simplified map $\hat{\pi}(\underline{x})$ based on the samples. Here a two-
steps approach is used; first the state-space is segmented into zones based on
the optimal discrete action $u_{2}$, then regression is used to compute the
$u_{1}$ map independently in each zone.

#### Segmentation

For the resulting optimal policy found at Fig. 6b, three main zones can be
identified. The boundaries are approximated by two lines at $\dot{x}=\pm\;20$,
and the samples are separated into two groups:

$\displaystyle\left\\{\underline{u}^{i},\underline{x}^{i}\right\\}\rightarrow\text{Group
1}\quad\text{if}\;|\dot{x}|\geq 20$ (19)
$\displaystyle\left\\{\underline{u}^{i},\underline{x}^{i}\right\\}\rightarrow\text{Group
2}\quad\text{if}\;|\dot{x}|<20$ (20)

For more complex segmentation, machine learning techniques could be used to
perform the classification.

#### Continuous maps

In each zones, the optimal continuous input $u_{1}^{*}$ is very close to a
planar surface with exceptions when the input saturates or near the
boundaries, see Fig. 6b. Here, it is proposed to forgo absolute optimality and
instead use the simplest control law that can globally approximate $u_{1}^{*}$
map. Hence, simple linear plane are fitted on the samples group by group. The
regression is defined as follow, the dependent variable to approximate is
$u_{1}$, the independent variables are the state coordinates
$\underline{x}^{T}=[x\;\dot{x}]$ and the parameter vector is
$\underline{\alpha}^{T}=[k_{p}\;k_{d}]$. Hence the continuous control variable
$\hat{u}_{1}$ is approximated with the following linear map:

$\displaystyle\hat{u}_{1}=\underline{\alpha}^{T}\underline{x}=k_{p}x+k_{d}\dot{x}$
(21)

which is essentially a simple proportional-derivative control law. Then the
parameters are estimated using a least-square criterion,
$\underline{\alpha}_{1}$ using the samples in group 1 and
$\underline{\alpha}_{2}$ samples in group 2. Then combining the segmentation
and both linear maps, the global control policy is approximated with the
following law:

$\displaystyle\underline{\hat{u}}$
$\displaystyle=\hat{\pi}(\underline{x})=\left\\{\begin{array}[]{c}\left[\begin{array}[]{c}\underline{\alpha}_{1}^{T}\underline{x}\\\
1\end{array}\right]\text{if}\;|\dot{x}|\geq 20\\\ \\\
\left[\begin{array}[]{c}\underline{\alpha}_{2}^{T}\underline{x}\\\
2\end{array}\right]\text{if}\;|\dot{x}|<20\end{array}\right.$ (29)

Fig. 9a shows the resulting continuous control policy $\hat{u}_{1}$ for the
full state-space, and is an approximation of Fig. 6b map.

(a) Map for $u_{1}$ (Nm)

(b) Phase plane behavior

Figure 9: Approximated control laws and resulting behavior

The resulting behavior with the control law of eq.(29) is illustrated in the
phase plane at Fig. 9b, where the arrows illustrate the state derivative
throughout the state space. The blue arrows illustrate the natural dynamic of
the system and the red arrows illustrate the closed behavior with the control
policy of eq.(29). At high speed with the small reduction ratio, the
controller only takes small actions illustrated by the fact that red arrows
only have small deviation compared to blue arrows. However at low speed, the
large reduction ratio is used to drastically change the natural behavior of
the system, illustrated by the large red arrows.

###  II-G Stability analysis

Here in this section, the global stability of the controller generated for the
quadratic cost is demonstrated. When the two-speed actuator is controlled with
the feedback law of eq.(29), the closed loop behavior is given by:

$\displaystyle\ddot{x}$
$\displaystyle=\left\\{\begin{array}[]{c}\frac{1}{m_{r,1}}\left[-(b_{r,1}+\frac{R_{1}}{L_{o}}k_{d1})\dot{x}-\frac{R_{1}}{L_{o}}k_{p1}x\right]\quad\text{if}\;|\dot{x}|\geq
20\\\
\frac{1}{m_{r,2}}\left[-(b_{r,2}+\frac{R_{2}}{L_{o}}k_{d2})\dot{x}-\frac{R_{2}}{L_{o}}k_{p2}x\right]\quad\text{if}\;|\dot{x}|<20\end{array}\right.$
(32)

It will be shown that the energy (including virtual energy in the controller)
of the system is always decreasing. The energy function is discontinuous
across the switching surfaces, equal to $E_{1}$ when $|\dot{x}|\geq 20$ and
$E_{2}$ when $|\dot{x}|<20$, with value and time derivative given by:

$\displaystyle E_{i}$
$\displaystyle=\frac{1}{2}\bigg{[}m_{r,i}\bigg{]}\dot{x}^{2}+\frac{1}{2}\left[\frac{R_{i}}{L_{o}}k_{p,i}\right]x^{2}$
$\displaystyle i=\\{1,2\\}$ (33) $\displaystyle\dot{E}_{i}$
$\displaystyle=-\left[b_{r,i}+\frac{R_{i}}{L_{o}}k_{d,i}\right]\dot{x}^{2}$
$\displaystyle i=\\{1,2\\}$ (34)

If $k_{p,i}$ and $(b_{r,i}+\frac{R_{i}}{L_{o}}k_{d,i})$ are positive values
then $E_{i}>0$ and $\dot{E}_{i}<0$ for any non-zero point in the state space.
Hence all trajectories staying in their respective zone have decreasing
energy. Fig. 10a show a general trajectory on the phase plane, and Fig. 10b
shows the corresponding time-trajectory of the energy. Taken independently the
energy is decreasing on trajectory segment A ($E_{b}<E_{a}$) and B
($E_{d}<E_{c}$). However, when crossing a switching surface the energy can
increase, for instance $E_{c}>E_{b}$. For a down shift the kinetic energy
increase as the output speed stay the same but the reflected inertia is
greater with the larger reduction ratio. Note that this is not an artifact of
using a simplified actuator model: the motor needs to be accelerated to match
the output speed during the gear shift process.

(a) Phase plane trajectory

(b) Time trajectory

Figure 10: Energy along a trajectory

To prove global stability, it must be shown that the energy cannot increase by
going back and forth into different operating modes. Hence, that $E_{e}<E_{b}$
for any trajectory. First a mapping $E_{2}=\textit{h}(E_{1})$ of the energy
change when crossing the switching surface is needed. On the crossing surface
the energy is only a function of the position:

$\displaystyle E_{i}(x)$
$\displaystyle=\frac{1}{2}\bigg{[}m_{r,i}\bigg{]}s^{2}+\frac{1}{2}\left[\frac{R_{i}}{L_{o}}k_{p,i}\right]x^{2}$
(35)

where $s$ is the speed on the switching surface. The energy correspondence is:

$\displaystyle E_{2}=\textit{h}(E_{1})=\frac{1}{2}\bigg{[}\scriptstyle
m_{r,2}-\frac{R_{2}k_{p,2}}{R_{1}k_{p,1}}\bigg{]}\displaystyle
s^{2}+\bigg{[}\scriptstyle\frac{R_{2}k_{p,2}}{R_{1}k_{p,1}}\bigg{]}\displaystyle
E_{1}$ (36)

Then it is possible to prove that $E_{e}<E_{b}$ for any trajectory starting
with the inequality $E_{d}<E_{c}$ and substituting with the mapping of
eq.(36):

$\displaystyle
E_{d}<E_{c}\;\Rightarrow\;\textit{h}(E_{e})<\textit{h}(E_{b})\;\Rightarrow\;E_{e}<E_{b}$
(37)

Hence, any trajectory going into the low-speed zone and coming back will have
a smaller $E_{1}$ value on its return, thus $E_{1}$ can only decrease.
Similarly it is possible to show that $E_{2}$ is also always decreasing, for
any scenario. All in all, even though the energy can jumps when crossing a
switching surface, it is globally decreasing since both $E_{1}$ and $E_{2}$
can only decrease. Hence the system is stable and will converge to the origin.
Note that this stability analysis is only applicable to inertial loads with
linear damping. With a different type of load, this controller could be
unstable or exhibiting undesirable behavior such as rapid back-and-forth
switching between the operating modes.

##  III Experiment results and discussion

A trajectory of the system using the approximated control law of eq.(29)
starting from rest at $x$ = -120 mm is simulated using the simplified model
and results are illustrated on Fig. 9b and Fig. 11. The same experiment is
conducted using the real actuator prototype and results are illustrated on
Fig. 12. On both experiments it is possible to see clearly that the controller
actions are more aggressive when using the large reduction ratio, which lead
the desired optimal non-linear behavior, as discussed in section II-E.

Figure 11: Time trajectory (simulation) Figure 12: Time trajectory (experiment
with prototype)

The experiment with the real prototype shows that even though the control
policy was generated with a simplified linear model, the desired global
behavior is still reached. The discrepancy between the simulation and the
experiment is mainly due to the fact that the model used in the simulation
neglect the transitions dynamic, and also because the linear damping used in
the model does not describes very well the friction in the lead screw. Diverse
experiments, shown in the attached video, demonstrate that the proposed
controller takes advantage of the mode selection for a wide range of
situations. Even though the controller is not the optimal policy anymore in
those other scenarios, optimizing for an archetype point-to-point task is
found to be a practical way to tune global feedback laws that take advantage
of the mode selection.

##  IV CONCLUSION AND OUTREACH

In this paper, a practical approach based on dynamic programming is proposed
to generate nearly optimal feedback laws for robotic systems with both
continuous and discrete input variables. Optimal control policies for point-
to-point motions with a two-speed actuator prototype are derived numerically
for minimizing the time to arrival, the wasted energy or a quadratic cost
function. The optimal feedback policy for the quadratic cost is then
approximated to a simple piece-wise analytical function, global stability is
proven and performance is demonstrated with the prototype. The main advantages
of the technique are: the possibility to include any kind of nonlinearities
and constraints (such as discrete input variables), and the generation of
analytical feedback laws. Limitations are that approximations must be used for
both the system dynamic and the feedback policy function. Next, to follow this
initial work on a single actuator, the control of robots using multiple two-
speed actuators will be studied. The main challenge will be to find good
approximations so that the problem is still tractable computationally.

## References

  * [1] S. Hirose, “A Study of Design and Control of a Quadruped Walking Vehicle,” _The International Journal of Robotics Research_ , vol. 3, no. 2, pp. 113–133, June 1984.
  * [2] A. Girard and H. Asada, “A two-speed actuator for robotics with fast seamless gear shifting,” in _IEEE/RSJ Int. Conf. on Int. Robots and Systems (IROS)_ , Sept 2015, pp. 4704–4711.
  * [3] W. Zoels, I. Vittorias, and G. Bachmaier, “A novel piezohydraulic actuator as artificial muscle in robotic applications,” in _IEEE/RSJ Int. Conf. on Int. Robots and Systems (IROS)_ , Sept 2015, pp. 4712–4717.
  * [4] D. Leach, F. Gunther, N. Maheshwari, and F. Iida, “Linear multi-modal actuation through discrete coupling,” in _IEEE/RSJ Int. Conf. on Int. Robots and Systems (IROS)_ , Oct. 2012, pp. 2437–2442.
  * [5] H. J. Lee, Y. Shin, K.-S. Kim, and S. Kim, “Finger module design based on dual-mode twisting actuation for explosive ordnance disposal gripper,” in _Int. Symposium on Robotics (ISR)_ , Oct. 2013, pp. 1–3.
  * [6] J. Betts, _Practical Methods for Optimal Control and Estimation Using Nonlinear Programming_. Society for Industrial and Applied Mathematics, Jan. 2010.
  * [7] X. Xu and P. Antsaklis, “Optimal control of switched systems based on parameterization of the switching instants,” _IEEE Transactions on Automatic Control_ , vol. 49, no. 1, pp. 2–16, Jan. 2004.
  * [8] N. Majdoub, A. Sakly, and M. Benrejeb, “Hybrid approach for optimal control problem of switched systems,” in _IEEE Int. Conf. on Systems Man and Cybernetics (SMC)_ , Oct. 2010, pp. 4161–4168.
  * [9] A. Richards, T. Schouwenaars, J. P. How, and E. Feron, “Spacecraft Trajectory Planning with Avoidance Constraints Using Mixed-Integer Linear Programming,” _Journal of Guidance, Control, and Dynamics_ , vol. 25, no. 4, pp. 755–764, 2002.
  * [10] M. Gerdts, “Solving mixed-integer optimal control problems by branch&bound,” _Optimal Control Applications and Methods_ , vol. 26, no. 1, pp. 1–18, Jan. 2005.
  * [11] F. Borrelli, M. Baotić, A. Bemporad, and M. Morari, “Dynamic programming for constrained optimal control of discrete-time linear hybrid systems,” _Automatica_ , vol. 41, no. 10, pp. 1709–1721, Oct. 2005.
  * [12] Donald E. Kirk, _Optimal Control Theory: An Introduction_. Mineola, N.Y: Dover Publications, Apr. 2004.
  * [13] D. P. Bertsekas, _Dynamic programming and optimal control_. Athena Scientific Belmont, MA, 2012, vol. 2.