# Vibrational Control of Cluster Synchronization:
Connections with Deep Brain Stimulation

Yuzhen Qin, Danielle S. Bassett, and Fabio Pasqualetti This work was supported
in part by awards ARO-W911NF1910360 and NSF-NCS-FO-1926829. Y. Qin and F.
Pasqualetti are with the Department of Mechanical Engineering, University of
California at Riverside, {yuzhenqin<EMAIL_ADDRESS>D. S. Bassett is
with the Department of Bioengineering, the Department of Electrical & Systems
Engineering, the Department of Physics & Astronomy, the Department of
Psychiatry, and the Department of Neurology, University of Pennsylvania, and
The Santa Fe Institute<EMAIL_ADDRESS>
###### Abstract

Cluster synchronization underlies various functions in the brain. Abnormal
patterns of cluster synchronization are often associated with neurological
disorders. Deep brain stimulation (DBS) is a neurosurgical technique used to
treat several brain diseases, which has been observed to regulate neuronal
synchrony patterns. Despite its widespread use, the mechanisms of DBS remain
largely unknown. In this paper, we hypothesize that DBS plays a role similar
to _vibrational control_ since they both highly rely on high-frequency
excitation to function. Under the framework of Kuramoto-oscillator networks,
we study how vibrations introduced to network connections can stabilize
cluster synchronization. We derive some sufficient conditions and also provide
an effective approach to design vibrational control. Also, a numerical example
is presented to demonstrate our theoretical findings.

##  I INTRODUCTION

Cluster synchronization describes a phenomenon where the behavior of the units
in a network evolves into different, yet synchronized, clusters. It has been
widely observed in the brain as correlated neural activity. Patterns of
cluster synchronization underlie various brain functions such as neuronal
communication, memory formation, and cognition [1, 2]. Pathological synchrony
patterns also characterize many brain disorders, e.g., Parkinson’s disease [3]
and epilepsy [4].

DBS is a therapeutic technique that uses intracranial electrodes and high-
frequency electrical stimulation (see Fig. 1 (a)) to treat such brain
disorders [5]. Despite its ubiquitous use, the underlying mechanisms of DBS
are still elusive. Selecting the proper dose of stimulation and the best
location in the brain, which is crucial for achieving the greatest clinical
benefits of DBS, still largely relies on trial and error [5]. It is thus
critical to further characterize the mechanisms that underpin DBS’s efficacy,
which, in turn, can inform the design of optimal and personalized brain
stimulation.

In this paper, we provide a hypothetical mechanism for DBS from the
perspective of control systems. We propose that high-frequency signals
delivered by DBS to brain regions work as vibrational control. Vibrational
control is an open-loop control strategy that has been used to stabilize
engineering systems, e.g., inverted pendulums, chemical reactors, and under-
actuated robots (see [6, 7, 8, 9] and the references therein). Similar to DBS,
vibrational control relies on high-frequency dithers. It is particularly
useful in situations where feedback control is not applicable since online
measurement of systems’ states or outputs is unavailable. Using the Kuramoto-
oscillator framework, we study how vibrational control can stabilize cluster
synchronization. Despite its simplicity, the Kuramoto model has been widely
used to study synchronization in neural systems [10].

Related work. Cluster synchronization has recently attracted extensive
attention. Some studies have found that cluster synchronization is closely
related to network symmetries [11, 12, 13] and equitable partitions [14].
Stability conditions for cluster synchronization are constructed in networks
with dyadic [15, 16, 17, 18] and hyper connections [19, 20]. Different control
strategies, such as pinning control [21], and intervention of network
connections or nodes’ dynamics [22, 23, 24], are proposed to control cluster
synchronization. By contrast, we employ vibrational control, a much more
realistic control strategy especially for neural oscillation regulation since
it does not need to modify structural brain networks or the intrinsic dynamics
of neurons.

Figure 1: Vibrational control of cluster synchronization inspired by DBS. (a)
DBS consists of an intracranial electrode, an extension wire, and a pulse
generator. (b) DBS predominantly affects axons and dendrites rather than the
soma. (c) Vibrations are introduced to network connections to stabilize
cluster synchronization (different colors represent different clusters).

Paper contribution. The contribution of this paper is threefold. First, we use
vibration control, a potential mechanism of DBS, to regulate cluster
synchronization. Different from the existing computational studies on DBS
(e.g., see [25]), we introduce control inputs to network connections rather
than nodes’ dynamics (see Fig. 1 (c)). This is inspired by the observation
that brain stimulation predominantly affects dendrites and axons near the
electrode, rather than the soma [26] (see Fig. 1 (b)). Second, we derive
sufficient conditions on the vibrational inputs such that they can stabilize
originally unstable cluster synchronization. The key idea is that vibrations
are required improve the robustness of synchronization within clusters to a
level that is sufficient to overcome perturbations resulting from inter-
cluster connections. Averaging methods and perturbation theory are used for
the analysis. Third, we propose an analytically tractable approach to design
vibrational control. Finally, we provide an example to demonstrate our
theoretical findings.

Notation. We denote the unit circle as $\mathbb{S}^{1}$, a point of which is a
phase. The $n$-dimensional torus is denoted by
$\mathbb{T}^{n}=\mathbb{S}\times\dots\times\mathbb{S}$. Given a set
$\mathcal{S}$, $|\mathcal{S}|$ stands for its cardinality. For two sets
$\mathcal{S}_{1}$ and $\mathcal{S}_{2}$,
$\mathcal{S}_{1}\backslash\mathcal{S}_{2}$ is their set-theoretic difference.
Given $A\in\mathbb{R}^{n\times m}$, $B:=[A]^{+}$ is the matrix where
$b_{ij}=a_{ij}$ if $a_{ij}>0$ and $b_{ij}=0$ otherwise. Given matrices
$A_{1},\dots,A_{r}$, ${\rm blkdiag}(A_{1},\dots,A_{r})$ is the block-diagonal
matrix of them.

##  II Problem Formulation and Preliminary

###  II-A Problem formulation

We consider $n$ oscillators coupled by an undirected network
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ associated with the weighted adjacency
matrix $A=A^{\top}=[a_{ij}]_{n\times n}$. The dynamics of the oscillators are
governed by the following Kuramoto model:

$\dot{\theta}_{i}=\omega_{i}+\sum_{j=1}^{n}a_{ij}\sin(\theta_{j}-\theta_{i}),$
(1)

where $\theta_{i}\in\mathbb{S}^{1}$ is the $i$th oscillator’s phase, and
$\omega_{i}\in\mathbb{R}$ is its natural frequency. Let
$\theta:=[\theta_{1},\dots,\theta_{n}]^{\top}$ and
$\omega:=[\omega_{i},\dots,\omega_{n}]^{\top}$.

Since we are interested in studying cluster synchronization, we define cluster
synchronization manifolds as follows.

######  Definition 1 (Synchronization manifold).

For the the network $\mathcal{G}=(\mathcal{V},\mathcal{E})$, consider the
partition
$\mathcal{C}:=\\{\mathcal{C}_{1},\mathcal{C}_{2},\dots,\mathcal{C}_{r}\\}$,
where $\mathcal{C}_{k}\subset\mathcal{V}$,
$\mathcal{C}_{k}\cap\mathcal{C}_{\ell}=\emptyset$ for any $k\neq\ell$, and
$\cup_{k=1}^{r}\mathcal{C}_{k}=\mathcal{V}$. The cluster synchronization
manifold associated with the partition $\mathcal{C}$ is defined as

$\mathcal{M}:=\\{\theta\in\mathbb{T}^{n}:\theta_{i}=\theta_{j},\forall
i,j\in\mathcal{C}_{k},k=1,\dots,r\\}.\hskip 12.0pt\hbox{}\hfill\triangle$

The cluster synchronization manifold is invariant along the system (1) if
starting from any $\theta(0)\in\mathcal{M}$, the solution to (1) satisfies
$\theta(t)\in\mathcal{M}$ for all $t\geq 0$. Following [15, 27], we make the
following assumption to ensure that $\mathcal{M}$ is invariant.

######  Assumption 1 (Invariance).

For $k=1,2,\dots,r$: i) the natural frequencies satisfy
$\omega_{i}=\omega_{j}$ for any $i,j\in\mathcal{C}_{k}$; and ii) the coupling
strengths satisfy that, for any $\ell\in\\{1,2,\dots,r\\}\backslash\\{k\\}$,
$\sum_{p\in\mathcal{C}_{\ell}}(a_{ip}-a_{jp})=0$ for any
$i,j\in\mathcal{C}_{k}$.  $\triangle$

In addition to Assumption 1, to guarantee that the cluster synchronization
represented by $\mathcal{M}$ can appear in the network, $\mathcal{M}$ is
required to be stable. Further conditions combining the connection weights and
natural frequencies need to be satisfied to ensure the stability of
$\mathcal{M}$ (e.g., some sufficient conditions are constructed in [15, 16,
14]).

Yet, changes on the anatomical connectivity or the intrinsic dynamics of
neuronal populations, caused by brain aging or disorders [28], may violate
such conditions so that cluster synchronization patterns needed for normal
brain functions are no longer stable. In this paper, we aim to investigate how
vibrational control, a control strategy that resembles DBS [5], can restore
the stability of such desired synchrony patterns.

Specifically, we consider the following controlled model

$\scalebox{1.0}{\mbox{$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sum_{j=1}^{n}\Big{(}a_{ij}+q(a_{ij})u_{ij}(t)\Big{)}\sin(\theta_{j}-\theta_{i})$}},$
(2)

where $q(a_{ij})=1$ if $a_{ij}>0$ and $q(a_{ij})=0$, otherwise111This setting
ensures that control inputs can only introduce vibrations to the existing
connections and do not unreasonably create new ones., and $u_{ij}(t)$ is the
corresponding vibrational control input. Although the couplings are symmetric,
satisfying $a_{ij}=a_{ji}$ for any $i$ and $j$, we allow the effect of the
vibrational control on each edge to be asymmetric, that is, it is allowed that
$u_{ij}(t)\neq u_{ji}(t)$. As typically done in the literature (e.g., [8, 7,
6]), we assume that $u_{ij}(t)$ is periodic with period $T$ and has high
frequency and zero mean. It can be seen that
$\int_{0}^{T}a_{ij}+q(a_{ij})u_{ij}(t)dt=a_{ij}$, which means that the effect
of vibrations on network weights is zero in average. Particularly, we consider
sinusoidal vibrations in this paper, i.e.,

$\displaystyle
u_{ij}(t)=\frac{u_{ij}}{\varepsilon}\sin\Big{(}\frac{t}{\varepsilon}\Big{)},$

where $\varepsilon>0$ determines the amplitude and frequency of the
vibrations. Then, the period becomes $T=2\varepsilon\pi$.

It remains to design the parameters of the vibrational control, i.e.,
$u_{ij}$’s and $\varepsilon$, to stabilize the cluster synchronization
manifold $\mathcal{M}$. To this end, it is necessary to ensure that the
vibrational inputs preserve the invariance of $\mathcal{M}$. In this paper, we
consider a special configuration of vibrations that satisfies

$\displaystyle{u_{ij}(t)=0,}$ $\displaystyle\forall
i\in\mathcal{C}_{k},j\in\mathcal{C}_{\ell},\forall k\neq\ell,$

so that the invariance of $\mathcal{M}$ is preserved. In other words,
vibrations are only introduced to the connections between oscillators within
the same clusters (such connections are referred to as intra-cluster
connections).

###  II-B Preliminary

To facilitate analysis in the reminder of this paper, we introduce some
notations (summarized in Fig. 2) and derive a compact model for the system
(2).

Figure 2: Graphic illustration of the main notations we employ in this paper:
different subgraphs and their corresponding incidence matrices.

For the network $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and a partition
$\mathcal{C}:=\\{\mathcal{C}_{1},\mathcal{C}_{2},\dots,\mathcal{C}_{r}\\}$,
define $\mathcal{G}_{k}=(\mathcal{C}_{k},\mathcal{E}_{k})$ where
$\mathcal{E}_{k}:=\\{(i,j)\in\mathcal{E}:i,j\in\mathcal{C}_{k}\\}$. For each
$k$, denote $n_{k}:=|\mathcal{C}_{k}|$ as the number of nodes in
$\mathcal{G}_{k}$. We assume that each $\mathcal{G}_{k}$ is connected and
contains at least 2 nodes. Let $\mathcal{G}_{\rm
intra}=(\mathcal{V},\mathcal{E}_{\rm intra})=\cup_{k=1}^{r}\mathcal{G}_{k}$
and we refer to it as the intra-cluster graph. Let $\mathcal{G}_{{\rm
inter}}=(\mathcal{V},\mathcal{E}_{\rm inter})$ be the inter-cluster subgraph,
where $\mathcal{E}_{\rm inter}:=\mathcal{E}\backslash\mathcal{E}_{\rm intra}$.
Define the directed graph associated with $\mathcal{G}_{\rm intra}$ by
$\mathcal{G}_{\rm dir}=(\mathcal{V},\mathcal{E}_{\rm dir})$, where
$\mathcal{E}_{\rm dir}$ is defined in a way such that every edge in
$\mathcal{G}_{\rm intra}$ becomes two directed edges in $\mathcal{E}_{\rm
dir}$. Let $\mathcal{G}_{{\rm dir}1}=(\mathcal{V},\mathcal{E}_{{\rm dir}1})$
be the graph where $\mathcal{E}_{{\rm dir}1}$ contains only one directed edge
between each pair of nodes in $\mathcal{E}_{\rm dir}$. Let $\mathcal{G}_{{\rm
dir}2}=(\mathcal{V},\mathcal{E}_{{\rm dir}2})$ with $\mathcal{E}_{{\rm
dir}2}:=\mathcal{E}_{{\rm dir}}\backslash\mathcal{E}_{{\rm dir}1}$. Let
$\hat{\mathcal{G}}=(\mathcal{V},\hat{\mathcal{E}})$ be any spanning tree of
$\mathcal{G}$, and thus $|\hat{\mathcal{E}}|=n-1$. For each $k$, Let
$\hat{\mathcal{G}}_{k}=(\mathcal{C}_{k},\hat{\mathcal{E}}_{k})$ with
$\hat{\mathcal{E}}_{k}:=\hat{\mathcal{E}}\cap\mathcal{E}_{k}$. The intra-
cluster subgraph of the spanning tree $\hat{\mathcal{G}}$ is
$\hat{\mathcal{G}}_{\rm intra}=(\mathcal{V},\hat{\mathcal{E}}_{\rm
intra}):=\hat{\cup}_{k=1}^{r}\hat{\mathcal{G}}_{k}$; the inter-cluster
subgraph of $\hat{\mathcal{G}}$ is $\hat{\mathcal{G}}_{\rm
inter}=(\mathcal{V},\hat{\mathcal{E}}_{\rm inter})$ with
$\hat{\mathcal{E}}_{\rm
inter}:=\hat{\mathcal{E}}\backslash\hat{\mathcal{E}}_{\rm intra}$.

We further define the corresponding (oriented) incidence matrices of the above
graphs and subgraphs as in Fig. 2. Without loss of generality, we order the
columns of these incidence matrices in a way such that

$\displaystyle\begin{matrix}[l]B=[B_{\rm intra},B_{\rm inter}],&B_{\rm
intra}={\rm blkdiag}(B^{(1)}_{{\rm intra}},\dots,B^{(r)}_{{\rm intra}})\\\
\hat{B}=[\hat{B}_{\rm intra},\hat{B}_{\rm inter}],&\hat{B}_{\rm intra}={\rm
blkdiag}(\hat{B}^{(1)}_{{\rm intra}},\dots,\hat{B}^{(r)}_{{\rm intra}}),\\\
B_{\rm dir}=[B_{{\rm dir}1},B_{{\rm dir}2}]&B_{{\rm dir}1}=B_{\rm
intra},\hskip 12.0ptB_{{\rm dir}2}=-B_{\rm intra}.\end{matrix}$

(3)

Let $W=\scalebox{0.7}{\mbox{$\displaystyle\begin{bmatrix}W_{\rm intra}&0\\\
0&W_{\rm inter}\end{bmatrix}$}}$, where  $\displaystyle W_{\rm intra}:={\rm
diag}\\{a_{ij},(i,j)\in\mathcal{G}_{\rm intra}\\}$ and $W_{\rm inter}:={\rm
diag}\\{a_{ij},(i,j)\in\mathcal{G}_{\rm inter}\\}$ are diagonal weight
matrices of $\mathcal{G}_{\rm intra}$ and $\mathcal{G}_{\rm inter}$,
respectively. Similarly, since vibrations are only introduced to the intra-
cluster connections and can have asymmetric influence on each of them, let
$U(t):={\rm diag}\\{u_{ij}\sin(t),(i,j)\in\mathcal{E}_{\rm dir}\\}$. It holds
that $U(t)={\rm blkdiag}(U_{1}(t),U_{2}(t))$ for $U_{1}(t):={\rm
diag}\\{u_{ij}\sin(t),(i,j)\in\mathcal{E}_{{\rm dir}1}\\}$ and $U_{2}(t):={\rm
diag}\\{u_{ij}\sin(t),(i,j)\in\mathcal{E}_{{\rm dir}2}\\}$. Further, $U_{1}$
and $U_{2}$ can be decomposed as $U_{1}(t)={\rm
blkdiag}(U^{(1)}_{1}(t),\dots,U^{(r)}_{1}(t))$ and $U_{2}(t)={\rm
blkdiag}(U^{(1)}_{2}(t),\dots,U^{(r)}_{2}(t))$, where $U^{(k)}_{1}(t))$ and
$U^{(k)}_{2}(t))$ represent the vibrations to the $k$th cluster.

Let $x:=\hat{B}_{\rm intra}^{\top}\theta$ and $y:=\hat{B}_{\rm
inter}^{\top}\theta$, and they define the intra- and inter-cluster phase
differences, respectively. The dynamics of $x$ and $y$ are (see Appendix -A
for the derivation and the expressions of the matrices $R_{1},R_{2}$, and
$R_{3}$)

$\displaystyle\dot{x}=f_{\rm intra}(x)+f_{\rm
inter}(x,y)+{\small\begin{matrix}\frac{1}{\varepsilon}\end{matrix}}f_{\rm
ctr}\big{(}U\big{(}{\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}\big{)},x\big{)},$
(4a)
$\displaystyle\dot{y}=g(x,y)+{\small\begin{matrix}\frac{1}{\varepsilon}\end{matrix}}g_{\rm
ctr}(U({\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}),x),$ (4b)

where

$\displaystyle\scalebox{0.8}{\mbox{$\displaystyle f_{\rm
intra}(x)=-\hat{B}_{\rm intra}^{\top}B_{\rm intra}W_{\rm
intra}\sin(R_{1}x)$}},$ (5) $\displaystyle\scalebox{0.8}{\mbox{$\displaystyle
f_{\rm inter}(x,y)=-\hat{B}_{\rm intra}^{\top}B_{\rm inter}W_{\rm
inter}\sin(R_{2}x+R_{3}y)$}},$
$\displaystyle\scalebox{0.8}{\mbox{$\displaystyle f_{\rm
ctr}(U,x)=-\hat{B}_{\rm intra}^{\top}\Big{(}[B_{{\rm
intra}}]^{+}U_{1}(t)-[-B_{{\rm intra}}]^{+}U_{2}(t)\Big{)}\sin(R_{1}x)$}}$
$\displaystyle\scalebox{0.8}{\mbox{$\displaystyle g(x,y)=\hat{B}_{\rm
inter}^{\top}\omega-\hat{B}_{\rm inter}^{\top}B_{\rm intra}W_{\rm
inter}\sin(R_{1}x)$}}$ $\displaystyle\hskip
56.9055pt\scalebox{0.8}{\mbox{$\displaystyle-\hat{B}_{\rm inter}^{\top}B_{\rm
inter}W_{\rm inter}\sin(R_{2}x+R_{3}y)$}},$
$\displaystyle\scalebox{0.8}{\mbox{$\displaystyle g_{\rm
ctr}(U,x)=-\hat{B}_{\rm inter}\Big{(}[B_{{\rm intra}}]^{+}U_{1}-[-B_{{\rm
intra}}]^{+}U_{2}\Big{)}\sin(R_{1}x)$}},$

Note that $f_{\rm intra}$, $f_{\rm inter}$, and $f_{\rm ctr}$ describe the
dynamics induced by the intra connections, inter connections, and vibrational
control inputs, respectively. Here, we single out the important properties of
these functions that will be used later: i) $f_{\rm intra}(0)=0$, ii) $f_{\rm
inter}(0,y)=0$ for any $y$, iii) $f_{\rm ctr}(0,x)=0$ for any $x$, and (iv)
$f_{\rm ctr}(U,0)=0$ for any $U$.

Notice that $x=0$ corresponds to the cluster synchronization manifold
$\mathcal{M}$. The manifold $\mathcal{M}$ is exponentially stable along the
system (2) if $x=0$ is exponentially stable uniformly222The equilibrium $x=0$
is said to be exponentially stable uniformly in $y$ if it is exponentially
stable starting from any $y(0)$ [29, Chap. 4]. in $y$ along (4). To stabilize
the cluster synchronization manifold $\mathcal{M}$, it suffices to design
vibration control inputs to ensure the uniform exponential stability of $x=0$
for the system (4).

##  III Vibrational Control: General Results

In this section, we derive some general results on vibrational stabilization
of $x=0$ of the system (4).

The term $f_{\rm inter}(x,y)$ in (4a) can be taken as a vanishing perturbation
dependent of $y$ to the controlled nominal system

$\displaystyle\dot{x}=f_{\rm
intra}(x)+{\small\begin{matrix}\frac{1}{\varepsilon}\end{matrix}}f_{\rm
ctr}(U({\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}),x).$ (6)

It can be decomposed as $f_{\rm inter}=[(f^{(1)}_{{\rm
inter}})^{\top},\dots,(f^{(r)}_{{\rm inter}})^{\top}]^{\top}$, where
$f^{(k)}_{{\rm inter}}=-(\hat{B}^{(k)}_{{\rm intra}})^{\top}B_{\rm
inter}W_{\rm inter}\sin(R_{2}x+R_{3}y)$ describes the inter-cluster dynamics
of the $k$th cluster.

Next, we show how the vibrational control term $f_{\rm ctr}(U,x)$ can
stabilize $x=0$ in the presence of the perturbation $f_{\rm inter}(x,y)$. To
this end, we consider a partially linearized model of the system (4a), which
is

$\displaystyle\scalebox{1.0}{\mbox{$\displaystyle\dot{x}=(J+{\small\begin{matrix}\frac{1}{\varepsilon}\end{matrix}}P({\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}))x+f_{\rm
inter}(x,y)$}},$ (7)

where

$\displaystyle{J=\frac{\partial f_{\rm intra}}{\partial x}(0)={\rm
blkdiag}(J^{(1)},\dots,J^{(r)})},\text{and}$ (8)
$\displaystyle{P(t)=\frac{\partial f_{\rm ctr}}{\partial x}(0)={\rm
blkdiag}(P^{(1)}(t),\dots,P^{(r)}(t))},$

with ${J^{(k)}=-(\hat{B}^{(k)}_{{\rm intra}})^{\top}B_{\rm
intra}^{(k)}W^{(k)}_{\rm intra}R_{1}}$ and

$\displaystyle\scalebox{0.8}{\mbox{$\displaystyle
P^{(k)}(t)=-(\hat{B}^{(k)}_{{\rm intra}})^{\top}\Big{(}[B^{(k)}_{{\rm
intra}}]^{+}U^{(k)}_{1}(t)-[-B^{(k)}_{{\rm
intra}}]^{+}U^{(k)}_{2}(t)\Big{)}R_{1}^{(k)}$}},$ (9)

for $k=1,\dots,r$. Note that only the first and third terms of Eq. (4a) are
linearized in Eq. (7). Observe that $P(t)$ is periodic and has the same period
$T$ as $U(t)$.

######  Lemma 1 (Connecting the stability of Systems (4) and (7)).

If the equilibrium $x=0$ is exponentially stable uniformly in $y$ for the
system (7), it is also exponentially stable uniformly in $y$ for the system
(4).

The proof can be found in Appendix -B. From this lemma, one can see that to
stabilize the cluster synchronization manifold $\mathcal{M}$, it suffices to
configure the vibrational control such that $P(t)$ stabilizes $x=0$ of the
system (7).

Let $s={t}/{\varepsilon}$. The system (7) can be rewritten as

$\displaystyle\frac{dx}{ds}=(\varepsilon J+P(s))x+\varepsilon f_{\rm
inter}(x,y).$ (10)

Next, we use averaging methods to analyze this system. However, the standard
first-order averaging is not applicable here. Recall that $P(s)$ has zero
mean. Then, applying the first-order averaging to (10) just eliminates the
$P(s)$ term and results in the uncontrolled system $\frac{dx}{ds}=\varepsilon
Jx+\varepsilon f_{\rm inter}(x,y)$.

To avoid that, we change the coordinates of (10) first before using averaging
method. To to that, we introduce an auxiliary system

$\displaystyle{\frac{d\hat{x}}{ds}=P(s)\hat{x}},$ (11)

and let $\Phi(s,s_{0})$ be its state transition matrix. Since $P$ is block-
diagonal, it holds that $\Phi={\rm blkdiag}(\Phi^{(1)},\dots,\Phi^{(r)}),$
where $\Phi^{(k)}$ is the transition matrix of the subsystem in the $k$th
cluster ${d\hat{x}_{k}}/{ds}=P^{(k)}(s)\hat{x}_{k}$.

Consider the change of coordinates $z(s)=\Phi^{-1}(s,s_{0})x(s)$. It follows
from the system (10) that

$\displaystyle{\frac{dz}{ds}=\varepsilon(\Phi^{-1}J\Phi z+\Phi^{-1}f_{\rm
inter}(\Phi z,y))}.$ (12)

Since $P(s)$ is $T$-periodic, $\Phi$ and $\Phi^{-1}$ are also $T$-periodic.
Then, we associate (12) with a partially averaged system

$\displaystyle{\frac{dz}{ds}=\varepsilon(\bar{J}z+\Phi^{-1}f_{\rm inter}(\Phi
z,y))},$ (13)

where

$\displaystyle{\bar{J}=\frac{1}{T}\int_{s_{0}}^{s_{0}+T}\Phi^{-1}(s,s_{0})J\Phi(s,s_{0})ds}.$
(14)

As both $J$ and $\Phi$ are block-diagonal, one can derive that $\bar{J}$ is
also block-diagonal satisfying $\bar{J}={\rm
blkdiag}(\bar{J}^{(1)},\dots,\bar{J}^{(r)})$ with

$\displaystyle\bar{J}^{(k)}=\frac{1}{T}\int_{s_{0}}^{s_{0}+T}\Big{(}\Phi^{(k)}(s,s_{0})\Big{)}^{-1}J^{(k)}\Phi^{(k)}(s,s_{0})ds.$

Recall that $\Phi$ and $\Phi^{-1}$ are periodic, $\|\Phi\|$ and
$\|\Phi^{-1}\|$ are both bounded. Then, it can be shown that there exist
$\bar{\gamma}_{k\ell}>0,k,\ell=1,\dots,r$, such that

$\displaystyle\|\Phi^{-1}f^{(k)}_{\rm inter}(\Phi
z,y)\|\leq\sum_{\ell=1}^{r}\bar{\gamma}_{k\ell}\|z_{\ell}\|$

for each $k$ (see Lemma 6 in Appendix -C for more details).

######  Theorem 1 (Sufficient condition for vibrational stabilization).

Assume that $\bar{J}={\rm blkdiag}(\bar{J}^{(1)},\dots,\bar{J}^{(r)})$ in Eq.
(14) is Hurwitz. Let $\bar{X}_{k}$ be the solution to the Lyapunov equation

$\displaystyle{(\bar{J}^{(k)})^{\top}\bar{X}_{k}+\bar{X}_{k}^{\top}\bar{J}^{(k)}=-I}.$
(15)

Define the matrix $S=[s_{k\ell}]_{r\times r}$ with

$\displaystyle
s_{k\ell}=\Big{\\{}\begin{matrix}[l]\lambda_{\max}^{-1}(\bar{X}_{k})-\bar{\gamma}_{kk},&\text{
if }k=\ell,\\\ s_{k\ell}=-\bar{\gamma}_{k\ell},&\text{ if
}k\neq\ell.\end{matrix}$ (16)

If $S$ is an $M$-matrix, then there exists $\varepsilon^{*}>0$ such that, for
any $\varepsilon<\varepsilon^{*}$:

(i) the equilibrium $x=0$ of the system (4) is exponentially stable uniformly
in $y$;

(ii) the cluster synchronization manifold $\mathcal{M}$ of the system (2) is
exponentially stable.

The proof can be found in Appendix -C. Theorem 1 provides a guideline to
design vibrational control. Any vibrational input that satisfies the following
three conditions stabilizes the cluster synchronization manifold
$\mathcal{M}$: i) $\bar{J}$ in (14) is Hurwitz, ii) $S$ defined in (16) is an
$M$-matrix, and iii) the frequency of the vibrations is sufficiently high,
i.e., $\varepsilon>0$ is sufficiently small.

We shall point out that the condition in Theorem 1 is still conservative. An
example will be presented in Sec. IV-B to show that a vibrational control can
effectively stabilize cluster synchronization without satisfying this
condition.

######  Remark 1 (Synchronization robustness).

For a stable linear system $\dot{x}=Ax$, previous studies propose that
$\lambda_{\max}(X)$, where $X$ is the solution to the Lyapunov equation
$A^{\top}X+XA=-2I$, can measure the robustness of this system [30, 31]. In our
case, from (7), the uncontrolled intra-dynamics around the manifold
$\mathcal{M}$ are described by

$\displaystyle\dot{x}_{k}=J^{(k)}x_{k}+f^{(k)}_{\rm inter}(x,y),$
$\displaystyle k=1,\dots,r,$ (17)

where $J^{(k)}$ is stable and $f^{(k)}_{\rm inter}(x,y)$ is taken as the
perturbation. Here, $x_{k}=0$ means synchronization of the oscillators in the
$k$th cluster. Let $X_{k}$ be the solution to

$\displaystyle\scalebox{0.85}{\mbox{$\displaystyle(J^{(k)})^{\top}X_{k}+X_{k}^{\top}J^{(k)}=-I$}}.$
(18)

Likewise, one can use $\lambda^{-1}_{\max}(X_{k})$ to measure the robustness
of synchronization in the $k$th cluster. If the intra-cluster synchronization
is sufficiently robust (i.e., large $\lambda^{-1}_{\max}(X_{k})$’s) to
dominate the perturbations resulted from inter-cluster connections, the
cluster synchronization is stable. A sufficient condition is constructed in
[15, Th. 3.2]. By contrast, if $\lambda_{\max}(X_{k})$’s are not large enough,
the cluster synchronization can lose its stability. Yet, the robustness of the
intra-cluster synchronization can be reshaped by introducing vibrations to the
local network connections. The new robustness is instead measured by
$\lambda^{-1}_{\max}(\bar{X}_{k})$, where $\bar{X}_{k}$ is the solution to Eq.
(15).  $\triangle$

The following example illustrates how vibrations can improve synchronization
robustness.

######  Example 1 (Improving robustness by vibrational control).

We restrict our attention to only one cluster in this example. Specifically,
consider that $J^{(k)}$ in (17) is

$\displaystyle\scalebox{0.85}{\mbox{$\displaystyle
J^{(k)}=\begin{bmatrix}-1&4\\\ 0&-2\end{bmatrix}$}}.$

Suppose that there is a vibrational control $U^{(k)}(t)$ to the $k$th cluster
that results in the state transition matrix

$\displaystyle{\Phi^{(k)}(t,t_{0})=\begin{bmatrix}1&0\\\
\cos(t_{0})-\cos(t)&1\end{bmatrix}}.$

It can be calculated that

$\displaystyle\scalebox{0.9}{\mbox{$\displaystyle\bar{J}^{(k)}=\frac{1}{2\pi}\int_{0}^{2\pi}\big{(}\Phi^{(k)}(t,t_{0})\big{)}^{-1}J^{(k)}\Phi^{(k)}(t,t_{0})dt=\begin{bmatrix}-1&4\\\
-2&-2\end{bmatrix}$}}.$

Solving the equations (15) and (18), one can obtain
$\lambda^{-1}_{\max}(X_{k})\approx 0.52$ and
$\lambda^{-1}_{\max}(\bar{X}_{k})=2$, which means that the synchronization in
this cluster becomes more robust due to the vibrational control.  $\triangle$

In fact, whether the synchronization robustness can be improved also depends
on the intra-cluster network. The following lemma provides a necessary
condition.

######  Lemma 2 (Necessary condition for robustness improvement).

Let $\lambda^{-1}_{\max}(\bar{X}_{k})$ and $\lambda^{-1}_{\max}(X_{k})$ be the
solutions to the equations (15) and (18), respectively. There exists a
vibrational control input such that
$\lambda^{-1}_{\max}(\bar{X}_{k})>\lambda^{-1}_{\max}(X_{k})$ _only if_
$\mathcal{G}_{k}$ is not a uniformly weighted complete graph333A uniformly
weighted complete graph is a complete graph in which every edge has the same
weight..

Figure 3: Uncontrollable cluster: a complete and uniformly weighted network.

When $\mathcal{G}_{k}$ is complete and uniformly weighted, it holds that
$J_{k}=-n_{k}I_{n_{k}-1}$, where $n_{k}$ is the number of nodes (see an
example in Fig. 3). As a consequence,
$\lambda^{-1}_{\max}(X_{k})=\lambda^{-1}_{\max}(\hat{X}_{k})$ for any
vibrational control since

$\displaystyle\scalebox{0.9}{\mbox{$\displaystyle\bar{J}^{(k)}=\int_{0}^{2\pi}\Phi^{-1}(t)J^{(k)}\Phi(t)=-n_{k}\int_{0}^{2\pi}\Phi^{-1}(t)\Phi(t)dt$}}=J^{(k)}.$

If all the clusters in a network are complete and uniformly weighted, one can
see that the synchronization robustness cannot be improved by any vibrational
control.

##  IV Configuring Vibrational Control

###  IV-A Configuring lower triangular transition matrices

Notice from Theorem 1 that, when designing a vibrational control, the state
transition matrix $\Phi$ of the periodic system (11) plays an important role.
However, deriving an explicit expression of $\Phi$ is often not
straightforward. In this subsection, we provide a special approach to
configure vibrational control so that the expression of $\Phi$ can be obtained
explicitly.

Inspired by [6], we consider vibration control inputs $U(t)$ such that
$P(t)={\rm blkdiag}(P^{(1)}(t),\dots,P^{(r)}(t))$ in Eq. (8) satisfies that,
for any $k$, $P^{(k)}$ has the following strictly lower triangular form (where
$\bar{n}_{k}=n_{k}-1$):

$\displaystyle\scalebox{0.75}{\mbox{$\displaystyle\begin{bmatrix}0&0&0&\cdots&0&0\\\
u_{21}^{(k)}\sin(t)&0&0&\cdots&0&0\\\
u_{31}^{(k)}\sin(t)&u_{32}^{(k)}\sin(t)&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
u_{\bar{n}_{k},1}^{(k)}\sin(t)&u_{\bar{n}_{k},2}^{(k)}\sin(t)&u_{\bar{n}_{k},3}^{(k)}\sin(t)&\cdots&u_{\bar{n}_{k},\bar{n}_{k}-1}^{(k)}\sin(t)&0\end{bmatrix}$}}.$
(19)

######  Lemma 3 (Lower triangular transition matrix).

Assume that the vibrational control is such that $P(t)$ in the system (8)
satisfies (IV-A). Then, the state transition matrix $\Phi(t,t_{0})={\rm
blkdiag}(\Phi_{1},\dots,\Phi_{r})$ satisfies

$\displaystyle\scalebox{0.9}{\mbox{$\displaystyle\Phi_{k}=\begin{bmatrix}1&0&0&\cdots&0&0\\\
\Phi_{21}^{(k)}&1&0&\cdots&0&0\\\
\Phi_{31}^{(k)}&\Phi_{32}^{(k)}&1&\cdots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
\Phi_{\bar{n}_{k},1}^{(k)}&\Phi_{\bar{n}_{k},2}^{(k)}&\Phi_{\bar{n}_{k},3}^{(k)}&\cdots&\Phi_{\bar{n}_{k},\bar{n}_{k}-1}^{(k)}&1\end{bmatrix}$}},$

where  $\displaystyle\Phi_{ij}^{(k)}(t)=-u_{ij}^{(k)}(\cos(t)-\cos(t_{0}))$ .

The next question is whether there exist vibrational control inputs such that
$P(t)$ has the form described by (IV-A).

######  Lemma 4 (Existence of vibrational control).

For each intra-cluster network $\mathcal{G}_{k},k=1,\dots,r$, there always
exist vibrational control inputs such that $P^{(k)}(t)$ is strictly lower
triangular and satisfies $P^{(k)}(t)\neq 0$.

The proof can be found in Appendix -D. From this lemma, one can always
carefully configure the vibrational control so that $P(t)$ is strictly lower
triangular. Subsequently, the explicit form of the transition matrix $\Phi$
can be derived, making the design of vibrations more analytically tractable.

It is worth noting that there certainly exist other configurations of
vibration control different from the one we present in this subsection, which
also allow for easy derivation of explicit transition matrices.

###  IV-B A numerical example

In this subsection, we use an example to show how the obtained results can be
used to design vibrational control.

We consider the network depicted in Fig. 4 (a), where the oscillators are
partitioned into $3$ clusters. The matrix $J={\rm
blkdiag}(J^{(1)},\dots,J^{(3)})$ as in the system (7) satisfies

$\displaystyle{J^{(1)}=\begin{bmatrix}-0.06&0.06\\\
0&-0.18\end{bmatrix},J^{(2)}=-3I_{3},\text{ and }J^{(3)}=-6I_{6}.}$

It can be seen in Fig. 4 (b) that the cluster synchronization manifold
$\mathcal{M}$ in this network is unstable.

Figure 4: Vibrational stabilization. (a) The network structure and some
connection weights, where $\alpha=0.03$ and $\beta=3$. All other unspecified
weights are $1$. (b) Phase differences without control, indicating that the
cluster synchronization is unstable. (c) Phase differences under vibrational
control to the cluster $\mathcal{C}_{1}$, showing that the cluster
synchronization has been stabilized by just local vibrations
($\varepsilon=0.02$). The natural frequencies in $\mathcal{C}_{1}$,
$\mathcal{C}_{2}$, and $\mathcal{C}_{3}$ are $\omega_{1}=1$, $\omega_{2}=10$,
and $\omega_{3}=6$, respectively.

To stabilize the cluster synchronization manifold, we only introduce
vibrations to the cluster $\mathcal{C}_{1}$. Specifically, we let

$\displaystyle\begin{matrix}[l]&u_{12}(t)=\scalebox{0.8}{\mbox{$\displaystyle\frac{u}{2\varepsilon}$}}\sin({\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}),&u_{13}(t)=-\scalebox{0.8}{\mbox{$\displaystyle\frac{u}{2\varepsilon}$}}\sin({\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}),\\\
&u_{32}(t)=-\scalebox{0.8}{\mbox{$\displaystyle\frac{u}{2\varepsilon}$}}\sin({\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}),&u_{31}(t)=\scalebox{0.8}{\mbox{$\displaystyle\frac{u}{2\varepsilon}$}}\sin({\small\begin{matrix}\frac{t}{\varepsilon}\end{matrix}}),\end{matrix}$

and $u_{ij}=0$ for any other $i$ and $j$. This configuration of vibrational
control ensures that the matrix $P^{(1)}(s)$ in Eq. (8) is strictly lower
triangular, i.e.,

$\displaystyle\scalebox{0.85}{\mbox{$\displaystyle
P^{(1)}(s)=\begin{bmatrix}0&0\\\ -u\sin(s)&0\end{bmatrix}$}}.$

Then, the transition matrix of the periodic system (11) can be calculated as

$\displaystyle\scalebox{0.85}{\mbox{$\displaystyle\Phi^{(1)}(s,s_{0})=\begin{bmatrix}1&0\\\
u(\cos(s_{0})-\cos(t))&1\end{bmatrix}$}}.$

Subsequently, the matrix $\bar{J}^{(1)}$ in Eq. (14) is

$\displaystyle\scalebox{0.85}{\mbox{$\displaystyle\bar{J}^{(1)}=\frac{1}{2\pi}\int_{0}^{2\pi}\big{(}\Phi^{(1)}(s,s_{0})\big{)}^{-1}J^{(1)}\Phi^{(1)}(s,s_{0})ds=\begin{bmatrix}-0.06&0.06\\\
\frac{-0.06u^{2}}{8}&-0.18\end{bmatrix}$}}.$

Let $u=2\sqrt{2}$. Solving the equations (15) and (18), respectively, one
obtains $\lambda^{-1}_{\max}(X_{1})\approx 0.109$ and
$\lambda^{-1}_{\max}(\bar{X}_{1})\approx 0.133$. This indicates that the
synchronization robustness in $\mathcal{C}_{1}$ is improved by the vibrational
control. From Fig. 4 (c), it can be observed that the local vibrations
introduced to $\mathcal{C}_{1}$, with frequency $50$ rad/s (i.e.,
$\varepsilon=0.02$), successfully stabilize the cluster synchronization in the
entire network.

It is worth mentioning that the condition in Theorem 1 is not satisfied. Yet
and remarkably, the slight improvement on synchronization robustness still
stabilizes the cluster synchronization. This suggests that vibrational control
can practically work in a much broader region than the theoretical one we
identified in Theorem 1.

##  V Conclusions

In this work, we hypothesize that the vibrational control is a mechanism of
DBS, a widely-used neurosurgical technique to treat some brain disorders,
based on the observation that they both employ high-frequency dithers and
regulate synchronization patterns. We study how the vibrational control can
stabilize cluster synchronization in Kuramoto-oscillator networks. Some
sufficient conditions on the vibrations introduced to the network connections
are obtained. We also provide an analytically tractable approach to design
vibrational inputs. Finally, we demonstrate our theoretical findings by a
numerical example. We believe that our work provides a new angle to understand
DBS and can inform the design of DBS algorithms.

###  -A Derivaton of Eq. (4)

To derive Eq. (4), we first present an instrumental lemma.

######  Lemma 5.

For the incidence matrices $B$ and $\hat{B}$, there exists
$R={\scriptsize\begin{bmatrix}R_{1}&0\\\ R_{2}&R_{3}\end{bmatrix}}$ such that
$B^{\top}=R\hat{B}^{\top}$, where

and ${\small\begin{matrix}R_{2}=B_{\rm inter}^{\top}(\hat{B}_{\rm
intra}^{\top}P_{\rm inter})^{\dagger}\end{matrix}}$,
${\small\begin{matrix}R_{3}=B_{\rm inter}^{\top}(\hat{B}_{\rm inter}P_{\rm
intra})^{\dagger}\end{matrix}}$ with ${\small\begin{matrix}P_{\rm
intra}=I_{n}-\hat{B}_{\rm intra}\hat{B}_{\rm intra}^{\dagger}\end{matrix}}$
and ${\small\begin{matrix}P_{\rm inter}=I_{n}-\hat{B}_{\rm inter}\hat{B}_{\rm
inter}^{\dagger}\end{matrix}}$.  $\triangle$

###### Proof.

The expressions of $R_{2}$ and $R_{3}$ follow from [18]. Recall that ${B_{\rm
intra}={\rm blkdiag}(B_{{\rm intra}}^{(1)},\dots,B_{{\rm intra}}^{(r))}}$and
${\hat{B}_{\rm intra}={\rm blkdiag}(\hat{B}^{(1)}_{{\rm
intra}},\dots,\hat{B}^{(r)}_{{\rm intra}})}$. Since ${(B^{(k)}_{{\rm
intra}})^{\top}}$ has full row rank, we have

$\displaystyle(B^{(k)}_{{\rm intra}})^{\top}=(B^{(k)}_{{\rm
intra}})^{\top}((\hat{B}^{(k)}_{{\rm
intra}})^{\top})^{\dagger}(\hat{B}^{(k)}_{{\rm intra}})^{\top},$

which completes the proof. ∎

One can rewrite Eq. (2) into a compact form

$\displaystyle\dot{\theta}=\omega-BW\sin(B^{\top}\theta)-B_{\rm
dir}^{+}U\sin(B_{\rm dir}^{\top}\theta).$

From Lemma 5, it holds that $B^{\top}\theta=R\hat{B}^{\top}\theta$ and $B_{\rm
dir}^{\top}\theta=R[\hat{B}_{\rm intra},-\hat{B}_{\rm intra}]^{\top}\theta$.
Subsequently, one can derive

$\displaystyle\hat{B}^{\top}\dot{\theta}=$
$\displaystyle\hat{B}^{\top}\omega-\hat{B}^{\top}BW\sin(R\hat{B}^{\top}\theta)$
$\displaystyle-B_{\rm dir}^{+}U\sin(R[\hat{B}_{\rm intra},-\hat{B}_{\rm
intra}]^{\top}\theta).$

Since $x=\hat{B}_{\rm intra}^{\top}\theta$ and $y=\hat{B}_{\rm
inter}^{\top}\theta$, one can obtain the compact dynamics (4) using the fact
that $[B_{{\rm dir}1}]^{+}=[B_{\rm intra}]^{+}$ and $[B_{{\rm
dir}1}]^{+}=[-B_{\rm intra}]^{+}$ and the relations between matrices in Eq.
(3) .

###  -B Proof of Lemma 1

Since $x=0$ is exponentially stable uniformly in $y$ for the system (7),
according to the converse Lyapunov theorem (see [29, Th. 4.4] and [32]) there
exists $\mathcal{D}=\\{x\in\mathbb{R}^{n-r}:\|x\|\leq\rho_{1}\\}$ and a
continuously differentiable function function
$V:[0,\infty]\times\mathcal{D}\times\mathbb{R}^{r}\to\mathbb{R}$ such that

$\displaystyle\frac{\partial V}{\partial t}+\frac{\partial V}{\partial
x}\big{(}(J+P(t))x+f_{\rm inter}(x,y)\big{)}\leq-c_{1}\|x\|^{2}$

and $\|\frac{\partial V}{\partial x}\|\leq c_{2}\|x\|$ for some constants
$c_{1},c_{2}>0$. Let $h(t,x)=f_{\rm intra}(x)+f_{\rm ctr}(U(t),x)$ and
$\Delta(t,x)=h(t,x)-(J+P(t))x.$ It can be checked that $\partial h/\partial x$
is bounded and Lipschitz on $\mathcal{D}$. Then, similar to the proof of [33,
Th. 4.13], one can show that $\|\Delta(t,x)\|\leq c_{3}\|x\|^{2}$ for some
$c_{3}>0$. The time derivative along the system (4) satisfies

$\displaystyle{\frac{\partial V}{\partial t}+\frac{\partial V}{\partial
x}\big{(}(J+P(t))x+\Delta(t,x)+f_{\rm inter}(x,y)\big{)}}$ $\displaystyle\leq-
c_{1}\|x\|^{2}+c_{2}c_{3}\|x\|^{3}$
$\displaystyle\leq-(c_{1}-c_{2}c_{3}\rho)\|x\|^{2},\forall\|x\|<\rho.$

Choosing $\rho=\min\\{\rho_{1},c_{1}/c_{2}c_{3}\\}$ completes the proof.

###  -C Proof of Theorem 1

One can observe that (i) implies (ii). Then, it suffices to prove the
exponential stability of $x=0$ for (10). To do that, we first present the
following lemma, whose proof follows similar lines as Lemma 3.1 of [15].

######  Lemma 6 (Growth bound of perturbations).

There exist some constants $\bar{\gamma}_{k\ell}>0$, $k,\ell=1,\dots,r$, such
that, for any $k$, it holds that

$\displaystyle\|\Phi^{-1}f^{(k)}_{\rm inter}(\Phi
z,y)\|\leq\sum\nolimits_{\ell=1}^{r}\bar{\gamma}_{k\ell}\|z_{\ell}\|.\hskip
20.0pt\hbox{}\hfill\triangle$

Let $V_{k}=z_{k}^{\top}X_{k}z_{k}$, and it holds that $\partial V_{k}/\partial
z_{k}\leq\lambda_{\max}(X_{k})$. Choose $V(z)=\sum_{k=1}^{r}d_{k}V_{k}$ as a
Lyapunov candidate. The time derivative of $V(z)$ satisfies

$\displaystyle\dot{V}(z)$

$\displaystyle=\scalebox{0.85}{\mbox{$\displaystyle\sum_{k=1}^{r}d_{k}[z_{k}^{\top}((\bar{J}^{(k)})^{\top}X_{k}+X_{k}\bar{J}^{(k)})z_{k}+\frac{\partial
V}{\partial z_{k}}\Phi^{-1}f_{\rm inter}(\Phi z,y)]$}}$
$\displaystyle\scalebox{0.85}{\mbox{$\displaystyle\leq\sum_{k=1}^{r}d_{k}[-\|z_{k}\|^{2}+\lambda_{\max}(X_{k})\sum_{k=1}^{r}\bar{\gamma}_{k\ell}\|z_{k}\|\|z_{\ell}\|]$}},$

where the second inequality has used Lemma 6.

Let ${D:={\rm diag}(d_{1},\dots,d_{r})}$ and ${\hat{S}=[\hat{s}_{ij}]_{r\times
r}}$ where

$\displaystyle\hat{s}_{k\ell}=\Big{\\{}\begin{matrix}[l]1-\lambda_{\max}(X_{k})\bar{\gamma}_{kk},&\text{
if }k=\ell,\\\
\hat{s}_{k\ell}=-\lambda_{\max}(X_{k})\bar{\gamma}_{k\ell},&\text{ if
}k\neq\ell.\end{matrix}$

Then, one can rewrite ${\dot{V}(z)\leq-\frac{1}{2}z^{\top}(DS+S^{\top}D)z}$.
By assumption, $S$ is an $M$-matrix, and so is $\hat{S}$ since
$\hat{S}=-\lambda_{\max}(X_{k})S$. It follows from [33, Th. 9.2] that the
system (13) is exponentially stable.

Following similar steps as in [33, Th. 10.4], one can prove that there exists
$\varepsilon^{*}>0$ such that for any $\varepsilon<\varepsilon^{*}$, $z=0$ is
exponentially stable uniformly in $y$ for the system (12). Since
$x(s)=\Phi(s,s_{0})z(s)$ and $\|\Phi\|$ is bounded, then $x=0$ is also
exponentially stable uniformly in $y$ for (10), which completes the proof.

###  -D Proof of Lemma 4

###### Proof.

We construct the proof by considering a specific configuration of vibrational
control.

Without loss of generality, we assume that each $B_{{\rm intra}}^{(k)}$ and
$\hat{B}_{{\rm intra}}^{(k)}$ are ordered in a way such that their first
$\bar{n}_{k}=n_{k}-1$ columns are the same, and the first two column are
$e_{2}-e_{1}$ and $e_{3}-e_{2}$, respectively, where $e_{i}$ is the $i$th
column of $I_{n_{k}}$. Then, there exists a matrix $\tilde{R}^{(k)}_{1}$ such
that $R_{1}^{(k)}={\begin{bmatrix}I_{\bar{n}_{k}}\\\
\tilde{R}^{(k)}_{1}\end{bmatrix}}.$

Let $U^{(k)}_{1}(t)=ue_{1}e_{1}^{\top}\sin(t)$ and
$U^{(k)}_{2}(t)=-ue_{1}e_{1}^{\top}\sin(t)$. It can be calculated that

$\displaystyle P^{(k)}(t)$ $\displaystyle=-(\hat{B}^{(k)}_{{\rm
intra}})^{\top}[B^{(k)}_{{\rm intra}}]^{+}U^{(k)}_{1}(t)-[-B^{(k)}_{{\rm
intra}}]^{+}U^{(k)}_{2}(t)R_{1}^{(k)}$
$\displaystyle\scalebox{1.0}{\mbox{$\displaystyle=-u(\hat{B}^{(k)}_{{\rm
intra}})^{\top}(e_{1}+e_{2})(e_{1}^{\prime})^{\top}\sin(t)$}},$

where the second equality follows from
$e_{1}^{\top}\scalebox{0.7}{\mbox{$\displaystyle\begin{bmatrix}I_{\bar{n}_{k}}\\\
\tilde{R}^{(k)}_{1}\end{bmatrix}$}}=(e_{1}^{\prime})^{\top}$ with
$e_{1}^{\prime}$ being the first column of $I_{\bar{n}_{k}}$. Then, it can be
seen that the first two rows of $P^{(k)}(t)$ are
$-u(e_{2}-e_{1})^{\top}(e_{1}+e_{2})(e_{1}^{\prime})^{\top}\sin(t)=0,$ and
$-u(e_{3}-e_{2})^{\top}(e_{1}+e_{2})(e_{1}^{\prime})^{\top}=(e^{\prime}_{1})^{\top}\sin(t),$
respectively. The rest rows of $P^{(k)}(t)$ can only have non-zero values at
the first column. Therefore, $P^{(k)}(t)$ have the strictly lower triangular
form and satisfies $P^{(k)}(t)\neq 0$, which completes the proof. ∎

## References

  * [1] G. Hahn, A. Ponce-Alvarez, G. Deco, A. Aertsen, and A. Kumar, “Portraits of communication in neuronal networks,” _Nature Reviews Neuroscience_ , vol. 20, no. 2, pp. 117–127, 2019.
  * [2] J. Fell and N. Axmacher, “The role of phase synchronization in memory processes,” _Nature Reviews Neuroscience_ , vol. 12, no. 2, pp. 105–118, 2011.
  * [3] C. Hammond, H. Bergman, and P. Brown, “Pathological synchronization in Parkinson’s disease: Networks, models and treatments,” _Trends in Neurosciences_ , vol. 30, no. 7, pp. 357–364, 2007.
  * [4] P. Jiruska, M. De Curtis, J. G. Jefferys, C. A. Schevon, S. J. Schiff, and K. Schindler, “Synchronization and desynchronization in epilepsy: controversies and hypotheses,” _The Journal of Physiology_ , vol. 591, no. 4, pp. 787–797, 2013.
  * [5] J. K. Krauss, N. Lipsman, T. Aziz, A. Boutet, P. Brown, J. W. Chang, B. Davidson, W. M. Grill, M. I. Hariz, A. Horn _et al._ , “Technology of deep brain stimulation: current status and future directions,” _Nature Reviews Neurology_ , vol. 17, no. 2, pp. 75–87, 2021.
  * [6] S. M. Meerkov, “Principle of vibrational control: Theory and applications,” _IEEE Transactions on Automatic Control_ , vol. 25, no. 4, 1980.
  * [7] R. E. Bellman, J. Bentsman, and S. M. Meerkov, “Vibrational control of nonlinear systems: Vibrational controllability and transient behavior,” _IEEE Transactions on Automatic Control_ , vol. 31, no. 8, pp. 717–724, 1986\.
  * [8] B. Shapiro and B. T. Zinn, “High-frequency nonlinear vibrational control,” _IEEE Transactions on Automatic Control_ , vol. 42, no. 1, pp. 83–90, 1997\.
  * [9] X. Cheng, Y. Tan, and I. Mareels, “On robustness analysis of linear vibrational control systems,” _Automatica_ , vol. 87, pp. 202–209, 2018\.
  * [10] M. Breakspear, S. Heitmann, and A. Daffertshofer, “Generative models of cortical oscillations: Neurobiological implications of the kuramoto model,” _Frontiers in Human Neuroscience_ , vol. 4, 2010.
  * [11] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy, “Cluster synchronization and isolated desynchronization in complex networks with symmetries,” _Nature Communications_ , vol. 5, no. 1, pp. 1–8, 2014\.
  * [12] Y. S. Cho, T. Nishikawa, and A. E. Motter, “Stable chimeras and independently synchronizable clusters,” _Physical Review Letters_ , vol. 119, no. 8, p. 084101, 2017.
  * [13] Y. Qin, M. Cao, B. D. O. Anderson, D. S. Bassett, and F. Pasqualetti, “Mediated remote synchronization: the number of mediators matters,” _IEEE Control Systems Letters_ , vol. 5, no. 3, pp. 767–772, 2020.
  * [14] M. T. Schaub, N. O’Clery, Y. N. Billeh, J.-C. Delvenne, R. Lambiotte, and M. Barahona, “Graph partitions and cluster synchronization in networks of oscillators,” _Chaos_ , vol. 26, no. 9, p. 094821, 2016.
  * [15] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti, “Stability conditions for cluster synchronization in networks of heterogeneous Kuramoto oscillators,” _IEEE Transactions on Control of Network Systems_ , vol. 7, no. 1, pp. 302 – 314, 2020.
  * [16] Y. Qin, Y. Kawano, O. Portoles, and M. Cao, “Partial phase cohesiveness in networks of networks of kuramoto oscillators,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 12, pp. 6100–6107, 2021.
  * [17] P. Feketa, A. Schaum, and T. Meurer, “Synchronization and multicluster capabilities of oscillatory networks with adaptive coupling,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 7, pp. 3084–3096, 2020.
  * [18] R. Kato and H. Ishii, “Averaging and cluster synchronization of Kuramoto oscillators,” in _European Control Conference_ , 2021, pp. 1497–1502.
  * [19] A. Salova and R. M. D’Souza, “Cluster synchronization on hypergraphs,” _arXiv preprint arXiv:2101.05464_ , 2021.
  * [20] ——, “Analyzing states beyond full synchronization on hypergraphs requires methods beyond projected networks,” _arXiv preprint arXiv:2107.13712_ , 2021\.
  * [21] W. Wu, W. Zhou, and T. Chen, “Cluster synchronization of linearly coupled complex networks under pinning control,” _IEEE Transactions on Circuits and Systems_ , vol. 56, no. 4, pp. 829–839, 2009.
  * [22] L. V. Gambuzza, M. Frasca, and V. Latora, “Distributed control of synchronization of a group of network nodes,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 1, pp. 365–372, 2018.
  * [23] D. Fiore, G. Russo, and M. di Bernardo, “Exploiting nodes symmetries to control synchronization and consensus patterns in multiagent systems,” _Control Systems Letters_ , vol. 1, no. 2, pp. 364–369, 2017.
  * [24] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti, “Functional control of oscillator networks,” _arXiv:2012.04217_ , 2021, submitted.
  * [25] G. Weerasinghe, B. Duchet, H. Cagnan, P. Brown, C. Bick, and R. Bogacz, “Predicting the effects of deep brain stimulation using a reduced coupled oscillator model,” _PLoS computational biology_ , vol. 15, no. 8, p. e1006575, 2019.
  * [26] M. H. Histed, V. Bonin, and R. C. Reid, “Direct activation of sparse, distributed populations of cortical neurons by electrical microstimulation,” _Neuron_ , vol. 63, no. 4, pp. 508–522, 2009.
  * [27] L. Tiberi, C. Favaretto, M. Innocenti, D. S. Bassett, and F. Pasqualetti, “Synchronization patterns in networks of Kuramoto oscillators: A geometric approach for analysis and control,” in _IEEE Conf. on Decision and Control_ , Melbourne, Australia, Dec. 2017, pp. 481–486.
  * [28] M. Wang, S. Jiang, Y. Yuan, L. Zhang, J. Ding, J. Wang, J. Zhang, K. Zhang, and J. Wang, “Alterations of functional and structural connectivity of freezing of gait in parkinson’s disease,” _Journal of neurology_ , vol. 263, no. 8, pp. 1583–1592, 2016.
  * [29] W. M. Haddad and V. Chellaboina, _Nonlinear Dynamical Systems and Control: A Lyapunov-Based Approach_. Princeton, NJ, USA: Princeton University Press, 2011.
  * [30] R. Patel and M. Toda, “Quantitative measures of robustness for multivariable systems,” in _Joint Automation Control Conference_ , no. 17, 1980, p. 35.
  * [31] R. Yedavalli, “Improved measures of stability robustness for linear state space models,” _IEEE Transactions on Automatic Control_ , vol. 30, no. 6, pp. 577–579, 1985.
  * [32] Y. Qin, Y. Kawano, B. D. Anderson, and M. Cao, “Partial exponential stability analysis of slow-fast systems via periodic averaging,” _IEEE Transactions on Automatic Control_ , 2021.
  * [33] H. K. Khalil, _Nonlinear Systems_. Prentice Hall, 2002.


#  Where to drive? Free Space Detection with one Fisheye Camera

Tobias Scheck  Adarsh Mallandur  Christian Wiede  Gangolf Hirtz Chemnitz
University of Technology

###### Abstract

The development in the field of autonomous driving goes hand in hand with ever
new developments in the field of image processing and machine learning
methods. In order to fully exploit the advantages of deep learning, it is
necessary to have sufficient labeled training data available. This is
especially not the case for omnidirectional fisheye cameras. As a solution, we
propose in this paper to use synthetic training data based on Unity3D. A five-
pass algorithm is used to create a virtual fisheye camera. This synthetic
training data is evaluated for the application of free space detection for
different deep learning network architectures. The results indicate that
synthetic fisheye images can be used in deep learning context.

###### keywords:

Free Space Detection, Fisheye Camera, Deep Learning, CNN, Synthetic Data
Creation.

SGD

    stochastic gradient descent
CNN

    convolutional neural network
IoU

    Intersection-over-Union

##  1 Introduction

In recent years, the idea of autonomous driving has gained tremendous
momentum. All automotive companies and suppliers have set up large research
departments to turn the dream of driver-less driving into reality. The
development is currently between levels 2 (Partial automation) and 3
(Conditional automation) of the level of automation. These includes systems
such as adaptive cruise control, parking assistance or driving on highway. In
order to reach levels 4 (High Automation) and 5 (Full automation), however,
further challenges have to be mastered. A fundamental challenge is the
detection of free spaces in which the car can generally move. These can be
streets, side walks or meadows but not parking vehicles, lanterns, trees or
people occupying the place. Therefore, it is essential for the navigation to
detect these free spaces. This is not sufficiently fulfilled by existing
systems. In addition, many external sensors are required.

In this paper, we propose a way to detect free spaces using a single
omnidirectional camera with a fisheye lens. Fisheye cameras can capture a wide
field of view. This is at the expense of strong tangential and radial
disturbances in the image. For a proper training with machine learning
algorithms a huge amount of training data is necessary. However, less or no
labeled training data exists for fisheye pictures. Especially the semantic
segmentation requires a pixel by pixel labeling. This cannot be done manually
without an extreme amount of resources. Instead, we suggest to use synthetic
data for the training. If these are as realistic as possible, these images can
be used for training and create a model working on real world data. In order
to create these images urban traffic scenes by Unity3D are used. This
environment can realize different backgrounds and dynamic driving scenarios.
Furthermore, we propose a five-pass algorithm in combination with synthetic
data to generate virtual fisheye images and a transfer learned convolutional
neural network (CNN) to generate a semantic segmentation.

This work is structured as follows: In Sect. 2, the related work of deep
learning based segmentation and synthetic data acquisition is outlined and the
research gap highlighted. The synthetic data generation by means of Unity3D is
described in Sect. 3. This is followed by the presentation of the deep
learning in Sect. 4. In Sect. 5, we present our experimental results, which is
accompanied by a discussion. Finally, we summarise our outcomes and outline
future work.

##  2 Related Work

Identifying drivable areas and detection of surrounding obstacles is a crucial
task for autonomous driving. A mixture of different sensors such as radar,
lidar, camera, high precision GPS and prior map information is used to
determine free space in current autonomous vehicles [1, 2]. Research in
autonomous driving started in Europe in 1986 with the PROMETHEUS project [3]
involving more than 13 vehicle manufacturers and several universities from
Europe. VaMP driverless car was a notable result of this project which covered
95 % of about a distance of 1,600 km fully automatically [4] and is considered
as one of the first autonomous cars [5]. Projects such as these were usually
computer vision systems for lateral and longitudinal vehicle guidance [6, 7],
lane departure warning [8] or collision avoidance [9]. The success of
PROMETHEUS and other similar projects [9, 10] influenced researchers to delve
into more complex urban environments from the much simpler highway scenarios.
The term urban environment means that the focus is towards the problem of
obstacle detection or free space detection rather than vehicle guidance.
Franke et al. [11] propose an ”Intelligent Stop & Go system”, using depth-
based obstacle detection and tracking from stereo images for urban
environment. In the 2000s, the problem of road detection was solved by using
low-level features such as color, shape [12] and texture [13] of the road.
Detection of road borders or road markings using laser [14], radar [15], Hough
transform [16] etc. is also used as a technique to detect drivable space.
However, such algorithms fail under severe lighting variations and depend
strongly on the structure of the roads. Hence, more and more algorithms relied
on geometric modeling using stereo images. Badino et al. introduce a method
for free space detection in complex scenarios using stochastic occupancy grids
[17]. These grids are built using data from stereo cameras and are integrated
over time with Kalman filters. Occupancy grids are in a way segmentation
between free and occupied regions. Apart from segmentation of pixels, other
representations such as Stixels [18], which groups objects with vertical
surfaces, have also been used for the task of urban traffic scene
understanding.

In the past decade, segmentation of road scenes using computer vision
techniques has been considered as a standard approach to detect free spaces
[19]. There are techniques ranging from Watershed Algorithms [20], Dense
Stereo Maps [21], Structure from Motion (SfM) [22] to global methods such as
Conditional Random Fields (CRF) [23] or Boosting Algorithms [24]. Inspired
from the success of CNNs in tasks such as classification, researchers have
found ways to apply CNNs in the task of road segmentation [19, 25, 26].
Alvarez et al. used CNN based algorithm to recover 3D scene layout of a road
image [19]. Hereby, the problems of overfitting and manual label generation
are solved by using noisy labels generated from a classifier, which is trained
on a general image dataset. The authors in [25] introduce convolutional patch
networks, which perform a pixel-wise labeling of input road scene images, by
first applying patch segmentation or classification of image patches. The
deconvolution networks used in [26] show a good performance on the KITTI
benchmark. However, this approach is computationally expensive and is not
suited for real-time road segmentation [27]. The most notable work in the
field of semantic segmentation are the fully convolutional networks (FCN)
[28], which are realized by replacing fully connected layers of a CNN by
convolutional layers and then up-sampling the resulting feature map into a
pixel-wise segmented image thus enabling end-to-end learning. The absence of
fully connected layers meant that the input can be of any size and that the
spatial information is preserved. Furthermore, the resulting architecture has
a smaller number of trainable parameters, thus achieving better performance
than patch networks.

Although state-of-the-art semantic segmentation networks achieve excellent
accuracy, maximum information about the surrounding area is necessary for safe
navigation. Fisheye cameras provide a wider field of view than narrow-angle
pin-hole cameras in complex urban traffic scenes and are becoming more popular
in vehicles since they are cheap and easy to handle [29, 30, 31]. However,
fisheye images are distorted due to strong perspective projections and are
consequently unwarped for practical usage. This unwarping process decreases
the image quality, especially at the image boundaries. Thus, CNNs trained on
pin-hole camera images do not perform well on fisheye images. Several research
works have been carried out to develop algorithms that directly apply on
distorted fisheye images [32, 33, 34, 35, 36]. The main problem of these
algorithms is the lack of labeled fisheye urban scene datasets that have good
quality and rich scenes. Therefore, Qian et al. [36] devise adversarial method
of training CNNs using images from pin-hole camera images dataset to detect
pedestrians in fisheye images. An adversarial network generates fisheye
pedestrian features, which are difficult for the detector to classify. However
adversarial networks are considered difficult to train because of various
reasons such as unstable parameters, sensitivity to hyperparameters, mode
collapses, vanishing gradients and overfitting [37]. In [33, 34, 35], the
authors synthetically generate fisheye images with labels from Cityscapes [38]
and SYNTHIA [39] datasets to train a CNN. The problem is that a transformation
from a perspective to a fisheye image, loses image details in the centre
during warping. Our work will solve this problem by devising a method to
generate realistic synthetic fisheye images along with the labels based on
5-pass rendering algorithm implemented in Unity3D. It solves the problem of
dataset generation as well as that of training a CNN directly on distorted
fisheye images by providing a way to generate fisheye datasets of different
urban traffic scenes.

##  3 Synthetic Data Creation

In this section, we describe the generation of fisheye images along with their
labels by a 5-pass rendering algorithm by using Unity3D. Currently, different
types of cameras such as perspective or omnidirectional cameras are used for
autonomous driving. Omnidirectional cameras provide a wider field of view
(FOV) than perspective cameras thus eliminating the need for more cameras or
mechanically rotatable cameras. An ideal omnidirectional camera can capture
light from all directions falling onto the focal point, thus covering a full
sphere ($360\degree$ FOV). There are various omnidirectional cameras ranging
from catadioptric cameras with curved mirrors to dioptric cameras with purely
dioptric fisheye lenses. The name fisheye lens was coined by R. W. Wood to
refer any lens capable of imaging the entire hemisphere in object space onto a
finite circle in the focal plane [40]. Fisheye lenses are designed to cover
the whole hemispherical field in front of the camera with a FOV of about
$180\degree$. Such fisheye lens cameras cannot be modeled by a linear
projection due to very high distortions in the images caused by refractions in
fisheye lenses. The distortions in fisheye images where straight lines appear
curved are known as radial symmetric distortions and are usually caused by the
shape of the lens [41]. Furthermore, these distortions are not uniformly
distributed over all spatial areas [32]. There are several models to describe
the behavior of fisheye lenses such as equidistant projection, stereographic
projection, equisolid angle projection and orthogonal projection. However,
real lenses do not exactly follow the designed projection model [42]. Using
such models, the authors in [33, 34, 35, 43] generated synthetic fisheye
images from existing datasets with standard perspective projection images
captured from a single camera. Instead of transforming existing perspective
projection images to fisheye projection images, our method directly
synthesizes fisheye projection images using virtual cameras within Unity3D
environment.

Figure 1:  Pipeline of synthetic data acquisition beginning with 4 cameras,
each camera projection represents a face of a cube. In the second step, each
projection result is used as a texture for a precomputed mesh to form a
fisheye image. The previous steps are performed for each type of synthetic
data (RGB image, segmentation mask and depth map) and exported as an image.

In order to generate a dataset for training, Microsoft Airsim simulator [44]
designed for Unity3D is used. Airsim provides a plugin with Python interface
and a vehicle asset to be used in any Unity3D environment. To help Artificial
Intelligence (AI) research, Unity3D offers free Windridge City environment
inside which Airsim asset vehicle can be imported. The vehicle can be
controlled either by keyboard or through the Airsim Python interface. The
generation of fisheye images is based on well-known technique of rendering to
cubic maps, where a $360\degree$ view of the entire scene is rendered onto a
cube. In our case, we render to four faces of a cube by using each a
perspective camera with a FOV of $90\degree$, as shown in Fig. 1. The image
captured by each camera is rendered into a texture, which is then used for the
corresponding precomputed mesh to form a fisheye image. The texture
coordinates of these four meshes are designed in such a way that they capture
the radial symmetric distortions and produce a fisheye projection. Further
information about the mesh generation can be obtained such as described in
[45]. This is a computationally inexpensive method compared to a single pass
rendering algorithm with a vertex shader that pre-distorts the geometry of the
world to produce a fisheye projection [45].

In the five-pass rendering algorithm, four orthogonal cameras each with
$90\degree$ FOV are created inside Unity3D and are placed above the virtual
car. The fifth camera is orthographic and is positioned to capture the fisheye
projection formed by the combination of the four meshes. During one frame
update procedure of Unity3D, for each orthogonal camera, three additional
orthogonal hidden cameras at the same position are created: one for the source
image rendering, one for the label image rendering and one for the depth image
rendering. All these three hidden cameras render to the same texture as their
parent cameras. This approach is based on [46] and we modified it to use it
with the five-pass rendering algorithm. In order to render the label image,
the layer id of each object is mapped to a color and is passed to the vertex
shader for rendering. The depth label is computed inside the vertex shader by
calculating the normalized distance value of each pixel from the near plane
and assigning it to that pixel intensity. Therefore, for one image-saving
cycle, three fisheye projections are obtained: source fisheye image, label
fisheye image and depth fisheye image. To preserve the floating-point values
in the depth image, the .exr file format is used.

Each object inside Unity3D is assigned with a layer. Thus, the total number of
classes in the training set is equal to the number of layers in the Unity3D
environment. Unity3D supports a maximum of 32 layers. These layer IDs along
with their names and mapped colors are exported as JSON file in order to
retain the mapping information. The files were saved at a rate of one frame
per second in order to reduce redundant images and the CPU load. In order to
train a classifier, 12028 images were generated by driving the virtual car
around the Windridge city. We call the generated dataset OmniCity and the
usage of a 80/20 training/validation split ratio results in 9623 training and
2405 validation images. Thus, the proposed method provides a flexible way to
generate a dataset in a urban driving environment.

##  4 Deep Learning

This section introduces the CNN free space detection architectures and the
hyperparameters used for trainings. The selected learning strategies and data
augmentation methods during the trainings are described as well.

The detection of free spaces in a scene requires a general image understanding
on pixel level. In our approach we are using CNNs designed for semantic
segmentation tasks to train and classify pixel-wise. As shown in Fig. 2, the
first model is based on the implementation of SegNet [47], which is an
encoder-decoder architecture. This architecture gradually reduces the spatial
dimension in the encoder part, while the decoder part behaves in the opposite
direction. In concrete terms, this means as the CNN progresses, object details
and spatial dimensions are restored.

Figure 2:  This overview visualizes the used encoder-decoder structure for
free space detection. While the encoder produces sparse feature maps using max
pooling, the decoder uses a transposed convolution for a learnable upsampling.
Visualization based on [48]

DeepLabV3 was chosen as additional CNN architecture [49] with ResNet101 v2
[50] as feature extractor. Pooling layers, used in encoder-decoder
architectures, are designed in such a way that they increase the field of view
and aggregate necessary information but at the same time discards the local
context. However, this information is fundamental to semantic segmentation and
should not be dispensed with. To address this problem, the implementation of
DeepLabV3 uses atrous convolutions. They are also known as dilated
convolutions and allow an enlargement of the field of view without reducing
the spatial dimension and this at the same time efficiently [51, 52].

Both architectures use RGB images as input with a resolution of $512\times
512$ pixel. By using ResNet101 v2 as feature extractor, weights, pre-trained
on ILSVRC2012-CLS image classification dataset [53], are used for DeepLabV3.
The encoder-decoder weights are initialized using the technique described in
Glorot et al. [54]. To train all architectures we use stochastic gradient
descent (SGD) with a momentum of 0.9. An initial learning rate of 0.0001 for
DeepLabV3 and 0.001 for the encoder-decoder architecture is selected. A
training process can benefit from a reduction in the learning rate during
training progress. To address this property, we use cosine decay [55] as our
learning rate strategy. We train each architecture for 50 epochs with a mini-
batch size of 4. In the training process, we shuffle the dataset after each
epoch to ensure that each training example is used only once during an epoch.
As objective function, cross-entropy loss is used for training which is summed
over all pixel values in a mini-batch. The last layer for each architecture
calculates a per-pixel propability using softmax. In this work, the CNNs
distinguish only between free space and background.

In order to avoid over-fitting and ensure better generalization, data
augmentation methods are applied during training [56, 57]. All methods depend
on a random factor to decide a method is applied or not. We have selected
horizontal flipping, brightness changing with a max_delta of 0.5 and random
Gaussian noise with a mean of 0.0 and a standard deviation of 8.0. For a
training with the generated synthetic data we also change the hue of a RGB
image and the saturation.

##  5 Results and Discussion

In this section, we outline our experiments and discuss the corresponding
result. In the first experiment, we use the training set of the Cityscapes
dataset [38], containing 2975 images, to train both CNN architectures on
perspective images of urban street scenes. Cityscapes contains segmentation
masks for 30 classes. However, for free space detection we map the classes:
road, sidewalk, parking, rail track and terrain to the class free space and
everything else to the class background. The second experiment aims to
simulate a fisheye distortion by mapping perspective images into fisheye
images. Accordingly all images of Cityscapes are distorted with a focal length
of 159 based on the approach described in Deng et al. [33], named Fisheye
Cityscapes. In the third experiment, we use our approach for synthetic fisheye
data generation to train all CNNs on OmniCity.

After training the architectures, we validate the performance using a sequence
of real fisheye images. For this purpose, a subset of the fisheye data set
DriveA sequence [58] was manually labeled, in a coarse manner, for our free
space scenario111https://gitlab.com/tschec/annotations. In order to evaluate
the 150 labeled images of this subset, the Intersection-over-Union (IoU) was
used. The IoU is described as the quotient of the area of overlap of the
ground truth and detected mask and area of union of both masks. During the
evaluation all pixels beyond the fisheye FOV are labeled as void and ignored.

Table 1:  IoU (%) for free space detection after training using Cityscapes, Fisheye Cityscapes and our proposed synthetic dataset OmniCity. | CNN
---|---
Dateset | Encoder-Decoder | Deeplab v3
Cityscapes | 93.1 | 84.6
Fisheye Cityscapes | 90.7 | 72.5
OmniCity | 86.7 | 78.3

As shown in Tab. 1, the results for the CNN architectures trained on
Cityscapes are those that achieve the highest IoU. To the contrary, the
results for Fisheye Cityscapes decrease by 14.3 % for Deeplab v3 and by 2.58 %
for the Encoder-Decoder structure. It seems that the distortions caused by a
fisheye lense do not affect the free space detection results. This fact is
supported by the results achieved by a training with OmniCity. The generated
synthetic data, following the equiangular projection, achieves comparable
results to Fisheye Cityscapes. The IoU decreases, compared to Cityscapes, by
7.45 % using Deeplab v3 and about 6.87 % with the Encoder-Decoder. However, it
should be noted that, during training with OmniCity, none of the CNN
architectures have ever seen real data.

An explanation for the poorer performance using fisheye distorted images can
be the fact that the class free space is mainly not affected by those. Unlike
objects, it is difficult to recognize a pattern in free space. Roads do not
follow clear lines or terrain like meadows and fields rarely have flat
surfaces, but these properties are most affected by distortion. It can be
assumed that for more complex classes such as cars, people or buildings, the
distortion has a greater impact than that studied in our scenario.

$Deeplab\ v3$ (a) $GT$ (b) $Cityscapes$ (c) $Fisheye\ Cityscapes$ (d)
$OmniCity$ $Encoder-Decoder$ $Deeplab\ v3$ $Encoder-Decoder$

Figure 3: Qualitative results of the evaluated CNN architectures, trained on
different datasets.

In Fig. 3, qualitative free space segmentation results for the experiments are
shown. In general, Deeplab v3 is found to be slightly worse than the Encoder-
Decoder architecture. However, this also reflects the results from Tab. 1. In
the first selected scenario (row 1 and 2), a recognizable gap of quality
between the results of all Encoder-Decoder architectures seems minimal. The
difference in segmentation quality only becomes apparent in the last scenario
(row 3 and 4). While the CNNs trained with Cityscapes and Fisheye Cityscapes
are able to segment roads, footpaths and even terrain between them, the model
trained with OmniCity is not detecting these accurately. This is particularly
evident in the grassy areas and footpaths. One possible cause may be a lack of
training examples with such constellations.

However, we could demonstrate that it is possible to learn deep learning
architectures based on synthetic fisheye data and are able to apply that model
to real world data. Nevertheless, further work as to be carried out in this
field.

##  6 Conclusion

In this work, we proposed an approach for free space detection using a single
fisheye camera to benefit from its wider FOV. Due the lack of omnidirectional
training data, we use a data augmentation method, introduced in [33], to
simulate a fisheye distortion on perspective images. Additionally, we
described a method to generate synthetic omnidirectional images, segmentation
masks and depth maps, using Unity3D. A new omnidirectional dataset, called
OmniCity, with urban synthetic city scenes is introduced. Furthermore, an
Encoder-Decoder and the Deeplab v3 CNN architecture in the context of free
space detection was evaluated, which is commonly used for semantic
segmentations tasks.

It can been shown that for free space detection on omnidirectional images, the
use of fisheye distorted images for training is not mandatory, but can be
advantageous in case no real world data is available. However, during our
evaluation we have shown that CNNs trained only on synthetic data can achieve
comparable results. In future research, the evaluation of more complex
semantic segmentation scenarios should be investigated to address the usage of
omnidirectional images. For this purpose, we will improve our OmniCity dataset
to include a higher variation in urban scenes, different light and weather
conditions, different vegetation and more classes like people and additional
vehicles.

## REFERENCES

  * [1] Montemerlo, M., Becker, J., Bhat, S., Dahlkamp, H., Dolgov, D., Ettinger, S., Haehnel, D., Hilden, T., Hoffmann, G., Huhnke, B., et al., “Junior: The stanford entry in the urban challenge,” Journal of field Robotics 25(9), 569–597 (2008).
  * [2] Urmson, C., Anhalt, J., Bagnell, D., Baker, C., Bittner, R., Clark, M., Dolan, J., Duggins, D., Galatali, T., Geyer, C., et al., “Autonomous driving in urban environments: Boss and the urban challenge,” Journal of Field Robotics 25(8), 425–466 (2008).
  * [3] Williams, M., “Prometheus-the european research programme for optimising the road transport system in europe,” in [IEE Colloquium on Driver Information ], 1–1, IET (1988).
  * [4] Dickmanns, E. D. et al., “Vehicles capable of dynamic vision,” in [IJCAI ], 97, 1577–1592 (1997).
  * [5] Dickmanns, E. D., [Dynamic vision for perception and control of motion ], Springer Science & Business Media (2007).
  * [6] Weber, J., Koller, D., Luong, Q.-T., and Malik, J., “New results in stereo-based automatic vehicle guidance,” in [Proceedings of the Intelligent Vehicles’ 95. Symposium ], 530–535, IEEE (1995).
  * [7] Pomerleau, D. and Jochem, T., “Rapidly adapting machine vision for automated vehicle steering,” IEEE expert 11(2), 19–27 (1996).
  * [8] Franke, U., Mehring, S., Suissa, A., and Hahn, S., “The daimler-benz steering assistant: a spin-off from autonomous driving,” in [Proceedings of the Intelligent Vehicles’ 94 Symposium ], 120–124, IEEE (1994).
  * [9] Thorpe, C., Hebert, M. H., Kanade, T., and Shafer, S. A., “Vision and navigation for the carnegie-mellon navlab,” IEEE Transactions on Pattern Analysis and Machine Intelligence 10(3), 362–373 (1988).
  * [10] Morimoto, H., Koizumi, M., Inoue, H., and Nitadori, K., “Ahs road-to-vehicle communication system,” in [Proceedings 199 IEEE/IEEJ/JSAI International Conference on Intelligent Transportation Systems (Cat. No. 99TH8383) ], 327–334, IEEE (1999).
  * [11] Franke, U., Gavrila, D., Görzig, S., Lindner, F., Paetzold, F., and Wöhler, C., “Autonomous driving goes downtown,” IEEE Intelligent systems (6), 40–48 (1998).
  * [12] Sotelo, M. A., Rodriguez, F. J., Magdalena, L., Bergasa, L. M., and Boquete, L., “A color vision-based lane tracking system for autonomous driving on unmarked roads,” Autonomous Robots 16(1), 95–116 (2004).
  * [13] Lombardi, P., Zanin, M., and Messelodi, S., “Switching models for vision-based on-board road detection,” in [IEEE Conference on Intelligent Transportation Systems ], 67–72 (2005).
  * [14] Sparbert, J., Dietmayer, K., and Streller, D., “Lane detection and street type classification using laser range images,” in [ITSC 2001. 2001 IEEE Intelligent Transportation Systems. Proceedings (Cat. No. 01TH8585) ], 454–459, IEEE (2001).
  * [15] Ma, B., Lakshmanan, S., and Hero, A. O., “Simultaneous detection of lane and pavement boundaries using model-based multisensor fusion,” IEEE Transactions on Intelligent Transportation Systems 1(3), 135–147 (2000).
  * [16] Yu, B. and Jain, A. K., “Lane boundary detection using a multiresolution hough transform,” in [Proceedings of International Conference on Image Processing ], 2, 748–751, IEEE (1997).
  * [17] Badino, H., Franke, U., and Mester, R., “Free space computation using stochastic occupancy grids and dynamic programming,” in [Workshop on Dynamical Vision, ICCV, Rio de Janeiro, Brazil ], 20, Citeseer (2007).
  * [18] Badino, H., Franke, U., and Pfeiffer, D., “The stixel world-a compact medium level representation of the 3d-world,” in [Joint Pattern Recognition Symposium ], 51–60, Springer (2009).
  * [19] Alvarez, J. M., Gevers, T., LeCun, Y., and Lopez, A. M., “Road scene segmentation from a single image,” in [European Conference on Computer Vision ], 376–389, Springer (2012).
  * [20] Beucher, S., Bilodeau, M., and Yu, X., “Road segmentation by watershed algorithms,” in [PROMETHEUS Workshop, Sophia Antipolis, France ], (1990).
  * [21] Ladickỳ, L., Sturgess, P., Russell, C., Sengupta, S., Bastanlar, Y., Clocksin, W., and Torr, P. H., “Joint optimization for object class segmentation and dense stereo reconstruction,” International Journal of Computer Vision 100(2), 122–133 (2012).
  * [22] Sturgess, P., Alahari, K., Ladicky, L., and Torr, P. H., “Combining appearance and structure from motion features for road scene understanding,” in [BMVC-British Machine Vision Conference ], BMVA (2009).
  * [23] Passani, M., Yebes, J. J., and Bergasa, L. M., “Crf-based semantic labeling in miniaturized road scenes,” in [17th International IEEE Conference on Intelligent Transportation Systems (ITSC) ], 1902–1903, IEEE (2014).
  * [24] Vitor, G. B., Victorino, A. C., and Ferreira, J. V., “A probabilistic distribution approach for the classification of urban roads in complex environments,” in [IEEE Workshop on International Conference on Robotics and Automation ], (2014).
  * [25] Brust, C.-A., Sickert, S., Simon, M., Rodner, E., and Denzler, J., “Convolutional patch networks with spatial prior for road detection and urban scene understanding,” arXiv preprint arXiv:1502.06344 (2015).
  * [26] Mohan, R., “Deep deconvolutional networks for scene parsing,” arXiv preprint arXiv:1411.4101 (2014).
  * [27] Oliveira, G. L., Burgard, W., and Brox, T., “Efficient deep models for monocular road segmentation,” in [2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) ], 4885–4891, IEEE (2016).
  * [28] Shelhamer, E., Long, J., and Darrell, T., “Fully convolutional networks for semantic segmentation,” IEEE Trans. Pattern Anal. Mach. Intell. 39, 640–651 (Apr. 2017).
  * [29] Wang, C., Zhang, H., Yang, M., Wang, X., Ye, L., and Guo, C., “Automatic parking based on a bird’s eye view vision system,” Advances in Mechanical Engineering 6, 847406 (2014).
  * [30] Liu, Y.-C., Lin, K.-Y., and Chen, Y.-S., “Bird’s-eye view vision system for vehicle surrounding monitoring,” in [International Workshop on Robot Vision ], 207–218, Springer (2008).
  * [31] Haltakov, V., Belzner, H., and Ilic, S., “Scene understanding from a moving camera for object detection and free space estimation,” in [2012 IEEE Intelligent Vehicles Symposium ], 105–110, IEEE (2012).
  * [32] Fremont, V., Bui, M., Boukerroui, D., and Letort, P., “Vision-based people detection system for heavy machine applications,” Sensors 16(1), 128 (2016).
  * [33] Deng, L., Yang, M., Qian, Y., Wang, C., and Wang, B., “Cnn based semantic segmentation for urban traffic scenes using fisheye camera,” in [2017 IEEE Intelligent Vehicles Symposium (IV) ], 231–236, IEEE (2017).
  * [34] Deng, L., Yang, M., Li, H., Li, T., Hu, B., and Wang, C., “Restricted deformable convolution based road scene semantic segmentation using surround view cameras,” arXiv preprint arXiv:1801.00708 (2018).
  * [35] Sáez, Á., Bergasa, L. M., López-Guillén, E., Romera, E., Tradacete, M., Gómez-Huélamo, C., and del Egido, J., “Real-time semantic segmentation for fisheye urban driving images based on erfnet,” Sensors 19(3), 503 (2019).
  * [36] Qian, Y., Yang, M., Wang, C., and Wang, B., “Pedestrian feature generation in fish-eye images via adversary,” in [2018 IEEE International Conference on Robotics and Automation (ICRA) ], 2007–2012, IEEE (2018).
  * [37] Arjovsky, M., Chintala, S., and Bottou, L., “Wasserstein generative adversarial networks,” in [International Conference on Machine Learning ], 214–223 (2017).
  * [38] Cordts, M., Omran, M., Ramos, S., Scharwächter, T., Enzweiler, M., Benenson, R., Franke, U., Roth, S., and Schiele, B., “The cityscapes dataset,” in [CVPR Workshop on the Future of Datasets in Vision ], 2 (2015).
  * [39] Ros, G., Sellart, L., Materzynska, J., Vazquez, D., and Lopez, A. M., “The synthia dataset: A large collection of synthetic images for semantic segmentation of urban scenes,” in [Proceedings of the IEEE conference on computer vision and pattern recognition ], 3234–3243 (2016).
  * [40] Kingslake, R., [A history of the photographic lens ], Elsevier (1989).
  * [41] Tardif, J.-P., Sturm, P., and Roy, S., “Self-calibration of a general radially symmetric distortion model,” in [European conference on computer vision ], 186–199, Springer (2006).
  * [42] Kannala, J. and Brandt, S. S., “A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses,” IEEE transactions on pattern analysis and machine intelligence 28(8), 1335–1340 (2006).
  * [43] Sáez, A., Bergasa, L. M., Romeral, E., López, E., Barea, R., and Sanz, R., “Cnn-based fisheye image real-time semantic segmentation,” in [2018 IEEE Intelligent Vehicles Symposium (IV) ], 1039–1044, IEEE (2018).
  * [44] Shah, S., Dey, D., Lovett, C., and Kapoor, A., “Airsim: High-fidelity visual and physical simulation for autonomous vehicles,” in [Field and service robotics ], 621–635, Springer (2018).
  * [45] Bourke, P., “idome: Immersive gaming with the unity3d game engine,” CGAT09 Computer Games, Multimedia and Allied Technology 9, 265–272 (2009).
  * [46] “Unity-technologies / ml-imagesynthesis.”  (Accessed on 06/22/2019).
  * [47] Badrinarayanan, V., Kendall, A., and Cipolla, R., “Segnet: A deep convolutional encoder-decoder architecture for image segmentation,” IEEE transactions on pattern analysis and machine intelligence 39(12), 2481–2495 (2017).
  * [48] Iqbal, H., “Harisiqbal88/plotneuralnet v1.0.0,” (Dec. 2018).
  * [49] Chen, L., Papandreou, G., Schroff, F., and Adam, H., “Rethinking atrous convolution for semantic image segmentation,” CoRR abs/1706.05587 (2017).
  * [50] He, K., Zhang, X., Ren, S., and Sun, J., “Identity mappings in deep residual networks,” in [European conference on computer vision ], 630–645, Springer (2016).
  * [51] Yu, F. and Koltun, V., “Multi-scale context aggregation by dilated convolutions,” arXiv preprint arXiv:1511.07122 (2015).
  * [52] Chen, L.-C., Papandreou, G., Kokkinos, I., Murphy, K., and Yuille, A. L., “Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs,” IEEE transactions on pattern analysis and machine intelligence 40(4), 834–848 (2017).
  * [53] Russakovsky, O., Deng, J., Su, H., Krause, J., Satheesh, S., Ma, S., Huang, Z., Karpathy, A., Khosla, A., Bernstein, M., Berg, A. C., and Fei-Fei, L., “ImageNet Large Scale Visual Recognition Challenge,” International Journal of Computer Vision (IJCV) 115(3), 211–252 (2015).
  * [54] Glorot, X. and Bengio, Y., “Understanding the difficulty of training deep feedforward neural networks,” in [Proceedings of the thirteenth international conference on artificial intelligence and statistics ], 249–256 (2010).
  * [55] Loshchilov, I. and Hutter, F., “Sgdr: Stochastic gradient descent with warm restarts,” arXiv preprint arXiv:1608.03983 (2016).
  * [56] Takahashi, R., Matsubara, T., and Uehara, K., “Data augmentation using random image cropping and patching for deep cnns,” arXiv preprint arXiv:1811.09030 (2018).
  * [57] Krizhevsky, A., Sutskever, I., and Hinton, G. E., “Imagenet classification with deep convolutional neural networks,” in [Advances in neural information processing systems ], 1097–1105 (2012).
  * [58] Eichenseer, A. and Kaup, A., “A data set providing synthetic and real-world fisheye video sequences,” in [2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) ], 1541–1545, IEEE (2016).

  *[CNN]: convolutional neural network
  *[CNNs]: convolutional neural network
  *[SGD]: stochastic gradient descent
  *[IoU]: Intersection-over-Union