# Implementation of the algorithm for constructing homogeneous approximations
of nonlinear control systems

Grigory SKLYAR  Faculty of Computer Science and Information Technology
West Pomeranian University of Technology in Szczecin, Żołnierska 49, 71-210
Szczecin, Poland
e-mail<EMAIL_ADDRESS>Institute of Mathematics
University of Szczecin, Wielkopolska 15, 70-451 Szczecin, Poland
e-mail<EMAIL_ADDRESS>Pavel BARKHAYEV  Department of Mathematical
Sciences
Norwegian University of Science and Technology, Høgskoleringen 1, 7491
Trondheim, Norway
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>B.Verkin Institute for Low
Temperature Physics and Engineering
Academy of Sciences of Ukraine, Nauki Ave. 47, 61103 Kharkiv, Ukraine
Svetlana IGNATOVICH  Department of Applied Mathematics
V. N. Karazin Kharkiv National University, Svobody Sq. 4, 61022 Kharkiv,
Ukraine
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>Viktor RUSAKOV  Department
of Applied Mathematics
V. N. Karazin Kharkiv National University, Svobody Sq. 4, 61022 Kharkiv,
Ukraine
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a “calculator” for constructing a homogeneous approximation of
nonlinear control systems, which is based on the algebraic approach developed
by the authors in their previous papers. This approach mainly uses linear
algebraic and combinatorial tools, so, it is perfectly adapted to computer
realization. We describe the algorithm and discuss its capabilities and
limitations. We present its implementation as a web application and show by
example how this app works.

Keywords:  nonlinear control systems, homogeneous approximation, nonlinear
power moments, free algebra, algorithm for constructing an approximating
system.

##  1 Introduction

When studying properties of dynamical systems or problems of control design,
the first step commonly used in mechanics and engineering is to linearize the
systems under consideration. Namely, the nonlinear system is substituted by a
linear one, which is nevertheless close enough to the original one, so that
the results obtained for the linear system are applicable to the original
problem. A brilliant example of such a method is the Lyapunov stability
theory.

However, in many cases it is impossible to choose an appropriate linear
system. As a simple example, let us consider a two-input driftless control
system

$\dot{x}=X_{1}(x)u_{1}+X_{2}(x)u_{2},$ (1)

where $X_{1}(x)$ and $X_{2}(x)$ are smooth vector fields defined in a domain
of ${\mathbb{R}}^{n}$ with $n\geq 3$. Suppose that this system is locally
controllable at some point $x^{0}\in{\mathbb{R}}^{n}$, i.e., any point from a
neighborhood $U(x^{0})$ can be reached from $x(0)=x^{0}$ by a piecewise
continuous control. However, its linearization, which should be of the form

$\dot{x}=b_{1}u_{1}+b_{2}u_{2},$

where $b_{1},b_{2}\in{\mathbb{R}}^{n}$ are constant vectors, necessarily is
uncontrollable. Thus, any linear system looses the most important property of
the original nonlinear system, therefore, hardly can be considered as its
approximation.

In such a case, one can try to use an appropriate _nonlinear_ system as an
approximation. It is natural to expect that an approximating system, being
nonlinear, has a simpler structure than the original one. Deep and long-term
studies in this direction were carried out by several researches, who analyzed
different problems. But despite the diversity of approaches, in the end all
these authors came to the same concept of a so-called _homogeneous
approximation_. We mention several important original publications [6, 4, 35,
10, 1, 3, 36]; this list is undoubtedly far from being complete. The results
obtained were summarized in [11] and [2]. In the above works, the differential
geometry language and tools were used.

Another way to introduce a nonlinear approximation was related to the
interpretation of nonlinear control systems as series of iterated integrals.
This idea was originated by M. Fliess [7, 9], who proposed to apply Chen
series [5] to control problems. The approach uses the tools of free algebras
[24] to study control systems [8, 18, 16, 14].

In the mid 90s of the last century, two of the authors of the present paper,
Grigory Sklyar and Svetlana Ignatovich, considered another approximation
problem, namely, approximation in the sense of time optimality [25, 26]. They
proposed a way for analysis of the time optimal problem for nonlinear systems
based on algebraic methods [28]. The main goal was to generalize the moment
approach, which proved to be an effective tool in the linear setting [19, 20].
As a result, they gave an explicit construction of approximating systems in
the sense of time optimality using the free algebras language. It turned out
that in essence, from the algebraic point of view, the approximation in the
sense of time optimality and the homogeneous approximation are similar
concepts though related to different free algebras [27, 30, 31]. In general,
algebraic approach proved to be very natural and allowed a deeper
understanding of nonlinear systems [34, 29, 12, 13, 31, 33].

It turns out that the proposed construction, in addition to its theoretical
significance, has the following advantage, which is essential for computer
implementation: it mainly consists of purely linear algebraic and
combinatorial steps and, therefore, promises to become a universal tool
adapted to arbitrary control-affine systems. However, in all their previous
papers, Grigory Sklyar and Svetlana Ignatovich did not bring the approach to
practical calculations, although they perfectly understood the possibilities
provided by the algorithmic nature of the method. Another author of the
present paper, Pavel Barkhayev, being deeply involved in the matter,
initialized this work and created the first version of the algorithm
realization. When the last author, Viktor Rusakov, joined, the computer
program was improved and, finally, a “calculator” of homogeneous
approximations was obtained.

The main goal of the paper is to present and discuss the algorithm for
constructing a homogeneous approximation based on the mentioned algebraic
approach, its implementation, capabilities and limitations.

Concluding the introduction section, we illustrate possible advantages of a
homogeneous approximation by a very simple example. Let us consider a
stabilization problem for a rotating rigid body

$\dot{x}_{1}=\alpha_{1}x_{2}x_{3}+u_{1},\
\dot{x}_{2}=\alpha_{2}x_{1}x_{3}+u_{2},\ \dot{x}_{3}=\alpha_{3}x_{1}x_{2},$
(2)

where $\alpha_{3}\neq 0$. In the paper [37], time-varying feedback controls
were proposed, which stabilize the equilibrium $x=0$ and can be easily found
explicitly. Using a homogeneous approximation, one can simplify the proposed
controls.

Figure 1: Norms of the trajectories $x(t)$ of the system (2).

The system (2) is not homogeneous, however, its homogeneous approximation can
be easily guessed; it has the form

$\dot{x}_{1}=u_{1},\ \ \dot{x}_{2}=u_{2},\ \
\dot{x}_{3}=\alpha_{3}x_{1}x_{2}.$ (3)

Now, following [37], we find time-varying feedback controls stabilizing the
homogeneous system (3)

$\begin{array}[]{l}u_{1}(t,x)=-\gamma
x_{1}+\frac{4\pi\sqrt{|a_{120}|}}{\varepsilon}\cos\left(\frac{2\pi
t}{\varepsilon}\right),\\\\[3.0pt] u_{2}(t,x)=-\gamma
x_{2}+\frac{4\pi\sqrt{|a_{120}|}}{\varepsilon}\mathrm{sign}(a_{120})\cos\left(\frac{2\pi
t}{\varepsilon}\right),\end{array}$ (4)

where $\varepsilon>0$ is rather small, $\gamma>0$, and

$\textstyle a_{120}=-\frac{\gamma}{2\alpha_{3}}x_{3}-\frac{1}{2}x_{1}x_{2}.$

It turns out that these controls also stabilize the original system (2). We
note that they are slightly simpler than those proposed in the mentioned paper
for the system (2), since they are designed for a simpler system (3). In Fig.
1 we compare two trajectories of the system (2): the first trajectory (dotted
line) corresponds to the controls proposed in [37] and the second one (solid
line) corresponds to the controls (4); in both cases $\alpha_{1}=3$,
$\alpha_{2}=2$, $\alpha_{3}=1$, $\gamma=5$, $\varepsilon=1$, and
$x(0)=(3,2,1)^{\top}$.

We note that the controls (4) do not depend on $\alpha_{1},\alpha_{2}$ and,
therefore, can be applied to systems with different $\alpha_{1},\alpha_{2}$ or
with $\alpha_{1},\alpha_{2}$ that vary within some bounds, which leads to
robust control design. Moreover, the controls (4) stabilize, along with (2),
other systems for which (3) is a homogeneous approximation.

##  2 Algebraic approach to the homogeneous approximation problem

By definition, a homogeneous approximation is of a local character, so, we
consider control systems in a neighborhood of a given point; without loss of
generality we assume this point to be the origin. The algorithm we present in
this paper is applicable in two situations: (i) the system is autonomous,
control-linear and the starting point is fixed,

$\dot{x}=\sum_{i=1}^{m}X_{i}(x)u_{i},\ x(0)=0;$ (5)

(ii) the system is non-autonomous, control-affine and the end point, which is
an equilibrium of the system, is fixed,

$\dot{x}=a(t,x)+\sum_{i=1}^{m}b_{i}(t,x)u_{i},\ a(t,0)\equiv 0,\ x(\theta)=0.$
(6)

In this paper we concentrate on the second case, where we assume $m=1$ for
simplicity, though, the main ideas for both cases are similar. A description
of the algorithm for systems (5) can be found in [32].

In the following subsections, we briefly outline our approach, which was first
described in [28].

###  2.1 Series representation of control systems

We consider a non-autonomous single-input system of the form

$\dot{x}=a(t,x)+b(t,x)u,\ \ a(t,0)\equiv 0,$ (7)

where $a(t,x)$, $b(t,x)$ are analytic in a neighborhood of the origin in
${\mathbb{R}}^{n+1}$, and we are interested in trajectories ending at the
origin, i.e., $x(\theta)=0$, for some time moment $\theta>0$. The condition
$a(t,0)\equiv 0$ means that the origin is an equilibrium of the system. The
set of admissible controls is defined as

$u\in B^{\theta}=\\{u(t)\in L_{\infty}[0,\theta]:|u(t)|\leq 1\mbox{ a.e.}\\}.$
(8)

Then there exists $T>0$ such that for any $\theta\in(0,T)$ and any $u\in
B^{\theta}$ the solution of the differential equation
$\dot{x}=a(t,x)+b(t,x)u(t)$ satisfying the end condition ${x(\theta)=0}$ is
well defined for $t\in[0,\theta]$. Therefore, the system (7) defines the
operator $S_{a,b}$, which maps a pair $(\theta,u)$ to the starting point,

$S_{a,b}(\theta,u)=x^{0},\ \ \theta\in(0,T),\ u\in B^{\theta}.$

The operator $S_{a,b}$ admits the following representation in the form of a
series

$S_{a,b}(\theta,u)=\sum_{k\geq 1,m_{i}\geq 0}v_{m_{1}\ldots
m_{k}}\xi_{m_{1}\ldots m_{k}}(\theta,u),$ (9)

where

$\begin{array}[]{l}\xi_{m_{1}\ldots m_{k}}(\theta,u)=\\\
\displaystyle=\int_{0}^{\theta}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{k-1}}\prod_{i=1}^{k}(\tau_{i}^{m_{i}}u(\tau_{i}))d\tau_{k}\cdots
d\tau_{1}\end{array}$ (10)

are _nonlinear power moments_ and $v_{m_{1}\ldots m_{k}}\in{\mathbb{R}}^{n}$
are constant vector coefficients. A precise formula for $v_{m_{1}\ldots
m_{k}}$ is given below, see (20).

Since we are interested in local analysis, we assume that $\theta$ is small.
Admissible controls are from the set (8), hence, the following natural
asymptotic equivalence can be adopted for nonlinear power moments (10):

$\xi_{m_{1}\ldots m_{k}}(\theta,u)\sim\theta^{m_{1}+\cdots+m_{k}+k},$ (11)

i.e., the number $m=m_{1}+\cdots+m_{k}+k$ can be considered as an _order_ of
the integral $\xi_{m_{1}\ldots m_{k}}(\theta,u)$.

This order naturally defines the homogeneity property, which allows us to
introduce a homogeneous approximation of the system (7) in terms of its
series.

###### Definition 2.1

Suppose that, after some coordinate transformation in the system (7), the
series $S_{a,b}(\theta,u)$ takes the form

$(S_{a,b}(\theta,u))_{k}=z_{k}(\theta,u)+\rho_{k}(\theta,u),\ k=1,\ldots,n,$
(12)

where

  * (i)

$z_{k}(\theta,u)$ is a sum of integrals (10) of order $w_{k}$,

  * (ii)

$\rho_{k}(\theta,u)$ contains integrals (10) of order greater than $w_{k}$,

  * (iii)

the system

$\dot{x}=\widehat{a}(t,x)+\widehat{b}(t,x)u,$ (13)

which corresponds to the series

$(S_{\widehat{a},\widehat{b}}(\theta,u))_{k}=z_{k}(\theta,u),\ k=1,\ldots,n,$

is controllable.

Then (13) is called a homogeneous approximation of the system (7).

It might seem that the problem is to find a coordinate transformation taking
the series to the form (12). However, an approximation itself is defined by
$z_{k}(\theta,u)$, which are linear combinations of integrals (10). The
elements $z_{k}(\theta,u)$ can be found directly, without finding an
appropriate transformation. Moreover, only basic linear algebraic and simple
combinatorial tools can be used, as is shown in the next subsections.

###  2.2 Algebra of nonlinear power moments

It turns out that the linear span of nonlinear power moments (10) can be
considered as a free associative algebra

${\mathcal{A}}=\mathrm{Lin}\\{\xi_{m_{1}\ldots m_{k}}:m_{1},\ldots,m_{k}\geq
0,\ k\geq 1\\},$

equipped with the (formal) multiplication defined on basis elements as

$\xi_{m_{1}\ldots m_{k}}\xi_{n_{1}\ldots n_{r}}=\xi_{m_{1}\ldots
m_{k}n_{1}\ldots n_{r}}.$

The series $S_{a,b}$ defines the linear map
$v:{\mathcal{A}}\to{\mathbb{R}}^{n}$ by

$v(\xi_{m_{1}\ldots m_{k}})=v_{m_{1}\ldots m_{k}}.$

The asymptotic property (11) leads to a _canonical grading structure_ in the
algebra ${\mathcal{A}}$ defined as

${\mathcal{A}}=\sum_{m=1}^{\infty}{\mathcal{A}}^{m},$ (14)

where

${\mathcal{A}}^{m}=\mathrm{Lin}\\{\xi_{m_{1}\ldots
m_{k}}:m_{1}+\cdots+m_{k}+k=m\\}$

are the subspaces of homogeneous elements of order $m$. Below we write
$\mathrm{ord}(z)=m$ if $z\in{\mathcal{A}}^{m}$. We emphasize that subspaces
${\mathcal{A}}^{m}$ are finite-dimensional.

The free associative algebra ${\mathcal{A}}$ is closely related to the free
Lie algebra generated by the elements $\xi_{m}$, $m\geq 0$, with the Lie
bracket operation $[\ell_{1},\ell_{2}]=\ell_{1}\ell_{2}-\ell_{2}\ell_{1}$; we
denote this free Lie algebra by ${\mathcal{L}}$. The Lie algebra
${\mathcal{L}}$ is graded as well,

${\mathcal{L}}=\sum_{m=1}^{\infty}{\mathcal{L}}^{m},\ \ \
{\mathcal{L}}^{m}={\mathcal{L}}\cap{\mathcal{A}}^{m}.$

Omitting details, which can be found in [31], we mention two properties of the
map $v$.

1) Below we assume that the system (7) is locally accessible, i.e., the
Rashevsky-Chow condition is satisfied. This can be expressed as

$v({\mathcal{L}})={\mathbb{R}}^{n}.$ (15)

2) The series (9) is “realizable”, i.e., it corresponds to a control system.
In our case this requirement reduces to the following property

$\mbox{if }\ell\in{\mathcal{L}}\mbox{ and }v(\ell)=0,\mbox{ then }v(\ell
z)=0\mbox{ for any }z\in{\mathcal{A}}.$ (16)

An analytic change of coordinates in the system (7) can be expressed as a
transformation over the series $S_{a,b}$. In order to find a result of such a
transformation, one needs to find products of integrals (10). The product of
two integrals corresponds to a _shuffle_ operation in the algebra
${\mathcal{A}}$, which can be defined recursively as

$\xi_{m_{1}}\shuffle\xi_{n_{1}}=\xi_{m_{1}n_{1}}+\xi_{n_{1}m_{1}},$
$\xi_{m_{1}}\shuffle\xi_{n_{1}\ldots n_{r}}=\xi_{m_{1}n_{1}\ldots
n_{r}}+\xi_{n_{1}}(\xi_{m_{1}}\shuffle\xi_{n_{2}\ldots n_{r}}),$
$\xi_{m_{1}\ldots m_{k}}\shuffle\xi_{n_{1}\ldots n_{r}}=$
$=\xi_{m_{1}}(\xi_{m_{2}\ldots m_{k}}\shuffle\xi_{n_{1}\ldots
n_{r}})+\xi_{n_{1}}(\xi_{m_{1}\ldots m_{k}}\shuffle\xi_{n_{2}\ldots n_{r}}).$

For example, the product

$\xi_{m_{1}}(u,\theta)\cdot\xi_{n_{1}}(u,\theta)=\int_{0}^{\theta}\tau_{1}^{m_{1}}u(\tau_{1})d\tau_{1}\cdot\int_{0}^{\theta}\tau_{2}^{n_{1}}u(\tau_{2})d\tau_{2}$

equals the integral of the function
$\tau_{1}^{m_{1}}\tau_{2}^{n_{1}}u(\tau_{1})u(\tau_{2})$ over the square
$\\{(\tau_{1},\tau_{2}):0\leq\tau_{1},\tau_{2}\leq\theta\\}$, which can be
written out as a sum of two integrals of the same function over the triangles
$\\{(\tau_{1},\tau_{2}):0\leq\tau_{1}\leq\tau_{2}\leq\theta\\}$ and
$\\{(\tau_{1},\tau_{2}):0\leq\tau_{2}\leq\tau_{1}\leq\theta\\}$. Therefore,

$\xi_{m_{1}}(\theta,u)\cdot\xi_{n_{1}}(\theta,u)=\xi_{m_{1}n_{1}}(\theta,u)+\xi_{n_{1}m_{1}}(\theta,u),$

which corresponds to the shuffle product in ${\mathcal{A}}$.

###  2.3 Finding a series representation for a homogeneous approximation

Now let us return to Definition 2.1. Justifying the terms “homogeneous” and
“approximation”, it also gives a clue to finding a homogeneous approximation.
Actually, one can start with the series $S_{a,b}$ and transform it, aiming to
separate the “main part” with respect to the order in ${\mathcal{A}}$. The
goal is to find a set of homogeneous elements $z_{1},\ldots,z_{n}$ such that
the “series” $S=(z_{1},\ldots,z_{n})^{\top}$ is realizable and accessible,
i.e., satisfied conditions 1) and 2) mentioned above (see (15) and (16)).
Desired elements can be obtained step by step; many examples can be found in
[34].

However, one can obtain elements $z_{1},\ldots,z_{n}$ explicitly, without such
a brutal force way. For a given series $S_{a,b}$ of the form (9) or, what is
the same, for the linear map $v$ satisfying conditions 1) and 2), we introduce
the linear subspaces

${\mathcal{P}}^{m}=\\{\ell\in{\mathcal{L}}^{m}:v(\ell)\in
v({\mathcal{L}}^{1}+\cdots+{\mathcal{L}}^{m-1})\\},\ m\geq 1,$

and denote

${\mathcal{L}}_{a,b}=\sum_{m=1}^{\infty}{\mathcal{P}}^{m}.$

It can be shown that ${\mathcal{L}}_{a,b}$ is a Lie subalgebra of
${\mathcal{L}}$; by definition, it is graded and has codimension $n$ in
${\mathcal{L}}$. We call it _a core Lie subalgebra_ of the system (7). Also we
introduce the following graded _right ideal_ generated by
${\mathcal{L}}_{a,b}$,

${\mathcal{J}}_{a,b}=\mathrm{Lin}\\{\ell
z:\ell\in{\mathcal{L}}_{a,b},z\in{\mathcal{A}}+{\mathbb{R}}\\}.$

Finally, we introduce the inner product $\langle\cdot,\cdot\rangle$ in the
algebra ${\mathcal{A}}$ declaring elements $\xi_{m_{1}\ldots m_{k}}$ to form
an orthonormal basis. Below we denote by $\widetilde{z}$ the orthogonal
projection of an element $z$ onto the subspace ${\mathcal{J}}_{a,b}^{\perp}$,
i.e., onto the orthogonal complement to the subspace ${\mathcal{J}}_{a,b}$.

###### Theorem 2.2

For a given system (7), let us choose $n$ homogeneous elements
$\ell_{1},\ldots,\ell_{n}\in{\mathcal{L}}$ such that

$\mathrm{Lin}\\{\ell_{1},\ldots,\ell_{n}\\}+{\mathcal{L}}_{a,b}={\mathcal{L}}.$

Then the elements $z_{1},\ldots,z_{n}$ from Definition 2.1 can be chosen as
$z_{k}=\widetilde{\ell}_{k}$, $k=1,\ldots,n$. Moreover, all possible sets of
such elements can be obtained from the mentioned elements by polynomial
transformations that preserve the order.

For the proof, examples and comments, see [28], [34].

Thus, the homogeneous approximation is, in essence, unique and is defined by
the core Lie subalgebra ${\mathcal{L}}_{a,b}$. Moreover, finding a series
$S_{\widehat{a},\widehat{b}}$ for the homogeneous approximation can be carried
out in a purely algebraic way.

###  2.4 Constructing an approximating system: general case

First, we notice that the system (13) is not unique. However, if
$\widehat{a}(t,x)\equiv 0$, then $\widehat{b}(t,x)$ is defined uniquely. We
can construct components $\widehat{b}_{i}(t,x)$ of the vector function
$\widehat{b}(t,x)=(\widehat{b}_{1}(t,x),\ldots,\widehat{b}_{n}(t,x))^{\top}$
one by one in the following way.

According to Definition 2.1, the element $z_{k}=\widetilde{\ell}_{k}$ is of
order $w_{k}$, $k=1,\ldots,n$. To start with, one can show that
$\widetilde{\ell}_{1}=\alpha_{1}\xi_{w_{1}-1}$, where
$\alpha_{1}\in{\mathbb{R}}$ is nonzero. Then we set
$\widehat{b}_{1}(t,x)=-\alpha_{1}t^{w_{1}-1}$. Further we act by induction.
For any $i=2,\ldots,n$, we write the element $z_{i}=\widetilde{\ell}_{i}$ in
the form

$\widetilde{\ell}_{i}=\sum_{j=0}^{w_{i}-2}y_{j}\xi_{j}+\alpha_{i}\xi_{w_{i}-1},\
y_{j}\in{\mathcal{A}}^{w_{i}-j-1},\ \alpha_{i}\in{\mathbb{R}}.$

One can show that $y_{j}$ are shuffle polynomials of the elements
$\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1}$. This follows form the
fact that $y_{j}\in{\mathcal{J}}_{a,b}^{\perp}$ and elements
$\widetilde{\ell}_{1}^{\shuffle
q_{1}}\shuffle\cdots\shuffle\widetilde{\ell}_{n}^{\shuffle q_{n}}$ form a
basis of the subspace ${\mathcal{J}}_{a,b}^{\perp}$. Here and below we use the
notation for the “shuffle power”: $\ell^{\shuffle
q}=\ell\shuffle\cdots\shuffle\ell$ ($q$ times). Therefore,

$y_{j}=p_{j}(\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1}),\ 0\leq j\leq
w_{i}-2,$

where $p_{j}$ denotes a shuffle polynomial,

$p_{j}(\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1})=\sum\beta^{j}_{q_{1}\ldots
q_{i-1}}\widetilde{\ell}_{1}^{\shuffle
q_{1}}\shuffle\cdots\shuffle\widetilde{\ell}_{i-1}^{\shuffle q_{i-1}}$

and the sum is taken over all integers $q_{1},\ldots,q_{i-1}\geq 0$ such that
$w_{1}q_{1}+\cdots+w_{i-1}q_{i-1}=w_{i}-j-1$. Then we set

$\widehat{b}_{i}(t,x)=-\sum_{j=0}^{w_{i}-2}P_{j}(x_{1},\ldots,x_{i-1})t^{j}-\alpha_{i}t^{w_{i}-1},$

where

$P_{j}(x_{1},\ldots,x_{i-1})=\sum\beta^{j}_{q_{1}\ldots
q_{i-1}}x_{1}^{q_{1}}\cdots x_{i-1}^{q_{i-1}}.$

So, the polynomial $P_{j}$ in real variables $x_{1},\ldots,x_{i-1}$ is
constructed as follows: for any monomial $x_{1}^{q_{1}}\cdots
x_{i-1}^{q_{i-1}}$, its coefficient is the same as the coefficient of the
monomial $\widetilde{\ell}_{1}^{\shuffle
q_{1}}\shuffle\cdots\shuffle\widetilde{\ell}_{i-1}^{\shuffle q_{i-1}}$ in the
shuffle polynomial $p_{j}$.

###  2.5 Constructing an approximating system: autonomous case

If the initial system is autonomous, i.e., has the form

$\dot{x}=a(x)+b(x)u,\ \ a(0)=0,$ (17)

then a homogeneous approximation can also be chosen as an autonomous system,
which is uniquely defined. It can be constructed as follows.

Let us introduce a _derivation_ $\varphi$ in the algebra ${\mathcal{A}}$
defined by

$\varphi(\xi_{0})=0,\ \varphi(\xi_{m})=m\xi_{m-1},\ m\geq 1.$

Recall that a derivation in the algebra is a linear map satisfying the Leibniz
product rule. Hence,

$\varphi(\xi_{m_{1}\ldots m_{k}})=\sum_{i:m_{i}\geq
1}m_{i}\xi_{m_{1}\ldots(m_{i}-1)\ldots m_{k}}.$

One can show that for the autonomous system (17), if
$z\in{\mathcal{J}}_{a,b}^{\perp}$, then
$\varphi(z)\in{\mathcal{J}}_{a,b}^{\perp}$.

Denote by $\psi$ a linear map in ${\mathcal{A}}$ defined as

$\psi(\xi_{0})=1,\ \ \psi(\xi_{m_{1}\ldots m_{k-1}0})=\xi_{m_{1}\ldots
m_{k-1}}$

and $\psi(\xi_{m_{1}\ldots m_{k}})=0$ if $m_{k}\geq 1$. Again, if
$z\in{\mathcal{J}}_{a,b}^{\perp}$ and $\mathrm{ord}(z)\geq 2$, then
$\psi(z)\in{\mathcal{J}}_{a,b}^{\perp}$.

We notice that if $z\in{\mathcal{A}}^{m}$, then
$\varphi(z)\in{\mathcal{A}}^{m-1}$ and $\psi(z)\in{\mathcal{A}}^{m-1}$ if
$m\geq 2$.

Now we are ready to construct the autonomous homogeneous approximation

$\dot{x}=\widehat{a}(x)+\widehat{b}(x)u,\ \ \widehat{a}(0)=0,$

for the system (17). First, one can show that
$\widetilde{\ell}_{1}=\alpha_{1}\xi_{0}$, where $\alpha_{1}\neq 0$; we set
$\widehat{a}_{1}(x)=0$, $\widehat{b}_{1}(x)=-\alpha_{1}$. Then we act by
induction. For any $i=2,\ldots,n$, we consider the elements
$\varphi(\widetilde{\ell}_{i})$ and $\psi(\widetilde{\ell}_{i})$; since they
belong to ${\mathcal{J}}_{a,b}^{\perp}$, they can be written in the form

$\varphi(\widetilde{\ell}_{i})=p_{1,i}(\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1}),\
\psi(\widetilde{\ell}_{i})=p_{2,i}(\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1}),$

where $p_{1,i}$ and $p_{2,i}$ are shuffle polynomials,

$p_{j,i}(\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1})=\sum\beta^{j}_{q_{1}\ldots
q_{i-1}}\widetilde{\ell}_{1}^{\shuffle
q_{1}}\shuffle\cdots\shuffle\widetilde{\ell}_{i-1}^{\shuffle q_{i-1}},$

$j=1,2$, and the sums are taken over all integers $q_{1},\ldots,q_{i-1}\geq 0$
such that $w_{1}q_{1}+\cdots+w_{i-1}q_{i-1}=w_{i}-1$. Then we set

$\begin{array}[]{c}\widehat{a}_{i}(x)=-P_{1,i}(x_{1},\ldots,x_{i-1}),\\\
\widehat{b}_{i}(x)=-P_{2,i}(x_{1},\ldots,x_{i-1}),\end{array}$

where

$P_{j,i}(x_{1},\ldots,x_{i-1})=\sum\beta^{j}_{q_{1}\ldots
q_{i-1}}x_{1}^{q_{1}}\cdots x_{i-1}^{q_{i-1}}.$

Thus, the polynomial $P_{j,i}$ in real variables $x_{1},\ldots,x_{i-1}$ is
constructed so that for any monomial $x_{1}^{q_{1}}\cdots x_{i-1}^{q_{i-1}}$,
its coefficient is the same as the coefficient of the monomial
$\widetilde{\ell}_{1}^{\shuffle
q_{1}}\shuffle\cdots\shuffle\widetilde{\ell}_{i-1}^{\shuffle q_{i-1}}$ in the
shuffle polynomial $p_{j,i}$.

##  3 The algorithm

The pipeline of the algorithm for constructing a homogeneous approximation for
a given affine control system may be described naturally in a step-by-step
manner. In what follows we give a detailed description of the algorithm
splitted into nine steps.

We illustrate all the steps by the following example:

$\left\\{\begin{array}[]{l}\dot{x}_{1}=-u\cos x_{1},\\\
\dot{x}_{2}=-\sin^{2}x_{1}+t^{2}u,\\\ \dot{x}_{3}=2x_{1}^{2}\sin
t-x_{2}u,\end{array}\right.$ (18)

i.e., for $n=3$, $a(t,x)=(0,-\sin^{2}x_{1},2x_{1}^{2}\sin t)^{T}$,
$b(t,x)=(-\cos x_{1},t^{2},-x_{2})^{T}$.

Step 1. Construction of the algebra ${\mathcal{A}}$.

We fix $N\in{\mathbb{N}}$ and for any $m=1,\ldots,N$ we construct the graded
subspace ${\mathcal{A}}^{m}$, i.e., we find all the sequences
$\\{m_{1},\ldots,m_{k}\\}$ of length $k\leq m$ satisfying the condition
$m_{1}+\cdots+m_{k}=m-k$ (admissible sequences). The corresponding nonlinear
power moments $\xi_{m_{1}\ldots m_{k}}$ belong to ${\mathcal{A}}^{m}$.

To find these sequences, one may consider all ${(m-k+1)}$-ary numbers with no
more than $k$ digits and check the sum of the digits. We note that
$\dim{\mathcal{A}}^{m}=2^{m-1}$.

Having in mind the system (18), in Table 1 we give the admissible sequences
for $N=4$.

Table 1: Basis of the algebra ${\mathcal{A}}$. $m$ | Admissible sequences
---|---
1 | {0}
2 | {1}; {0, 0}
3 | {2}; {0, 1}; {1, 0}; {0, 0, 0}
4 | {3}; {0, 2}; {2, 0}; {1, 1}; {0, 0, 1}; {0, 1, 0};
| {1, 0, 0}; {0, 0, 0, 0}

Step 2. Construction of the Lie algebra ${\mathcal{L}}$.

Lie elements of the form

$[\xi_{m_{1}},[\ldots[\xi_{m_{k-1}},\xi_{m_{k}}]\ldots]],\quad k\geq
1,m_{i}\geq 0$ (19)

are called right normed elements [24, 17]. It is known that the set of all
right normed elements is complete in the Lie algebra ${\mathcal{L}}$ and does
not form a basis. However, discarding linearly dependent elements from the set
of right normed elements, one can construct a basis of ${\mathcal{L}}$.

Along with the canonical grading (14) in ${\mathcal{A}}$, we consider the
_length grading_

${\mathcal{A}}=\sum_{k=1}^{\infty}{\mathcal{B}}^{k},$

where ${\mathcal{B}}^{k}=\mathrm{Lin}\\{\xi_{m_{1}\ldots
m_{k}}:m_{1},\ldots,m_{k}\geq 0\\}$. We note that every right normed element
is homogeneous in the sense that it belongs to some subspace
${\mathcal{A}}^{m}$ as well as to some subspace ${\mathcal{B}}^{k}$. This
means that a basis of ${\mathcal{L}}$ may be constructed independently in
every single subspace ${{\mathcal{A}}^{m}\cap{\mathcal{B}}^{k}}$, $m,k\geq 1$.

For any $m=1,\ldots,N$ (the number $N$ has been chosen at Step 1) and $k\leq
m$ we consider admissible sequences $\\{m_{1},\ldots,m_{k}\\}$ such that
$\xi_{m_{1}\ldots m_{k}}\in{\mathcal{A}}^{m}\cap{\mathcal{B}}^{k}$ and
construct the corresponding right normed elements (19).

To find linearly independent elements, we consider a realization of the
subspaces ${\mathcal{A}}^{m}\cap{\mathcal{B}}^{k}$ in ${\mathbb{R}}^{q}$,
$q=\dim({\mathcal{A}}^{m}\cap{\mathcal{B}}^{k})$. Namely, let $F_{m,k}$ be a
bijection from the set of all elements $\xi_{m_{1}\ldots
m_{k}}\in{\mathcal{A}}^{m}\cap{\mathcal{B}}^{k}$ to the set of unit vectors
$\\{e_{j}\\}_{j=1}^{q}\in{\mathbb{R}}^{q}$, and we extend it by linearity to
the mapping from ${\mathcal{A}}^{m}\cap{\mathcal{B}}^{k}$ to
${\mathbb{R}}^{q}$. Further, we represent right normed elements from
${\mathcal{A}}^{m}\cap{\mathcal{B}}^{k}$ as sums of $\xi_{m_{1}\ldots m_{k}}$
and apply $F_{m,k}$. The obtained vectors in ${\mathbb{R}}^{q}$ are being
checked for linear independence (adding one by one), and we discard linearly
dependent vectors.

We store the basis of the Lie algebra ${\mathcal{L}}$ as the set of sequences
$\\{m_{1},\ldots,m_{k}\\}$ and as the corresponding vectors in
${\mathbb{R}}^{q}$.

Further we refer to the obtained set of linearly independent right normed
elements as $\\{g_{i}\\}_{i\geq 1}$, assuming that
$\mathrm{ord}(g_{i})\leq\mathrm{ord}(g_{j})$ as $i<j$.

Table 2 gives elements $g_{i}$ for $N=4$.

Table 2: Basis of the Lie algebra ${\mathcal{L}}$. $m$ | $g_{i}$
---|---
1 | $g_{1}=\xi_{0}$
2 | $g_{2}=\xi_{1}$
3 |  $g_{3}=\xi_{2}$, $g_{4}=[\xi_{0},\xi_{1}]$
4 |  $g_{5}=\xi_{3}$, $g_{6}=[\xi_{0},\xi_{2}]$, $g_{7}=[\xi_{0},[\xi_{1},\xi_{0}]]$

###### Remark 3.1

We note that there exist purely algebraic methods of constructing a basis of
the Lie algebra: Hall, Lyndon, Shirshov bases. All these algorithms are based
on the remarkable Lazard elimination theorem, see, e.g., [24].

###### Remark 3.2

Dimensions of the subspaces ${\mathcal{L}}^{m}$ are be given explicitly by
Witt’s formula

$\dim{\mathcal{L}}^{m}=\frac{1}{m}\mathop{\sum}_{d|m}\mu(d)2^{m/d}\ \mbox{ for
}m>1,$

where $\mu:{\mathbb{N}}\rightarrow\\{-1,0,1\\}$ is a Möbius function,

$\mu(d)=\left\\{\begin{array}[]{ll}1,&d\ \mbox{ is a square-free positive
integer with}\\\ &\mbox{an even number of prime factors;}\\\ -1,&d\ \mbox{ is
a square-free positive integer with}\\\ &\mbox{an odd number of prime
factors;}\\\ 0,&d\ \mbox{ has a squared prime factor.}\end{array}\right.$

In Table 3 we give comparison between $\dim{\mathcal{A}}^{m}$ and
$\dim{\mathcal{L}}^{m}$ for several first values of $m$.

Table 3: Dimensions of ${\mathcal{A}}^{m}$ and ${\mathcal{L}}^{m}$. $m$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$\dim{\mathcal{A}}^{m}$ | 1 | 2 | 4 | 8 | 16 | 32 | 64 | 128 | 256 | 512
$\dim{\mathcal{L}}^{m}$ | 1 | 1 | 2 | 3 | 6 | 9 | 18 | 30 | 56 | 99

###### Remark 3.3

Steps 1 and 2 do not depend on a specific control system; here we construct
and store basis elements of the free associative algebra ${\mathcal{A}}$ and
the corresponding Lie algebra ${\mathcal{L}}$, which can be used later for any
specific system.

We also note that both steps allow parallelizing in computing.

Step 3. Computation of vector coefficients $v_{m_{1}\ldots m_{k}}$ of the
series (9).

Let us introduce differential operators $R_{a}$, $R_{b}$ acting on a function
$\varphi=\varphi(t,x)$ as

$R_{a}\varphi=\varphi_{t}+\varphi_{x}a,\quad R_{b}\varphi=\varphi_{x}b,$

where $\varphi_{t}=\frac{\partial\varphi}{\partial t}$ and
$\varphi_{x}=(\frac{\partial\varphi}{\partial
x_{1}},\ldots,\frac{\partial\varphi}{\partial x_{n}})$. Below we apply
$R_{a}$, $R_{b}$ to vector functions assuming that they act componentwise.
Also, let us introduce the iterative operators $\mathrm{ad}_{R_{a}}^{j}$
acting on $R_{b}$ as

$\mathrm{ad}_{R_{a}}^{0}R_{b}=R_{b},\quad\mathrm{ad}_{R_{a}}^{j}R_{b}=[R_{a},\mathrm{ad}_{R_{a}}^{j-1}R_{b}],\
j\geq 1,$

where the bracket means the operator commutator. Then coefficients
$v_{m_{1}\ldots m_{k}}$ of the series (9) equal

$v_{m_{1}\ldots m_{k}}=\frac{(-1)^{k}}{m_{1}!\cdots
m_{k}!}\mathrm{ad}_{R_{a}}^{m_{1}}R_{b}\cdots\mathrm{ad}_{R_{a}}^{m_{k}}R_{b}E(x){\phantom{*}\atop\Bigl{|}_{{t=0\atop
x=0}\atop\phantom{x}}},$ (20)

where $E(x)\equiv x$ (see, e.g., [28]).

We implement the differential operators $R_{a}$, $R_{b}$ by using symbolic
differentiation by $x$ and $t$, and the operators ${\rm ad}^{j}_{R_{a}}$ are
implemented recursively. We note that the usage of recursion may be avoided,
however in this case we would have to store a set of functions.

For any $m=1,\ldots,N$ and any $\xi_{m_{1}\ldots m_{k}}\in{\mathcal{A}}^{m}$
we apply the operators $\mathrm{ad}^{m_{i}}_{R_{a}}R_{b}$ iteratively to
obtain $v_{m_{1}\ldots m_{k}}$.

For example, for the system (18), finding coefficients $v_{m_{1}\ldots m_{k}}$
for all admissible sequences $\\{m_{1},\ldots,m_{k}\\}$ (see Table 1), we
obtain the first terms of the (infinite) series (9):

$S_{a,b}=\begin{pmatrix}1\\\ 0\\\ 0\end{pmatrix}\xi_{0}-\begin{pmatrix}0\\\
1\\\ 0\end{pmatrix}\xi_{2}+\begin{pmatrix}0\\\ 2\\\
0\end{pmatrix}\xi_{01}-\begin{pmatrix}1\\\ 0\\\ 0\end{pmatrix}\xi_{000}$
$-\begin{pmatrix}0\\\ 0\\\ 1\end{pmatrix}\xi_{20}-\begin{pmatrix}0\\\ 0\\\
2\end{pmatrix}\xi_{02}+\begin{pmatrix}0\\\ 0\\\
2\end{pmatrix}\xi_{010}-\begin{pmatrix}0\\\ 0\\\
2\end{pmatrix}\xi_{001}+\ldots$

Step 4. Computation of coefficients of basis Lie elements $v(g_{i})$.

Here we use the representation of the basis Lie elements in ${\mathbb{R}}^{q}$
as linear combinations of basis elements $\xi_{m_{1}\ldots m_{k}}$ and then
use the values of $v_{m_{1}\ldots m_{k}}$ computed at the previous step taking
into account the linearity of the map $v$.

For our example (18), using Table 2 and taking into account that

$\begin{array}[]{l}g_{4}=[\xi_{0},\xi_{1}]=\xi_{01}-\xi_{10},\ \
g_{6}=[\xi_{0},\xi_{2}]=\xi_{02}-\xi_{20},\\\
g_{7}=[\xi_{0},[\xi_{1},\xi_{0}]]=2\xi_{010}-\xi_{001}-\xi_{100},\end{array}$

we have

$v(g_{1})=\begin{pmatrix}1\\\ 0\\\ 0\end{pmatrix},\ v(g_{2})=v(g_{5})=0,\
v(g_{3})=\begin{pmatrix}0\\\ -1\\\ 0\end{pmatrix},$
$v(g_{4})=\begin{pmatrix}0\\\ 2\\\ 0\end{pmatrix},\
v(g_{6})=\begin{pmatrix}0\\\ 0\\\ -1\end{pmatrix},\
v(g_{7})=\begin{pmatrix}0\\\ 0\\\ 6\end{pmatrix}.$

Step 5. Finding the first (by order) $n$ Lie elements
$\\{\ell_{j}\\}_{j=1}^{n}$ with linearly independent vector coefficients and
the generating set of the right ideal ${\mathcal{J}}_{a,b}$.

Successively for $m=1,\ldots,N$ we divide the basis Lie elements
$\\{g_{i}\\}\subset{\mathcal{L}}^{m}$ into two groups as follows.

At the beginning, both groups are empty and we put $g_{1}$ to the first group
if $v(g_{1})\not=0$ and to the second one otherwise. In the latter case we
proceed until we find $g_{i}$ such that $v(g_{i})\not=0$ and set it to the
first group. The first group is non-empty now.

General step: suppose there are $k$ ($k<n$) elements in the first group
denoted by $\ell_{1},\ldots,\ell_{k}$, and we examine the basis Lie element
$g_{i}$ with $\mathrm{ord}(g_{i})=m$.

If $v(g_{i})\not\in\mathrm{Lin}\\{v(\ell_{1}),\ldots,v(\ell_{k})\\}$, then we
add $g_{i}$ to the first group and denote $\ell_{k+1}=g_{i}$.

If not, then, if $v(g_{i})\in
v({\mathcal{L}}^{1}+\cdots+{\mathcal{L}}^{m-1})$, we add $g_{i}$ to the second
group.

If not again, we find a linear combination of elements
$\ell_{j},\ldots,\ell_{k}\in{\mathcal{L}}^{m}$ such that
$v(g_{i}+c_{j}\ell_{j}+\cdots+c_{k}\ell_{k})\in
v({\mathcal{L}}^{1}+\cdots+{\mathcal{L}}^{m-1})$ and add this linear
combination to the second group.

We perform this splitting until the first group consists of $n$ elements and
then until all $g_{i}$ such that $\mathrm{ord}(g_{i})=\mathrm{ord}(\ell_{n})$
are being checked.

We denote the elements of the second group by $\\{d_{j}\\}$; below we use them
to construct the right ideal ${\mathcal{J}}_{a,b}$.

For our example (18) we get

$\begin{array}[]{lll}\ell_{1}=g_{1},&\ell_{2}=g_{3},&\ell_{3}=g_{6}\end{array}$

and

$\begin{array}[]{llll}d_{1}=g_{2},&d_{2}=g_{4}+2g_{3},&d_{3}=g_{5},&d_{4}=g_{7}+6g_{6}.\end{array}$

Let us explain the last step, i.e., the choice of $d_{4}$. For
$g_{7}\in{\mathcal{L}}^{4}$ we see that
$v(g_{7})\in\mathrm{Lin}\\{v(\ell_{1}),v(\ell_{2}),v(\ell_{3})\\}$, however,
$v(g_{7})\not\in
v({\mathcal{L}}^{1}+{\mathcal{L}}^{2}+{\mathcal{L}}^{3})=\mathrm{Lin}\\{v(\ell_{1}),v(\ell_{2})\\}$.
Then we find that
$v(g_{7}+6\ell_{3})\in\mathrm{Lin}\\{v(\ell_{1}),v(\ell_{2})\\}$ (more
specifically, $v(g_{7}+6\ell_{3})=0$), therefore, we set
$d_{4}=g_{7}+6\ell_{3}=g_{7}+6g_{6}$.

Step 6. Construction of the right ideal ${\mathcal{J}}_{a,b}$.

For every $m\in\\{\mathrm{ord}(\ell_{i}):i=1,\ldots,n\\}$ we construct the
intersection ${\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{m}$ of the right ideal
with the subspace of elements of order $m$. To this end, we take every $d_{j}$
such that $\mathrm{ord}(d_{j})\leq m$ and construct the set

$D_{j}^{m}=\\{d_{j}\xi_{s_{1}\ldots s_{k}}:\;\xi_{s_{1}\ldots
s_{k}}\in{\mathcal{A}}^{r},\ r=m-\mathrm{ord}(d_{j})\\}$

if $\mathrm{ord}(d_{j})<m$ and $D_{j}^{m}=\\{d_{j}\\}$ if
$\mathrm{ord}(d_{j})=m$. Due to construction one has

${\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{m}=\mathrm{Lin}\\{D_{j}^{m}:\mathrm{ord}(d_{j})\leq
m\\}.$

Further we find the representation of the obtained elements
$\\{D_{j}^{m}:\mathrm{ord}(d_{j})\leq m\\}$ in ${\mathbb{R}}^{q}$,
$q=\dim{\mathcal{A}}^{m}$, considering them as linear combinations of basis
elements $\xi_{m_{1}\ldots m_{k}}$. Similarly to Step 1, we fix a bijection
$F_{m}$ from the set of all elements $\xi_{m_{1}\ldots
m_{k}}\in{\mathcal{A}}^{m}$ to the set of unit vectors
$\\{e_{i}\\}_{i=1}^{q}\in{\mathbb{R}}^{q}$, $q=\dim{\mathcal{A}}^{m}$. Then,
for any element $d_{j}\xi_{s_{1}\ldots s_{k}}\in D_{j}^{m}$, we express
$F_{m}(d_{j}\xi_{s_{1}\ldots s_{k}})=c_{1}e_{1}+\cdots+c_{q}e_{q}$ and compose
a matrix $J_{m}$ with rows $(c_{1},\ldots,c_{q})$.

It may happen that the rank of the matrix $J_{m}$ is not full. E.g., suppose
that $\xi_{0},\xi_{1}\in{\mathcal{J}}_{a,b}$, then necessarily
$[\xi_{0},\xi_{1}]\in{\mathcal{J}}_{a,b}$. If at the previous step we denoted
$d_{1}=\xi_{0}$, $d_{2}=\xi_{1}$, $d_{3}=[\xi_{0},\xi_{1}]$, then $d_{3}\in
D_{3}^{3}$ is a linear combination of elements from $D_{1}^{3}$ and
$D_{2}^{3}$ (namely, of $d_{1}\xi_{1}$ and $d_{2}\xi_{0}$). In this case the
row of $J_{3}$ corresponding to $d_{3}$ linearly depends on the rows
corresponding to $d_{1}\xi_{1}$ and $d_{2}\xi_{0}$.

It is crucial for the next step that $J_{m}$ is of full rank, so we discard
linearly dependent rows from it.

For the system (18), $\\{\mathrm{ord}(\ell_{i}):i=1,\ldots,n\\}=\\{1,3,4\\}$,
and we get

$\begin{array}[]{rcl}{\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{1}&=&\\{0\\},\\\
{\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{3}&=&\mathrm{Lin}\\{D_{1}^{3},D_{2}^{3}\\},\\\
{\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{4}&=&\mathrm{Lin}\\{D_{1}^{4},D_{2}^{4},D_{3}^{4},D_{4}^{4}\\},\end{array}$

where

$\begin{array}[]{rcl}D_{1}^{3}&=&\\{d_{1}\xi_{0}\\}=\\{\xi_{10}\\},\\\
D_{2}^{3}&=&\\{d_{2}\\}=\\{\xi_{01}-\xi_{10}+2\xi_{2}\\},\\\
D_{1}^{4}&=&\\{d_{1}\xi_{00},d_{1}\xi_{1}\\}=\\{\xi_{100},\xi_{11}\\},\\\
D_{2}^{4}&=&\\{d_{2}\xi_{0}\\}=\\{\xi_{010}-\xi_{100}+2\xi_{20}\\},\\\
D_{3}^{4}&=&\\{d_{3}\\}=\\{\xi_{3}\\},\\\
D_{4}^{4}&=&\\{d_{4}\\}=\\{2\xi_{010}-\xi_{001}-\xi_{100}+6\xi_{02}-6\xi_{20}\\}.\end{array}$

For $m=3$, considering the bijection
$F_{3}:{\mathcal{A}}^{3}\rightarrow{\mathbb{R}}^{4}$ that maps
$\\{\xi_{2},\xi_{01},\xi_{10},\xi_{000}\\}$ to $\\{e_{i}\\}_{i=1}^{4}$, we get

$J_{3}=\left(\begin{array}[]{rrrr}0&0&1&0\\\ 2&1&-1&0\end{array}\right).$

For $m=4$, considering the bijection
$F_{4}:{\mathcal{A}}^{4}\rightarrow{\mathbb{R}}^{8}$ that maps
$\\{\xi_{3},\xi_{02},\xi_{20},\xi_{11},\xi_{001},\xi_{010},\xi_{100},\xi_{0000}\\}$
to $\\{e_{i}\\}_{i=1}^{8}$, we get

$J_{4}=\left(\begin{array}[]{rrrrrrrr}0&0&0&0&0&0&1&0\\\ 0&0&0&1&0&0&0&0\\\
0&0&2&0&0&1&-1&0\\\ 1&0&0&0&0&0&0&0\\\ 0&6&-6&0&-1&2&-1&0\end{array}\right).$

Step 7. Construction of the projections of the elements
$\\{\ell_{j}\\}_{j=1}^{n}$ onto the orthogonal complement to the right ideal
${\mathcal{J}}_{a,b}$ (the main part of the series).

For every $i=1,\ldots,n$ we construct the projection of $\ell_{i}$ onto
${\mathcal{J}}_{a,b}^{\perp}$, which we denote by $\widetilde{\ell}_{i}$.

Let $\mathrm{ord}(\ell_{i})=m$. Then, actually, we construct the projection of
$\ell_{i}\in{\mathcal{A}}^{m}$ onto
${\mathcal{J}}_{a,b}^{\perp}\cap{\mathcal{A}}^{m}$. Let us consider the
bijection $F_{m}:{\mathcal{A}}^{m}\to{\mathbb{R}}^{q}$ used at Step 6. Since
$\widetilde{\ell}_{i}\in{\mathcal{J}}_{a,b}^{\perp}$, we get
$J_{m}F_{m}(\widetilde{\ell}_{i})=0$.

Let us denote
$y_{i}=\ell_{i}-\widetilde{\ell}_{i}\in{\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{m}$.
Then $F_{m}(y_{i})$ equals a linear combination of the (transposed) rows of
the matrix $J_{m}$, i.e., $F_{m}(y_{i})=J_{m}^{T}x_{i}$ for some column vector
$x_{i}$. On the other hand,
$J_{m}F_{m}(y_{i})=J_{m}F_{m}(\ell_{i}-\widetilde{\ell}_{i})=J_{m}F_{m}(\ell_{i})$.

Therefore, we get the following system of linear equations for $x_{i}$

$J_{m}J_{m}^{T}x_{i}=J_{m}F_{m}(\ell_{i}).$ (21)

The matrix $J_{m}$ is of full rank, hence, the system (21) has a unique
solution. Having found it, we then calculate $F_{m}(y_{i})=J_{m}^{T}x_{i}$,
i.e., express $y_{i}$ as a linear combination of basis elements of
${\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{m}$ (collected in $D_{j}^{m}$). Then we
find $\widetilde{\ell}_{i}=\ell_{i}-y_{i}$.

For the system (18) we have $\widetilde{\ell}_{1}=\ell_{1}=\xi_{0}$ since
${\mathcal{J}}_{a,b}\cap{\mathcal{A}}^{1}=\\{0\\}$. Further, for $\ell_{2}$,
using the bijection $F_{3}$ and the matrix $J_{3}$ from the previous step, we
find $x_{2}$ as a solution of the the equation (21), which takes the form

$\left(\begin{array}[]{rr}1&-1\\\
-1&6\end{array}\right)x_{2}=\left(\begin{array}[]{r}0\\\ 2\end{array}\right).$

Thus, $x_{2}=\frac{2}{5}(1,1)^{T}$ and therefore $y_{2}$ equals $\frac{2}{5}$
times the sum of the elements from $D_{1}^{3}$ and $D_{2}^{3}$, i.e.,
$y_{2}=\frac{2}{5}(\xi_{10}+2\xi_{2}+\xi_{01}-\xi_{10})=\frac{2}{5}(\xi_{01}+2\xi_{2})$.
Hence,

$\widetilde{\ell}_{2}=\ell_{2}-y_{2}=\xi_{2}-{\textstyle\frac{2}{5}}(\xi_{01}+2\xi_{2})={\textstyle\frac{1}{5}}\xi_{2}-{\textstyle\frac{2}{5}}\xi_{01}.$

Applying the same procedure to $\ell_{3}$ we get

$\widetilde{\ell}_{3}={\textstyle\frac{3}{19}}\xi_{02}+{\textstyle\frac{23}{285}}\xi_{20}+{\textstyle\frac{8}{57}}\xi_{001}-{\textstyle\frac{46}{285}}\xi_{010}.$

Step 8. Construction of the non-autonomous approximating system.

We follow the procedure described in Subsection 2.4. First we implement
recursively the shuffle product for basis elements:

$\xi_{m_{1}\ldots m_{k}}\shuffle\xi_{n_{1}\ldots n_{r}}.$

The result is a homogeneous sum of basis elements and we represent it as a
vector in the corresponding space ${\mathbb{R}}^{q}$, $q=2^{\sum m_{i}+\sum
n_{j}+k+r-1}$. Then we extend by linearity the shuffle operation on arbitrary
homogeneous elements of the algebra ${\mathcal{A}}$.

Further for every $i=1,\ldots,n$ we construct the function
$\widehat{b}_{i}(t,x)$ based on the representation of
$\widetilde{\ell}_{i}\in{\mathcal{A}}_{m}$:

$\widetilde{\ell}_{i}=\sum_{j=0}^{m-2}y_{j}\xi_{j}+\alpha_{i}\xi_{m-1}.$

We represent each $y_{j}$ as a shuffle polynomial with respect to
$\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1}$. To find all shuffle
products of $\\{\widetilde{\ell}_{j}\\}$ of the same order $m-j-1$, we use the
method similar to one described at the Step 1.

Let us apply this procedure to the system (18). Since
$\widetilde{\ell}_{1}=\xi_{0}$, we have $\widehat{b}_{1}(t,x)=-1$. Further,

$\widetilde{\ell}_{2}={\textstyle\frac{1}{5}}\xi_{2}-{\textstyle\frac{2}{5}}\xi_{0}\xi_{1}={\textstyle\frac{1}{5}}\xi_{2}-{\textstyle\frac{2}{5}}\widetilde{\ell}_{1}\xi_{1},$

which gives $\widehat{b}_{2}(t,x)=-\frac{1}{5}t^{2}+\frac{2}{5}tx_{1}$.
Finally,

$\widetilde{\ell}_{3}={\textstyle\frac{3}{19}}\xi_{0}\xi_{2}+{\textstyle\frac{8}{57}}\underbrace{\xi_{00}}_{=\frac{1}{2}\widetilde{\ell}_{1}^{\shuffle
2}}\xi_{1}+{\textstyle\frac{23}{285}}(\underbrace{\xi_{2}-2\xi_{01}}_{=5\widetilde{\ell}_{2}})\xi_{0}$
$={\textstyle\frac{3}{19}}\widetilde{\ell}_{1}\xi_{2}+{\textstyle\frac{4}{57}}\widetilde{\ell}_{1}^{\shuffle
2}\xi_{1}+{\textstyle\frac{23}{57}}\widetilde{\ell}_{2}\xi_{0},$

which gives
$\widehat{b}_{3}(t,x)=-\frac{3}{19}t^{2}x_{1}-\frac{4}{57}tx_{1}^{2}-\frac{23}{57}x_{2}$,
and we obtain the following non-stationary approximation:

$\left\\{\begin{array}[]{l}\dot{x}_{1}=-u,\\\
\dot{x}_{2}=(-\frac{1}{5}t^{2}+\frac{2}{5}tx_{1})u,\\\
\dot{x}_{3}=(-\frac{3}{19}t^{2}x_{1}-\frac{4}{57}tx_{1}^{2}-\frac{23}{57}x_{2})u.\end{array}\right.$
(22)

Step 9. Construction of the autonomous approximating system.

We follow the procedure described in Subsection 2.5. To this end, we implement
the linear maps $\varphi$ and $\psi$ in the algebra ${\mathcal{A}}$ and use
the shuffle product implemented at the previous step.

For every $i=1,\ldots,n$ we find $\varphi(\widetilde{\ell}_{i})$,
$\psi(\widetilde{\ell}_{i})$ and represent them as shuffle polynomials with
respect to $\widetilde{\ell}_{1},\ldots,\widetilde{\ell}_{i-1}$ if possible.
It may happen that the corresponding system for
$\varphi(\widetilde{\ell}_{i})$ is unsolvable, which means that there is no
autonomous approximating system. Again, to find all shuffle products of
$\\{\widetilde{\ell}_{j}\\}$ of a given order, we use the method similar to
one described at Step 1.

Applying this procedure to the system (18), we see that the autonomous
approximating system does not exist since
$\varphi(\widetilde{\ell}_{2})=\frac{2}{5}(\xi_{1}-\xi_{00})$ cannot be
represented as a shuffle polynomial with respect to
$\widetilde{\ell}_{1}=\xi_{0}$.

Let us change the system (18) a little bit by adding the term $-2tx_{1}$ to
the second equation:

$\left\\{\begin{array}[]{l}\dot{x}_{1}=-u\cos x_{1},\\\
\dot{x}_{2}=-\sin^{2}x_{1}-2tx_{1}+t^{2}u,\\\ \dot{x}_{3}=2x_{1}^{2}\sin
t-x_{2}u.\end{array}\right.$ (23)

For this system, $v(\xi_{2})=0$ and

$\ell_{1}=\xi_{0},\quad\ell_{2}=[\xi_{0},\xi_{1}],\quad\ell_{3}=[\xi_{0},\xi_{2}].$

Applying the described steps to the system (23), we construct the following
autonomous approximation:

$\left\\{\begin{array}[]{l}\dot{x}_{1}=-u,\\\
\dot{x}_{2}=-\frac{1}{2}x_{1}^{2},\\\
\dot{x}_{3}=\frac{1}{27}x_{1}^{3}-\frac{10}{9}x_{2}+\frac{4}{9}x_{2}u.\end{array}\right.$

##  4 Implementation of the algorithm and perspectives

The described above functionality and the pipeline of the algorithm were
implemented as a package and a web application by using Python programming
language. Here are the links:

the web application: https://happycontrol.pythonanywhere.com/,

the package: https://github.com/ViktorRusakov/happy-control.

The web application was built with the Django framework and has the following
functionality. First, one has to choose system’s dimension and type of the
problem to be solved. Then one enters coefficients of the system by using
LaTeX text formatting (and possibly the point in which neighborhood the system
is considered). The output is a pdf file with detailed calculations similar to
the ones from the example given in the previous section. In Fig. 2 we give a
screenshot of the application running.

Figure 2: Screenshot of the web application: setting a system for
approximation

The functionality on which the application is built, was implemented as an
open source package, so it is possible to use it both as an end-to-end
solution or as particular functions.

The core packages used for calculations are SymPy and NumPy:

  * •

SymPy is used for symbolic calculations when constructing series
representation of control system (Step 3 requires differentiation operation)
and when printing out nice formatted output pdf files;

  * •

NumPy is used for numerical calculations: finding linearly independent Lie
elements, their projections, etc. Basically all matrix operations were done
with NumPy.

The most computationally demanding are Steps 7, 8 and 9 since they require
solving high dimensional systems of linear equations. Since values of
$\dim{\mathcal{A}}^{m}$ grow as $2^{m-1}$, one can conclude that the “depth”
of a system, i.e., the order of the last linearly independent Lie element
$\mathrm{ord}(\ell_{n})$, determines the complexity and time of computations
rather than dimension $n$ of a system.

The present-day implementation of the algorithm possesses the following
restrictions on the depth of systems: the graded subspaces and the basis of
the Lie algebra were constructed for orders up to $15$; the web application
can handle in real time systems with ${\mathrm{ord}(\ell_{n})\leq 9}$; by
using the functionality of the package one can construct approximations for
systems with ${\mathrm{ord}(\ell_{n})\leq 13}$ on an “average” local computer
in a “reasonable” time.

As possible developments of the proposed algorithm and its implementation we
consider the following:

  1. 1)

optimization of the present implementation to be able to construct
approximations for “deeper” systems;

  2. 2)

implementation of the algorithm for the case of multi-dimensional control (for
that, optimization is even more crucial since dimensions of graded subspaces
for a system with $k$ controls grows as $(k+1)^{m-1}$);

  3. 3)

applying the algorithm as an effective tool for analysis of certain specific
classes of systems, in particular, Goursat systems, see, e.g., [23].

##  5 Physical example: Goursat systems and their homogeneous approximations

Let us consider a driftless control system (1) defined by two vector fields
and introduce the following notation for the “derived flag”

$\begin{array}[]{l}D^{1}=\mathrm{Lin}\\{X_{1},X_{2}\\},\\\
D^{k+1}=D^{k}+[D^{k},D^{k}],\ k\geq 1,\end{array}$

where $[D^{k},D^{k}]=\mathrm{Lin}\\{[Y,Z]:Y,Z\in D^{k}\\}$ and $[\cdot,\cdot]$
stands for the Lie bracket of vector fields. By definition, (1) is called a
Goursat system if $\dim(D^{k})=k+1$ for $k=1,\ldots,n-1$. A convenient
coordinate description was found in [21]: it turned out that, in appropriate
local coordinates, vector fields $X_{1}(x)$ and $X_{2}(x)$ are polynomial and
can include constant parameters.

These systems have an important application discovered by F. Jean [15]: it
relates to a kinematic model of a car pulling $n$ trailers. Roughly speaking,
variables here are orientation angles of the car and the trailers. If the car
is not perpendicular to the first trailer and none of the consecutive trailers
are perpendicular to each other, then the position is non-singular. In this
case, constant parameters vanish, and all non-singular positions are described
by a unique system (of a chained form). Otherwise, different singularities
arise; a thorough research was made in [15].

The study of possible homogeneous approximations of Goursat systems, which
describe a possible local behavior of the mentioned kinematic model, is very
difficult in general; the main question here is if constant parameters
mentioned above are included in a homogeneous approximation. A deep
investigation was carried out by P. Mormul [22, 23], who performed a lot of
sophisticated calculations applying very special properties of Goursat
systems. Actually, Goursat systems were our main motivation in developing an
implementation of the algorithm presented above.

The described above restrictions on the depth of systems for the current
implementation of the algorithm applies to Goursat systems as well. One of the
“deepest” system for which we were able to calculate a homogeneous
approximation is GGSGSG (written in a standard coding for Goursat systems, we
refer to [15], [23] for more details): this is a system of dimension 8 and Lie
brackets up to length 12 are involved. However, we believe that based on the
special structure of Goursat systems, the implementation of the algorithm can
be significantly optimized, and we will be able to find homogeneous
approximations for Goursat systems of higher dimension.

### Acknowledgments

The work was financially supported by Polish National Science Centre grant no.
2017/25/B /ST1/01892. Pavel Barkhayev was supported by the Norwegian Research
Council project ”COMAN” No. 275113.

## References

  * [1] A. A. Agrachev, R. V. Gamkrelidze, and A. V. Sarychev.  Local invariants of smooth control systems.  Acta Appl. Math., 14:191–237, 1989.  DOI: 10.1007/BF01307214.
  * [2] A. Bellaïche.  The tangent space in sub-Riemannian geometry.  In A. Bellaïche and J. J. Risler, editors, Sub-Riemannian Geometry, volume 144 of Progress in Mathematics, pages 1–78. Birkhäuser Basel, 1996.  DOI: 10.1007/978-3-0348-9210-0_1.
  * [3] R. M. Bianchini and G. Stefani.  Graded approximation and controllability along a trajectory.  SIAM J. Control Optimiz., 28:903–924, 1990.  DOI: 10.1137/0328050.
  * [4] A. Bressan.  Local asymptotic approximation of nonlinear control systems.  Internat. J. Control, 41:1331–1336, 1985.
  * [5] K. T. Chen.  Integration of paths – a faithful representation of parths by noncommutative formal power series.  Trans. Amer. Math. Soc, 89:395–407, 1958.
  * [6] P. E. Crouch.  Solvable approximations to control systems.  SIAM J. Control Optimiz., 22:40–54, 1984.  DOI: 10.1137/0322004.
  * [7] M. Fliess.  Développements fonctionnels en indéterminées non commutatives des solutions d’équations différentielles non linéaires forcées.  C. R. Acad. Sci. Paris Ser. A-B, 287:1133–1135, 1978.
  * [8] M. Fliess.  Realizations of nonlinear systems and abstract transitive Lie algebras.  Bull. Amer. Math. Soc., 2:444–446, 1980.  DOI: 10.1090/S0273-0979-1980-14760-6.
  * [9] M. Fliess.  Fonctionnelles causales non linéaires et indéterminées non commutatives.  Bull. Soc. Math. France, 109:3–40, 1981.
  * [10] H. Hermes.  Nilpotent approximations of control systems and distributions.  SIAM J. Control Optimiz., 24:731–736, 1986.  DOI: 10.1137/0324045.
  * [11] H. Hermes.  Nilpotent and high-order approximations of vector field systems.  SIAM Rev., 33:238–264, 1991.  DOI: 10.1137/1033050.
  * [12] S. Yu. Ignatovich.  Realizable growth vectors of affine control systems.  J. Dyn. Control Syst., 15:557–585, 2009.  10.1007/s10883-009-9075-y.
  * [13] S. Yu. Ignatovich.  Normalization of homogeneous approximations of symmetric affine control systems with two controls.  J. Dyn. Control Syst., 17:1–48, 2011.  10.1007/s10883-011-9109-0.
  * [14] A. Isidori.  Nonlinear Control Systems.  Springer-Verlag London, 3rd edition, 1995.
  * [15] F. Jean.  The car with $n$ trailers: Characterization of the singular configurations.  ESAIM Control Optim. Calc. Var., 1:241–266, 1996.  10.1051/cocv:1996108.
  * [16] M. Kawski.  Nonlinear control and combinatorics of words.  In Geometry of Feedback and Optimal Control, pages 305–346. Dekker, 1997.
  * [17] M. Kawski.  The combinatorics of nonlinear controllability and noncommuting flows.  volume 8 of Abdus Salam ICTP Lecture Notes Series, pages 223–312. 2002.
  * [18] M. Kawski and H. J. Sussmann.  Noncommutative power series and formal Lie-algebraic techniques in nonlinear control theory.  In Operators, Systems and Linear Algebra. European Consortium for Mathematics in Industry, pages 111–128. Teubner, 1997.  DOI: 10.1007/978-3-663-09823-2_10.
  * [19] V. I. Korobov and G. M. Sklyar.  Time-optimality and the power moment problem (Russian).  Mat. Sb. (N.S.), 134(176)(2):186–206, 1987.  Translation: Math. USSR-Sb., vol. 62, no. 1, pp. 185-206, 1989, DOI: 10.1070/SM1989v062n01ABEH003235.
  * [20] V. I. Korobov and G. M. Sklyar.  The Markov moment min-problem and time optimality (Russian).  Sibirsk. Mat. Zh., 32(1):60–71, 1991.  Translation: Siberian Math. J., vol. 32, no. 1, pp. 46-55, 1991, DOI: 10.1007/BF00970159.
  * [21] A. Kumpera and C. Ruiz.  Sur l’équilalence locale des systèmes de Pfaff en drapeau.  In Monge-Ampère Equations and Related Topics, Florence 1980, Inst. Alta Math. F. Severi, Roma, pages 201–248. 1982.
  * [22] P. Mormul.  Real moduli in local classification of Goursat flags.  Hokkaido Math J., 34:1–35, 2005.  DOI: 10.14492/hokmj/1285766199.
  * [23] P. Mormul.  Do moduli of Goursat distributions appear on the level of nilpotent approximations?  In J.P. Brasselet and M.A.S. Ruas, editors, Real and Complex Singularities, volume 39 of Trends in Mathematics, pages 229–246. Birkhäuser Basel, 2006.  DOI: 10.1007/978-3-7643-7776-2_17.
  * [24] C. Reutenauer.  Free Lie Algebras.  Clarendon Press, Oxford, 1993.
  * [25] G. M. Sklyar and S. Yu. Ignatovich.  Local equivalence of time-optimal control problems and Markov moment min-problem.  Z. Angew. Math. Mech., 76(S.3):561–562, 1996.
  * [26] G. M. Sklyar and S. Yu. Ignatovich.  Moment approach to nonlinear time optimality.  SIAM J. Control Optimiz., 38(6):1707–1728, 2000.  DOI: 10.1137/S0363012997329767.
  * [27] G. M. Sklyar and S. Yu. Ignatovich.  Representations of control systems in the Fliess algebra and in the algebra of nonlinear power moments.  Systems Control Lett., 47(3):227–235, 2002.  DOI: 10.1016/S0167-6911(02)00201-3.
  * [28] G. M. Sklyar and S. Yu. Ignatovich.  Approximation of time-optimal control problems via nonlinear power moment min-problems.  SIAM J. Control Optimiz., 42:1325–1346, 2003.  DOI: 10.1137/S0363012901398253.
  * [29] G. M. Sklyar and S. Yu. Ignatovich.  Description of all privileged coordinates in the homogeneous approximation problem for nonlinear control systems.  C. R. Math. Acad. Sci. Paris, 344:109–114, 2007.  DOI: 10.1016/j.crma.2006.11.016.
  * [30] G. M. Sklyar and S. Yu. Ignatovich.  Fliess series, a generalization of the Ree’s theorem, and an algebraic approach to a homogeneous approximation problem.  Int. J. Control, 81:369–378, 2008.  DOI: 10.1080/00207170701561427.
  * [31] G. M. Sklyar and S. Yu. Ignatovich.  Free algebras and noncommutative power series in the analysis of nonlinear control systems: an application to approximation problems.  Dissertationes Math. (Rozprawy Mat.), 504:1–88, 2014.  DOI: 10.4064/dm504-0-1.
  * [32] G. M. Sklyar and S. Yu. Ignatovich.  Constructing of a homogeneous approximation.  In A. Bartoszewicz, J. Kabziński, and J. Kacprzyk, editors, Advanced, Contemporary Control, volume 1196 of Advances in Intelligent Systems and Computing, pages 611–624. Springer, Cham., 2020.  DOI: 10.1007/978-3-030-50936-1_52.
  * [33] G. M. Sklyar and S. Yu. Ignatovich.  Subspaces of maximal singularity for homogeneous control systems.  J. Dyn. Control Syst., 27:585–616, 2021.  DOI: 10.1007/s10883-020-09518-x.
  * [34] G. M. Sklyar, S. Yu. Ignatovich, and P. Yu. Barkhayev.  Algebraic classification of nonlinear steering problems with constraints on control.  In Advances in mathematics research. Vol. 6, pages 37–96. Nova Sci. Publ., New York, 2005.
  * [35] G. Stefani.  Polynomial approximations to control systems and local controllability.  In 1985 24th IEEE Conference on Decision and Control, pages 33–38, 1985.  DOI: 10.1109/CDC.1985.268467.
  * [36] H. J. Sussmann.  New differential geometric methods in nonholonomic path finding.  In Systems, Models and Feedback: Theory and Applications, volume 12 of Progress in Systems and Control Theory, pages 365–384. Birkhäuser, Boston, MA, 1992.  DOI: 10.1007/978-1-4757-2204-8_24.
  * [37] A. Zuyev and V. Grushkovskaya.  On stabilization of nonlinear systems with drift by time-varying feedback laws.  In 2019 12th International Workshop on Robot Motion and Control (RoMoCo), Poznan, Poland, pages 9–14, 2019.  doi: 10.1109/RoMoCo.2019.8787353.


$d(\rho_{\textup{ColdMW}}^{1},\rho_{\textup{ColdMW}}^{\textup{no corr}})\leq
k$ and that trace distance decreases under partial trace. By (P.ii), Eq. (356)
tells us that $\forall i$,

$[\rho^{1}_{\textup{ColdMW}}]_{i}=[\underline{\rho^{1}_{\textup{ColdMW}}}]_{i}+o(k).$
(357)

With Eq.(357), we may relate the $F_{\infty}$ quantities for the states
$\rho_{\textup{ColdMW}}^{1}$ and $\underline{\rho_{\textup{ColdMW}}^{1}}$.
From Eq. (37), we have

$\displaystyle F_{\infty}\big{(}\rho^{1}_{\textup{ColdMW}}(k,\rho^{\textup{{no
corr}}}_{\textup{Cold}\textup{M}\textup{W}}(g)$
$\displaystyle,\rho^{\textup{{corr}}}_{\textup{Cold}\textup{M}\textup{W}}),\tau^{h}_{\textup{Cold}\textup{M}\textup{W}}\big{)}$
(358)
$\displaystyle=\ln\max_{i}\left\\{\frac{[\rho^{1}_{\textup{Cold}\textup{M}\textup{W}}]_{i}}{\tau_{i}}\right\\},$
(359)
$\displaystyle=\ln\max_{i}\left\\{\frac{[\underline{\rho^{1}_{\textup{ColdMW}}}]_{i}}{\tau_{i}}\right\\}+o(k),$
(360)
$\displaystyle=F_{\infty}\left(\tau(g)\otimes\rho^{0}_{\textup{M}}(g)\otimes\rho^{1}_{\textup{W}},\tau^{h}_{\textup{Cold}\textup{M}\textup{W}}\right)+o(k),$
(361)

where we used Eq. (37) in the last line.

The next step is to evaluate the restriction on $W_{\textup{ext}}$ that
satisfies

$\displaystyle
F_{\infty}(\tau_{\textup{Cold}}^{0}\otimes\rho^{0}_{\textup{M}}\otimes\rho^{0}_{\textup{W}},\tau^{h}_{\textup{Cold}\textup{M}\textup{W}})$
$\displaystyle\geq
F_{\infty}\left(\rho^{1}_{\textup{Cold}\textup{M}\textup{W}}(k,\rho^{\textup{{no
corr}}}_{\textup{Cold}\textup{M}\textup{W}}(g),\rho^{\textup{{corr}}}_{\textup{Cold}\textup{M}\textup{W}}),\tau^{h}_{\textup{Cold}\textup{M}\textup{W}}\right)$
(362)
$\displaystyle=F_{\infty}\left(\tau(g)\otimes\rho^{0}_{\textup{M}}(g)\otimes\rho^{1}_{\textup{W}},\tau^{h}_{\textup{Cold}\textup{M}\textup{W}}\right)+o(k),$
(363)

for $W_{\textup{ext}}$ up to order $o(k)$. Taking into account the additivity
of $F_{\infty}$ under tensor product, we can rearrange Eq. (363) to provide an
upper bound on $W_{\textup{ext}}$,

$W_{\textup{ext}}\leq\frac{ng}{\beta_{h}}\left[\gamma(\infty)+{\Theta}(f(g))+o(k/g)\right],$
(364)

where $\displaystyle\lim_{g\rightarrow 0^{+}}f(g)=0$, $\gamma(\infty)$ is
given by Eq. (252). The bound in (364) is achievable since the $F_{\infty}$
conditions imposed by Eq. (362) are achievable with equality.

Lastly, by using the expression for efficiency in Eqs. (291), and substituting
$W_{\textup{ext}}$ from Eq. (364) (with equality for the maximum possible
$W_{\textup{ext}}$) followed by $\Delta C$ from Eq. (255), we have

$\displaystyle\sup_{W_{\textup{ext}}>0}\eta(\rho^{1}_{\textup{Cold}},W_{\textup{ext}})$
$\displaystyle=\sup_{W_{\textup{ext}}>0}\left(1-\varepsilon+\frac{\Delta
C}{W_{\textup{ext}}}\right)^{-1}$ (365)
$\displaystyle=\left[1+\frac{\beta_{h}}{(\beta_{c}-\beta_{h})}\frac{\gamma(1)}{\gamma(\infty)}\right]^{-1}+{\Theta}(f(g))+o(k/g)+{\Theta}(g)+{\Theta}(\varepsilon).$
(366)

Hence using Eqs. (292), (366) we find Eq. (349). Note that in Eq (349) we have
an inequality, this is due to the fact that in the optimisation problem Eq.
(292), there is an additional constraint (namely mean energy conservation)
which is not taken into account in the derivation of Eq. (366). ∎

Finally, the above lemmas allow us to conclude that allowing further
correlations in the final state cannot allow us to achieve the Carnot
efficiency when $\Omega>1$.

######  Theorem 3.

[Correlations do not improve efficiency] Suppose that $\Omega>1$.
Parametrizing the final state of the heat engine by Eq. (326)-(329), the
quantum efficiency $\eta^{\textup{qm}}$ defined in Eq. (287) is strictly upper
bounded by the Carnot efficiency,

$\sup_{k\in[0,1],~{}\rho^{\textup{no
corr}}_{\textup{ColdMW}}}\eta^{\textup{qm}}\left(\rho^{1}_{\textup{ColdMW}}(k,\rho^{\textup{no
corr}}_{\textup{ColdMW}},\,\rho^{\textup{corr}}_{\textup{ColdMW}})\right)<1-\frac{\beta_{h}}{\beta_{c}}.$
(367)

###### Proof.

From Lemma 17, we have that both
$\eta^{\textup{qm}}\leq\eta^{\textup{qm}}_{1}$ and
$\eta^{\textup{qm}}\leq\eta^{\textup{qm}}_{\infty}$ hold. Thus a necessary
condition to achieve the Carnot efficiency for a particular
$\rho^{1}_{\textup{ColdMW}}$, is that both $\eta^{\textup{qm}}_{1}$ and
$\eta^{\textup{qm}}_{\infty}$ are equal to or greater than the Carnot
efficiency.

Lemma 20 proves that Eqs. (331) and (332) are necessary conditions for
$\eta^{\textup{qm}}_{1}$ to achieve the Carnot efficiency. However, when Eqs.
(331), (332) are satisfied, then Lemma 21 provides an upper bound on the
efficiency $\eta^{\textup{qm}}_{\infty}$ in Eq. (349).

Now, suppose $\Omega>1$. Since it is shown in Eq. (253) that
$\gamma(1)/\gamma(\infty)=\Omega$, plugging this into the leading term
appearing in Eq. (349)

$\left[1+\frac{\beta_{h}}{(\beta_{c}-\beta_{h})}\frac{\gamma(1)}{\gamma(\infty)}\right]^{-1},$
(368)

we have that the quantity $\eta_{\infty}^{\textup{qm}}$ (and therefore also
$\eta^{\textup{qm}}$) is strictly less than the Carnot efficiency
$1-\beta_{h}/\beta_{c}$. ∎

###  F.2 A more general final battery state

For the simplicity of our analysis, we have assumed that the battery is left
in the specific final state described in Eq. (23), i.e. an amount of work
$W_{\textup{ext}}=E_{k}-E_{j}$ is extracted, except with failure probability
$\varepsilon$ that the battery remains in the initial state
$\left|E_{j}\right\rangle\\!\\!\left\langle E_{j}\right|_{\textup{W}}$. In
this section, we show that this is a simplification which can be removed in
general, i.e. the final battery state is allowed to be any state within the
$\varepsilon$-ball of $\left|E_{k}\right\rangle\\!\\!\left\langle
E_{k}\right|_{\textup{W}}$. In particular, our result that the Carnot
Efficiency cannot be achieved when $\Omega>1$ still holds.

In Lemma 22, we show that for any final state of the cold bath
$\rho_{\textup{Cold}}^{1}$, allowing a more general final battery state does
not affect the amount of work bounded by the $F_{\infty}$ condition. We then
use this to prove in Theorem 4 that when $\Omega>1$, Carnot cannot be achieved
even if we allow a more general battery final state.

######  Lemma 22.

For any given $\rho_{\textup{Cold}}^{0},\rho_{\textup{Cold}}^{1}$, with
$\rho_{\textup{W}}^{0}=\left|E_{j}\right\rangle\\!\\!\left\langle
E_{j}\right|_{\textup{W}}$, consider the maximum
$W_{\infty}^{1}:=E_{k_{1}}-E_{j}$ such that
$\rho_{\textup{Cold}}^{0}\otimes\rho_{\textup{W}}^{0}\rightarrow\rho_{\textup{Cold}}^{1}\otimes\rho_{\textup{W}}^{1}$
is allowed by the non-increasing $F_{\infty}$ condition (Eq. (36)) i.e.

$D_{\infty}(\rho_{\textup{Cold}}^{0}\|\tau_{\textup{Cold}}^{\beta_{h}})+D_{\infty}(\rho_{\textup{W}}^{0}\|\tau_{\textup{W}}^{\beta_{h}})\geq
D_{\infty}(\rho_{\textup{Cold}}^{1}\|\tau_{\textup{Cold}}^{\beta_{h}})+D_{\infty}(\rho_{\textup{W}}^{1}\|\tau_{\textup{W}}^{\beta_{h}}),$
(369)

with

$\rho_{\textup{W}}^{1}=(1-\varepsilon)\left|E_{k_{1}}\right\rangle\\!\\!\left\langle
E_{k_{1}}\right|_{\textup{W}}+\varepsilon\left|E_{j}\right\rangle\\!\\!\left\langle
E_{j}\right|_{\textup{W}}.$ (370)

On the other hand, consider any battery final state

$\rho_{\textup{W}}^{2}=(1-\varepsilon)\left|E_{k_{2}}\right\rangle\\!\\!\left\langle
E_{k_{2}}\right|_{\textup{W}}+\varepsilon\rho_{\textup{W}}^{\textup{junk}},$
(371)

where $\rho_{\textup{W}}^{\textup{junk}}$ is an energy-diagonal state
orthogonal to $\left|E_{k_{2}}\right\rangle\\!\\!\left\langle
E_{k_{2}}\right|_{\textup{W}}$ which may depend on $\varepsilon$, i.e.
$\rho_{\textup{W}}^{\textup{junk}}=\sum_{i}p_{i}\left|E_{i}\right\rangle\\!\\!\left\langle
E_{i}\right|_{\textup{W}}$ with $p_{k_{2}}=0$ and $\sum_{i}p_{i}=1$. Define
$W_{\infty}^{2}:=E_{k_{2}}-E_{j}$ such that
$\rho_{\textup{Cold}}^{0}\otimes\rho_{\textup{W}}^{0}\rightarrow\rho_{\textup{Cold}}^{1}\otimes\rho_{\textup{W}}^{2}$
is allowed by the non-increasing $F_{\infty}$ condition, i.e.

$D_{\infty}(\rho_{\textup{Cold}}^{0}\|\tau_{\textup{Cold}}^{\beta_{h}})+D_{\infty}(\rho_{\textup{W}}^{0}\|\tau_{\textup{W}}^{\beta_{h}})\geq
D_{\infty}(\rho_{\textup{Cold}}^{1}\|\tau_{\textup{Cold}}^{\beta_{h}})+D_{\infty}(\rho_{\textup{W}}^{2}\|\tau_{\textup{W}}^{\beta_{h}}).$
(372)

Then for all $0<\varepsilon\leq\hat{\varepsilon}=\left[1+e^{\beta_{h}(E_{\rm
max}-E_{j})}\right]^{-1}$, we have $W_{\infty}^{1}=W_{\infty}^{2}$.

###### Proof.

Firstly, note that any energy-diagonal state $\rho_{\textup{W}}^{2}$ with
trace distance
$d(\rho_{\textup{W}}^{2},\left|E_{k_{2}}\right\rangle\\!\\!\left\langle
E_{k_{2}}\right|_{\textup{W}})=\varepsilon$ can be written in the form of Eq.
(371). Rearranging the terms in Eq. (369),

$\displaystyle
D_{\infty}(\rho_{\textup{W}}^{1}\|\tau_{\textup{W}}^{\beta_{h}})$
$\displaystyle\leq
D_{\infty}(\rho_{\textup{W}}^{0}\|\tau_{\textup{W}}^{\beta_{h}})+D_{\infty}(\rho_{\textup{Cold}}^{0}\|\tau_{\textup{Cold}}^{\beta_{h}})-D_{\infty}(\rho_{\textup{Cold}}^{1}\|\tau_{\textup{Cold}}^{\beta_{h}})=:A.$
(373)

One can use the definition of $D_{\infty}$ in Eq. (41) to expand the L.H.S. of
Eq. (373), obtaining

$\displaystyle\log\max\\{(1-\varepsilon)e^{\beta_{h}E_{k_{1}}},\varepsilon
e^{\beta_{h}E_{j}}\\}\leq A-\log Z_{\textup{W}}^{\beta_{h}}.$ (374)

We know that since near perfect work is extracted, $\varepsilon$ is
arbitrarily small. This implies that for $\varepsilon$ small enough,
$\displaystyle\max\\{(1-\varepsilon)e^{\beta_{h}E_{k_{1}}},\varepsilon
e^{\beta_{h}E_{j}}\\}=(1-\varepsilon)e^{\beta_{h}E_{k_{1}}}$.

Similarly, one can evaluate Eq. (369) to obtain

$\displaystyle\log\max\\{(1-\varepsilon)e^{\beta_{h}E_{k_{2}}},\\{\varepsilon
p_{i}e^{\beta_{h}E_{i}}\\}_{i\neq k_{2}}\\}$ $\displaystyle\leq A-\log
Z_{\textup{W}}^{\beta_{h}}.$ (375)

Note that the maximization in Eq. (375) only picks out the maximum value. In
particular, denoting $E_{\rm max}$ to be the largest energy eigenvalue of the
battery, then whenever

$(1-\varepsilon)e^{\beta_{h}E_{k_{2}}}\geq\varepsilon e^{\beta_{h}E_{\rm
max}},$ (376)

or equivalently

$\varepsilon\leq\left[1+e^{\beta_{h}(E_{\rm max}-E_{k_{2}})}\right]^{-1},$
(377)

then $\displaystyle\max\\{(1-\varepsilon)e^{\beta_{h}E_{k_{2}}},\\{\varepsilon
p_{i}e^{\beta_{h}E_{i}}\\}_{i\neq
k_{2}}\\}=(1-\varepsilon)e^{\beta_{h}E_{k_{2}}}$. In other words, as long as
$\varepsilon$ is upper bounded by Eq. (377), we know which terms attains the
maximization in Eq. (374). However, we also want an upper bound that is
independent of any limit involving the final state
$\rho^{1}_{\textup{Cold}\textup{M}\textup{W}}$ we wish to take, or any amount
of work extracted (and therefore, we want the bound to be independent of
$E_{k_{2}}$). As such, let us construct the following upper bound
$\varepsilon\leq\hat{\varepsilon}$ where,

$\hat{\varepsilon}:=\inf_{\begin{subarray}{c}E_{k_{2}}\\\
W^{2}_{\infty}>0\end{subarray}}\left[1+e^{\beta_{h}(E_{\rm
max}-E_{k_{2}})}\right]^{-1}=\left[1+e^{\beta_{h}(E_{\rm
max}-E_{j})}\right]^{-1}$ (378)

Now, we see that $E_{k_{1}}$ and $E_{k_{2}}$ correspond to the solutions for
Eq. (374) and Eq. (375), which for $\varepsilon\leq\hat{\varepsilon}$ reduce
to exactly the same equation. Therefore, $E_{k_{1}}=E_{k_{2}}$ and hence
$W_{\infty}^{1}=W_{\infty}^{2}$. ∎

We will use Lemma 22 to prove Theorem 4. But before we proceed, let us fix
some notation: we define the efficiency as a function of $\alpha\geq 0$ :

$\displaystyle\eta_{\alpha}^{J}(\rho_{\textup{Cold}}^{1})=\sup_{E_{k_{J}}-E_{j}>0}\eta(\rho_{\textup{Cold}}^{1})\quad\textup{
subject to }$ $\displaystyle\quad
F_{\alpha}(\rho^{0}_{\textup{W}}\otimes\tau_{\textup{Cold}}^{0},\tau_{\textup{ColdW}}^{h})\geq
F_{\alpha}(\rho^{J}_{\textup{W}}\otimes\rho_{\textup{Cold}}^{1},\tau_{\textup{ColdW}}^{h}),$
(379) and
$\displaystyle\quad\textup{tr}(\hat{H}_{t}\rho_{\textup{ColdHotMW}}^{0})=\textup{tr}(\hat{H}_{t}\rho_{\textup{ColdHotMW}}^{1,J}).$
(380)

with $J=1,2$ denoting the final battery state $\rho^{J}_{\textup{W}}$. We also
define an $\alpha$ independent efficiency:

$\displaystyle\eta^{J}(\rho_{\textup{Cold}}^{1})=\sup_{E_{k_{J}}-E_{j}>0}\eta(\rho_{\textup{Cold}}^{1})\quad\textup{
subject to }$ (381) $\displaystyle F_{\alpha}(\rho^{0}_{\textup{W}}\otimes$
$\displaystyle\tau_{\textup{Cold}}^{0},\tau_{\textup{ColdW}}^{h})\geq
F_{\alpha}(\rho^{J}_{\textup{W}}\otimes\rho_{\textup{Cold}}^{1},\tau_{\textup{ColdW}}^{h})\,\forall\alpha\geq
0.$ (382)

For any $\alpha\geq 0$, and any state $\rho_{\textup{Cold}}^{1}$,
$\eta_{\alpha}^{J}(\rho_{\textup{Cold}}^{1})\geq\eta^{J}(\rho_{\textup{Cold}}^{1})$
holds.

We already know that when $\Omega>1$, for any final cold bath state
$\rho_{\textup{Cold}}^{1}$, the efficiency
$\eta^{1}(\rho_{\textup{Cold}}^{1})$ is strictly less than the Carnot
efficiency. Theorem 4 shows that this is also true for
$\eta^{2}(\rho_{\textup{Cold}}^{1})$, i.e. when allowing a more general
battery final state.

######  Theorem 4.

[General battery states do not improve efficiency] Consider a heat engine with
a cold bath consisting of $n$ qubits, and consider the case where $\Omega>1$
(recall the definition of $\Omega$ in Eq. (245)). Then for any final cold bath
state $\rho_{\textup{Cold}}^{1}$, the efficiency
$\eta^{2}(\rho_{\textup{Cold}}^{1})$ is strictly less than the Carnot
efficiency.

###### Proof.

Firstly, suppose that $\Omega>1$. By Lemma 15 we know that the infimum is
obtained at $\alpha=\infty$, and by Lemma 16 we know that the efficiency for
quasi-static heat engine is strictly less than the Carnot value:

$\lim_{g\rightarrow 0}\eta^{1}(\tau_{\beta_{f}})=\lim_{g\rightarrow
0}\eta^{1}_{\infty}(\tau_{\beta_{f}})<\eta_{C}.$ (383)

On the other hand, we also know from Lemma 6 that
$\eta^{2}(\rho_{\textup{Cold}}^{1})$ can only possibly achieve Carnot
efficiency in the quasi-static limit. In other words, for all other final
states $\rho_{\textup{Cold}}^{1}$ we know that Carnot efficiency cannot be
achieved. Therefore, it suffices to see that in the quasi-static limit,

$\lim_{g\rightarrow 0}\eta^{2}(\tau_{\beta_{f}})\leq\lim_{g\rightarrow
0}\eta^{2}_{\infty}(\tau_{\beta_{f}})=\lim_{g\rightarrow
0}\eta^{1}_{\infty}(\tau_{\beta_{f}})=\lim_{g\rightarrow
0}\eta^{1}(\tau_{\beta_{f}})<\eta_{C}.$ (384)

The second equality is obtained by noting that for any state
$\tilde{\rho}_{\textup{Cold}}^{1}$ (and therefore for $\tau_{\beta_{f}}$):

  1. 1.

$\Delta C$ is the same for both expressions of efficiency
$\eta^{1}_{\infty}(\rho_{\textup{Cold}}^{1})$ and
$\eta^{2}_{\infty}(\rho_{\textup{Cold}}^{1})$.

  2. 2.

By Lemma 22, for all $0<\varepsilon<\left[1+e^{\beta_{h}(E_{\rm
max}-E_{j})}\right]^{-1}$,
$W_{\infty}^{1}(\tilde{\rho}_{\textup{Cold}}^{1})=W_{\infty}^{2}(\tilde{\rho}_{\textup{Cold}}^{1})$.

Hence, from Items 1 and 2, one concludes that
$\eta^{1}_{\infty}(\tilde{\rho}_{\textup{Cold}}^{1})=\eta^{2}_{\infty}(\tilde{\rho}_{\textup{Cold}}^{1})$.
The third equality in Eq. (384) comes directly from Eq. (383). ∎