# A Finite Test for the Linearizability of Two-Input Systems by a Two-
Dimensional Endogenous Dynamic Feedback

Conrad Gstöttner1 Bernd Kolar1 Markus Schöberl1 *The first author and the
second author have been supported by the Austrian Science Fund (FWF) under
grant number P 32151 and P 29964.1Institute of Automatic Control and Control
Systems Technology, Johannes Kepler University, Linz, Austria
{conrad.gstoettner,bernd.kolar,
<EMAIL_ADDRESS>
###### Abstract

We propose an algorithmic test to check whether a two-input system is
linearizable by an endogenous dynamic feedback with a dimension of at most
two. This test furthermore provides a procedure for systematically deriving
flat outputs for this class of systems.

##  I INTRODUCTION

Flatness was introduced in control theory by Fliess, Lévine, Martin and
Rouchon, see e. g. [1, 2]. The flatness property allows an elegant systematic
solution of feed-forward and feedback problems, see e. g. [2]. Roughly
speaking, a nonlinear control system

$\displaystyle\begin{aligned} \dot{x}&=f(x,u)\end{aligned}$

with $\mathrm{dim}(x)=n$ states and $\mathrm{dim}(u)=m$ inputs is flat, if
there exist $m$ differentially independent functions
$y^{j}=\varphi^{j}(x,u,u_{1},\ldots,u_{q})$, $u_{k}$ denoting the $k$-th time
derivative of $u$, such that $x$ and $u$ can locally be parameterized by $y$
and its time derivatives. For this parameterization we write

$\displaystyle\begin{aligned} x&=F_{x}(y,y_{1},\ldots,y_{r-1})\\\
u&=F_{u}(y,y_{1},\ldots,y_{r})\end{aligned}$

and refer to it as the parameterizing map with respect to the flat output $y$.
For a given flat output, $F_{x}$ and $F_{u}$ are unique. If the parameterizing
map is invertible, i. e. $y$ and all the time derivatives of $y$ present in
the map can be expressed solely as functions of $x$ and $u$, the system is
static feedback linearizable. In this case we call $y$ a linearizing output of
the static feedback linearizable system. The static feedback linearization
problem has been completely solved, see [3, 4]. However, there do not exist
verifiable necessary and sufficient conditions for flatness, except for
certain classes of systems, including driftless systems with two inputs [5]
and systems which are linearizable by a one-fold prolongation of a suitably
chosen control [6]. Necessary and sufficient conditions for $(x,u)$-flatness
of two-input control affine systems with four states can be found in [7].

It is well known that every flat system can be rendered static feedback
linearizable by an endogenous dynamic feedback. If a flat output is known,
such a linearizing feedback can be constructed systematically, see e. g. [8].
In this contribution we propose an algorithmic test to check whether a two
input system is linearizable by an endogenous dynamic feedback with a
dimension of at most two. This test furthermore provides a procedure for
systematically deriving flat outputs of such systems. The main idea of our
algorithmic test is to successively split off or add endogenous dynamic
feedbacks to the system in order to obtain a static feedback linearizable
system. If the algorithm does not yield a static feedback linearizable system,
i. e. if the test fails, the considered system is not linearizable by an
endogenous dynamic feedback with a dimension of at most two. The linearizing
outputs of the eventually obtained static feedback linearizable system are
then flat outputs of the original system. In the literature, constructive
approaches for deriving flat outputs can be found, see e. g. [9], [10] or
[11]. All these procedures are applicable to a broad variety of systems and
are not restricted to the two-input case. However, there usually occur degrees
of freedom, which may lead to an infinite number of branches. In contrast to
that, in every step of our procedure, a certain action has to be applied
depending on easily verifiable properties of the system. We will see that
systems which are linearizable by an endogenous dynamic feedback with a
dimension of at most two are actually linearizable by a special subclass of
endogenous dynamic feedbacks, namely prolongations of a suitably chosen
control. In [12], a complete solution for the flatness problem of two-input
systems linearizable via a one-fold prolongation of a suitably chosen control
is provided. In [13], two-input systems linearizable via a two-fold
prolongation of a suitably chosen control are considered. The present
contribution is greatly influenced by these results. However, in [13] a
complete solution for the flatness problem of two-input systems linearizable
via a two-fold prolongation is not provided, due to Assumption 2 therein. Our
procedure also applies to the cases which are not covered by the results in
[13]. In Section V, we apply our procedure to three examples, to none of which
the results in [13] are applicable. In [14] we proposed a structurally flat
triangular form based on the extended chained form. Our first two examples
turn out to be static feedback equivalent to the proposed triangular form,
which provides another way for systematically deriving flat outputs for our
first two examples. However, the results in [14] are not applicable to our
third example.

##  II NOTATION

Let $\mathcal{X}$ be an $n$-dimensional smooth manifold, equipped with local
coordinates $x^{i}$, $i=1,\ldots,n$. Its tangent bundle is denoted by
$(\mathcal{T}(\mathcal{X}),\tau_{\mathcal{X}},\mathcal{X})$, for which we have
the induced local coordinates $(x^{i},\dot{x}^{i})$ with respect to the basis
$\\{\partial_{x^{i}}\\}$. We make use of the Einstein summation convention. By
$\partial_{x}h$ we denote the $m\times n$ Jacobian matrix of
$h=(h^{1},\ldots,h^{m})$ with respect to $x=(x^{1},\ldots,x^{n})$. The
$k$-fold Lie derivative of a function $\varphi$ along a vector field $v$ is
denoted by $\mathrm{L}_{v}^{k}\varphi$. Let $v$ and $w$ be two vector fields.
Their Lie bracket is denoted by $[v,w]$. Let furthermore $D_{1}$ and $D_{2}$
be two distributions. By $[v,D_{1}]$ we denote the distribution spanned by the
Lie bracket of $v$ with all basis vector fields of $D_{1}$, and by
$[D_{1},D_{2}]$ the distribution spanned by the Lie brackets of all possible
pairs of basis vector fields of $D_{1}$ and $D_{2}$. The $i$-th derived flag
of a distribution $D$ is denoted by $D^{(i)}$ and defined by $D^{(0)}=D$ and
$D^{(i+1)}=D^{(i)}+[D^{(i)},D^{(i)}]$ for $i\geq 0$. We denote the Cauchy
characteristic distribution of $D$ by $\mathcal{C}(D)$. It is spanned by all
vector fields $c$ which belong to $D$ and satisfy $[c,D]\subset D$. We make
use of multi-indices, in particular, by $R=(r_{1},r_{2})$ we denote the unique
multi-index associated to a flat output of a system with two inputs, where
$r_{j}$ denotes the order of the highest derivative of $y^{j}$ needed to
parameterize $x$ and $u$ by this flat output, i. e.
$y_{[R]}=(y^{1},y^{1}_{1},\ldots,y^{1}_{r_{1}},y^{2},y^{2}_{1},\ldots,y^{2}_{r_{2}})$.
Furthermore, we define $R\pm c=(r_{1}\pm c,r_{2}\pm c)$ with an integer $c$,
and $\\#R=r_{1}+r_{2}$.

##  III PRELIMINARIES

Throughout, we assume all distributions to have locally constant dimension, we
consider generic points only. Consider a nonlinear two-input system of the
form

$\displaystyle\begin{aligned} \dot{x}&=f(x,u)\,,\end{aligned}$ (1)

with $\mathrm{dim}(x)=n$, $\mathrm{dim}(u)=2$ and
$\mathrm{rank}(\partial_{u}f)=2$. Let
$y^{j}=\varphi^{j}(x,u,u_{1},\ldots,u_{q})$, $j=1,2$ be a flat output of (1)
and

$\displaystyle\begin{aligned} x&=F_{x}(y_{[R-1]})\\\
u&=F_{u}(y_{[R]})\end{aligned}$ (2)

the parameterizing map with respect to this flat output. The parameterizing
map (2) is a submersion. We denote the difference of the dimensions of the
domain and the codomain of (2) by $d$, i. e. $d=\\#R+2-(n+2)=\\#R-n$. In [12]
and [13], the number $\\#R+2$ is called the differential weight of the flat
output. (The differential weight of a flat output with difference $d$ is thus
given by $d+n+2$.) The difference $d$ is the minimal dimension of an
endogenous dynamic feedback needed to render (1) static feedback linearizable
such that $y$ forms a linearizing output of the closed loop system. The
endogenous feedback can be constructed systematically, see e. g. [8]. If
$d=0$, the map (2) degenerates to a diffeomorphism and the system is static
feedback linearizable with $y$ being a linearizing output. A flat output $y$
is called a minimal flat output if the difference $d$ is minimal compared to
all other possible flat outputs of the system. We define the difference $d$ of
a flat system to be the difference of a minimal flat output. The difference
$d$ of a system (1) therefore measures its distance from static feedback
linearizability, i. e. $d$ is the minimal possible dimension of an endogenous
dynamic feedback needed to render the system static feedback linearizable. For
(1), we define the distributions
$D_{0}=\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}}\\}$ and
$D_{i}=D_{i-1}+[f,D_{i-1}]$, $i\geq 1$ on the state and input manifold
$\mathcal{X}\times\mathcal{U}$, where $f=f^{i}(x,u)\partial_{x^{i}}$.

###### Theorem 1

The two-input system (1) is linearizable by static feedback if and only if all
the distributions $D_{i}$ are involutive and $\mathrm{dim}(D_{n-1})=n+2$.

For a proof of this theorem, we refer to [4]. The involutivity of $D_{i-1}$
implies the invariance of $D_{i}$ with respect to regular input
transformations $\bar{u}=\Phi_{u}(x,u)$. In particular, since
$D_{0}=\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}}\\}$ is always
involutive, $D_{1}=D_{0}+[f,D_{0}]$ is always feedback invariant, also when
the system is not static feedback linearizable. A special case of (1) are
affine input systems (AI-systems)

$\displaystyle\begin{aligned}
\dot{x}&=a(x)+b_{1}(x)u^{1}+b_{2}(x)u^{2}\,.\end{aligned}$ (3)

###### Lemma 1

The system (1) allows an AI representation (3) if and only if
$D_{0}\subset\mathcal{C}(D_{1})$.

It follows from the proof of this lemma, which can be found in Section VII,
that the input transformation $\bar{u}^{j}=f^{j}(x,u)$, $j=1,2$, which may
require a renumbering of the state variables, yields an AI representation if
an AI representation is indeed possible. For AI-systems, we can omit
$\partial_{u^{j}}$ in $D_{1}$, i. e. for AI-systems we define
$D_{1}=\mathrm{span}\\{b_{1},b_{2}\\}$, where
$b_{j}=b_{j}^{i}(x)\partial_{x^{i}}$ are vector fields on the state manifold
$\mathcal{X}$. A central role in this contribution plays the more general form

$\displaystyle\begin{aligned}
\dot{x}&=a(x,u^{1})+b(x,u^{1})u^{2}\,,\end{aligned}$ (4)

in which at least the input $u^{2}$ occurs affine. This form was introduced in
[9] and we refer to it as partial affine input form (PAI form)111In [15] it is
shown that the existence of a PAI representation is actually a necessary
condition for flatness and this necessary condition is not restricted to
systems with two inputs. See also [16]..

###### Lemma 2

The system (1) allows a PAI representation (4) if and only if the condition

$\displaystyle\begin{aligned}
(\alpha^{1})^{2}\partial_{u^{1}}^{2}f+2\alpha^{1}\alpha^{2}\partial_{u^{1}}\partial_{u^{2}}f+(\alpha^{2})^{2}\partial_{u^{2}}^{2}f&\overset{!}{\in}D_{1}\,,\end{aligned}$
(5)

admits nontrivial solutions $\alpha^{j}(x,u)$ and at least one of the vector
fields $v_{c}=\alpha^{1}\partial_{u^{1}}+\alpha^{2}\partial_{u^{2}}$ obtained
from a solution of (5) meets $v_{c}\in\mathcal{C}(\Delta_{1})$ where
$\Delta_{1}=\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}},[v_{c},f]\\}$.

A proof of this lemma is provided in Section VII, where we also discuss a
crucial property of the condition (5), namely that for systems which do not
allow an AI representation, i. e. $D_{0}\not\subset\mathcal{C}(D_{1})$, the
condition (5) admits at most two independent nontrivial solutions and thus, it
yields at most two non-collinear candidates for a vector field $v_{c}$. To
each vector field $v_{c}$ which actually fulfills
$v_{c}\in\mathcal{C}(\Delta_{1})$, there exists a corresponding PAI
representation. To derive it, apply the following two steps:

  1. 1.

Straighten out the distribution $\mathrm{span}\\{v_{c}\\}$, i. e. apply a
suitable input transformation $\bar{u}=\Phi_{u}(x,u)$ such that
$\mathrm{span}\\{v_{c}\\}=\mathrm{span}\\{\partial_{\bar{u}^{2}}\\}$.

  2. 2.

Normalize one equation which explicitly depends on $\bar{u}^{2}$ (i. e.
introduce e. g. $\tilde{u}^{2}=\bar{f}^{1}(x,\bar{u}^{1},\bar{u}^{2})$, which
may require a renumbering of the states of the system).

Based on the fact that at most two non-collinear vector fields $v_{c}$ exist
and that to every PAI representation of the system, a vector field $v_{c}$
belongs, it can be shown that there exist at most two fundamentally different
PAI representations of the system, all others are equivalent to one of those
two by a transformation of the form

$\displaystyle\begin{aligned} \bar{x}&=\Phi_{x}(x)\\\
\bar{u}^{1}&=g^{1}(x,u^{1})\\\
\bar{u}^{2}&=g^{2}(x,u^{1})+m(x,u^{1})u^{2}\,.\end{aligned}$ (6)

##  IV FINITE ALGORITHMIC TEST

The main idea of our test for flatness with $d\leq 2$ is to successively split
off or add endogenous dynamic feedbacks to the system, in order to eventually
obtain a static feedback linearizable system. The linearizing outputs of the
eventually obtained static feedback linearizable system are then flat outputs
with $d\leq 2$ of the original system. As we will see in the next section, for
systems with $d\leq 2$, there always applies one of the following cases.

  1. 1.

The system allows an AI representation (3) with an involutive input
distribution $D_{1}=\mathrm{span}\\{b_{1},b_{2}\\}$. In this case, straighten
out $D_{1}$ by a suitable state transformation $\bar{x}=\Phi_{x}(x)$, in order
to obtain a decomposition of (3) into the form

$\displaystyle\begin{aligned}
\Sigma_{1}:&&\dot{\bar{x}}_{1}^{i_{1}}&=\bar{f}_{1}^{i_{1}}(\bar{x}_{1},\bar{x}_{2})\,,&i_{1}&=1,\ldots,n-2\\\
\Sigma_{2}:&&\dot{\bar{x}}_{2}^{i_{2}}&=\bar{f}_{2}^{i_{1}}(\bar{x}_{1},\bar{x}_{2},u)\,,&i_{2}&=1,2\,.\end{aligned}$
(7)

Continue the analysis with the subsystem $\Sigma_{1}$ with the state
$\bar{x}_{1}$ and the input $\bar{x}_{2}$.

  2. 2.

The system allows an AI representation (3) but the input distribution
$D_{1}=\mathrm{span}\\{b_{1},b_{2}\\}$ is non-involutive. Two subcases are
possible. If $\mathrm{dim}(\overline{D}_{1})=3$, construct a nontrivial linear
combination $b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ such that
$[a,b_{c}]\in\overline{D}_{1}$. If $\mathrm{dim}(\overline{D}_{1})=4$,
construct a nontrivial linear combination
$b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ such that
$b_{c}\in\mathcal{C}(D_{1}^{(1)})$. In either subcase, introduce the new input
$\bar{u}^{1}=\alpha^{2}u^{1}-\alpha^{1}u^{2}$ and one-fold prolong it, i. e.
add the equation $\dot{\bar{u}}^{1}=\bar{u}^{1}_{1}$. Continue the analysis
with the prolonged system

$\displaystyle\begin{aligned}
\dot{x}&=a(x)+\bar{b}_{1}(x)\bar{u}^{1}+\bar{b}_{2}(x)\bar{u}^{2}\,,\end{aligned}$
$\displaystyle\begin{aligned} \dot{\bar{u}}^{1}&=\bar{u}^{1}_{1}\end{aligned}$
(8)

with the state $(x,\bar{u}^{1})$ and the input
$(\bar{u}^{1}_{1},\bar{u}^{2})$.

  3. 3.

The system only allows a PAI representation (4). Transform the system into PAI
form and one-fold prolong the non-affine occurring input $\bar{u}^{1}$, i. e.
add the equation $\dot{\bar{u}}^{1}=\bar{u}^{1}_{1}$. Continue the analysis
with the prolonged system

$\displaystyle\begin{aligned}
\dot{x}&=a(x,\bar{u}^{1})+b(x,\bar{u}^{1})\bar{u}^{2}\,,\end{aligned}$
$\displaystyle\begin{aligned} \dot{\bar{u}}^{1}&=\bar{u}^{1}_{1}\end{aligned}$
(9)

with the state $(x,\bar{u}^{1})$ and the input
$(\bar{u}^{1}_{1},\bar{u}^{2})$. If two non-equivalent PAI representation
exist, we have to continue with both of them (branching point).

The procedure succeeds if a static feedback linearizable system is obtained.
In the next section, we explain in detail the effects of the individual steps
which are applied depending on which case applies. As we will see, the
subsystem $\Sigma_{1}$ in (7) which is obtained if case 1 applies, has the
same difference $d$ as the complete system. Such a step therefore only reduces
the dimension of the state, but not the distance from static feedback
linearizability. In contrast to that, if for a system with $d\leq 2$ the cases
2 or 3 apply, we obtain a prolonged system (8) or (9) with a difference of
$d-1$ compared to the non-prolonged system. If case 3 applies, we may
encounter a branching point since there may exist two non-equivalent PAI
representations which both have to be considered. Nevertheless, at least one
of the AI-systems obtained by prolonging the non-affine input has a difference
of $d-1$ compared to the non-prolonged system. Therefore, applying the
procedure to a system with $d\leq 2$, after at most two steps which actually
decrease the difference (cases 2 and 3), we must obtain a static feedback
linearizable system. Otherwise, the original system has a difference of $d\geq
3$ or is not flat.

###  IV-A Necessity and sufficiency for systems with $d\leq 2$

In the following, we explain why applying the procedure to a system with a
difference of $d\leq 2$ necessarily yields a static feedback linearizable
system, and why the linearizing outputs of the obtained static feedback
linearizable system are flat outputs with $d\leq 2$ of the original system. We
have the following three results, which are proven in Section VII.

###### Lemma 3

Every flat output with difference $d$ of the complete system (7) is also a
flat output with the same difference $d$ for the subsystem $\Sigma_{1}$ in
(7). Conversely, every flat output of $\Sigma_{1}$ is also a flat output of
the complete system and the differences again coincide.

###### Lemma 4

For an AI-system with a difference of $d\leq 2$ which has a non-involutive
input distributions $D_{1}=\mathrm{span}\\{b_{1},b_{2}\\}$, we either have
$\mathrm{dim}(\overline{D}_{1})=3$ and there exists a non-trivial vector field
$b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ satisfying
$[a,b_{c}]\in\overline{D}_{1}$, or we have $\mathrm{dim}(\overline{D}_{1})=4$
and there exists a non-trivial vector field
$b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ satisfying
$b_{c}\in\mathcal{C}(D_{1}^{(1)})$. Furthermore, the prolonged system (8) has
a difference of $d-1$ (i. e. it is either static feedback linearizable or has
a difference of $1$). The linearizing outputs or flat outputs with difference
$1$ of the prolonged system are flat outputs with difference $1$ or $2$ of the
original system.

###### Lemma 5

A system with $d\leq 2$ which does not allow an AI representation allows a PAI
representation. There exist at most two non-equivalent PAI representations and
at least one of the prolonged systems (9) obtained from these possibly two PAI
representations has a difference of $d-1$ (i. e. it is either static feedback
linearizable or has a difference of $1$). The linearizing outputs or flat
outputs with difference $1$ of the prolonged system are flat outputs with
difference $1$ or $2$ of the original system.

Based on these results, we first explain why the procedure necessarily
succeeds when it is applied to a system with a difference of $d\leq 2$. Assume
we have a system (1) with $d\leq 2$. If the system allows an AI
representation, its input distribution $D_{1}=\mathrm{span}\\{b_{1},b_{2}\\}$
can either be involutive or not. If it is involutive, case 1 applies, which
according to Lemma 3 yields a system with the same difference $d$ as the
original system. We can exclude the case $\mathrm{dim}(D_{1})\leq 1$, since
then we effectively have a single input system which is either static feedback
linearizable or not flat at all, see [17]. Furthermore, we can exclude the
case that the subsystem $\Sigma_{1}$ in (7) has redundant inputs222If
$\Sigma_{1}$ has redundant inputs, i. e. if $\Sigma_{1}$ effectively reduces
to a single input system, it can either be static feedback linearizable or not
flat at all. However, if $\Sigma_{1}$ is static feedback linearizable, then
also the complete system is static feedback linearizable and the procedure
already succeeded before this decomposition. On the other hand, if
$\Sigma_{1}$ is not flat then due to Lemma 3, also (7) is not flat. The system
(7) need not be the original system, it may be the outcome of previously
applied decompositions or prolongations according to the cases 1 to 3.
Nevertheless, due to the Lemmas 3 to 5, if $\Sigma_{1}$ is not flat, also the
original system cannot be flat with a difference of $d\leq 2$. In fact, it
even follows that the original system cannot be flat at all.. If the procedure
succeeds, it necessarily terminates with case 2 or 3 since case 1 always
yields a subsystem $\Sigma_{1}$ with the same difference as the complete
system (7).

If the system allows an AI representation with a non-involutive input
distribution $D_{1}$, case 2 applies and due to Lemma 4, the prolonged system
(8) has a difference of $d\leq 1$.

If the system does not allow an AI representation, according to Lemma 5 it
necessarily allows a PAI representation and thus, case 3 applies. Since there
exist at most two non-equivalent PAI representations and in case 3 we require
to consider both of them, one of the possibly two prolonged systems (9) has a
difference of $d\leq 1$.

Each step therefore either yields a system with a lower state dimension or a
lower difference (for at least one of the possibly two branches if case 3
applies). Since by assumption the original system has a difference of $d\leq
2$, after at most two steps which actually yield a system with a lower
difference (cases 2 and 3) we necessarily have a static feedback linearizable
system.

On the other hand, due to Lemma 3 to 5 it follows that if the procedure
succeeds, i. e. if after at most two steps which actually yield a system with
a lower difference (cases 2 and 3) we obtain a static feedback linearizable
system, the linearizing outputs of this static feedback linearizable system
are flat outputs with $d=1$ or $d=2$ of the original system.

##  V EXAMPLES

As already mentioned in the introduction, our first two examples can also be
handled with the structurally flat triangular form which we proposed in [14].
However, neither [14] nor [13], nor both of them together completely cover the
class of systems addressed in this paper. The results in [14] and [13] do not
apply to our third example.

###  V-A Planar VTOL aircraft

Consider the planar VTOL aircraft, also treated in e. g. [8], [18], [19] or
[14] and given by

$\displaystyle\begin{aligned}
\dot{x}&=v_{x}&&&\dot{v}_{x}&=\epsilon\cos(\theta)u^{2}-\sin(\theta)u^{1}\\\
\dot{z}&=v_{z}&&&\dot{v}_{z}&=\cos(\theta)u^{1}+\epsilon\sin(\theta)u^{2}-1\\\
\dot{\theta}&=\omega&&&\dot{\omega}&=u^{2}\,.\end{aligned}$ (10)

The input vector fields of this system are given by
$b_{1}=-\sin(\theta)\partial_{v_{x}}+\cos(\theta)\partial_{v_{z}}$ and
$b_{2}=\epsilon\cos(\theta)\partial_{v_{x}}+\epsilon\sin(\theta)\partial_{v_{z}}+\partial_{\omega}$,
which span the involutive input distribution
$D_{1}=\mathrm{span}\\{b_{1},b_{2}\\}$. Therefore, case 1 of the procedure
applies. In order to decompose the system as described in case 1, we
straighten out $D_{1}$ by applying the state transformation
$\bar{v}_{x}=\cos(\theta)v_{x}+\sin(\theta)v_{z}-\epsilon\omega$, with the
rest of the coordinates left unchanged. This results in the decomposition

$\displaystyle\begin{aligned} &\Sigma_{1}:\quad\begin{aligned}
\dot{x}&=\tfrac{1}{\cos(\theta)}(\bar{v}_{x}-\sin(\theta)v_{z}+\epsilon\omega)\\\
\dot{z}&=v_{z}\\\ \dot{\theta}&=\omega\\\
\dot{\bar{v}}_{x}&=\tfrac{\omega}{\cos(\theta)}(v_{z}-\sin(\theta)(\bar{v}_{x}+\epsilon\omega))-\sin(\theta)\\\
\end{aligned}\\\\[2.15277pt] &\Sigma_{2}:\quad\begin{aligned}
\dot{v}_{z}&=\cos(\theta)u^{1}+\epsilon\sin(\theta)u^{2}-1\\\
\dot{\omega}&=u^{2}\,.\end{aligned}\end{aligned}$ (11)

We proceed with the subsystem $\Sigma_{1}$, for which $v_{z}$ and $\omega$ act
as inputs. This system does not allow an AI representation, so we are in case
3 of the procedure. In order to derive a PAI representation for $\Sigma_{1}$,
we first solve (5), which yields
$\alpha^{2}(\alpha^{1}-\alpha^{2}\epsilon\sin(\theta))\overset{!}{=}0$, and
has the two independent solutions $\alpha^{1}=\lambda$, $\alpha^{2}=0$ and
$\alpha^{1}=\lambda\epsilon\sin(\theta)$, $\alpha^{2}=\lambda$, with an
arbitrary function $\lambda\neq 0$ (we can choose e. g. $\lambda=1$). Both
vector fields $v_{c,1}=\partial_{v_{z}}$ and
$v_{c,2}=\epsilon\sin(\theta)\partial_{v_{z}}+\partial_{\omega}$ obtained from
these solutions fulfill
$v_{c,i}\in\mathcal{C}(\mathrm{span}\\{\partial_{v_{z}},\partial_{\omega},[v_{c,i},f]\\})$,
i. e. $\Sigma_{1}$ allows two non-equivalent PAI representations. The vector
field $v_{c,1}=\partial_{v_{z}}$ is already straightened out and $\Sigma_{1}$
in (11) is actually already in PAI form with the input $v_{z}$ occurring
affine. However, the AI-system obtained by one-fold prolonging the non-affine
input $\omega$ cannot have a difference of $d=1$, since none of the cases of
the procedure applies to this AI-system. Straightening out
$v_{c,2}=\epsilon\sin(\theta)\partial_{v_{z}}+\partial_{\omega}$, by the input
transformation $\bar{v}_{z}=v_{z}-\epsilon\sin(\theta)\omega$, results in
$v_{c,2}=\partial_{\omega}$ and the PAI representation333Note that in these
coordinates in fact both vector fields $v_{c,1}$ and $v_{c,2}$ are
straightened out simultaneously and that in (12) actually either input could
be interpreted as the affine entering input.

$\displaystyle\begin{aligned}
\dot{x}&=\tfrac{1}{\cos(\theta)}(\bar{v}_{x}-\sin(\theta)\bar{v}_{z})+\epsilon\cos(\theta)\omega\\\
\dot{z}&=\bar{v}_{z}+\epsilon\sin(\theta)\omega\\\ \dot{\theta}&=\omega\\\
\dot{\bar{v}}_{x}&=\tfrac{\omega}{\cos(\theta)}(\bar{v}_{z}-\sin(\theta)\bar{v}_{x})-\sin(\theta)\,.\end{aligned}$
(12)

For the AI-system obtained by one-fold prolonging $\bar{v}_{z}$, i. e. adding
the equation $\dot{\bar{v}}_{z}=\bar{v}_{z,1}$, we have
$\mathrm{dim}(\overline{D})=3$, where $D_{1}=\mathrm{span}\\{b_{1},b_{2}\\}$
is the input distribution spanned by its input vector fields
$b_{1}=\partial_{\bar{v}_{z}}$ and
$b_{2}=\epsilon\cos(\theta)\partial_{x}+\epsilon\sin(\theta)\partial_{z}+\partial_{\theta}+\tfrac{\bar{v}_{z}-\sin(\theta)\bar{v}_{x}}{\cos(\theta)}\partial_{\bar{v}_{x}}$.
The drift of this AI-system is given by
$a=\tfrac{1}{\cos(\theta)}(\bar{v}_{x}-\sin(\theta)\bar{v}_{z})\partial_{x}+\bar{v}_{z}\partial_{z}-\sin(\theta)\partial_{\bar{v}_{x}}$
and there indeed exists a linear combination
$b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ of the input vector fields which
satisfies $[a,b_{c}]\in\overline{D}_{1}$, namely $b_{c}=\lambda b_{2}$ with an
arbitrary function $\lambda\neq 0$ (e. g. $\lambda=1$ and thus $\alpha^{1}=0$
and $\alpha^{2}=1$). So the conditions for case 2 are met and we prolong the
input $\tilde{u}^{1}=\alpha^{2}\bar{v}_{z,1}-\alpha^{1}\omega=\bar{v}_{z,1}$,
i. e. we simply have to prolong the input $\bar{v}_{z,1}$, no actual input
transformation is required in this case, and we obtain the static feedback
linearizable system

$\displaystyle\begin{aligned}
\dot{x}&=\tfrac{1}{\cos(\theta)}(\bar{v}_{x}-\sin(\theta)\bar{v}_{z})+\epsilon\cos(\theta)\omega\\\
\dot{z}&=\bar{v}_{z}+\epsilon\sin(\theta)\omega\\\ \dot{\theta}&=\omega\\\
\dot{\bar{v}}_{x}&=\tfrac{\omega}{\cos(\theta)}(\bar{v}_{z}-\sin(\theta)\bar{v}_{x})-\sin(\theta)\end{aligned}$
$\displaystyle\begin{aligned} \dot{\bar{v}}_{z}&=\bar{v}_{z,1}\\\
\dot{\bar{v}}_{z,1}&=\bar{v}_{z,2}\,,\end{aligned}$

with the input $(\omega,\bar{v}_{z,2})$. We have
$y=(x-\epsilon\sin(\theta),z+\epsilon\cos(\theta))$ as a possible linearizing
output of this system, which in turn is a flat output with a difference of
$d=2$ of the original system (10).

###  V-B Academic example I

Consider the system

$\displaystyle\begin{aligned}
\dot{x}^{1}&=u^{1}\,,&\dot{x}^{2}&=u^{2}\,,&\dot{x}^{3}&=\sin(\tfrac{u^{1}}{u^{2}})\,,\end{aligned}$
(13)

also considered in [20], [21] and [14]. This system does not allow an AI
representation, so we are in case 3. In order to derive a PAI representation,
we solve (5), which for this system yields

$\displaystyle\begin{aligned}
&(\alpha^{1})^{2}\sin(\tfrac{u^{1}}{u^{2}})(u^{2})^{2}+2\alpha^{1}\alpha^{2}(\cos(\tfrac{u^{1}}{u^{2}})u^{2}-\sin(\tfrac{u^{1}}{u^{2}})u^{1})u^{2}+\\\
&\hskip
60.00009pt(\alpha^{2})^{2}(\sin(\tfrac{u^{1}}{u^{2}})u^{1}-2\cos(\tfrac{u^{1}}{u^{2}})u^{2})u^{1}\overset{!}{=}0\end{aligned}$

and has the two independent non-trivial solutions $\alpha^{1}=\lambda u^{1}$,
$\alpha^{2}=\lambda u^{2}$ and
$\alpha^{1}=\lambda(u^{1}\tan(\tfrac{u^{1}}{u^{2}})-2u^{2})$,
$\alpha^{2}=\lambda u^{2}\tan(\tfrac{u^{1}}{u^{2}})$ with an arbitrary
function $\lambda\neq 0$ (e. g. $\lambda=1$). Only the vector field
$v_{c}=u^{1}\partial_{u^{1}}+u^{2}\partial_{u^{2}}$, which is obtained from
the first solution, meets
$v_{c}\in\mathcal{C}(\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}},[v_{c},f]\\})$.
To derive the corresponding PAI representation, we straighten out
$\mathrm{span}\\{v_{c}\\}$ by the input transformation
$\bar{u}^{1}=\tfrac{u^{1}}{u^{2}}$, $\bar{u}^{2}=u^{2}$, resulting in

$\displaystyle\begin{aligned}
\dot{x}^{1}&=\bar{u}^{1}\bar{u}^{2}\,,&\dot{x}^{2}&=\bar{u}^{2}\,,&\dot{x}^{3}&=\sin(\bar{u}^{1})\end{aligned}$

and proceed with the AI-system obtained by one-fold prolonging $\bar{u}^{1}$,
i. e. adding the equation $\dot{\bar{u}}^{1}=\bar{u}^{1}_{1}$. Its input
vector fields are given by $b_{1}=\partial_{\bar{u}^{1}}$,
$b_{2}=\bar{u}^{1}\partial_{x^{1}}+\partial_{x^{2}}$ and its drift reads
$a=\sin(\bar{u}^{1})\partial_{x^{3}}$. We have
$\mathrm{dim}(\overline{D}_{1})=3$ and there exists a linear combination
$b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ which meets
$[a,b_{c}]\in\overline{D}_{1}$, namely $b_{c}=\lambda b_{2}$, with an
arbitrary function $\lambda\neq 0$ (e. g. $\lambda=1$ and thus $\alpha^{1}=0$
and $\alpha^{2}=1$). So the conditions for case 2 are met and we prolong the
input
$\tilde{u}^{1}=\alpha^{2}\bar{u}^{1}_{1}-\alpha^{1}\bar{u}^{2}=\bar{u}^{1}_{1}$,
i. e. we simply have to prolong the input $\bar{u}^{1}_{1}$, no actual input
transformation is required, and we obtain the static feedback linearizable
system

$\displaystyle\begin{aligned} \dot{x}^{1}&=\bar{u}^{1}\bar{u}^{2}\\\
\dot{x}^{2}&=\bar{u}^{2}\\\ \dot{x}^{3}&=\sin(\bar{u}^{1})\\\ \end{aligned}$
$\displaystyle\begin{aligned} \dot{\bar{u}}^{1}&=\bar{u}^{1}_{1}\\\
\dot{\bar{u}}^{1}_{1}&=\bar{u}^{1}_{2}\end{aligned}$ (14)

with the input $(\bar{u}^{1}_{2},\bar{u}^{2})$. We have
$y=(x^{3},x^{1}-x^{2}\bar{u}^{1})=(x^{3},x^{1}-x^{2}\tfrac{u^{1}}{u^{2}})$ as
a possible linearizing output of (14), which in turn is a flat output with a
difference of $d=2$ of the original system (13).

###  V-C Academic example II

As a final example, consider the system

$\displaystyle\begin{aligned}
\dot{x}^{1}&=\arcsin(\tfrac{u^{1}+u^{2}}{x^{2}})-x^{4}\\\
\dot{x}^{2}&=x^{4}\end{aligned}$ $\displaystyle\begin{aligned}
\dot{x}^{3}&=u^{1}\\\ \dot{x}^{4}&=u^{2}\,.\end{aligned}$

With the algorithmic test proposed in this paper, it can be shown that this
system has a difference of $d=2$ and $y=(x^{1}+x^{2},x^{3}+x^{4})$ follows as
a corresponding flat output with $d=2$. In the first step case 3 applies, in
the second step case 1 applies, and then case 3 applies again.

##  VI CONCLUSIONS

We have proposed a finite algorithmic test for linearizability by an
endogenous dynamic feedback with a dimension of at most two, together with a
procedure for systematically deriving flat outputs for such systems. The
conditions are in principle verifiable, however, applying the test requires
straightening out involutive distributions, which from a computational point
of view is unfavorable. Further research will be devoted to deriving necessary
and sufficient conditions based on the proposed test which overcome this
drawback.

##  VII PROOFS

### Proof of Lemma 1

Necessity. For an AI-system (3), we have
$D_{1}=\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}},b_{1},b_{2}\\}$ and
since the vector fields $b_{j}$ do not depend on $u^{1}$ and $u^{2}$, we
obviously have
$D_{0}=\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}}\\}\subset\mathcal{C}(D_{1})$.

Sufficiency. Given a system $\dot{x}=f(x,u)$ with
$\mathrm{rank}(\partial_{u}f)=2$, we can always apply an input transformation
$\bar{u}^{j}=f^{j}(x,u)$, $j=1,2$ (this transformation may require a
renumbering of the sate variables) in order to obtain

$\displaystyle\begin{aligned} \dot{x}^{j}&=\bar{u}^{j}\,,&j&=1,2\\\
\dot{x}^{i}&=\bar{f}^{i}(x,\bar{u})\,,&i&=3,\ldots,n\,.\end{aligned}$ (15)

Calculating $D_{1}$ in these coordinates yields
$D_{1}=\mathrm{span}\\{\partial_{\bar{u}^{1}},\partial_{\bar{u}^{2}},v_{1},v_{2}\\}$
with
$v_{1}=\partial_{x^{1}}+\partial_{\bar{u}^{1}}\bar{f}^{i}\partial_{x^{i}}$,
$i=3,\ldots,n$ and
$v_{2}=\partial_{x^{2}}+\partial_{\bar{u}^{2}}\bar{f}^{i}\partial_{x^{i}}$,
$i=3,\ldots,n$. The condition $D_{0}\subset\mathcal{C}(D_{1})$ implies
$[\partial_{\bar{u}^{j}},v_{k}]\in D_{1}$, $j,k\in\\{1,2\\}$, which can only
hold if $[\partial_{\bar{u}^{j}},v_{k}]=0$, implying that (15) is actually in
AI form. ∎

### Proof of Lemma 2

For a system in PAI form (4), we have
$D_{1}=\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}},(\partial_{u^{1}}a^{i}+u^{2}\partial_{u^{1}}b^{i})\partial_{x^{i}},b^{i}\partial_{x^{i}}\\}$.
Consider the subdistribution
$\Delta_{1}=\mathrm{span}\\{\partial_{u^{1}},\partial_{u^{2}},b^{i}\partial_{x^{i}}\\}$
of $D_{1}$. We have $\partial_{u^{2}}\in\mathcal{C}(\Delta_{1})$ and
$\Delta_{1}=D_{0}+\mathrm{span}\\{[\partial_{u^{2}},f]\\}$. Therefore, if a
two input system (1) allows a PAI representation, then there also must exist a
vector field $v_{c}\in D_{0}$, i. e.
$v_{c}=\alpha^{1}\partial_{u^{1}}+\alpha^{2}\partial_{u^{2}}$ with functions
$\alpha^{j}=\alpha^{j}(x,u)$, such that with
$\Delta_{1}=D_{0}+\mathrm{span}\\{[v_{c},f]\\}$, we have
$v_{c}\in\mathcal{C}(\Delta_{1})$. Since $v_{c}\in D_{0}$ and $D_{0}$ is
involutive, this is equivalent to $[v_{c},[v_{c},f]]\in\Delta_{1}$. Since we
have $\Delta_{1}\subset D_{1}$, this also implies $[v_{c},[v_{c},f]]\in
D_{1}$. Inserting
$v_{c}=\alpha^{1}\partial_{u^{1}}+\alpha^{2}\partial_{u^{2}}$ into the latter
relation yields (5). On the other hand, given a vector field
$v_{c}=\alpha^{1}\partial_{u^{1}}+\alpha^{2}\partial_{u^{2}}$ which meets
$[v_{c},[v_{c},f]]\in\Delta_{1}$, a PAI representation can be derived
following the two steps in Section III. To show that this procedure indeed
leads to a PAI representation, evaluate the condition
$[v_{c},[v_{c},f]]\in\Delta_{1}$ after applying the two steps of the procedure
in Section III.

We do not show in detail here why the condition (5) admits at most two
independent nontrivial solutions for systems which do not allow an AI
representation, i. e. $D_{0}\not\subset\mathcal{C}(D_{1})$. The proof is based
on the fact that $D_{0}\not\subset\mathcal{C}(D_{1})$ implies that at least
one of the vector fields $\partial_{u^{1}}^{2}f$,
$\partial_{u^{1}}\partial_{u^{2}}f$ or $\partial_{u^{2}}^{2}f$ in (5) is not
contained in $D_{1}$. If non of them is contained in $D_{1}$, it follows that
(5) admits no non-trivial solution. If one of them is$\mathrm{~{}mod~{}}D_{1}$
linearly dependent of the other two, (5) admits at most one independent non-
trivial solution. If two of them are$\mathrm{~{}mod~{}}D_{1}$ linearly
dependent of one of them, (5) admits at most two independent non-trivial
solutions. ∎

### Proof of Lemma 3

Given a flat output $y$ of the complete system (7), we have a certain
difference $d=\\#R-n$ (which by assumption is at most two) of the dimensions
of the domain and the codomain of the corresponding parameterizing map

$\displaystyle(x_{1},x_{2})$ $\displaystyle=F_{x}(y_{[R-1]})$ (16a)
$\displaystyle u$ $\displaystyle=F_{u}(y_{[R]})\,.$ (16b)

The parameterizing map of the subsystem $\Sigma_{1}$ in (7) only consists of
(16a), i. e. $y$ is also a flat output of the subsystem $\Sigma_{1}$. For the
parameterizing map (16a) of $\Sigma_{1}$ we have the same difference of the
dimensions of the domain and the codomain, i. e. the same value for $d$.
Indeed, in (16a), we have a domain of dimension $\\#(R-1)+2=\\#R$ and a
codomain of dimension $n$ and thus again a difference of $d=\\#R-n$.
Similarly, it can be shown that every flat output of $\Sigma_{1}$ with a
certain difference $d$ is also a flat output for the complete system with the
same difference $d$. ∎

The proofs of Lemma 4 and 5 are based on the following Theorem, which is
proven at the end of this section.

###### Theorem 2

A system (1) with $d\leq 2$ can be rendered static feedback linearizable by
$d$-fold prolonging a suitably chosen (new) input after a suitable input
transformation has been applied. If the system under consideration is an AI-
system, there always exists an affine input transformation which generates the
required (new) input.

### Proof of Lemma 4

We prove the cases $d=1$ and $d=2$ separately. Let us first assume that the
AI-system has a difference of $d=1$. According to Theorem 2, it becomes static
feedback linearizable after applying a suitable affine input transformation
$\bar{u}^{j}=g(x)^{j}+m^{j}_{k}(x)u^{k}$ and subsequently one-fold prolonging
the new input $\bar{u}^{1}$, i. e. the prolonged system

$\displaystyle\begin{aligned}
\dot{x}&=\bar{a}(x)+\bar{b}_{1}(x)\bar{u}^{1}+\bar{b}_{2}(x)\bar{u}^{2}\,,\end{aligned}$
$\displaystyle\begin{aligned} \dot{\bar{u}}^{1}&=\bar{u}^{1}_{1}\end{aligned}$

with the state $(x,\bar{u}^{1})$ and the input $(\bar{u}^{1}_{1},\bar{u}^{2})$
is static feedback linearizable. Thus, the distributions
$\Delta_{1}=\mathrm{span}\\{\partial_{\bar{u}^{1}},\bar{b}_{2}\\}$,
$\Delta_{2}=\mathrm{span}\\{\partial_{\bar{u}^{1}},\bar{b}_{1},\bar{b}_{2},[\bar{a},\bar{b}_{2}]+\bar{u}^{1}[\bar{b}_{1},\bar{b}_{2}]\\}$,
$\Delta_{3},\ldots$ are all involutive. The involutivity of $\Delta_{2}$
implies that
$[\bar{a},\bar{b}_{2}]\in\mathrm{span}\\{\bar{b}_{1},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}$
and thus, we actually have
$\Delta_{2}=\mathrm{span}\\{\partial_{\bar{u}^{1}},\bar{b}_{1},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}=\mathrm{span}\\{\partial_{\bar{u}^{1}}\\}+D_{1}^{(1)}$,
which in turn implies that $D_{1}^{(1)}$ is involutive, i. e.
$\overline{D}_{1}=D_{1}^{(1)}$. Based on these considerations, in the
following we explain how to find an input transformation required to generate
an input which needs to be prolonged in order to render the system static
feedback linearizable. If $n=3$, it can be shown that we can prolong any input
$\tilde{u}^{1}=g(x)^{1}+m^{1}_{k}(x)u^{k}$ of the system and obtain a static
feedback linearizable prolonged system. For $n\geq 4$, the direction of the
vector field $\bar{b}_{2}$ is uniquely determined by the condition
$[a,\bar{b}_{2}]\in D_{1}^{(1)}$ (i. e. it can be shown that
$[a,\bar{b}_{1}]\not\in D_{1}^{(1)}$). The construction is as follows.
Calculate $D_{1}^{(1)}=\mathrm{span}\\{b_{1},b_{2},[b_{1},b_{2}]\\}$ and find
functions $\alpha^{1}$ and $\alpha^{2}$ such that
$\alpha^{1}[a,b_{1}]+\alpha^{2}[a,b_{2}]\in D_{1}^{(1)}$. The functions
$\alpha^{1}$ and $\alpha^{2}$ are of course only unique up to a multiplicative
factor. The vector field $b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ is then
collinear with the vector field $\bar{b}_{2}$ from above. By applying the
input transformation $\tilde{u}^{1}=\alpha^{2}u^{1}-\alpha^{1}u^{2}$,
$\tilde{u}^{2}=\beta^{1}u^{1}+\beta^{2}u^{2}$, with $\beta^{1}$ and
$\beta^{2}$ chosen such that the transformation is invertible, i. e.
$det=\alpha^{1}\beta^{1}+\alpha^{2}\beta^{2}\neq 0$, to the original system,
we obtain

$\displaystyle\begin{aligned}
\dot{x}&=a(x)+\tfrac{1}{det}((\beta^{2}b_{1}-\beta^{1}b_{2})\tilde{u}^{1}+\overbrace{(\alpha^{1}b_{1}+\alpha^{2}b_{2})}^{b_{c}}\tilde{u}^{2})\,.\end{aligned}$

By prolonging the input $\tilde{u}^{1}$, i. e. adding the equation
$\dot{\tilde{u}}^{1}=\tilde{u}^{1}_{1}$, we obtain

$\displaystyle\begin{aligned}
\Sigma_{1}:&&\dot{x}&=a(x)+\tfrac{1}{det}((\beta^{2}b_{1}-\beta^{1}b_{2})\tilde{u}^{1}+b_{c}\tilde{u}^{2})\\\
\Sigma_{2}:&&\dot{\tilde{u}}^{1}&=\tilde{u}^{1}_{1}\,,\end{aligned}$ (17)

with the state $x_{p}=(x,\tilde{u}^{1})$ and the input
$u_{p}=(\tilde{u}^{1}_{1},\tilde{u}^{2})$. Calculating the distributions
involved in the test for static feedback linearizability of (17), we obtain
the distributions $\mathrm{span}\\{\partial_{\tilde{u}^{1}},b_{c}\\}$,
$\mathrm{span}\\{\partial_{\tilde{u}^{1}},b_{1},b_{2},[b_{1},b_{2}]\\}\,,\ldots$
which all can be shown to be involutive based on the involutivity of the
distributions $\Delta_{i}$ from above, and thus, (17) is static feedback
linearizable. The linearizing outputs of (17) are flat outputs with a
difference of $d=1$ of the original system $\Sigma_{1}$. Indeed, let $y$ be a
linearizing output of (17). For the parameterization of the states and inputs
of (17) with respect to $y$, we have the diffeomorphism

$\displaystyle\begin{aligned} (x,\tilde{u}^{1})&=F_{x_{p}}(y_{[R-1]})\\\
(\tilde{u}^{1}_{1},\tilde{u}^{2})&=F_{u_{p}}(y_{[R]})\end{aligned}$

with $\\#R=n+1$. It contains the map $(x,\tilde{u}^{1})=F_{x_{p}}(y_{[R-1]})$,
$\tilde{u}^{2}=F_{u_{p}}^{1}(y_{[R]})$, which is the parameterization of the
state and the input $\tilde{u}$ of the original system with respect to the
flat output $y$. It is a submersion with a domain of dimension $\\#R+2=n+3$
and a codomain of dimension $n+2$, i. e. $d=n+3-(n+2)=1$. (It is immediate
that the parameterization of $\tilde{u}^{2}$ involves the highest derivatives
$y_{R}=(y^{1}_{r_{1}},y^{2}_{r_{2}})$, otherwise, the original system would be
static feedback linearizable, i. e. it would have a difference of $d=0$, which
contradicts with our assumption $d=1$.)

For $d=2$, it follows from Theorem 2 that there exists an affine input
transformation $\bar{u}^{j}=g(x)^{j}+m^{j}_{k}(x)u^{k}$ such that the
prolonged system

$\displaystyle\begin{aligned}
\dot{x}&=\bar{a}(x)+\bar{b}_{1}(x)\bar{u}^{1}+\bar{b}_{2}(x)\bar{u}^{2}\\\
\dot{\bar{u}}^{1}&=\bar{u}^{1}_{1}\\\
\dot{\bar{u}}^{1}_{1}&=\bar{u}^{1}_{2}\end{aligned}$ (18)

withe the state $(x,\bar{u}^{1},\bar{u}^{1}_{1})$ and the input
$(\bar{u}^{1}_{2},\bar{u}^{2})$ is static feedback linearizable. Thus, the
distributions
$\Delta_{1}=\mathrm{span}\\{\partial_{\bar{u}^{1}_{1}},\bar{b}_{2}\\}$,
$\Delta_{2}=\mathrm{span}\\{\partial_{\bar{u}^{1}_{1}},\partial_{\bar{u}^{1}},\bar{b}_{2},[\bar{a},\bar{b}_{2}]+\bar{u}^{1}[\bar{b}_{1},\bar{b}_{2}]\\}$,
$\Delta_{3}=\mathrm{span}\\{\partial_{\bar{u}^{1}_{1}},\partial_{\bar{u}^{1}},\bar{b}_{1},\bar{b}_{2},[\bar{a},\bar{b}_{2}]+\bar{u}^{1}[\bar{b}_{1},\bar{b}_{2}],[\bar{a},[\bar{a},\bar{b}_{2}]]+\bar{u}^{1}[\bar{a},[\bar{b}_{1},\bar{b}_{2}]]+\bar{u}^{1}[\bar{b}_{1},[\bar{a},\bar{b}_{2}]]+(\bar{u}^{1})^{2}[\bar{b}_{1},[\bar{b}_{1},\bar{b}_{2}]]+\bar{u}^{2}[\bar{b}_{2},[\bar{a},\bar{b}_{2}]]+\bar{u}^{1}\bar{u}^{2}[\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]]\\}$,
$\Delta_{4},\ldots$ are all involutive. The involutivity of $\Delta_{2}$
implies that
$[\bar{a},\bar{b}_{2}]\in\mathrm{span}\\{\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}$,
i. e. we actually have
$\Delta_{2}=\mathrm{span}\\{\partial_{\bar{u}^{1}_{1}},\partial_{\bar{u}^{1}},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}$,
which in turn implies that $\bar{b}_{2}\in\mathcal{C}(D_{1}^{(1)})$ (we have
$D_{1}=\mathrm{span}\\{\bar{b}_{1},\bar{b}_{2}\\}$ and
$D_{1}^{(1)}=\mathrm{span}\\{\bar{b}_{1},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}$).
The distribution $\Delta_{3}$ simplifies to
$\Delta_{3}=\mathrm{span}\\{\partial_{\bar{u}^{1}_{1}},\partial_{\bar{u}^{1}},\bar{b}_{1},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}],[\bar{a},[\bar{b}_{1},\bar{b}_{2}]],[\bar{b}_{1},[\bar{b}_{1},\bar{b}_{2}]]\\}$,
where the vector fields $[\bar{a},[\bar{b}_{1},\bar{b}_{2}]]$,
$[\bar{b}_{1},[\bar{b}_{1},\bar{b}_{2}]]$ are linearly
dependent$\mathrm{~{}mod~{}}D_{1}^{(1)}=\mathrm{span}\\{\bar{b}_{1},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}$.
We have to distinguish between two subcases, namely between $D_{1}^{(1)}$
being involutive or not. If $D_{1}^{(1)}$ is non-involutive, we necessarily
have $[\bar{b}_{1},[\bar{b}_{1},\bar{b}_{2}]]\notin D_{1}^{(1)}$ and thus
$\mathrm{dim}(D_{1}^{(2)})=4$ and $\mathrm{dim}(\mathcal{C}(D_{1}^{(1)}))=1$.
Furthermore, we then have
$\Delta_{3}=\mathrm{span}\\{\partial_{\bar{u}^{1}_{1}},\partial_{\bar{u}^{1}}\\}+D_{1}^{(2)}$,
implying that $\overline{D}_{1}=D_{1}^{(2)}$. So in this case, the direction
of $\bar{b}_{2}$ is uniquely determined by the condition
$\bar{b}_{2}\in\mathcal{C}(D_{1}^{(1)})$. The construction is as follows.
Calculate the Cauchy characteristic distribution $\mathcal{C}(D_{1}^{(1)})$
and find functions $\alpha^{1}$ and $\alpha^{2}$ such that
$b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}\in\mathcal{C}(D_{1}^{(1)})$. The
functions $\alpha^{1}$ and $\alpha^{2}$ are again only unique up to a
multiplicative factor. The vector field
$b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ is then collinear with the vector
field $\bar{b}_{2}$ from above. If $D_{1}^{(1)}$ is involutive, i. e.
$\overline{D}_{1}=D_{1}^{(1)}$, we obviously have
$[\bar{b}_{1},[\bar{b}_{1},\bar{b}_{2}]]\in
D_{1}^{(1)}=\mathrm{span}\\{\bar{b}_{1},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}$
and thus
$\Delta_{3}=\mathrm{span}\\{\partial_{\bar{u}^{1}_{1}},\partial_{\bar{u}^{1}},\bar{b}_{1},\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}],[\bar{a},[\bar{b}_{1},\bar{b}_{2}]]\\}$.
In this case, the condition $\bar{b}_{2}\in\mathcal{C}(D_{1}^{(1)})$ is not
useful since $\mathcal{C}(D_{1}^{(1)})=D_{1}^{(1)}$, but then, analogous to
the previously considered case $d=1$, the condition $[\bar{a},\bar{b}_{2}]\in
D_{1}^{(1)}$ (which follows from
$[\bar{a},\bar{b}_{2}]\in\mathrm{span}\\{\bar{b}_{2},[\bar{b}_{1},\bar{b}_{2}]\\}\subset
D_{1}^{(1)}$) uniquely determines the direction of $\bar{b}_{2}$. The
construction of the vector field $b_{c}=\alpha^{1}b_{1}+\alpha^{2}b_{2}$ which
is collinear with $\bar{b}_{2}$ is then analogous to the previously considered
case $d=1$.

###### Remark 1

In the following, we explain how to derive the input transformation required
to generate the input which needs to be prolonged twice in order to render the
system static feedback linearizable. In contrast to the previously considered
case $d=1$, in the case $d=2$, additional functions $\gamma^{1},\gamma^{2}$
and $\delta$ of the state $x$ of the system are involved in the input
transformation. We do not discuss their construction. For the proof, only the
existence of these functions is of importance and for applying the procedure,
we do not have to construct these functions either.

In any of the two subcases, i. e. independently of $D_{1}^{(1)}$ being
involutive or not, by applying the input transformation
$\tilde{u}^{1}=(\alpha^{2}u^{1}-\alpha^{1}u^{2})\delta+\gamma^{1}$,
$\tilde{u}^{2}=\beta^{1}u^{1}+\beta^{2}u^{2}+\gamma^{2}$, with $\beta^{1}$ and
$\beta^{2}$ chosen such that the transformation is invertible, i. e.
$det=\alpha^{1}\beta^{1}+\alpha^{2}\beta^{2}\neq 0$, and functions
$\gamma^{1},\gamma^{2}$ and $\delta$ of the state $x$ of the system, to the
original system, we obtain

$\displaystyle\begin{aligned}
\dot{x}&=\underbrace{a(x)+\tfrac{\gamma^{1}}{\delta
det}(\beta^{2}b_{1}-\beta^{1}b_{2})-\tfrac{\gamma^{2}}{det}\overbrace{(\alpha^{1}b_{1}+\alpha^{2}b_{2})}^{b_{c}}}_{\tilde{a}}+\\\
&\hskip 30.00005pt\underbrace{\tfrac{1}{\delta
det}(\beta^{2}b_{1}-\beta^{1}b_{2})}_{\tilde{b}_{1}}\tilde{u}^{1}+\tfrac{1}{det}\underbrace{(\alpha^{1}b_{1}+\alpha^{2}b_{2})}_{b_{c}}\tilde{u}^{2}\,.\end{aligned}$

Assume that the functions $\gamma^{1}$, $\gamma^{2}$ and $\delta$ are chosen
such that $\tilde{a}$ coincides$\mathrm{~{}mod~{}}\mathrm{span}\\{b_{c}\\}$
with $\bar{a}$ from above and $\tilde{b}_{1}$
coincides$\mathrm{~{}mod~{}}\mathrm{span}\\{b_{c}\\}$ with $\bar{b}_{1}$ from
above, which is indeed always possible (in fact, we can always choose
$\gamma^{2}=0$). Then, based on the properties of the vector fields $\bar{a}$,
$\bar{b}_{1}$, $\bar{b}_{2}$ and the distributions $\Delta_{i}$ form above, it
can be shown that we have
$[\tilde{a},b_{c}]\in\mathrm{span}\\{b_{c},[\tilde{b}_{1},b_{c}]\\}$ with
$\mathrm{span}\\{b_{c},[\tilde{b}_{1},b_{c}]\\}$ being involutive and that
$[\tilde{a},[\tilde{b}_{1},b_{c}]]$ and
$[\tilde{b}_{1},[\tilde{b}_{1},b_{c}]]$ are
collinear$\mathrm{~{}mod~{}}D_{1}^{(1)}$. Prolonging the input $\tilde{u}^{1}$
two-fold, i. e. adding the equations $\dot{\tilde{u}}^{1}=\bar{u}^{1}_{1}$ and
$\dot{\tilde{u}}^{1}_{1}=\tilde{u}^{1}_{2}$ we obtain

$\displaystyle\begin{aligned} &\Sigma_{1}:&&\begin{aligned}
\dot{x}&=\tilde{a}+\tilde{b}_{1}\tilde{u}^{1}+\tfrac{1}{det}b_{c}\tilde{u}^{2}\end{aligned}\\\\[4.30554pt]
&\Sigma_{2}:&&\begin{aligned} \dot{\tilde{u}}^{1}&=\tilde{u}^{1}_{1}\\\
\dot{\tilde{u}}^{1}_{1}&=\tilde{u}^{1}_{2}\,,\end{aligned}\end{aligned}$ (19)

with the state $x_{p}=(x,\tilde{u}^{1},\tilde{u}^{1}_{1})$ and the input
$u_{p}=(\tilde{u}^{1}_{2},\tilde{u}^{2})$. Calculating the distributions
involved in the test for static feedback linearizability of (19), we then
obtain $\mathrm{span}\\{\partial_{\tilde{u}^{1}_{1}},b_{c}\\}$,
$\mathrm{span}\\{\partial_{\tilde{u}^{1}_{1}},\partial_{\tilde{u}^{1}},b_{c},[\tilde{b}_{1},b_{c}]\\}$,
$\mathrm{span}\\{\partial_{\tilde{u}^{1}_{1}},\partial_{\tilde{u}^{1}},\tilde{b}_{1},b_{c},[\tilde{b}_{1},b_{c}],[\tilde{a},[\tilde{b}_{1},b_{c}]],[\tilde{b}_{1},[\tilde{b}_{1},b_{c}]]\\}\,,\ldots$
which all can be shown to be involutive based on the involutivity of the
distributions $\Delta_{i}$ from above. The linearizing outputs of (19) are
flat outputs with a difference of $d=2$ of the original system $\Sigma_{1}$.
Indeed, let $y$ be a linearizing output of (19). For the parameterization of
the states and inputs of (19) with respect to $y$, we have the diffeomorphism

$\displaystyle\begin{aligned}
(x,\tilde{u}^{1},\tilde{u}^{1}_{1})&=F_{x_{p}}(y_{[R-1]})\\\
(\tilde{u}^{1}_{2},\tilde{u}^{2})&=F_{u_{p}}(y_{[R]})\end{aligned}$

with $\\#R=n+2$ (in fact the parameterization of $\tilde{u}^{1}$ can only
involve derivatives of $y$ up to the order $R-2$ due to the structure of
(19)). It contains the map

$\displaystyle\begin{aligned}
(x,\tilde{u}^{1})&=F_{x_{p}}^{1,\ldots,n+1}(y_{[R-1]})\\\
\tilde{u}^{2}&=F_{u_{p}}^{1}(y_{[R]})\,,\end{aligned}$

which is the parameterization of the state and the input $\tilde{u}$ of the
original system with respect to the flat output $y$. It is a submersion with a
domain of dimension $\\#R+2=n+4$ and a codomain of dimension $n+2$, i. e.
$d=n+4-(n+2)=2$. (It is immediate that parameterization of $\tilde{u}^{2}$
involves the highest derivatives $y_{R}=(y^{1}_{r_{1}},y^{2}_{r_{2}})$,
otherwise, the system would have a difference of $d\leq 1$, which contradicts
with our assumption $d=2$.) As already mentioned in Remark 1, we do not need
to find the appropriate functions $\gamma^{1}$ and $\delta$. Since
$\tilde{u}^{1}$ with the correctly chosen functions $\gamma^{1}$ and $\delta$
only depends on $y_{[R-1]}$ (actually $y_{[R-2]}$ due to the structure of
(19)) and the functions $\gamma^{1}$ and $\delta$ only depend on the state $x$
of the original system, which on its own only depends on $y_{[R-1]}$, an
arbitrary choice for these functions, e. g. $\gamma^{1}=0$ and $\delta=1$
which would yield the input transformation
$\hat{u}^{1}=\alpha^{2}u^{1}-\alpha^{1}u^{2}$,
$\hat{u}^{2}=\beta^{1}u^{1}+\beta^{2}u^{2}$, still yields an input
$\hat{u}^{1}$ which only depends on $y_{[R-1]}$. The flat output $y$ with
$d=2$ of the original system is thus a flat output with $d=1$ for the
prolonged system

$\displaystyle\begin{aligned}
\dot{x}&=a(x)+\hat{b}_{1}(x)\hat{u}^{1}+\hat{b}_{2}(x)\hat{u}^{2}\,,\end{aligned}$
$\displaystyle\begin{aligned}
\dot{\hat{u}}^{1}&=\hat{u}^{1}_{1}\,,\end{aligned}$ (20)

i. e. (20) has a difference of $d=1$. Conversely, by an analogues reasoning as
before based on the parameterizing map of (20) with respect to any flat output
with $d=1$ of (20), it can be shown that any flat output with $d=1$ of (20) is
a flat output with $d=2$ of the original system. ∎

### Proof of Lemma 5

According to Theorem 2, the system becomes static feedback linearizable after
applying a suitable input transformation $\bar{u}=\Phi_{u}(x,u)$ and
subsequently $d$-fold prolonging the new input $\bar{u}^{1}$, i. e. the
prolonged system

$\displaystyle\begin{aligned} \dot{x}&=\bar{f}(x,\bar{u}^{1},\bar{u}^{2})\\\
\dot{\bar{u}}^{1}&=\bar{u}^{1}_{1}\end{aligned}$ or
$\displaystyle\begin{aligned} \dot{x}&=\bar{f}(x,\bar{u}^{1},\bar{u}^{2})\\\
\dot{\bar{u}}^{1}&=\bar{u}^{1}_{1}\\\
\dot{\bar{u}}^{1}_{1}&=\bar{u}^{1}_{2}\,,\end{aligned}$ (21)

depending on the actual value of $d$, is static feedback linearizable.
Analogous to the proof of Lemma 4, it can be shown that the linearizing
outputs $y$ of (21) are flat outputs with a difference of $d=1$ or $d=2$ of
the original system and that the parameterization of the input $\bar{u}^{1}$
by $y$ involves derivatives up to order $R-d$ only. Independent of the actual
value of $d$, the involutivity of
$\mathrm{span}\\{\partial_{\bar{u}^{1}_{d}},\partial_{\bar{u}^{2}},\partial_{\bar{u}^{1}_{d-1}},\partial_{\bar{u}^{2}}\bar{f}^{i}\partial_{x^{i}}\\}$
implies that $\dot{x}=\bar{f}(x,\bar{u}^{1},\bar{u}^{2})$ is actually of the
form
$\dot{x}=\bar{a}(x,\bar{u}^{1})+b(x,\bar{u}^{1})h(x,\bar{u}^{1},\bar{u}^{2})$
and introducing $\tilde{u}^{2}=h(x,\bar{u}^{1},\bar{u}^{2})$ thus results in
the PAI representation

$\displaystyle\begin{aligned}
\dot{x}&=a(x,\bar{u}^{1})+b(x,\bar{u}^{1})\tilde{u}^{2}\,.\end{aligned}$ (22)

This PAI representation is special, since the parameterization of the non-
affine occurring input $\bar{u}^{1}$ involves derivatives of $y$ up to order
$R-d$ only. For every PAI representation which is equivalent to (22) via a PAI
form preserving transformation (6), we still have that the non-affine input
only depends on $y_{[R-1]}$, since in such transformations the non-affine
input is only combined with the state $x$ of the system and we have
$x=F_{x}(y_{[R-1]})$. Recall that there exist at most two fundamentally
different PAI representations of a system, all others are equivalent to one of
those two by such a PAI form preserving transformation (6). Therefore, given
any two non-equivalent PAI representations of the system (which can be derived
systematically, provided they exist), one of them is equivalent to the special
PAI representation (22) via a transformation of the form (6). The non-affine
inputs of any two non-equivalent PAI representations are thus the candidates
for an input whose parameterization involves derivatives of $y$ up to order
$R-1$ only. Similar as in the proof of Lemma 4, it can then be shown that at
least one of the possibly two AI-systems obtained by prolonging the
corresponding non-affine inputs has a difference of $d\leq 1$ and that
furthermore the linearizing outputs or flat outputs with $d=1$ of the AI-
systems are flat outputs with $d=1$ or $d=2$ of the original system. ∎

### Proof of Theorem 2

The proof of Theorem 2 is based on the following results.

###### Lemma 6

A system (1) with $d\leq 2$ is $(x,u)$-flat, an AI-system (3) with $d\leq 2$
is $x$-flat.

Due to space limitations, we do not provide a proof of this lemma here. The
second part of the lemma, i. e. that AI-systems with $d\leq 2$ are $x$-flat,
can also be found in [13], Proposition 1.

###### Theorem 3

Every $(x,u)$-flat system with two inputs can be rendered static feedback
linearizable by $d=\\#R-n$ fold prolonging a suitably chosen (new) input.

A proof of this theorem can be found in [22]. Since according to Lemma 6,
systems with $d\leq 2$ are $(x,u)$-flat, Theorem 3 always applies to these
systems and in turn, such systems can always be rendered static feedback
linearizable by prolongations of a suitably chosen (new) input, which shows
the first part of Theorem 2. From the proof of Theorem 3 in [22], it follows
that the input transformation which has to be applied in order to generate the
required input which needs to be prolonged, is given by
$\bar{u}^{1}=\mathrm{L}_{f}^{k_{j}}\varphi^{j}(x,u)$, $j=1$ or $j=2$ where
$k_{j}\geq 0$ denotes the relative degree of the component
$y^{j}=\varphi^{j}(x,u)$ of the flat output, i. e.
$\mathrm{L}_{f}^{k_{j}-1}\varphi^{j}=\mathrm{L}_{f}^{k_{j}-1}\varphi^{j}(x)$
and
$\mathrm{L}_{f}^{k_{j}}\varphi^{j}=\mathrm{L}_{f}^{k_{j}}\varphi^{j}(x,u)$.
Since according to Lemma 6 an AI-system with $d\leq 2$ is $x$-flat, i. e.
$y^{j}=\varphi^{j}(x)$, there always exists an affine input transformation
which generates the required input since
$\bar{u}^{1}=\mathrm{L}_{(a^{i}+u^{1}b_{1}^{i}+u^{2}b_{2}^{i})\partial_{x^{i}}}^{k_{j}}\varphi^{j}(x)$
is an affine input transformation. ∎

## References

  * [1] M. Fliess, J. Lévine, P. Martin, and P. Rouchon, “Sur les systèmes non linéaires différentiellement plats,” _Comptes rendus de l’Académie des sciences. Série I, Mathématique_ , vol. 315, pp. 619–624, 1992.
  * [2] ——, “Flatness and defect of non-linear systems: introductory theory and examples,” _International Journal of Control_ , vol. 61, no. 6, pp. 1327–1361, 1995.
  * [3] B. Jakubczyk and W. Respondek, “On linearization of control systems,” _Bull. Acad. Polonaise Sci. Ser. Sci. Math._ , vol. 28, pp. 517–522, 1980\.
  * [4] H. Nijmeijer and A. van der Schaft, _Nonlinear Dynamical Control Systems_. New York: Springer, 1990.
  * [5] P. Martin and P. Rouchon, “Feedback linearization and driftless systems,” _Mathematics of Control, Signals and Systems_ , vol. 7, no. 3, pp. 235–254, 1994.
  * [6] F. Nicolau and W. Respondek, “Flatness of multi-input control-affine systems linearizable via one-fold prolongation,” _SIAM J. Control and Optimization_ , vol. 55, pp. 3171–3203, 2017.
  * [7] J. Pomet, “On dynamic feedback linearization of four-dimensional affine control systems with two inputs,” _ESAIM: Control, Optimisation and Calculus of Variations_ , vol. 2, pp. 151–230, 1997.
  * [8] M. Fliess, J. Lévine, P. Martin, and P. Rouchon, “A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems,” _IEEE Transactions on Automatic Control_ , vol. 44, no. 5, pp. 922–937, 1999.
  * [9] K. Schlacher and M. Schöberl, “A jet space approach to check Pfaffian systems for flatness,” in _Proceedings 52nd IEEE Conference on Decision and Control (CDC)_ , 2013, pp. 2576–2581.
  * [10] M. Schöberl and K. Schlacher, “On an implicit triangular decomposition of nonlinear control systems that are 1-flat - a constructive approach,” _Automatica_ , vol. 50, pp. 1649–1655, 2014.
  * [11] K. Schlacher, M. Schöberl, and B. Kolar, “A jet space approach to derive flat outputs,” in _Proceedings 1st IFAC Conference on Modelling, Identification and Control of Nonlinear Systems (MICNON)_ , 2015, iFAC-PapersOnLine, volume 48, issue 11, pages 131–136.
  * [12] F. Nicolau and W. Respondek, “Two-input control-affine systems linearizable via one-fold prolongation and their flatness,” _European Journal of Control_ , vol. 28, pp. 20–37, 2016.
  * [13] ——, “Flatness of two-input control-affine systems linearizable via a two-fold prolongation,” in _2016 IEEE 55th Conference on Decision and Control (CDC)_ , Dec 2016, pp. 3862–3867.
  * [14] C. Gstöttner, B. Kolar, and M. Schöberl, “A structurally flat triangular form based on the extended chained form,” _International Journal of Control_ , 2020.
  * [15] B. Kolar, _Contributions to the Differential Geometric Analysis and Control of Flat Systems_. Aachen: Shaker Verlag, 2017.
  * [16] W. Sluis, “A necessary condition for dynamic feedback linearization,” _Systems & Control Letters_, vol. 21, pp. 277–283, 1993.
  * [17] B. Charlet, J. Lévine, and R. Marino, “Sufficient conditions for dynamic state feedback linearization,” _SIAM J. Control Optim._ , vol. 29, no. 1, pp. 38–57, 1991.
  * [18] M. Schöberl, K. Rieger, and K. Schlacher, “System parametrization using affine derivative systems,” in _Proceedings 19th International Symposium on Mathematical Theory of Networks and Systems (MTNS)_ , 2010, pp. 1737–1743.
  * [19] M. Schöberl and K. Schlacher, “On calculating flat outputs for pfaffian systems by a reduction procedure - demonstrated by means of the vtol example,” in _9th IEEE International Conference on Control $\&$ Automation (ICCA11)_, 2011, pp. 477–482.
  * [20] J. Lévine, _Analysis and Control of Nonlinear Systems: A Flatness-based Approach_. Berlin: Springer, 2009.
  * [21] M. Schöberl, _Contributions to the Analysis of Structural Properties of Dynamical Systems in Control and Systems Theory - A Geometric Approach_. Aachen: Shaker Verlag, 2014\.
  * [22] C. Gstöttner, B. Kolar, and M. Schöberl, “On the linearization of flat two-input systems by prolongations and applications to control design,” in _IFAC-PapersOnLine_ , 2020, pp. 5553–5560, 21th IFAC World Congress (Virtual).


# A variational frequency-dependent stabilization for the Helmholtz equation
with noisy Cauchy data

Vo Anh Khoa<EMAIL_ADDRESS>Nguyen Dat Thuc  Ajith Gunaratne  Department
of Mathematics, Florida A&M University, Tallahassee, FL 32307, USA  Department
of Mathematics and Computer Science, VNUHCM-University of Science, Ho Chi Minh
City, Vietnam

###### Abstract

This article considers a Cauchy problem of Helmholtz equations whose solution
is well known to be exponentially unstable with respect to the inputs. In the
framework of variational quasi-reversibility method, a Fourier truncation is
applied to appropriately perturb the underlying problem, which allows us to
obtain a stable approximate solution. Unlike the original version of the
truncation, our perturbation is driven by the frequencies, rather than the
noise level. This deals with real-world circumstances that one only measures
data once or even does not know the noise level in advance, but needs to
truncate high Fourier frequencies appropriately. The corresponding approximate
problem is of a hyperbolic equation, which is also a crucial aspect of this
approach. Error estimate between the approximate and true solutions is derived
with respect to the noise level and to the frequencies is derived. From this
analysis, the Lipschitz stability with respect to the noise level follows.
Some numerical examples are provided to see how our numerical algorithm works
well.

###### keywords:

Stabilization, Helmholtz equation , Ill-posed problems , Convergence , Error
estimate

###### MSC:

[2010] 65J05 , 65J20 , 35K92

††journal: Applied Mathematics Letters

textheight=25cm textwidth=17cm

mytitlenotemytitlenotefootnotetext: This work was supported by the Faculty
Research Awards Program (FRAP) at Florida A&M University, under the project
“Approximation of a time-reversed reaction-diffusion system in cancer cell
population dynamics” (007633).

##  1 Statement of the Cauchy problem

In this work, we are concerned with the reconstruction problem of
electromagnetic field from its knowledge on a part of boundary of the physical
region $\Omega$. Here, $\Omega=(0,1)\times(0,1)$ is our computational domain
of interest, but it can be extended easily to $(0,a_{1})\times(0,a_{2})$,
where $a_{1},a_{2}$ are two positive numbers. Often, the propagation of the
electromagnetic wave field is governed by the system of the Maxwell’s
equations for the electric field $\mathbf{E}=\mathbf{E}(x,y,t)$ and the
magnetic field $\mathbf{B}=\mathbf{B}(x,y,t)$. Considering $\Omega$ as a
homogeneous medium in a region with free currents and charges, this system can
be reduced to the classical wave equations, cf. [1],

$\displaystyle\frac{\partial^{2}\mathbf{E}}{\partial
t^{2}}-c^{2}\Delta\mathbf{E}=0,\quad\frac{\partial^{2}\mathbf{B}}{\partial
t^{2}}-c^{2}\Delta\mathbf{B}=0\quad\text{in }\Omega\times(0,T),$ (1)

where $c>0$ is the speed of light and $T>0$ is the travel time. Consider the
frequency $\omega>0$. For $\text{i}=\sqrt{-1}$, we take
$\mathbf{E}(x,y,t)=e^{\text{i}\omega t}\mathbf{E}(x,y)$ and
$\mathbf{B}(x,y,t)=e^{\text{i}\omega t}\mathbf{B}(x,y)$. Then, setting
$k=\omega/c>0$, it follows from system (1) that

$\displaystyle\Delta\mathbf{E}(x,y)+k^{2}\mathbf{E}(x,y)=0,\quad\Delta\mathbf{B}(x,y)+k^{2}\mathbf{B}(x,y)=0,$
(2)

which form a system of Helmholtz equations. Since system (2) is uncoupled and
linear with respect to each component of $\mathbf{E}(x,y)$ and
$\mathbf{B}(x,y)$, it is pertinent to solve the following model:

$\displaystyle\Delta u\left(x,y\right)+k^{2}u\left(x,y\right)=0\quad\text{in
}\Omega.$ (3)

Note that in (2), $\mathbf{E}(x,y)$ and $\mathbf{B}(x,y)$ are complex-valued
components, but it is sufficient to find a real-valued function $u=u(x,y)$ in
(3). Physically, fields vanish far from the axes and thus, we can assume that
the electromagnetic field vanishes on the sides $\left\\{y=0\right\\}$ and
$\left\\{y=1\right\\}$ of the computational domain $\Omega$. Mathematically,
we consider

$\displaystyle u(x,0)=u(x,1)=0\quad\text{for }x\in(0,1).$ (4)

On the other hand, we assume to measure the electromagnetic Cauchy data at
$\left\\{x=0\right\\}$,

$\displaystyle u(0,y)=u_{0}(y),\quad u_{x}(0,y)=u_{1}(y)\quad\text{for
}y\in(0,1).$ (5)

Cf. [2], we remark that the second data (i.e. the Neumann data at $x=0$) in
(5) can be reduced to the zero boundary condition. In fact, let $U=U(x,y)$ be
a solution to the following system:

$\displaystyle\begin{cases}\Delta
U\left(x,y\right)+k^{2}U\left(x,y\right)=0&\text{in }\Omega,\\\
U\left(x,0\right)=U\left(x,1\right)=0&\text{for }x\in\left(0,1\right),\\\
U_{x}\left(0,y\right)=u_{1}\left(y\right),U\left(1,y\right)=0&\text{for
}y\in\left(0,1\right).\end{cases}$ (6)

Next, consider $V=V(x,y)$ as a solution to the following system:

$\displaystyle\begin{cases}\Delta
V\left(x,y\right)+k^{2}V\left(x,y\right)=0&\text{in }\Omega,\\\
V\left(x,0\right)=V\left(x,1\right)=0&\text{for }x\in\left(0,1\right),\\\
V\left(0,y\right)=u_{0}\left(y\right)-U\left(0,y\right),V_{x}\left(0,y\right)=0&\text{for
}y\in\left(0,1\right).\end{cases}$ (7)

With (6) and (7), it is clear that the solution $u$ to system (3)–(5) can be
computed via $u=U+V$. By [2, Lemma 1], we know that system (6) is well-posed
with $U$ in $H^{2}(\Omega)$ when $u_{1}\in L^{2}(0,1)$ and thus, $U(0,y)$
exists in $H^{2}(0,1)$ by the embedding $H^{2}(\Omega)\subset
C([0,1];H^{2}(0,1))$. Henceforth, from (7), instead of working on the Cauchy
data (5) we can assume that $u_{1}=0$ in (5) in our analysis below.

Combining (3), (4) and (5) with $u_{1}=0$ forms our Cauchy problem for the
Helmholtz equation. In this scenario, we want to reconstruct the whole wave
field in $\Omega$ and especially, the field at the boundary $x=1$.

######  Remark 1.

Cf. the appendix of [3], if the incident electric wave field has only one non-
zero component, then the propagation of this component in a heterogeneous
medium is governed equally well by a single Helmholtz equation. In other
words, the Helmholtz equation may play an equal role as the Maxwell’s system
when the medium is no longer homogeneous as assumed above.

##  2 Preliminary analysis and frequency-dependent stabilization

Let $A$ be either a Banach space or a Hilbert space. We call $A^{\prime}$ the
dual space of $A$. For a certain Banach space $A$, $\left\|\cdot\right\|_{A}$
stands for the $A$-norm. When $A$ is Hilbert, we define the $A$-norm of $u$ as
$\left\|u\right\|_{A}^{2}=\left\langle u,u\right\rangle_{A}$, where
$\left\langle\cdot,\cdot\right\rangle_{A}$ is the corresponding inner product.
Throughout the paper, we will use $\left\langle\cdot,\cdot\right\rangle$ to
indicate either the scalar product in $L^{2}(0,1)$ or the dual pairing of a
continuous linear functional and an element of a function space. We thereby
denote by $\left\|\cdot\right\|$ the $L^{2}(0,1)$-norm.

The Cauchy problem of Helmholtz equation is well known to be unstable with
respect to any small perturbation of the data. Based on the zero Dirichlet
boundary condition (4), the Laplace operator $-\partial^{2}/\partial y^{2}$ is
non-negative. According to the standard result for the Dirichlet eigenvalue
problem, there exists an orthonormal basis $\left\\{\phi_{j}\right\\}$ of
$L^{2}\left(0,1\right)$ such that $\phi_{j}\in H^{1}\left(0,1\right)\cap
C^{\infty}\left[0,1\right]$ and
$-d^{2}/dy^{2}\phi_{j}\left(y\right)=\mu_{j}\phi_{j}\left(y\right)$. The
Dirichlet eigenvalues $\mu_{j}$ in this case form an infinite sequence such
that $0\leq\mu_{0}<\mu_{1}<\mu_{2}<\ldots,\text{and
}\lim_{j\to\infty}\mu_{j}=\infty$. It follows from (3) and (5) that we obtain
the following initial-value differential system:

$\displaystyle\begin{cases}\frac{d^{2}}{dx^{2}}\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle-\lambda_{j,k}\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle=0,\\\ \left\langle
u\left(0,\cdot\right),\phi_{j}\right\rangle=\left\langle
u_{0},\phi_{j}\right\rangle,\quad\frac{d}{dx}\left\langle
u\left(0\right),\phi_{j}\right\rangle=0.\end{cases}$ (8)

In (8), $\lambda_{j,k}=\mu_{j}-k^{2}$. By this way, we solve system (8) in
each of the following set of Fourier frequencies:

$A_{1}:=\left\\{j\in\mathbb{N}:\lambda_{j,k}>0\right\\},\quad
A_{2}:=\left\\{j\in\mathbb{N}:\lambda_{j,k}=0\right\\},\quad
A_{3}:=\left\\{j\in\mathbb{N}:\lambda_{j,k}<0\right\\}.$

It is also straightforward to see that $\phi_{j}(y)=\sqrt{2}\sin\left(j\pi
y\right),\;\mu_{j}=j^{2}\pi^{2}$. In addition,
$\left\\{\phi_{j}^{\prime}/\sqrt{\mu_{j}}\right\\}_{j\in\mathbb{N}^{*}}$ is an
orthonormal basis of $L^{2}(0,1)$. Therefore, it holds that

$\displaystyle\left\|u_{y}\right\|^{2}=\sum_{j\in\mathbb{N}^{*}}\left|\left\langle
u_{y},\frac{\phi_{j}^{\prime}}{\sqrt{\mu_{j}}}\right\rangle\right|^{2}=\sum_{j\in\mathbb{N}^{*}}\left|\left\langle
u,\frac{\phi_{j}^{\prime\prime}}{\sqrt{\mu_{j}}}\right\rangle\right|^{2}=\sum_{j\in\mathbb{N}^{*}}\mu_{j}\left|\left\langle
u,\phi_{j}\right\rangle\right|^{2}.$ (9)

######  Theorem 1.

The Fourier coefficient $\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle$ has the form:

$\displaystyle\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle=\begin{cases}\cosh\left(\sqrt{\lambda_{j,k}}x\right)\left\langle
u_{0},\phi_{j}\right\rangle&\text{in }A_{1},\\\
\cos\left(\sqrt{-\lambda_{j,k}}x\right)\left\langle
u_{0},\phi_{j}\right\rangle&\text{in }A_{3}.\end{cases}$ (10)

###### Proof.

Proof of the theorem can be proceeded as in [2]. ∎

Now, we show a very important relation of these Fourier frequencies in the
following theorem.

######  Theorem 2.

Taking into account set $A_{1}$, the Fourier coefficient of $u$ satisfies the
following relation:

$\displaystyle\lambda_{j,k}e^{\left(1-x\right)\sqrt{\lambda_{j,k}}}\left(\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle+\frac{1}{\sqrt{\lambda_{j,k}}}\left\langle
u_{x}\left(x,\cdot\right),\phi_{j}\right\rangle\right)=\lambda_{j,k}\left\langle
u\left(1,\cdot\right),\phi_{j}\right\rangle+\sqrt{\lambda_{j,k}}\left\langle
u_{x}\left(1,\cdot\right),\phi_{j}\right\rangle.$ (11)

###### Proof.

From (10), we can compute that

$\displaystyle\frac{1}{\sqrt{\lambda_{j,k}}}\left\langle
u_{x}\left(x,\cdot\right),\phi_{j}\right\rangle=\sinh\left(\sqrt{\lambda_{j,k}}x\right)\left\langle
u_{0},\phi_{j}\right\rangle.$ (12)

Thus, we take $x=1$ in (12) and in (10) and then combine the resulting
formulations to obtain (11). We complete the proof of the theorem. ∎

Practically, the data $u_{0}$ in (5) always contain noise of measurement.
Therefore, we assume to have $u_{0}^{\varepsilon}\in H^{1}(0,1)$ as the noisy
data such that for $\varepsilon\in(0,1)$,

$\displaystyle\left\|u_{0}^{\varepsilon}-u_{0}\right\|_{H^{1}\left(0,1\right)}\leq\varepsilon.$
(13)

By Theorem 1, our Cauchy problem is exponentially unstable in $A_{1}$ due to
the natural growth of the hyperbolic cosine function. Any small perturbation
of the initial data $u_{0}$ may cause a huge error when computing solution $u$
of the Cauchy problem. In this work, we then adapt our recent modified quasi-
reversibility method (cf. [4] for elliptic operators and [5] for parabolic
operators) to solve our system (3)–(5). To do so, we rewrite (3) as

$u_{xx}-u_{yy}+2u_{yy}+k^{2}u=0.$ (14)

We then perturb (14) by a linear mapping $\mathbf{Q}$ and take
$\mathbf{P}=\mathbf{Q}+2\partial^{2}/\partial y^{2}$. Henceforth, we arrive at

$\displaystyle u_{xx}-u_{yy}+\mathbf{P}u+k^{2}u=0\quad\text{in }\Omega.$ (15)

It is worth mentioning that together with the boundary condition (4) and the
Cauchy data (5) with measurement $u_{0}^{\varepsilon}$, (15) forms a system of
wave equation. Herewith, $x$ becomes a parametric time variable. Since the
noise level $\varepsilon$ is involved, we then seek a sequence of
$\left\\{u^{\varepsilon}\right\\}_{\varepsilon>0}$ satisfying the following
system:

$\displaystyle\begin{cases}u_{xx}^{\varepsilon}-u_{yy}^{\varepsilon}+\mathbf{P}u^{\varepsilon}+k^{2}u^{\varepsilon}=0&\text{in
}\Omega,\\\
u^{\varepsilon}\left(x,0\right)=u^{\varepsilon}\left(x,1\right)=0&\text{for
}x\in\left(0,1\right),\\\
u^{\varepsilon}\left(0,y\right)=u_{0}^{\varepsilon}\left(y\right),\quad
u_{x}^{\varepsilon}\left(0,y\right)=0&\text{for
}y\in\left(0,1\right).\end{cases}$ (16)

Cf. [4, 5], $\mathbf{Q}$ is called perturbation as it is to “absorb” high
Fourier frequencies in the Laplace operator, and $\mathbf{P}$ is called
stabilized operator as it only contains large enough Fourier frequencies
serving for the convergence of the scheme. Let $\gamma>1$. Consider
$B:=\left\\{j\in A_{1}:\lambda_{j,k}>\log^{2}(\gamma)\right\\}$. We choose the
following truncation operator:

$\displaystyle\mathbf{Q}u(x,\cdot)=2\sum_{j\in B}\lambda_{j,k}\left\langle
u(x,\cdot),\phi_{j}\right\rangle\phi_{j}+2\sum_{j\in
A_{3}}\lambda_{j,k}\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle\phi_{j}:=\mathbf{Q_{1}}u(x,\cdot)+\mathbf{Q_{2}}u(x,\cdot).$
(17)

As to the corresponding stabilized operator $\mathbf{P}$, we find that

$\displaystyle\mathbf{P}u(x,\cdot)$ $\displaystyle=2\sum_{j\in B\cup
A_{3}}\lambda_{j,k}\left\langle
u(x,\cdot),\phi_{j}\right\rangle\phi_{j}-2\sum_{j\in\mathbb{N}}\mu_{j}\left\langle
u(x,\cdot),\phi_{j}\right\rangle\phi_{j}$ $\displaystyle=2\sum_{j\in B\cup
A_{3}}\left(\lambda_{j,k}-\mu_{j}\right)\left\langle
u(x,\cdot),\phi_{j}\right\rangle\phi_{j}-2\sum_{j\in\mathbb{N}\backslash(B\cup
A_{3})}\mu_{j}\left\langle u(x,\cdot),\phi_{j}\right\rangle\phi_{j}$
$\displaystyle=-2k^{2}\sum_{j\in\mathbb{N}}\left\langle
u(x,\cdot),\phi_{j}\right\rangle\phi_{j}-2\sum_{j\in\mathbb{N}\backslash(B\cup
A_{3})}\lambda_{j,k}\left\langle u(x,\cdot),\phi_{j}\right\rangle\phi_{j}.$

######  Remark 2.

Many works concerning the conventional truncation method (see, e.g., [4, 6, 7,
8, 9]) require the noise dependence to prove the strong convergence of the
truncation projection in different physical models. However, practically we
cannot know exactly the noise level or, if that noise is known, it is only
measured once. Thus, it is reasonable to get another way to truncate high
Fourier frequencies appropriately. We cannot truncate them arbitrarily due to
the natural severe ill-posedness of the problem. In the present work, the
frequency-dependent truncation is proposed, which is our novelty here. Since
the wavenumber $k$ is large in several physical applications; see our brief
discussion after Theorem 4, we can, in fact, exploit this largeness to prove
the error estimate of the scheme. By that way, we can obtain a stable
approximation for the Cauchy problem.

In view of Parseval’s identity, we now estimate that

$\displaystyle\left\|\mathbf{Q}_{1}u(x,\cdot)\right\|^{2}=4\sum_{j\in
B}e^{-2\sqrt{\lambda_{j,k}}}\lambda_{j,k}^{2}e^{2\sqrt{\lambda_{j,k}}}\left|\left\langle
u(x,\cdot),\phi_{j}\right\rangle\right|^{2}\leq 4\gamma^{-2}\sum_{j\in
B}\lambda_{j,k}^{2}e^{2\sqrt{\lambda_{j,k}}}\left|\left\langle
u(x,\cdot),\phi_{j}\right\rangle\right|^{2}.$ (18)

By using (11) obtained in Theorem 2, we have

$\displaystyle\sup_{x\in\left[0,1\right]}\sum_{j\in
B}\lambda_{j,k}^{2}e^{2\sqrt{\lambda_{j,k}}}\left|\left\langle
u(x,\cdot),\phi_{j}\right\rangle\right|^{2}$
$\displaystyle\leq\sup_{x\in\left[0,1\right]}\left[\sum_{j\in
B}\lambda_{j,k}^{2}e^{2\left(1-x\right)\sqrt{\lambda_{j,k}}}\left(\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle+\frac{1}{\sqrt{\lambda_{j,k}}}\left\langle
u_{x}\left(x,\cdot\right),\phi_{j}\right\rangle\right)^{2}\right]$
$\displaystyle\leq\sum_{j\in B}\left(\lambda_{j,k}\left\langle
u\left(1,\cdot\right),\phi_{j}\right\rangle+\sqrt{\lambda_{j,k}}\left\langle
u_{x}\left(1,\cdot\right),\phi_{j}\right\rangle\right)^{2}\leq
2\left\|u\left(1,\cdot\right)\right\|_{H^{2}\left(0,1\right)}^{2}+2\left\|u_{x}\left(1,\cdot\right)\right\|_{H^{1}\left(0,1\right)}^{2}.$

by means of $\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle\left\langle
u_{x}\left(x,\cdot\right),\phi_{j}\right\rangle\geq 0$; cf. (10) and (12). Now
we estimate $\mathbf{Q}_{2}u$ as follows. Observe that if
$\log\left(\gamma\right)\geq k$, then $\mu_{i}-k^{2}\geq k^{2}-\mu_{j}>0$ for
$i\in B$ and $j\in A_{3}$. This means that
$\lambda_{i,k}\geq\left|\lambda_{j,k}\right|$ for $i\in B$ and $j\in A_{3}$.
Therefore, we estimate that

$\displaystyle\left\|\mathbf{Q}_{2}u(x,\cdot)\right\|^{2}=4\sum_{j\in
A_{3}}\left|\lambda_{j,k}\right|^{2}\left|\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle\right|^{2}\leq 4\sum_{j\in
B}\left|\lambda_{j,k}\right|^{2}\left|\left\langle
u\left(x,\cdot\right),\phi_{j}\right\rangle\right|^{2}=\left\|\mathbf{Q}_{1}u(x,\cdot)\right\|^{2}.$
(19)

Henceforth, we can assume that the true solution satisfies $u(1,\cdot)\in
H^{2}(0,1)$ and $u_{x}(1,\cdot)\in H^{1}(0,1)$ to gain the strong convergence
of the scheme. Note now that we can rewrite $\mathbf{P}$ as

$\displaystyle\mathbf{P}u(x,\cdot)=-2k^{2}u(x,\cdot)-2\sum_{j\in\mathbb{N}\backslash(B\cup
A_{3})}\lambda_{j,k}\left\langle u(x,\cdot),\phi_{j}\right\rangle\phi_{j}$
(20)

which is computable for our numerical simulation. Moreover, since in
$\mathbb{N}\backslash(B\cup A_{3})$ it holds that
$0\leq\lambda_{j,k}\leq\log^{2}(\gamma)$, we, according to Parseval’s identity
and using (9), get that

$\displaystyle\left\|\mathbf{P}u(x,\cdot)\right\|^{2}\leq
4k^{4}\left\|u(x,\cdot)\right\|^{2}+4\log^{2}(\gamma)\left\|u_{y}(x,\cdot)\right\|^{2}.$
(21)

##  3 Convergence analysis

We now formulate theorems for our convergence analysis. Existence and
uniqueness theorems can be proven similarly to those theoretical findings in
[4]. Therefore, we only state those theorems below without giving any proof.

######  Definition 1 (Weak solution).

For each $\varepsilon>0$, a function $u^{\varepsilon}:[0,1]\to H^{1}_{0}(0,1)$
is said to be a weak solution to system (16) if

  * 1.

$u^{\varepsilon}\in C([0,1];H^{1}_{0}(0,1)),\partial_{x}u^{\varepsilon}\in
C([0,1];L^{2}(0,1)),\partial^{2}_{x^{2}}u^{\varepsilon}\in
L^{2}(0,1;(H^{1}(0,1))^{\prime})$;

  * 2.

For every test function $\psi\in H_{0}^{1}(0,1)$, it holds that

$\displaystyle\left\langle\frac{\partial^{2}u^{\varepsilon}}{\partial
x^{2}},\psi\right\rangle+\left\langle\frac{\partial u^{\varepsilon}}{\partial
y},\frac{\partial\psi}{\partial
y}\right\rangle+\left\langle\mathbf{P}u^{\varepsilon},\psi\right\rangle+k^{2}\left\langle
u^{\varepsilon},\psi\right\rangle=0\quad\text{for a.e. in }(0,1);$ (22)

  * 3.

$u^{\varepsilon}(0)=u_{0}^{\varepsilon}\in
H^{1}(0,1),\partial_{x}u^{\varepsilon}(0)=0$.

######  Theorem 3 (Existence and uniqueness of a weak solution).

For any $\varepsilon>0$, system (16) admits a unique weak solution in the
sense of Definition 1. Moreover, it holds that $u^{\varepsilon}\in
C([0,1];H_{0}^{1}(0,1))$ and $\partial_{x}u^{\varepsilon}\in
C([0,1];L^{2}(0,1))$.

######  Theorem 4 (Error estimate).

Let $u\in C([0,1];H^{2}(0,1))\cap C^{1}([0,1];H^{1}(0,1))$ be a unique
solution of the Cauchy problem of (3)–(5). Let $M>0$ independent of
$\varepsilon$ and $k$ be such that $\left\|u\right\|_{C([0,1];H^{2}(0,1))\cap
C^{1}([0,1];H^{1}(0,1))}\leq M$. Let $u^{\varepsilon}$ be a unique weak
solution of system (16) as analyzed in Theorem 3. Then, by choosing
$\log(\gamma)=2k-\eta$ for $0<\eta\leq k$, the following $L^{2}$ error
estimate holds:

$\displaystyle\left\|u^{\varepsilon}\left(x,\cdot\right)-u\left(x,\cdot\right)\right\|^{2}\leq\left(\frac{1}{k^{2}}+1\right)\varepsilon^{2}e^{\left(4k+1\right)x}+\frac{16M^{2}e^{4k(x-1)+2\eta+x}x}{k^{2}}.$
(23)

###### Proof.

Consider the difference function $w=u^{\varepsilon}-u$. It is straightforward
that this function satisfies the corresponding difference equation:

$\displaystyle w_{xx}-w_{yy}+k^{2}w=-\textbf{P}w+\textbf{Q}u\quad\text{in
}\Omega.$ (24)

Besides, this equation is associated with the Dirichlet boundary condition and
the initial conditions

$\displaystyle\begin{cases}w(x,0)=w(x,1)=0&\quad\text{for }x\in(0,1),\\\
w(0,y)=u_{0}^{\varepsilon}(y)-u_{0}(y),\;w_{x}(0,y)=0&\quad\text{for
}y\in(0,1).\end{cases}$ (25)

Multiply both sides of (24) by $w_{x}$ and integrate the resulting equation
with respect to $y$ from 0 to $1$. After that, dividing both sides of the
resulting equation by $k^{2}$ gives

$\displaystyle\frac{1}{k^{2}}\frac{d}{dx}\left\|w_{x}\left(x,\cdot\right)\right\|^{2}+\frac{1}{k^{2}}\frac{d}{dx}\left\|w_{y}\left(x,\cdot\right)\right\|^{2}+\frac{d}{dx}\left\|w\left(x,\cdot\right)\right\|^{2}=\frac{2}{k^{2}}\left\langle\mathbf{Q}u\left(x,\cdot\right),w_{x}\left(x,\cdot\right)\right\rangle-\frac{2}{k^{2}}\left\langle\mathbf{P}w\left(x,\cdot\right),w_{x}\left(x,\cdot\right)\right\rangle.$

Using the energy estimate (21) of the stabilized operator in combination of
the Cauchy–Schwarz inequality, we estimate that

$\displaystyle-\frac{2}{k^{2}}\left\langle\mathbf{P}w\left(x,\cdot\right),w_{x}\left(x,\cdot\right)\right\rangle$
$\displaystyle\leq\frac{1}{k^{2}}\left[\frac{1}{4k}\left\|\mathbf{P}w\left(x,\cdot\right)\right\|^{2}+4k\left\|w_{x}\left(x,\cdot\right)\right\|^{2}\right]$
$\displaystyle\leq\frac{k^{4}}{k^{3}}\left\|w\left(x,\cdot\right)\right\|^{2}+\frac{\log^{2}\left(\gamma\right)}{k^{3}}\left\|w_{y}\left(x,\cdot\right)\right\|^{2}+\frac{4}{k}\left\|w_{x}\left(x,\cdot\right)\right\|^{2}$

In the same vein, using (17), (18) and (19), we obtain the following estimate:

$\frac{2}{k^{2}}\left\langle\mathbf{Q}_{\varepsilon}u\left(x,\cdot\right),w_{x}\left(x,\cdot\right)\right\rangle\leq\frac{16\gamma^{-2}}{k^{2}}\left\|u\right\|_{C\left(\left[0,1\right];H^{2}(0,1)\right)\cap
C^{1}\left([0,1];H^{1}(0,1)\right)}^{2}+\frac{1}{k^{2}}\left\|w_{x}\left(x,\cdot\right)\right\|^{2}$

Therefore, $\log(\gamma)=2k-\eta$ where $0<\eta\leq k$, we have

$\displaystyle\frac{1}{k^{2}}\frac{d}{dx}\left\|w_{x}\left(x,\cdot\right)\right\|^{2}+\frac{1}{k^{2}}\frac{d}{dx}\left\|w_{y}\left(x,\cdot\right)\right\|^{2}+\frac{d}{dx}\left\|w\left(x,\cdot\right)\right\|^{2}$
$\displaystyle\leq\frac{16\gamma^{-2}}{k^{2}}M^{2}+\frac{4k+1}{k^{2}}\left\|w_{x}\left(x,\cdot\right)\right\|^{2}+k\left\|w\left(x,\cdot\right)\right\|^{2}+\frac{\log^{2}\left(\gamma\right)}{k^{3}}\left\|w_{y}\left(x,\cdot\right)\right\|^{2}$
$\displaystyle\leq\frac{16e^{-4k+2\eta}}{k^{2}}M^{2}+\frac{4k+1}{k^{2}}\left\|w_{x}\left(x,\cdot\right)\right\|^{2}+k\left\|w\left(x,\cdot\right)\right\|^{2}+\frac{(2k-\eta)^{2}}{k^{3}}\left\|w_{y}\left(x,\cdot\right)\right\|^{2}$

By Gronwall’s inequality and using (13) and (25), we, after dropping the
gradient norms on the left-hand side, find that

$\displaystyle\left\|w\left(x,\cdot\right)\right\|^{2}$
$\displaystyle\leq\left(\frac{1}{k^{2}}\left\|w_{x}\left(0,\cdot\right)\right\|^{2}+\frac{1}{k^{2}}\left\|w_{y}\left(0,\cdot\right)\right\|^{2}+\left\|w\left(0,\cdot\right)\right\|^{2}+\frac{16M^{2}e^{-4k+2\eta}x}{k^{2}}\right)e^{\left(4k+1\right)x}$
$\displaystyle\leq\left[\left(\frac{1}{k^{2}}+1\right)\varepsilon^{2}+\frac{16M^{2}e^{-4k+2\eta}x}{k^{2}}\right]e^{\left(4k+1\right)x}.$

This completes the proof of the theorem. ∎

From (23) in Theorem 4, if we only measure the data once, i.e. only one
$\varepsilon\in(0,1)$ is attained, then we can still obtain a fine
approximation to $u$ by collecting data for a set of multiple high
frequencies. In this case, we consider $\underline{k}\leq k\leq\overline{k}$,
where $\underline{k}>1$, and thus, $u\left(x,y\right)$ solutions to the Cauchy
problem is understood as $u\left(x,y,k\right)$. The perspective of multiple
frequencies is relevant to the inverse scattering problem for detection of
suspicious explosive-like objects; cf. e.g. [10], when the near-field
scattering is approximated by the far-field measurement. Moreover, as
practically shown in [10, 11], the wavenumber $k$ is large. Now, assume
$\varepsilon$ fixed such that $\overline{k}^{2}e^{4\overline{k}}\leq
K\varepsilon^{-2}$ for $K>0$, then $u^{\varepsilon}$ is close to $u$ in
$L^{2}$ sense with an error of $\mathcal{O}(k^{-2})$.

Similarly, if we fix a set of $k$ such that $4\underline{k}\geq
K\log\left(\varepsilon^{-\alpha}\right)$ for $\alpha>0$ and $K>0$, then
$u^{\varepsilon}$ is close to $u$ in $L^{2}$ sense with an error of
$\mathcal{O}\left(\varepsilon^{\min\left\\{\alpha,2\right\\}}\right)$. It is
also easy to deduce from the above proof that we can obtain a Lipschitz
stability of the scheme in terms of $\varepsilon$, as $\mathbf{Q}$ is not
involved.

##  4 Numerical examples

Given $M,N\in\mathbb{N}$, we consider uniform grids of mesh-points
$x_{m}=m\Delta x,y_{n}=n\Delta y$ for $m\leq M,n\leq N$ with $\Delta x,\Delta
y$ being the mesh-widths in $x$ and $y$, respectively. For any function
$u(x,y)$, we denote by $u_{m,n}\approx u(x_{m},y_{n})$ the corresponding
discrete function. To generate the data, we apply the central finite
difference method (FDM) to solve the Helmholtz equation (3) with the Dirichlet
boundary conditions imposed on four sides of $\Omega=[0,1]^{2}$, viz.

$\displaystyle u(0,y)=u_{0}(y),\quad u(1,y)=g(y),\quad u(x,0)=0,\quad
u(x,1)=0.$ (26)

In our numerical performance of the stabilization scheme below, we do not
choose the true solution of the Helmholtz equation (3). Instead, we choose its
boundary data $u_{0},g$ in (26) so that our choice is more flexible. This is
relevant because (3) with full data (26) is a well-posed problem and the
central FDM is well known to be stable and convergent with respect to the
refinement of $x$ and $y$. In this circumstance, one can consider the discrete
function $u_{m,n}$ obtained from that well-posed problem as a reliable true
solution. The Neumann data $u_{1}$ in (5) can be generated using the fact that

$u_{0}(y_{n})\approx u_{1}(y_{n})\Delta x+u(x_{1},y_{n}).$

The same FDM is applied when we solve $U(x,y)$ of system (6). Thereby, a
stable approximation of $U(0,y)$ can be obtained for the Dirichlet data in
(7). To solve for $V(x,y)$ in (7) numerically, we remark that the stabilized
operator $\mathbf{P}$ (cf. (20)) needs to be linearized, albeit it is a linear
mapping. Cf. [4] and based on (16), the linearized stabilized problem for (7)
is established as follows. We construct a sequence of
$\left\\{V^{\varepsilon,q}\right\\}_{q\in\mathbb{N}}$ satisfying

$\displaystyle\begin{cases}V_{xx}^{\varepsilon,q+1}-V_{yy}^{\varepsilon,q+1}+\mathbf{P}_{1}V^{\varepsilon,q}-k^{2}V^{\varepsilon,q+1}=0&\text{in
}\Omega,\\\
V^{\varepsilon,q+1}\left(x,0\right)=V^{\varepsilon,q+1}(x,1)=0&\text{for
}x\in\left(0,1\right),\\\
V^{\varepsilon,q+1}\left(0,y\right)=u_{0}^{\varepsilon}\left(y\right)-U\left(0,y\right),\quad
V_{x}^{\varepsilon,q+1}\left(0,y\right)=0&\text{for
}y\in\left(0,1\right).\end{cases}$ (27)

In (27), we use

$\displaystyle\mathbf{P}_{1}u(x,\cdot)=-2\sum_{j\in\mathbb{N}\backslash(B\cup
A_{3})}\lambda_{j,k}\left\langle u(x,\cdot),\phi_{j}\right\rangle\phi_{j},$

and the initial guess $V^{\varepsilon,0}$ (i.e. $q=0$) is chosen to be
$u_{0}^{\varepsilon}(y)-U(0,y)$. We choose this initial guess because it is a
unique function that contains many information of our sought $V^{\varepsilon}$
under stabilization. As to the measured data $u_{0}^{\varepsilon}$, we apply
the additive noise in the following sense:
$u_{0}^{\varepsilon}(y)=u_{0}(y)+\varepsilon\text{rand}(y)$, where rand is a
uniformly distributed random number such that
$\max_{y\in[0,1]}\left|\text{rand}(y)\right|\leq 1$. At the discretization
level, the gradient of $u_{0}^{\varepsilon}$ is then approximated by

$\partial_{y}u_{0}^{\varepsilon}(y_{n})\approx\frac{u_{0}^{\varepsilon}(y_{n+1})-u_{0}^{\varepsilon}(y_{n})}{\Delta
y}\approx\partial_{y}u_{0}(y_{n})+\varepsilon
N\left(\text{rand}(y_{n+1})-\text{rand}(y_{n})\right).$

Therefore, we can see that to fulfill assumption (13), having the measured
data $u_{0}^{\varepsilon}$ is practically sufficient if we choose $N$ such
that either $N\leq\varepsilon^{-\beta}$ or $N\leq k^{\beta}$ holds for
$\beta\in(0,1)$. Subsequently, a coarse mesh for the domain of $y$ should be
used.

For the (local) convergence of the linearization scheme for (27), we refer the
reader to the pioneering work [4], where the same convergence analysis is
studied. Since the convergence is local in $x$, we use a fine mesh for this
variable. We do not detail how small $\Delta x$ should be as it can be done as
in [4]. Generally, it requires that $\log(\gamma)\Delta x\gamma^{\Delta x}<1$.

For each iteration $q$, we apply the same central FDM to approximate
$V^{\varepsilon,q}$. The inner product in $\mathbf{P}_{1}$ can be approximated
by the standard Riemann sum. Henceforth, we can obtain the fully discrete
version of (27) in matrix form. We refer the reader to [4] for the standard
derivation of this form. Note that after having $V^{\varepsilon,q}$, we obtain
an approximation of $u^{\varepsilon}$ via
$u^{\varepsilon}_{m,n}=U_{m,n}+V^{\varepsilon,q}_{m,n}$.

Now, in all examples below, we choose $\eta=k,M=400,N=80$ and $q=1$. To show
how good our proposed scheme is, all numerical results are reported with
$\varepsilon=0.99,0.01$. Often, the noise level is not too small in practice.
Moreover, it can be seen that our scheme is not much expensive as $q$ is
relatively small. Below, we consider the following relative error:

$E=\frac{\sqrt{\sum_{m=0}^{M}\sum_{n=0}^{N}\left|u_{m,n}^{\varepsilon}-u_{m,n}\right|^{2}}}{\sqrt{\sum_{m=0}^{M}\sum_{n=0}^{N}\left|u_{m,n}\right|^{2}}}\times
100\%.$

### Example 1: Low frequency

In this test, we choose $k=5$ and

$u_{0}(y)=-e^{-2\left(0.5^{4}+(y-0.5)^{4}\right)}+0.5^{4}+(y-0.5)^{4},\quad
g(y)=u_{0}(y).$

Such $k$ is suitable in the context of landmine detection; cf. e.g. [10].
Depicted in Figure 1 are the graphical illustrations of the true solution and
its reconstructed with large noise ($\varepsilon=0.99$) and intermediate noise
($\varepsilon=0.1$). When $\varepsilon$ is smaller, the computed solution is
very close to the true one in terms of the value and, furthermore, the shape
and location of the yellow circular protrusion; see Figures 1a and 1c. On the
other hand, the relative error reduces from 30.00% for $\varepsilon=0.99$ to
1.91% for $\varepsilon=0.1$, which shows the efficiency of the scheme. In this
example, it is concrete that choosing a frequency-dependent stabilization
works very well in the reconstruction process, compared to the noise-dependent
choice in the previous works.

(a) True

(b) Computed ($\varepsilon=0.99$)

(c) Computed ($\varepsilon=0.1$)

Figure 1: Numerical results of Example 1 (low-frequency problem). (a)
Graphical illustration of the true solution. (b) and (c) Illustrations of the
reconstructed solution with $\varepsilon=0.99$ and $\varepsilon=0.1$,
respectively.

### Example 2: High frequency

In this test, we consider $k=50$ and

$u_{0}(y)=0,\quad
g(y)=\frac{-\sin\left(7\sqrt{0.001+(y-0.5)^{2}}\right)}{7\sqrt{1+(y-0.5)^{2}}}.$

Similar to Example 1, we can see the reconstruction becomes better when
$\varepsilon$ decreases from 0.99 to 0.1; see Figure 2. Especially, when
$\varepsilon=0.1$, the computed solution, cf. Figure 2c shows exactly the same
shape and location of all yellow bands in the true solution (Figure 2). In
Figure 2b for $\varepsilon=0.99$, those bands do not have the correct shape.
Moreover, the illustration when $\varepsilon=0.1$ is much clearer than that of
$\varepsilon=0.99$. This observation is exactly the same as in Example 1; see
Figures 1b and 1c. We also report that the relative error $E$ in this case
reduces significantly from 183.38% to 16.18%, when $\varepsilon$ decreases
from 0.99 to 0.1.

(a) True

(b) Computed ($\varepsilon=0.99$)

(c) Computed ($\varepsilon=0.1$)

Figure 2: Numerical results of Example 1 (high-frequency problem). (a)
Graphical illustration of the true solution. (b) and (c) Illustrations of the
reconstructed solution with $\varepsilon=0.99$ and $\varepsilon=0.1$,
respectively.

## Acknowledgment

V.A.K. thanks Prof. Dr. Roselyn Williams (Tallahassee, USA) for her support
during the time V.A.K working at Florida A&M University.

## References

  * [1] M. Born, E. Wolf, Principles of Optics, Cambridge University Press, 2019.
  * [2] N. H. Tuan, V. A. Khoa, M. N. Minh, T. Tran, Reconstruction of the electric field of the Helmholtz equation in three dimensions, Journal of Computational and Applied Mathematics 309 (2017) 56–78.  doi:10.1016/j.cam.2016.05.021.
  * [3] M. V. Klibanov, D.-L. Nguyen, L. H. Nguyen, A coefficient inverse problem with a single measurement of phaseless scattering data, SIAM Journal on Applied Mathematics 79 (1) (2019) 1–27.  doi:10.1137/18m1168303.
  * [4] V. A. Khoa, P. T. H. Nhan, Constructing a variational quasi-reversibility method for a Cauchy problem for elliptic equations, Mathematical Methods in the Applied Sciences 44 (5) (2020) 3334–3355.  doi:10.1002/mma.6945.
  * [5] H. T. Nguyen, V. A. Khoa, V. A. Vo, Analysis of a quasi-reversibility method for a terminal value quasi-linear parabolic problem with measurements, SIAM Journal on Mathematical Analysis 51 (1) (2019) 60–85.  doi:10.1137/18m1174064.
  * [6] P. T. Nam, D. D. Trong, N. H. Tuan, The truncation method for a two-dimensional nonhomogeneous backward heat problem, Applied Mathematics and Computation 216 (2010) 3423–3432.
  * [7] T. Regińska, K. Regiński, Approximate solution of a Cauchy problem for the Helmholtz equation, Inverse Problems 22 (3) (2006) 975–989.  doi:10.1088/0266-5611/22/3/015.
  * [8] N. H. Tuan, L. D. Thang, V. A. Khoa, T. Tran, On an inverse boundary value problem of a nonlinear elliptic equation in three dimensions, Journal of Mathematical Analysis and Applications 426 (2015) 1232–1261.
  * [9] N. H. Tuan, V. A. Khoa, P. T. K. Van, V. V. Au, An improved quasi-reversibility method for a terminal-boundary value multi-species model with white Gaussian noise, Journal of Computational and Applied Mathematics 384 (2021) 113176\.  doi:10.1016/j.cam.2020.113176.
  * [10] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. H. Nguyen, A. J. Sullivan, V. N. Astratov, Convexification and experimental data for a 3D inverse scattering problem with the moving point source, Inverse Problems 36 (8) (2020) 085007.  doi:10.1088/1361-6420/ab95aa.
  * [11] V. L. R. Murthy, S. V. J. Lakshman, Electronic absorption spectrum of cobalt antipyrine complex, Solid State Communications 38 (7) (1981) 651–652.  doi:10.1016/0038-1098(81)90960-1.