single = true UAV short=UAV, long = unmanned aerial vehicle, short-indefinite
= a, long-indefinite = an MAV short=MAV, long = micro aerial vehicle, OMAV
short=OMAV, long = omnidirectional micro aerial vehicle, short-indefinite =
an, long-indefinite = an IMU short=IMU, long = inertial measurement unit,
short-indefinite = an, long-indefinite = an RC short=RC, long = remote
control, short-indefinite = an, long-indefinite = a FPV short=FPV, long =
first person view, short-indefinite = an, long-indefinite = a MOCAP
short=MOCAP, long = motion capture system, short-indefinite = a, long-
indefinite = a

# Towards 6DoF Bilateral Teleoperation of an Omnidirectional Aerial Vehicle
for Aerial Physical Interaction

Mike Allenspach1, Nicholas Lawrance1, Marco Tognon1, and Roland Siegwart1 1
Authors are with the Autonomous Systems Lab, ETH Zürich, Zürich 8092,
Switzerland<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>research was in part supported by the
National Center of Competence in Research (NCCR) on Digital Fabrication, in
part by the NCCR Robotics, and in part by Armasuisse Science and Technology.

###### Abstract

Bilateral teleoperation offers an intriguing solution towards shared autonomy
with aerial vehicles in contact-based inspection and manipulation tasks.
Omnidirectional aerial robots allow for full pose operations, making them
particularly attractive in such tasks. Naturally, the question arises whether
standard bilateral teleoperation methodologies are suitable for use with these
vehicles. In this work, a fully decoupled 6DoF bilateral teleoperation
framework for aerial physical interaction is designed and tested for the first
time. The method is based on the well established rate control, recentering
and interaction force feedback policy. However, practical experiments evince
the difficulty of performing decoupled motions in a single axis only. As such,
this work shows that the trivial extension of standard methods is insufficient
for omnidirectional teleoperation, due to the operator’s physical inability to
properly decouple all input DoFs. This suggests that further studies on
enhanced haptic feedback are necessary.

##  I INTRODUCTION

The physical interaction between flying robots and the environment has gained
increasing interest in the robotics community in recent years. Aerial robots
with manipulation capabilities have already been successfully deployed in a
variety of interaction tasks [1, 2, 3]. The low-cost, high maneuverability and
nearly unlimited workspace of these aerial manipulators allow deployment in
hard-to-reach or remote places, as well as when contact-based inspection and
maintenance are too dangerous for human operators [4, 5]. In this regard, OMAV
offer a particularly compelling solution. Their capability to generate thrust
in any direction allows hovering at arbitrary orientations, as well as
independently controlling position and attitude. Thus, these platforms are
capable of precise motion and interaction force control, while simultaneously
rejecting disturbances [6, 7, 8].
Despite recent advances in autonomous control of these vehicles (e.g. [9, 10,
11]), existing regulations and safety requirements often still require a human
operator in the loop. Real-time inclusion of the human expert’s knowledge is
especially important when complex tasks must be performed in uncertain or
a-priori unknown environments, considering the limited decisional autonomy of
modern robots. Transferring the operator supervision and decision making
skills to the remote site requires careful design of teleoperation systems. On
the one hand, taking full control of the robot and exploiting its capabilities
is only possible if every DoF is individually controllable by the operator. In
view of the recent trend towards omnidirectional aerial manipulators,
teleoperation frameworks must therefore simultaneously provide fully decoupled
commands for all three translational and rotational axes. On the other hand,
meaningful system information must be reflected back to the operator through
haptic and/or visual feedback. This concept of bilateral teleoperation
improves the user situational awareness, which in turn supports their decision
making process [12].

State-of-the-art bilateral teleoperation approaches are almost exclusively
focused on underactuated platforms for which the operator can control the
position and yaw angle only. Naturally, the question arises if an extension to
omnidirectional vehicles is straightforward or if additional considerations
and/or problems must be addressed. Specifically, the goal of this work is to
evaluate whether it is possible to teleoperate an OMAV using standard
methodologies, in both contact-less and contact-based conditions.

###  I-A Related Work

Early works on teleoperation methods for aerial robots primarily considered
contact-free flight rather than interaction [13, 14], mostly focusing on
direct control of the vehicles and the obstacle avoidance problem. An
alternative teleoperation strategy was suggested in [15], where the operator
indirectly steers the MAV by modifying the parameters of a dynamic path.
Hereby, user feedback includes information about tracking performance or the
presence of obstacles.

Only recently, bilateral teleoperation of MAV has been extended to aerial
physical interaction as well. One of the earliest work to apply bilateral
teleoperation to aerial physical interaction has been presented in [16] but
was restricted to simulation. The authors proposed the use of a haptic device
with three actuated translational DoFs to command both the motion of the
vehicle and the interaction force when in contact, mapping the position of the
input device into the desired acceleration of the MAV. Simultaneously, the
device renders a feedback force aimed to recenter the input device, as well as
to provide an indication of the measured interaction force. A similar
reference and feedback generation scheme is used in [17], although
environmental forces are estimated using a risk field interaction model.
Relying on standard rate control and interaction force feedback, the framework
in [18, 19] makes use of passivity theory to ensure stability even under
communication delays.

The methods presented thus far focus exclusively on underactuated platforms
and are naturally limited to the position control of the vehicle. In fact, [7]
is the only work where the topic of omnidirectional bilateral teleoperation
for aerial robots is addressed. However, since it is not the main contribution
of the article, specifics regarding the implementation of reference and
feedback generation schemes, as well as detailed discussions are missing.
Furthermore, although the framework seems to support decoupled 6DoF bilateral
teleoperation in theory, the experimental verification is limited to
translational motion only.

Even though the extension of well established methods may appear trivial, a
detailed evaluation is still missing. It is important to note that standard
solutions developed for ground manipulators (e.g. Leader-Follower-
Configuration) are unsuitable for aerial robotics, since the limited input
workspace must be mapped to a virtually unlimited robot workspace. Additional
issues arise due to the complexity of the $\mathsf{SE}(3)$ control space,
containing an open vector space for translation and a closed isometry for
orientation. In summary, effective bilateral teleoperation for OMAV remains an
unsolved problem, let alone its application in physical aerial interaction
tasks.

###  I-B Contributions

In summary, the main contributions of this work are:

  * •

Design of a fully decoupled 6DoF teleoperation framework by extending the
established rate control, recentering and interaction force feedback policy to
$\mathsf{SE}(3)$.

  * •

Evaluation of the proposed policy in real-world flight experiments, including
free-flight omnidirectional reference generation, as well as push-and-slide
operation during physical interaction.

  * •

Discussion about limitations of standard policies, namely the operator’s
physical inability to properly decouple all input DoFs.

As such, this study serves as a first step towards remote controlled
omnidirectional aerial physical interaction, by identifying working features
and potential issues when adapting standard methods used for underactuated
MAV.

(a) Local environment: human and haptic device.

(b) Remote work environment: aerial robot and task objects.

Figure 1: Representation the 6DoF bilateral teleoperation setup. Local
Environment(1), (2)Reference Generation(8), (9)Feedback Generation(10)Remote
Environment(3)$\begin{bmatrix}\bm{p}_{H}\\\
\bm{q}^{M}_{H}\end{bmatrix}$$\begin{bmatrix}\bm{p}_{S,ref}\\\
\bm{v}_{S,ref}\\\ {\bm{R}^{W}_{S,ref}}\\\
\bm{\omega}_{{S,ref}}\end{bmatrix}$$\bm{\tau}_{fb,total}$$\bm{\tau}_{fb,ext}$
Figure 2: Detailed interactions between the different components of the
proposed teleoperation framework.

##  II Modeling

The system considered in a bilateral teleoperation framework consists of a
human operator, a haptic device, and the robot (in our case OMAV). Human and
haptic device constitute the local environment and are connected through a
virtual communication link to the robot in the remote work environment (see
Fig. 1). The human operator is in constant contact with the handle of the
haptic device which allows them to perform small-scale translational and
rotational perturbations. Hereby, the mechanical construction and actuation of
the input device must support decoupled 6DoF motion and force/torque feedback
rendering.

To describe the local configuration, we define the inertial frame
$\mathcal{F}_{M}=\\{O_{M},\bm{x}_{M},\bm{y}_{M},\bm{z}_{M}\\}$ with origin
$O_{M}$ and unit axes $\\{\bm{x}_{M},\bm{y}_{M},\bm{z}_{M}\\}$ corresponding
to the idle pose of the haptic device’s handle. Its current position and
orientation are captured by the frame
$\mathcal{F}_{H}=\\{O_{H},\bm{x}_{H},\bm{y}_{H},\bm{z}_{H}\\}$ with origin
$O_{H}$ rigidly attached to the handle. In particular, it will become clear
that when $\mathcal{F}_{H}$ and $\mathcal{F}_{M}$ coincide, the desired rate
commanded to the aerial robot is zero. In a similar fashion, two additional
frames are defined for the remote environment. The inertial world frame,
$\mathcal{F}_{W}=\\{O_{W},\bm{x}_{W},\bm{y}_{W},\bm{z}_{W}\\}$ is located at
an arbitrary origin point $O_{W}$, such that $\bm{z}_{W}$ is opposite to
gravity. Finally, the state of the aerial vehicle is described by the body
frame $\mathcal{F}_{S}=\\{O_{S},\bm{x}_{S},\bm{y}_{S},\bm{z}_{S}\\}$ whose
origin $O_{S}$ coincides with the OMAV’s center-of-mass and $\bm{x}_{S}$
points along the end-effector used during physical contact. A right-hand
superscript, e.g. $\bm{\star}^{S}$ , is used if a vector is not represented in
its original frame.

###  II-A Human Operator

Let $\bm{p}_{H}\in\mathbb{R}^{3}$ and $\bm{v}_{H}\in\mathbb{R}^{3}$, both
expressed in $\mathcal{F}_{M}$, denote the position and velocity of $O_{H}$
with respect to $O_{M}$, i.e. the haptic device end-effector position and
velocity with respect to the idle configuration. Similarly, we define the
attitude and angular rate of $\mathcal{F}_{H}$ with respect to
$\mathcal{F}_{M}$ as ${\bm{R}^{M}_{H}}\in\mathsf{SO}(3)$ and
$\bm{\omega}_{H}\in\mathbb{R}^{3}$, the latter expressed in $\mathcal{F}_{H}$.
Since the human is in constant contact with the end-effector, its pose and
twist values are similar to the ones of the end-effector. The dynamic relation
between these values is modeled in $\mathcal{F}_{H}$ as:

$\displaystyle\bm{M}_{H}\begin{bmatrix}{\dot{\bm{v}}_{H}}^{H}\\\
{\dot{\bm{\omega}}_{H}}\end{bmatrix}+\bm{D}_{H}\begin{bmatrix}\bm{v}_{H}^{H}\\\
\bm{\omega}_{H}\end{bmatrix}=-\bm{\tau}_{H}+\bm{\tau}_{H,act},$ (1)

where $\bm{\tau}_{H,act}\in\mathbb{R}^{6}$ are wrenches (stacked forces and
torques) from the muscles and $\bm{\tau}_{H}\in\mathbb{R}^{6}$ are the
interaction wrenches with the haptic device. The human inherent inertia and
damping are denoted as $\bm{M}_{H}\in\mathbb{R}^{6\times 6}$ and
$\bm{D}_{H}\in\mathbb{R}^{6\times 6}$.

Following standard practice in both RC and manned rotary-wing vehicle
piloting, it is assumed that the visual frame of the human is identical to
$\mathcal{F}_{S}$. In fact, the robot is commonly equipped with an onboard
camera that streams images back to the operator during remote operations.
However, the human’s point of view can be changed as desired by modifying the
coordinate frame conventions used in Section III accordingly.

###  II-B Haptic Device

Since the inertia and damping of the human dynamics in (1) are generally
unknown, the haptic device manipulator is interacting with an unknown
environment with low impedance. As explained in detail in [20], this is
generally undesired for torque-controlled systems, since the contact
constraints can not be accurately described. To still ensure compliant
interaction and haptic transparency, an admittance filter is introduced and
combined with a low-level joint position controller. Assuming perfect
tracking, the closed-loop robot arm dynamics can then be approximated as:

$\displaystyle{\bm{M}_{adm}}\begin{bmatrix}{\dot{\bm{v}}_{H}}^{H}\\\
{\dot{\bm{\omega}}_{H}}\end{bmatrix}+{\bm{D}_{adm}}\begin{bmatrix}\bm{v}_{H}^{H}\\\
\bm{\omega}_{H}\end{bmatrix}=\bm{\tau}_{H}+\bm{\tau}_{fb,total},$ (2)

with inertia
${\bm{M}_{adm}}=diag({\bm{M}_{adm,t}},{\bm{M}_{adm,r}})\in\mathbb{R}^{6\times
6}$ and damping
${\bm{D}_{adm}}=diag({\bm{D}_{adm,t}},{\bm{D}_{adm,r}})\in\mathbb{R}^{6\times
6}$ tuning matrices and where $\bm{\tau}_{fb,total}\in\mathbb{R}^{6}$ are the
desired feedback wrenches to be applied to the user (see Section III-B).

###  II-C OMAV

We denote the position and velocity of the OMAV’s center-of-mass $O_{S}$ with
respect to $\mathcal{F}_{W}$ with $\bm{p}_{S}\in\mathbb{R}^{3}$ and
$\bm{v}_{S}\in\mathbb{R}^{3}$. The orientation and angular rate of
$\mathcal{F}_{S}$ with respect to $\mathcal{F}_{W}$ is given as
${\bm{R}^{W}_{S}}\in\mathsf{SO}(3)$ and $\bm{\omega}_{S}\in\mathbb{R}^{3}$,
the latter expressed in $\mathcal{F}_{S}$. To allow for compliant interaction,
we assume that the robot is controlled by an impedance controller similar to
the one presented in [11]. Thus, the rendered closed-loop dynamics in
$\mathcal{F}_{S}$ are

$\displaystyle\bm{M}_{v}\begin{bmatrix}\dot{\bm{v}}_{S}^{S}\\\
\dot{{\bm{\omega}}}_{S}\end{bmatrix}+\bm{D}_{v}\begin{bmatrix}\bm{e}_{v}\\\
\bm{e}_{\omega}\end{bmatrix}+\bm{K}_{v}\begin{bmatrix}\bm{e}_{p}\\\
\bm{e}_{R}\end{bmatrix}=\hat{\bm{\tau}}_{ext}.$ (3)

The virtual inertia $\bm{M}_{v}\in\mathbb{R}^{6\times 6}$, damping
$\bm{D}_{v}\in\mathbb{R}^{6\times 6}$ and stiffness
$\bm{K}_{v}\in\mathbb{R}^{6\times 6}$ are tuning parameters of the on-board
controller and $\hat{\bm{\tau}}_{ext}\in\mathbb{R}^{6}$ describes external
disturbances acting on the platform. In the context of this work, we assume
that such external disturbances originate solely from physical interaction of
the robot with the environment.

Given a desired position $\bm{p}_{S,ref}\in\mathbb{R}^{3}$ and velocity
$\bm{v}_{S,ref}\in\mathbb{R}^{3}$ in $\mathcal{F}_{W}$, attitude
${\bm{R}^{W}_{S,ref}}\in\mathsf{SO}(3)$ and angular rate
$\bm{\omega}_{{S,ref}}\in\mathbb{R}^{3}$ in $\mathcal{F}_{S,ref}$ of the OMAV,
the tracking errors are defined in $\mathcal{F}_{S}$ as

$\displaystyle\bm{e}_{p}$
$\displaystyle={\bm{R}^{W}_{S}}^{\top}\left({\bm{p}_{S}}-{\bm{p}_{S,ref}}\right)$
(4) $\displaystyle\bm{e}_{R}$
$\displaystyle=\frac{1}{2}({\bm{R}^{W}_{S,ref}}^{\top}{\bm{R}^{W}_{S}}-{\bm{R}^{W}_{S}}^{\top}{\bm{R}^{W}_{S,ref}})^{\vee}$
(5) $\displaystyle\bm{e}_{v}$
$\displaystyle={\bm{R}^{W}_{S}}^{\top}\left(\bm{v}_{S}-\bm{v}_{S,ref}\right)$
(6) $\displaystyle\bm{e}_{\omega}$
$\displaystyle=\bm{\omega}_{S}-{\bm{R}^{W}_{S}}^{\top}{\bm{R}^{W}_{S,ref}}\bm{\omega}_{{S,ref}},$
(7)

where the vee-map $(\cdot)^{\vee}:\mathfrak{so}{(3)}\rightarrow\mathbb{R}^{3}$
is the inverse of the skew-symmetric operator
$[\cdot]_{\times}:\mathbb{R}^{3}\rightarrow\mathfrak{so}{(3)}$.

##  III TELEOPERATION

An overview of the omnidirectional bilateral teleoperation framework proposed
in this work is presented in Fig. 2. The employed reference and feedback
generation policies are explained in detail in the following sections.

###  III-A Rate Control Reference Generation

Rate control is a well established method to teleoperate aerial vehicles,
since it provides an intuitive mapping between the restricted input workspace
and the potentially infinite robot workspace [21]. However, in the state-of-
the-art literature, it is often limited to translation only due to the
underactuated nature of standard fixed-rotor MAV. This restriction does not
apply to OMAV, which is why the concept is extended in this work to include
rotational rate control as well. Essentially, any translational or rotational
deviation between $\mathcal{F}_{H}$ and $\mathcal{F}_{M}$ is translated into a
corresponding translational velocity or angular rate reference for the robot.
Based on this, the translational references are computed as

$\displaystyle\bm{v}_{S,ref}^{S}$ $\displaystyle=v_{max}\bm{p}_{H}$ (8a)
$\displaystyle\bm{p}_{S,ref}^{S}$
$\displaystyle=\int_{0}^{t}\bm{v}_{S,ref}^{S}(s)ds,$ (8b)

and similarly for rotation

$\displaystyle\bm{\omega}_{{S,ref}}^{S}$
$\displaystyle=\frac{\omega_{max}}{2}\left({\bm{R}^{M}_{H}}-{\bm{R}^{M}_{H}}^{\top}\right)^{\vee}$
(9a) $\displaystyle{\bm{R}^{W}_{S,ref}}$
$\displaystyle=\int_{0}^{t}{\bm{R}^{W}_{S,ref}}(s)[\bm{\omega}_{{S,ref}}(s)]_{\times}ds,$
(9b)

where $v_{max}$ and $\omega_{max}$ are used to tune how fast the vehicle
moves. In accordance with the assumption of the human’s point of view from an
onboard camera, all references are provided in the body frame
$\mathcal{F}_{S}$.

###  III-B Feedback Generation

As stated in Section I, the design of adequate feedback wrenches is crucial to
ensure ease of operation and situational awareness for the human operator.
Hereby, the overall feedback wrench $\bm{\tau}_{fb,total}\in\mathbb{R}^{6}$
expressed in $\mathcal{F}_{M}$ is often a combination of multiple
contributions, each one representing a different aspect of the current task
(e.g. object avoidance, guiding towards a waypoint). In the context of this
work, we restrict our analysis to the well established recentering
$\bm{\tau}_{fb,rec}\in\mathbb{R}^{6}$ and interaction wrench
$\bm{\tau}_{fb,ext}\in\mathbb{R}^{6}$ feedback. Eventually, the total feedback
wrench takes the form

$\displaystyle\bm{\tau}_{fb,total}=\bm{\tau}_{fb,rec}+\bm{\tau}_{fb,ext},$
(10)

with $\bm{\tau}_{fb,rec}$ and $\bm{\tau}_{fb,ext}$ computed as explained
below.

####  III-B1 Recentering Wrench $\bm{\tau}_{fb,rec}$

When using rate control reference generation, the most essential type of
feedback is the recentering. The recentering wrench $\bm{\tau}_{fb,rec}$ in
$\mathcal{F}_{M}$ aims to move the haptic device’s end-effector back to its
idle pose, in other words make $\mathcal{F}_{H}$ identical to
$\mathcal{F}_{M}$:

$\displaystyle\bm{\tau}_{fb,rec}=-\bm{K}_{rec}\begin{bmatrix}\bm{p}_{H}\\\
\frac{1}{2}\left({\bm{R}^{M}_{H}}-{{\bm{R}^{M}_{H}}}^{\top}\right)^{\vee}\end{bmatrix},$
(11)

where the stiffness
$\bm{K}_{rec}=diag(\bm{K}_{rec,t},\bm{K}_{rec,r})\in\mathbb{R}^{6\times 6}$ is
a tuning parameter. Under rate control, recentering the end-effector will
cause the robot to slow down and eventually hold position and attitude. In
that sense, adding a virtual spring on the human side translates to the
addition of a virtual damper on the robot side. Without recentering feedback,
it would be almost impossible to manually zero the haptic device in all six
directions and achieve static hover. Additionally, it allows the human
operator to let go of the handle at any time and the robot will automatically
stabilize at its current pose.

When targeting applications involving physical contact with the environment
however, the recentering wrench is no longer sufficient to ensure situational
awareness. Especially in cases where the camera view might not be conclusive
about whether the robot is in contact or not, additional interaction-specific
information must be provided.

####  III-B2 Interaction Wrench $\bm{\tau}_{fb,ext}$

In state-of-the-art literature, interaction specific feedback involves
reflecting the measured or estimated forces at the contact point
$\bm{\tau}_{c,1:3}\in\mathbb{R}^{3}$ back to the operator, i.e. the forces
being applied by the environment to the aerial robot expressed in
$\mathcal{F}_{S}$. This work proposes an extension for omnidirectional
vehicles, whereby the torques $\bm{\tau}_{c,4:6}\in\mathbb{R}^{3}$ in
$\mathcal{F}_{S}$ acting on the vehicle at the contact point are also
included. Considering the offset between the vehicle’s center-of-mass and the
part of it that is in contact (i.e. the tool)
$\bm{r}_{O_{S}T}\in\mathbb{R}^{3}$ in $\mathcal{F}_{S}$, the interaction
wrench feedback is then given as:

$\displaystyle\bm{\tau}_{fb,ext}=\begin{bmatrix}\bm{\tau}_{c,1:3}\\\
\bm{\tau}_{c,4:6}\end{bmatrix}+\begin{bmatrix}\bm{0}\\\
\bm{r}_{O_{S}T}\times\bm{\tau}_{c,1:3}\end{bmatrix},$ (12)

where $\bm{\tau}_{fb,ext}$ is expressed in $\mathcal{F}_{M}$, again using the
assumption that $\mathcal{F}_{M}$ is aligned with $\mathcal{F}_{S}$. Notice
that this corresponds to the external wrench $\hat{\bm{\tau}}_{ext}$ acting on
the vehicle’s center-of-mass during interaction, effectively representing the
same wrench a human being would feel when holding the tool.

###  III-C Stability Considerations

In this paper, no formal proof of the stability of the teleoperation system is
provided. While this will be the focus of future work, some stability-related
aspects are briefly discussed here. The structure of our proposed framework,
namely rate control in combination with recentering and environment force
feedback, has some similarity with the scheme presented in [14]. In that
paper, stability of the teleoperation loop was proven when subject to bounded
operator and external forces. Although an underactuated system is considered,
that proof suggests the existence of similar formal guarantees for the case of
omnidirectional vehicles. Furthermore, the flight experiments presented in the
following section already verify the practical stability of the developed
teleoperation policy.

##  IV RESULTS

Figure 3: Decoupled translational (top) and rotational (bottom) reference
generation. Solid lines and corresponding shaded areas indicate the mean and
standard deviation of the operator inputs, respectively and dashed lines
indicate the mean total feedback wrench over five repeated trials. Coupling
effects related to unintended input agitation despite recentering wrench can
be observed, especially $\bm{y}_{S},\bm{z}_{S}$ coupling during translation
inputs and $\bm{x}_{S},\bm{y}_{S}$ translation during rotation inputs.

###  IV-A Experimental Setup

Indoor flight experiments are conducted to evaluate the capabilities and
performance of the proposed bilateral teleoperation setup. The employed system
is shown in Fig. 1, consisting of a 7DoF Franka Emika Panda arm with a handle
attached to the end-effector for a haptic device and the OMAV based on the
tiltrotor aerial platform introduced in [22]. A measure of the interaction
force between the human and the robot arm is obtained at
$800\text{\,}\mathrm{Hz}$ with a 6-axis Rokubi force-torque sensor mounted
between the end-effector and the handle. The Panda arm is running the default
cartesian velocity controller in combination with an admittance filter,
effectively rendering the closed-loop dynamics in (2). The admittance gains
are set as low as possible to improve haptic transparency, while still
ensuring a minimum dissipation to maintain system stability. Reference
generation and recentering parameters are tuned to render robot velocities and
feedback wrenches comfortable for the user. The specific values used during
experiments are listed in Table I. It should be noted that the framework can
easily be combined with other types of haptic devices, given they support 6DoF
motion input and wrench feedback rendering.

Parameter |  | Value |
---|---|---|---
${\bm{M}_{adm,t}}$ | ${\bm{M}_{adm,r}}$ | $10\bm{I}_{3\times 3}[$\mathrm{kg}$]$ | $\bm{I}_{3\times 3}[$\mathrm{k}\mathrm{g}\mathrm{m}^{2}$]$
${\bm{D}_{adm,t}}$ | ${\bm{D}_{adm,r}}$ | $5\bm{I}_{3\times 3}[$\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}$]$ | $\bm{I}_{3\times 3}[$\mathrm{k}\mathrm{g}\mathrm{s}^{-1}\mathrm{m}^{2}$]$
$v_{max}$ | $\omega_{max}$ | $1[${\mathrm{s}}^{-1}$]$ | $1[${\mathrm{s}}^{-1}$]$
$\bm{K}_{rec,t}$ | $\bm{K}_{rec,r}$ | $50\bm{I}_{3\times 3}[$\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$]$ | $2\bm{I}_{3\times 3}[$\mathrm{N}\text{\,}\mathrm{m}$]$
TABLE I: System parameters for experiments.

The remote environment is equipped with a MOCAP, providing pose measurements
for the OMAV at $100\text{\,}\mathrm{Hz}$. An EKF-based state estimator
provides the full state estimate
$\bm{p}_{S},\bm{v}_{S},{\bm{R}^{W}_{S}},\bm{\omega}_{S}$ required by the
impedance controller, fusing MOCAP with onboard IMU data (accelerometer and
gyroscope). The controller is tuned according to [11]. A safety tether is
connected to the OMAV but kept slack to limit external disturbances.
Additionally, a rigid rod of approximately $0.6\text{\,}\mathrm{m}$ length
with a soft ball at the end is attached to the robot and acts as the
interaction tool for contact-based tasks. The rod provides sufficient
clearance for the propellers to allow successful interaction without the risk
of collision. The ball acts as a mechanical damper to soften hard impacts.
Interaction forces are measured by a 6-axis Rokubi force-torque sensor mounted
between the OMAV and the rod. Hereby, the sensor data is transformed
accordingly to represent forces and torques acting on the vehicle center-of-
mass (see Section III-B).

Note that the human is directly looking at the robot, instead of viewing from
the robot perspective as mentioned in earlier sections. However, their visual
frame is still well aligned with $\mathcal{F}_{S}$, since the haptic device is
placed directly behind the robot and the experiments only involve small
attitude changes ($<$10\text{\,}\mathrm{\SIUnitSymbolDegree}$$).

###  IV-B Translational and Rotational Reference Generation

Recall that the aim of this paper is to analyze the suitability of standard
bilateral teleoperation methodologies when extended to 6DoF for
omnidirectional vehicles. As a first criteria, the operator must be able to
generate decoupled motion commands in all translational and rotational DoF, in
order to fully exploit the omnidirectional capabilities of the OMAV. This
requires the human to accurately render decoupled forces and torques at a
single interaction point, namely the handle of the haptic device. Thus,
evaluating whether they are physically and cognitively capable of performing
such manipulation is crucial for the feasibility of the proposed framework.
This is tested experimentally by repeatedly tasking an operator with
sequentially actuating each individual input axis without introducing motion
in other directions. The resulting translational and rotational reference
velocity statistics from five trials with a single operator over the
experiment duration $T$, as well as the total feedback wrench acting on the
handle are shown in Fig. 3. Note that this test only involves free flight
operation, meaning that the displayed feedback wrench only consists of
recentering actions, i.e. $\bm{\tau}_{fb,total}=\bm{\tau}_{fb,rec}$. The
effect of this recentering term is clearly visible, shown by the constant
wrench opposing the twist commands, aiming to restore the handle’s idle pose
in the local environment.

The results clearly show an unintended coupling between the different axes on
the input device. During the translational reference generation along
$\bm{y}_{S}$ ($\tilde{t}\in[0.15,0.3]$) for example, non-zero velocity
references in $\bm{z}_{S}$ can be observed. Similarly, when trying to move
along $\bm{z}_{S}$ only ($\tilde{t}\in[0.3,0.5]$), additional rotation along
$\bm{y}_{S}$ is accidentally introduced. A similar phenomenon is observed when
the user is tasked with performing rotations only. These coupling effects
worsen the performance when precise maneuvering is required, such as for high-
accuracy tasks or when operating in confined spaces.

In summary, it appears that producing decoupled reference inputs, especially a
pure rotation at the handle, are physically challenging for the operator,
despite the supporting recentering wrench. Making the recentering gain
adaptive could help to constrain the user input to a single axis at a time, by
making the remaining axes more stiff. Hereby, the adaptive solution must still
allow full exploitation of the omnidirectional capabilities. A detailed study
of such methodologies is left for future work.

###  IV-C Push-and-Slide Interaction

Apart from omnidirectional reference generation in contact-free conditions, a
bilateral teleoperation framework for an OMAV must allow the operator to
perform physical interaction tasks as well. We evaluate this requirement by
performing a push-and-slide operation with a whiteboard, as shown in Fig. 4.
This is a common task in contact-based inspection applications. During the
first phase of the experiment, the operator is asked to approach in a
direction normal to the whiteboard surface and push against it. In a second
phase, once contact with the board is established, the user is tasked with
sliding along $\bm{z}_{W}$. The resulting interaction wrench being fed back to
the user is shown in Fig. 5. Notice that the recentering wrench is not
included here, since it is not the primary focus of the experiment.
Additionally, the position of the OMAV in $\mathcal{F}_{W}$ is visualized to
highlight the motion with respect to the whiteboard surface.

When pushing against the board (highlighted in blue), a clear spike in
$-\bm{x}_{S}$ feedback force can be observed, indicating the presence and
intensity of the contact to the operator. The non-zero torque around
$-\bm{y}_{S}$ originates from the second term in (12) and is caused by a
misalignment between the surface normal $\bm{n}$ and $\bm{x}_{S}$ (see Fig. 4
with $\alpha<0$), resulting in an external pitching torque. During vertical
sliding, friction effects acting at the tool tip cause the same behavior.
While moving upwards (highlighted in orange), the tool lags behind due to the
high friction force. Since the tool and the connecting rod are rigidly
attached to the OMAV, this causes the vehicle to pitch down slightly,
producing the observed positive feedback torque around $\bm{y}_{S}$. The
opposite behavior occurs when sliding downwards (highlighted in yellow).

Compared to the omnidirectional reference generation, no immediate limitation
was detected when using the proposed extension of the standard interaction
force feedback. That being said, different experiments showed a degradation in
the performance of the flight controller, such as oscillations or tool
skipping, in the presence of large magnitude interaction and friction forces
(see also complementary video). In this regard, the effectiveness of the
provided feedback could be improved further by including information about the
robot state and its limitations.

Figure 4: Experimental setup for vertical sliding experiment along
$\bm{z}_{W}$. Uncompensated friction opposing the motion results in an angular
offset $\alpha$ between the surface norm $\bm{n}$ and the body x-axis
$\bm{x}_{S}$. Figure 5: Push-and-slide contact experiment. Contact with the
wall (the dashed line) and sliding along $\pm\bm{z}_{W}$ are shown. External
torques due to friction effects (see Fig. 4) can be observed.

##  V CONCLUSION

This work investigates the suitability of standard bilateral teleoperation
methods in the context of OMAV. Based on a straightforward extension of the
well established rate control, recentering and interaction force feedback
policy, a bilateral teleoperation framework for an omnidirectional aerial
robot has been designed and evaluated. The human is in contact with the handle
of a haptic device and performs small-scale deviations from the idle pose,
thereby generating twist commands for the robot. Wrench feedback is provided
to the operator, on the one hand recentering the handle to restore its idle
pose and on the other hand reflecting external forces acting on the vehicle
when in contact with the environment.

Practical experiments including contact-free flight, as well as push-and-slide
operation during physical interaction are conducted to evaluate the potential
of the proposed approach. Although the operator is able to control all six
axes of the OMAV, performing decoupled motion in a single DoF only is
practically challenging. This is a fundamental issue of the straightforward
extension of standard methodologies on this new types of platforms, which
shows the need for additional measures to suppress unintended inputs. Being
able to prevent coupling effects is absolutely mandatory, especially when
precise maneuvering is required for high-accuracy tasks and in confined
spaces. Future work will focus on addressing this problem through the use of
adaptive axes stiffness and human intention detection.

In summary, effective bilateral teleoperation for OMAV cannot be achieved by
simple extension of standard teleoperation methodologies but rather requires
more sophisticated policies to fully exploit the capabilities of these
vehicles.

## References

  * [1] M. Tognon, H. A. Tello Chávez, E. Gasparin, Q. Sablé, D. Bicego, A. Mallet, M. Lany, G. Santi, B. Revaz, J. Cortés, and A. Franchi, “A truly redundant aerial manipulator system with application to push-and-slide inspection in industrial plants,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 1846–1851, 2019.
  * [2] K. Bodie, M. Brunner, M. Pantic, S. Walser, P. Pfändler, U. Angst, R. Siegwart, and J. Nieto, “An omnidirectional aerial manipulation platform for contact-based inspection,” in _Proceedings of Robotics: Science and Systems_ , FreiburgimBreisgau, Germany, June 2019.
  * [3] A. Ollero, M. Tognon, A. Suarez, D. Lee, and A. Franchi, “Past, Present, and Future of Aerial Robotic Manipulators,” _IEEE Transactions on Robotics_ , vol. 38, no. 1, pp. 626–645, 2022.
  * [4] F. Ruggiero, V. Lippiello, and A. Ollero, “Aerial Manipulation: A Literature Review,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1957–1964, 2018.
  * [5] X. Meng, Y. He, and J. Han, “Survey on Aerial Manipulator: System, Modeling, and Control,” _Robotica_ , vol. 38, no. 7, p. 1288–1317, 2020.
  * [6] K. Bodie, M. Tognon, and R. Siegwart, “Dynamic end effector tracking with an omnidirectional parallel aerial manipulator,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 8165–8172, 2021.
  * [7] S. Park, J. Lee, J. Ahn, M. Kim, J. Her, G. Yang, and D. Lee, “Odar: Aerial manipulation platform enabling omnidirectional wrench generation,” _IEEE/ASME Transactions on Mechatronics_ , vol. 23, no. 4, pp. 1907–1918, 2018.
  * [8] D. Brescianini and R. D’Andrea, “An omni-directional multirotor vehicle,” _Mechatronics_ , vol. 55, pp. 76–93, 2018.
  * [9] M. Ryll, G. Muscio, F. Pierri, E. Cataldi, G. Antonelli, F. Caccavale, D. Bicego, and A. Franchi, “6D interaction control with aerial robots: The flying end-effector paradigm,” _The International Journal of Robotics Research_ , vol. 38, no. 9, pp. 1045–1062, 2019.
  * [10] M. A. Trujillo, J. R. Martínez-de Dios, C. Martín, A. Viguria, and A. Ollero, “Novel Aerial Manipulator for Accurate and Robust Industrial NDT Contact Inspection: A New Tool for the Oil and Gas Inspection Industry,” _Sensors_ , vol. 19, no. 6, 2019.
  * [11] K. Bodie, M. Brunner, M. Pantic, S. Walser, P. Pfandler, U. Angst, R. Siegwart, and J. Nieto, “Active Interaction Force Control for Contact-Based Inspection With a Fully Actuated Aerial Vehicle,” _IEEE Transactions on Robotics_ , pp. 1–14, 2020.
  * [12] A. Y. Mersha, S. Stramigioli, and R. Carloni, “On bilateral teleoperation of aerial robots,” _IEEE Transactions on Robotics_ , vol. 30, no. 1, pp. 258–274, 2014.
  * [13] F. Schill, X. Hou, and R. Mahony, “Admittance mode framework for haptic teleoperation of hovering vehicles with unlimited workspace,” in _2010 Australasian Conf. on Robotics & Automation, (Brisbane, Australia)_, 12 2010.
  * [14] H. Rifaï, M.-D. Hua, T. Hamel, and P. Morin, “Haptic-based bilateral teleoperation of underactuated Unmanned Aerial Vehicles,” _IFAC Proceedings Volumes_ , vol. 44, no. 1, pp. 13 782–13 788, 2011, 18th IFAC World Congress.
  * [15] C. Masone, M. Mohammadi, P. R. Giordano, and A. Franchi, “Shared planning and control for mobile robots with integral haptic feedback,” _The International Journal of Robotics Research_ , vol. 37, no. 11, pp. 1395–1420, 2018.
  * [16] G. Gioioso, M. Mohammadi, A. Franchi, and D. Prattichizzo, “A force-based bilateral teleoperation framework for aerial robots in contact with the environment,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , 2015, pp. 318–324.
  * [17] S. Islam, R. Ashour, and A. Sunda-Meya, “Haptic and Virtual Reality Based Shared Control for MAV,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 55, no. 5, pp. 2337–2346, 2019.
  * [18] J. Lee, R. Balachandran, Y. S. Sarkisov, M. D. Stefano, A. Coelho, K. Shinde, M. J. Kim, R. Triebel, and K. Kondak, “Visual-Inertial Telepresence for Aerial Manipulation,” 2020.
  * [19] A. Coelho, H. Singh, K. Kondak, and C. Ott, “Whole-body bilateral teleoperation of a redundant aerial manipulator,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 9150–9156.
  * [20] Y. Zimmermann, E. B. Küçüktabak, F. Farshidian, R. Riener, and M. Hutter, “Towards Dynamic Transparency: Robust Interaction Force Tracking Using Multi-Sensory Control on an Arm Exoskeleton,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020, pp. 7417–7424.
  * [21] F. Conti and O. Khatib, “Spanning large workspaces using small haptic devices,” in _First Joint Eurohaptics Conference and Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems. World Haptics Conference_ , 2005, pp. 183–188.
  * [22] M. Allenspach, K. Bodie, M. Brunner, L. Rinsoz, Z. Taylor, M. Kamel, R. Siegwart, and J. Nieto, “Design and optimal control of a tiltrotor micro-aerial vehicle for efficient omnidirectional flight,” _The International Journal of Robotics Research_ , vol. 39, no. 10-11, pp. 1305–1325, 2020.


# Independent Control and Path Planning of Microswimmers with a Uniform
Magnetic Field

###### Abstract

Artificial bacteria flagella (ABFs) are magnetic helical micro-swimmers that
can be remotely controlled via a uniform, rotating magnetic field. Previous
studies have used the heterogeneous response of microswimmers to external
magnetic fields for achieving independent control. Here we introduce
analytical and reinforcement learning control strategies for path planning to
a target by multiple swimmers using a uniform magnetic field. The comparison
of the two algorithms shows the superiority of reinforcement learning in
achieving minimal travel time to a target. The results demonstrate, for the
first time, the effective independent navigation of realistic micro-swimmers
with a uniform magnetic field in a viscous flow field.

###### keywords:

micro-swimmers, reinforcement learning, magnetically driven

Lucas Amoudruz Petros Koumoutsakos*

L. Amoudruz, Prof. P. Koumoutsakos
Computational Science and Engineering Laboratory, ETH Zürich, CH-8092,
Switzerland.
Email<EMAIL_ADDRESS>

L. Amoudruz, Prof. P. Koumoutsakos
John A. Paulson School of Engineering and Applied Sciences, Harvard
University, Cambridge, MA, USA.

##  1 Introduction

The magnetic control of micro-swimming devices [1, 2, 3, 4, 5] through micro-
manipulation [6, 7], targeted drug delivery [8, 9] or convection-enhanced
transport [10], has created new frontiers for bio-medicine. A particularly
promising technology involves corkscrew-shaped magnetic micro-swimmers
(artificial bacterial flagella) that propel themselves when subjected to a
rotating magnetic field [11]. Rotating magnetic fields can form propulsive
gradients and they are arguably preferable to alternatives, such as electric
fields, for in-vivo conditions [12, 13]. However, the independent, yet
coordinated, control of individual ABFs is challenging as it requires
balancing between the magnetic forces and the hydrodynamic interactions
between the swimmers while the employed magnetic fields are practically
uniform over lengths of few micrometers. We note that independent navigation
of mm-sized micro-swimmers has been shown in [14] through experiments and
simulations, while in [15] a reinforcement learning (RL) algorithm was applied
to adjust the velocity of an idealized swimmer in simulations with one way
coupling with a complex flow field. Control of swimmers using two-way coupling
and RL have been demonstrated with linked-spheres at low Reynolds numbers [16]
and for artificial fish in macro-scales [17]. Similarly, genetic algorithms
have been used to navigate micro-swimmers towards high concentrations of
chemicals [18].

The problem of heterogeneous micro-robots navigation via a uniform input has
been studied in two dimensions on surfaces [19] and in a fluid at rest [20].
The steering of two micro-propellers along two distinct paths in 3 dimensions
has been accomplished with the help of magnetic fields gradients [21]. These
advances exploited the heterogeneous response of micro-swimmers to a uniform
input to achieve independent trajectories along a prescribed path. These
control methods are based on short horizon objectives (stay on the prescribed
path) and do not provide the trajectory that minimizes the travel time to a
target position, particularly in the presence of a background flow. In
addition, strong background flows restrict the set of feasible paths for given
micro-swimmers. To the best of our knowledge the steering of multiple micron-
sized swimmers towards a target in a minimal time under a background flow and
a uniform magnetic field has not been reported before.

In this work, we present two methods to independently guide two micro-ABFs
towards a single target in the presence of a uniform magnetic field. The two
methods rely on simulations of swimming ABFs using an ordinary differential
equation (ODE) model. The model is calibrated with the method of dissipative
particle dynamics (DPD) [22, 23], taking into account the particular geometry
of the swimmers and their interactions with the viscous fluid. We first
present a semi-analytical solution for the simple yet instructive setup of
multiple, geometrically distinct ABFs in free space, with zero background
flow. This result enables understanding of the design constraints for the ABFs
necessary for independent control and how their geometric characteristics
relate to their travel time. We then employ RL to control multiple ABFs
trajectories in a broad range of flow conditions including a non-zero
background flow.

##  2 Artificial bacterial flagella

The ABFs are modeled as microscopic rigid bodies of length $l$ with position
$\mathbf{x}$ and orientation $q$ (represented by a quaternion), immersed in a
viscous fluid and subjected to a rotating, uniform, magnetic field. We
estimate that the magnetic and hydrodynamic interactions between ABFs are
orders of magnitude smaller than those due to the magnetic field for dilute
systems (see supplementary material) and we ignore inertial effects due to
their low Reynolds number ($\mathrm{Re}\approx${10}^{-3}$$). Following this
approximation, the system is fully described by the position and orientation
of the ABFs.

Additionally, the linear and angular velocities of the ABF, $\mathbf{V}^{b}$
and $\bm{\Omega}^{b}$, are directly linked to the external force and torque,
$\mathbf{F}^{b}$ and $\mathbf{T}^{b}$, via the mobility matrix [24],

$\begin{bmatrix}\mathbf{V}^{b}\\\
\bm{\Omega}^{b}\end{bmatrix}=\begin{bmatrix}\Delta&Z\\\
Z^{\text{T}}&\Gamma\end{bmatrix}\begin{bmatrix}\mathbf{F}^{b}\\\
\mathbf{T}^{b}\end{bmatrix},$ (1)

where the superscript b indicates that the quantity is expressed in the ABF
frame of reference, for which $\Delta$, $Z$ and $\Gamma$ are diagonal. The
matrices $\Gamma$ and $Z$ represent the application of torque to changing the
angular and linear velocity, respectively. The ABFs are propelled by torque
applied through a magnetic field and we assume that it can swim only in the
direction of its main axis so that $Z$ has only one non-zero entry ($Z_{11})$.
The coefficients in the mobility matrix are often estimated by empirical
formulas for low Reynolds number flows [25]. Here we estimate the components
of the mobility matrix for the specific ABF by conducting flow simulations
using DPD [22, 23], which we validate against experimental data of [8] (see
supplementary material). We remark that the shape (pitch, diameter, length,
thickness) of the ABF influence the elements of these matrices and the present
approach allows to account for these geometries.

The ABF with a magnetic moment $\mathbf{m}$ is subjected to a uniform magnetic
field $\mathbf{B}$ and hence experiences a torque

$\mathbf{T}=\mathbf{m}\times\mathbf{B}.$ (2)

No other external force is applied to the ABF, hence $\mathbf{F}=\mathbf{0}$.
Combining eq. 1 with the kinematic equations for a rigid body gives the
following system of ODEs:

$\displaystyle\dot{\mathbf{x}}$ $\displaystyle=\mathbf{V},$ (3a)
$\displaystyle\dot{q}$ $\displaystyle=\frac{1}{2}q\otimes\hat{\bm{\Omega}},$
(3b) $\displaystyle\mathbf{V}^{b}$ $\displaystyle=Z\mathbf{T}^{b},$ (3c)
$\displaystyle\bm{\Omega}^{b}$ $\displaystyle=\Gamma\mathbf{T}^{b},$ (3d)

where $\otimes$ denotes the quaternion product, and $\hat{\bm{\Omega}}$ the
pure quaternion formed by the vector $\bm{\Omega}$. The transformations
between the laboratory frame of reference and that of the ABF are given by:

$\displaystyle\mathbf{T}^{b}$ $\displaystyle=R(q)\mathbf{T},$ (4a)
$\displaystyle\mathbf{m}$ $\displaystyle=R(q^{\star})\mathbf{m}^{b},$ (4b)
$\displaystyle\mathbf{V}$ $\displaystyle=R(q^{\star})\mathbf{V}^{b},$ (4c)
$\displaystyle\bm{\Omega}$ $\displaystyle=R(q^{\star})\bm{\Omega}^{b},$ (4d)

where $q^{\star}$ is the conjugate of $q$ and $R(q)$ is the rotation matrix
that corresponds to the rotation by a quaternion $q$ [26]. The system of
differential equations (2) and (3) is advanced in time with a fourth order
Runge-Kutta integrator.

We note that when simulating multiple non-interacting ABFs in free space, we
use the above ODE system for each swimmer with the common magnetic field but
different mobility coefficients and magnetic moments.

##  3 Forward velocity

ABF were designed to swim under a rotating, uniform magnetic field [27, 11].
We first study this scenario by applying the field
$\mathbf{B}(t)=B\left(0,\cos{\omega t},\sin{\omega t}\right)$ to ABFs
initially aligned with the $x$ axis of the laboratory frame. Note that in the
later sections, the magnetic field is able to rotate in any direction so that
the swimmers can navigate in three dimensions. We consider two ABFs with the
same length but different pitch and magnetic moments, as shown in fig. 1. In
both cases, the magnetic moment is perpendicular to the helical axis of the
ABF. Under these conditions, by symmetry of the problem, the swimmers swim
along the $x$ axis. The difference in pitch results in different coefficients
of the mobility matrix and along with the different magnetic moments results
in distinct propulsion velocities for the two ABFs. For each ABF velocity we
distinguish a linear and a non-linear variation with respect to the frequency
of the magnetic field. First, the ABF rotates at the same frequency as the
magnetic field and its forward velocity increases linearly with the frequency
of the magnetic field, consistent with the low Reynolds approximation [28, 8,
29, 3]. In the non-linear regime, the magnetic torque is no longer able to
sustain the same frequency of rotation as the magnetic field. The onset of
non-linearity depends on the geometry and magnetic moment of the ABF as well
as the imposed magnetic field. Indeed, the magnitude of the magnetic torque is
bounded while that of the hydrodynamic torque increases linearly with the ABF
angular velocity $\Omega$. The torque imbalance at high rotation frequencies
causes the ABF to slip, resulting in an alternating forward and backward
motion (see supplementary material). Increasing the frequency further
increases the effective slip and accordingly decreases the forward velocity.
The two regimes are distinguished by the step-out frequency $\omega_{c}$
corresponding to the maximum forward velocity of the ABF.

Figure 1:  Left: Dimensionless time averaged forward velocity of two ABFs,
differing in shape and magnetic moment, against the field rotation frequency
(in units of the step-out frequency of the first swimmer, $\omega_{c,1}$).
Right: The ABFs geometries. The arrows represent the magnetic moment of the
ABFs.

The differences in propulsion velocities for the ABFs can be exploited to
control independently their trajectories. The slope $V/\omega$ in the linear
regime depends only on the shape of the ABF. The step-out frequency depends on
both the shape and the magnetic moment (it can also be changed by varying the
surface wetability of the ABF [30]). These two properties can be chosen such
that the forward velocities of two ABFs react differently to the magnetic
rotation frequency (fig. 1). By changing $\omega$, it is then possible to
control the relative velocities of the two ABFs: one is faster that the other
in one regime while the opposite occurs in an other regime. This simple
observation constitutes the key idea for independent control of several ABFs
even with a uniform magnetic field. We remark that, while this potential has
been previously identified [28, 30, 31, 32, 33], the control of similar
systems have been performed in the simple case of free space, non interacting
propellers and no background flow [21, 20]. To the best of our knowledge, this
is the first time that such independent controlled navigation of multiple
micro-swimmers is materialised in three dimensions with a complex background
flow. In the following sections we propose two methods to tackle the problem
of steering ABFs towards a target in a minimal amount of time.

##  4 Independent control I: semi-analytical solution

In the absence of an external flow field, we derive a semi-analytical strategy
for the navigation of $N$ ABFs towards a particular target. Each ABF has a
distinct magnetic moment and without any loss of generality, we set the target
position of all swimmers to the origin and define the initial position of the
$i^{\text{th}}$ ABF as $\mathbf{x}^{(i)}$. We assume that the time required by
one ABF to align with the rotation direction of the field is much smaller than
$|\mathbf{x}^{(i)}|/v$, where $v$ is the typical forward velocity of the ABF.
The proposed strategy consists in gathering all ABFs along one direction
$\mathbf{n}_{k}$ at a time, such that $\mathbf{x}^{(i)}\cdot\mathbf{n}_{k}=0$,
$i=1,2,\dots,N$ after phase $k$. We choose a sequence of orthogonal
directions, $\mathbf{n}_{k}\cdot\mathbf{n}_{k^{\prime}}=\delta_{kk^{\prime}}$.
The choice of the orientations of $\mathbf{n}_{k}$ is not restricted to the
basis vectors of the laboratory frame and is described at the end of this
section. In three dimensions, the strategy consists of three phases,
$k=1,2,3$, until all ABFs have reached their target: they first gather on a
plane, then on a line and finally to the target.

All ABFs are gathered along a given direction $\mathbf{n}_{k}$ by exploiting
the different forward responses of the ABFs when we alternate the frequency of
rotation of the magnetic field. More specifically, for $N$ ABFs, the field
rotates in the direction $\mathbf{n}_{k}$ for $t_{j}$ time units at frequency
$\omega_{c,j}$, $j=1,2,\dots,N$, where $\omega_{c,j}$ is the step-out
frequency of the $j^{\text{th}}$ swimmer. We define the velocity matrix with
elements $U_{ij}=V_{i}\left(\omega_{c,j}\right)$, denoting the velocity of
swimmer $i$ when the field rotates with the step out frequency of swimmer $j$.
We can relate the above quantities to the (signed) distances $d_{j}$ covered
by the ABFs as

$d_{i}=\sum\limits_{j=1}^{N}s_{j}t_{j}U_{ij},$

where $s_{j}\in\\{-1,1\\}$ determines if the field rotates
clockwise/counterclockwise. Equivalently, the vector form of the above is
$\mathbf{d}=U\mathbf{\beta}$, where $\beta_{j}=t_{j}s_{j}$. Setting
$d_{i}=\mathbf{x}^{(i)}\cdot\mathbf{n}_{k}$, we can invert this linear system
of equations for each phase $k$ and obtain the times spent at each step-out
frequency $\mathbf{\beta}=U^{-1}X\mathbf{n}_{k}$, where we have set
$X_{ij}=x^{(i)}_{j}$. We emphasize that this result holds only if the velocity
matrix is invertible, restricting the design of the ABFs to achieve
independent control. The total time spent at phase $k$ is then given by

$T(\mathbf{n}_{k})=\sum\limits_{i=1}^{N}t_{i}=\sum\limits_{i=1}^{N}|\beta_{i}|=\left\lVert
U^{-1}X\mathbf{n}_{k}\right\rVert_{1}.$

The yet unknown directions $\mathbf{n}_{k}$, $k=1,2,3$, are chosen to minimize
the total travel time. The directions are parameterized as
$\mathbf{n}_{k}=R(\phi,\theta,\psi)\mathbf{e}_{k}$, $k=1,2,3$, where
$R(\phi,\theta,\psi)$ is the rotation matrix given by the three Euler angles
$\phi$, $\theta$ and $\psi$. Note that this choice of handedness of the three
directions does not influence the final result. The optimal angles satisfy

$\phi^{\star},\theta^{\star},\psi^{\star}=\operatorname*{arg\,min}\limits_{\phi,\theta,\psi}{\sum\limits_{k=1}^{3}T\left(R(\phi,\theta,\psi)\mathbf{e}_{k}\right)}.$

We solve the above minimization problem numerically with derandomised
evolution strategy with covariance matrix adaptation (CMA-ES) [34] (see
supplementary material for the configuration of the optimizer).

##  5 Independent control II: Reinforcement Learning

We now employ a RL approach to solve the problem introduced in section 4. Each
of the $N$ ABFs is initially placed at a random position
$\mathbf{x}_{i}\sim\mathcal{N}\left(\mathbf{x}_{i}^{0},\sigma\right)$,
$i=1,2,\dots,N$. The RL agent controls the magnetic field frequency of
rotation and direction, and has the goal of bringing all ABFs within a small
distance (here two body lengths, $d=2l$) from the target origin. This small
distance is justified by the assumption of non-interacting ABFs. The agent
sets the direction and magnitude of the magnetic field frequency every fixed
time interval. An episode is terminated if either of the two conditions occur:
(a) all ABFs reached the target within a small distance $d$, or (b) the
simulation time exceeds a maximum time $T_{\text{max}}$. The positions
$\mathbf{x}_{i}$ and orientations $q_{i}$ of the ABFs describe the state $s$
of the environment in the RL framework. The action performed by the agent
every $\Delta t$ time encodes the magnetic field rotation frequency and
orientation for the next time interval. The reward of the system is designed
so that all ABFs reach the target and the travel time is minimized.
Additionally, a shaping reward term [35] is added to improve the learning
process. The training is performed using VRACER, the off-policy actor critic
RL method described in [36]. More details on the method can be found in the
supplementary material.

##  6 Reaching the targets

In this section, we demonstrate the effectiveness of the two methods
introduced in sections 4 and 5. We first consider 2 ABFs in free space with
zero background flow. Figure 2 shows the distance of the ABFs to their target
over time, and the corresponding magnetic field rotation frequency for both
methods. In both cases, the ABFs successfully reach their target.
Interestingly, the rotation frequencies chosen by the RL agent correspond to
the step-out frequencies of the ABFs. Indeed, these frequencies allow the
fastest absolute velocity difference between the ABFs, so it is consistent
that they are part of the fastest solution found by the RL method.
Furthermore, the RL trained swimmer was about $25\%$ faster than the semi-
analytical swimmer. We remark that the RL solution amounts to first blocking
the forward motion of one swimmer while the other continues swimming (see fig.
3). The blocked swimmer is first reoriented such that its magnetic moment is
orthogonal to the plane of the magnetic field rotation, thus the resulting
magnetic torque applied to this swimmer is zero. On the other hand, the method
presented in section 4 makes both ABFs swim at all time, even if one of them
must go further from its target. In such situations, the “blocking” method
found by RL is advantageous over the other method.

Figure 2: Distance to target of the two controlled ABFs (in units of body
length $l$) against dimensionless time (and ) in free space, zero background
flow, and corresponding magnetic field rotation frequency (), where
$\omega_{c,1}$ is the step-out frequency of the first swimmer. Figure 3:
Trajectories of the ABFs from their initial positions (ABF representations) to
the target area (sphere) obtained with the two methods in three dimensions:
semi-analytical (dotted lines) and RL (solid lines). The arrows show the
successive axes of rotation of the magnetic field. The size of the ABFs has
been scaled up by a factor of $7$, for visualization purpose.

We now employ the RL method in the case of 2 ABFs swimming in a background
flow with non zero velocity. The assumptions required for deriving the semi-
analytical approach are violated and therefore we do not use this approach in
this case.

In the presence of a background flow $\mathbf{u}_{\infty}$, eqs. 4c and 4d
become

$\displaystyle\mathbf{V}$
$\displaystyle=R(q^{\star})\mathbf{V}^{b}+\mathbf{u}_{\infty}(\mathbf{x}),$
$\displaystyle\bm{\Omega}$
$\displaystyle=R(q^{\star})\bm{\Omega}^{b}+\frac{1}{2}\nabla\times\mathbf{u}_{\infty}(\mathbf{x})+\frac{\lambda^{2}-1}{\lambda^{2}+1}\mathbf{p}\times\left(E(\mathbf{x})\mathbf{p}\right),$

where we approximated the rotation component by the effect of the flow on an
axisymmetric ellipsoid of aspect ratio $\lambda$ (Jeffery orbits). Here
$E(\mathbf{x})=\left(\nabla\mathbf{u}_{\infty}(\mathbf{x})+\nabla\mathbf{u}_{\infty}^{T}(\mathbf{x})\right)/2$
is the deformation rate tensor of the background flow evaluated at the
swimmer’s position and $\mathbf{p}=R(q^{\star})\mathbf{e}_{x}$ is the
orientation of the ellipsoid. We used $\lambda=2$ in the subsequent
simulations. The background flow is set to the initial conditions of the
Taylor-Green vortex,

$\mathbf{u}_{\infty}(\mathbf{r})=\begin{bmatrix}A\cos{ax}\sin{by}\sin{cz}\\\
B\sin{ax}\cos{by}\sin{cz}\\\ C\sin{ax}\sin{by}\cos{cz}\end{bmatrix},$ (5)

with $A=B=C/2=V_{1}(\omega_{c,1})$ and $a=b=-c=2\pi/50l$. With these
parameters, the maximum velocity of the background flow is larger than the
maximum swimming speed of the ABFs.

Figure 4: Left: Distance to target of the two controlled ABFs (in units of
body length $l$) against dimensionless time (and ) in free space with the
background flow described by eq. 5 and corresponding magnetic field rotation
frequency (), where $\omega_{c,1}$ is the step-out frequency of the first
swimmer. Right: Trajectories of the ABFs with non-zero flow obtained with the
RL method. The arrows represent the velocity field and the colors represent
the magnitude of the vorticity field. The flow field is only shown for a
distance less than $4l$ from the trajectories, where $l$ is the length of the
swimmers.

The distances between the swimmers and the target over time are shown for the
RL method on fig. 4. Despite the background flow perturbation, the RL method
successfully navigates the ABFs to their target. The magnetic action space
exhibits a similar behavior as in the free space case: the rotation frequency
of the magnetic field oscillates between the step out frequencies of both
swimmers and never exceeds the highest of these frequencies, where the
swimming performance would degrade considerably. The trajectories of the ABFs
seem to make use of the velocity field to achieve a lower travel time: fig. 4
shows that the trajectories tend to be parallel to the velocity field. The RL
method not only found a solution, but also made use of its environment to
reduce the travel time.

##  7 Robustness of the RL policy

The robustness of the RL method is tested against two external perturbations,
unseen during the training phase. In both cases, the robustness of the method
is measured in terms of success rate (expected ratio between the number of
successful trajectories and the number of attempts).

First, a flow perturbation $\mathbf{\delta
u}(\mathbf{r})=\varepsilon\mathbf{u}_{\infty}(\mathbf{r}/p)$ is added to the
background flow described in the previous section, where $\varepsilon$
controls the strength of the perturbation and $p$ controls the wave length of
the perturbation with respect to the original one. Figure 5 shows the success
rate of the RL approach against the perturbation strength $\varepsilon$ for
different $p$. For large wave lengths ($p=2$), the RL agent is able to
successfully steer the ABFs to their target in more than $90\%$ of the cases
when the perturbation strengths of less than $20\%$ of the original flow. In
contrast, the success rate degrades more sharply for smaller wave lengths
($p=1/2$), suggesting that the method is less robust for short wave length
perturbations. The RL policy seems more robust to perturbations with the same
wavelengths as the original flow ($p=1$) for large perturbation strengths: the
success rate is above $30\%$ even for large perturbations.

Figure 5: Left: Success rate of the RL method to guide swimmers to their
target against the flow perturbation strength $\varepsilon$, for different
wave numbers, $p=1/2$ (), $p=1$ (), $p=2$ (). Right: Success rate of the RL
method to steer swimmers to their target against thermal fluctuation
$k_{B}T/k_{B}T_{0}$.

At small length scales, micro-swimmers are subjected to thermal fluctuations.
We investigate the robustness of the RL policy (trained with the background
flow, eq. 5) on swimmers subjected to thermal noise and background flow (eq.
5). The thermal fluctuations are modeled as an additive stochastic term to the
linear and angular velocities of each swimmer, following the Einstein relation
with the mobility tensor given by eq. 1. Defining the generalized undisturbed
velocity $\bar{\mathcal{V}}=(\mathbf{V},\bm{\Omega})$, the resulting
stochastic generalized velocities satisfy

$\displaystyle\mathcal{V}$
$\displaystyle=\bar{\mathcal{V}}+\delta\mathcal{V},$
$\displaystyle\langle\delta\mathcal{V}_{i}\rangle$
$\displaystyle=0,\;\;i=1,2,\dots,6,$
$\displaystyle\langle\delta\mathcal{V}_{i},\delta\mathcal{V}_{j}\rangle$
$\displaystyle=k_{B}TM_{ij},\;\;i,j=1,2,\dots,6,$

where $M$ is the mobility tensor and $k_{B}T$ is the temperature of the
system, in energy units. The above property is achieved by adding a scaled
Gaussian random noise with zero mean to the velocities at every time step of
the simulation.

The success rate of the policy is shown in fig. 5 for various temperatures
$k_{B}T$, in units of the room temperature $k_{B}T_{0}$. As expected, a large
thermal noise causes the policy to fail at its task. Nevertheless, this
failure only occurs at relatively high temperatures: the success rate falls
below $50\%$ for $k_{B}T>25k_{B}T_{0}$, which is well above the normal
operating conditions of ABFs. With temperatures below $2k_{B}T_{0}$, the RL
method sustains a success rate above $99\%$. We remark that this robustness is
achieved successfully even with a policy trained with $k_{B}T=0$.

##  8 Conclusion

We have presented two methods to guide multiple ABFs individually towards
targets with a uniform magnetic field. The semi-analytical method allows to
understand the basic mechanisms that allow independent control and we derive
the necessary condition for the independent control of multiple ABFs: their
velocity matrix must be invertible, a condition that can be accommodated by
suitably choosing the geometries of the swimmers. This result may help to
optimize the shapes of ABFs to further reduce the travel time.

The RL approach can control multiple ABFs in quiescent flow as well as in the
presence of a complex background flow. Additionally, this approach is
resilient to small flow perturbations and to thermal noise. When the
background flow vanishes, the RL method recovers a very similar behavior as
the first method: the rotation frequency is alternating between the step-out
frequencies of the swimmers. Furthermore, the RL method could reach lower
travel time than the first method by blocking one swimmer while the other is
swimming. Steering an increased number of swimmers requires longer travel
times according to the first method. We thus expect that applying the RL
approach to more than two swimmers requires longer training times and might
become prohibitively expensive as the number of swimmers increases. Possible
solutions to this problem might include pre-training the RL agent with the
policy found by the semi-analytical method.

The current work focused on the simplified case of non-interacting swimmers.
In practice, the ABFs may interact hydrodynamically and magnetically with each
other, encounter obstacles, evolve in confined geometries or experience time
varying flows. Nevertheless, we expect the RL method to be a good candidate to
overcome these variants, in the same way as it naturally handled the addition
of a background flow.

Acknowledgments

We acknowledge insightful discussions with Guido Novati (ETHZ) and his
technical support for the usage of smarties. We acknowledge support by the
European Research Council (ERC Advanced Grant 341117).

Conflicts of Interest

The authors declare no financial or commercial conflicts of interest.

## References

  * [1] Q. Cao, X. Han, L. Li,  _Lab Chip_ 2014, _14_ 2762\.
  * [2] L. Yang, L. Zhang,  _Annual Review of Control, Robotics, and Autonomous Systems_ 2021, _4_ , 1 null.
  * [3] P. Tierno, R. Golestanian, I. Pagonabarraga, F. Sagués,  _The Journal of Physical Chemistry B_ 2008, _112_ , 51 16525.
  * [4] Y. Liu, D. Ge, J. Cong, H.-G. Piao, X. Huang, Y. Xu, G. Lu, L. Pan, M. Liu,  _Small_ 2018, _14_ , 17 1704546.
  * [5] P. Liao, L. Xing, S. Zhang, D. Sun,  _Small_ 2019, _15_ , 36 1901197.
  * [6] L. Zhang, K. E. Peyer, B. J. Nelson,  _Lab on a Chip_ 2010, _10_ , 17 2203.
  * [7] Y. Yu, J. Guo, Y. Wang, C. Shao, Y. Wang, Y. Zhao,  _ACS applied materials & interfaces_ 2020, _12_ , 14 16097.
  * [8] R. Mhanna, F. Qiu, L. Zhang, Y. Ding, K. Sugihara, M. Zenobi-Wong, B. J. Nelson,  _Small_ 2014, _10_ , 10 1953.
  * [9] P. Sharan, A. Nsamela, S. C. Lesher-Pérez, J. Simmchen,  _Small_ 2021, 2007403.
  * [10] S. Schuerle, A. P. Soleimany, T. Yeh, G. M. Anand, M. Häberli, H. E. Fleming, N. Mirkhani, F. Qiu, S. Hauert, X. Wang, B. J. Nelson, S. N. Bhatia,  _Science Advances_ 2019, _5_ , 4 1.
  * [11] L. Zhang, J. J. Abbott, L. Dong, B. E. Kratochvil, D. Bell, B. J. Nelson,  _Applied Physics Letters_ 2009, _94_ , 6 064107.
  * [12] H. Gu, Q. Boehler, D. Ahmed, B. J. Nelson,  _Science Robotics_ 2019, _4_ , 35.
  * [13] K. Bente, A. Codutti, F. Bachmann, D. Faivre,  _Small_ 2018, _14_ , 29 1704374.
  * [14] D. Wong, E. B. Steager, V. Kumar,  _IEEE Robotics and Automation Letters_ 2016, _1_ , 1 554.
  * [15] S. Colabrese, K. Gustavsson, A. Celani, L. Biferale,  _Physical Review Letters_ 2017, _118_ , 15 1.
  * [16] A. C. H. Tsang, P. W. Tong, S. Nallan, O. S. Pak,  _Phys. Rev. Fluids_ 2020, _5_ 074101\.
  * [17] S. Verma, G. Novati, P. Koumoutsakos,  _Proceedings of the National Academy of Sciences of the United States of America_ 2018, _115_ , 23 5849.
  * [18] B. Hartl, M. Hübl, G. Kahl, A. Zöttl,  _Proceedings of the National Academy of Sciences_ 2021, _118_ , 19.
  * [19] S. Floyd, E. Diller, C. Pawashe, M. Sitti,  _The International Journal of Robotics Research_ 2011, _30_ , 13 1553.
  * [20] P. J. Vach, S. Klumpp, D. Faivre,  _Journal of Physics D: Applied Physics_ 2015, _49_ , 6 065003.
  * [21] E. Diller, J. Giltinan, M. Sitti,  _The International Journal of Robotics Research_ 2013, _32_ , 5 614.
  * [22] P. Espanol, P. Warren,  _EPL (Europhysics Letters)_ 1995, _30_ , 4 191.
  * [23] D. Alexeev, L. Amoudruz, S. Litvinov, P. Koumoutsakos,  _Comput. Phys. Commun._ 2020, 107298.
  * [24] J. Happel, H. Brenner,  _Low Reynolds number hydrodynamics: with special applications to particulate media_ , volume 1,  Springer Science & Business Media, 1981.
  * [25] S. Kim, S. J. Karrila,  _Microhydrodynamics: principles and selected applications_ ,  Courier Corporation, 2013.
  * [26] B. Graf,  _arXiv preprint arXiv:0811.2889_ 2008.
  * [27] A. Ghosh, P. Fischer,  _Nano letters_ 2009, _9_ , 6 2243.
  * [28] K. E. Peyer, L. Zhang, B. J. Nelson,  _Nanoscale_ 2013, _5_ , 4 1259.
  * [29] D. Li, M. Jeong, E. Oren, T. Yu, T. Qiu,  _Robotics_ 2019, _8_ , 4 87.
  * [30] X. Wang, C. Hu, L. Schurz, C. De Marco, X. Chen, S. Pané, B. J. Nelson,  _ACS Nano_ 2018, _12_ , 6 6210.
  * [31] F. Bachmann, K. Bente, A. Codutti, D. Faivre,  _Physical Review Applied_ 2019, _11_ , 3 034039.
  * [32] I. S. Khalil, A. F. Tabak, Y. Hamed, M. E. Mitwally, M. Tawakol, A. Klingner, M. Sitti,  _Advanced Science_ 2018, _5_ , 2 1700461.
  * [33] I. S. Khalil, A. F. Tabak, Y. Hamed, M. Tawakol, A. Klingner, N. El Gohary, B. Mizaikoff, M. Sitti,  _IEEE Robotics and Automation Letters_ 2018, _3_ , 3 1703.
  * [34] N. Hansen, S. D. Müller, P. Koumoutsakos,  _Evolutionary computation_ 2003, _11_ , 1 1.
  * [35] A. Y. Ng, D. Harada, S. Russell,  In _ICML_ , volume 99. 1999 278–287.
  * [36] G. Novati, P. Koumoutsakos,  In _Proceedings of the 36 th International Conference on Machine Learning_. 2019 .

  *[ABFs]: artificial bacterial flagellum
  *[RL]: reinforcement learning
  *[ODE]: ordinary differential equation
  *[DPD]: dissipative particle dynamics
  *[ABF]: artificial bacterial flagellum
  *[ODEs]: ordinary differential equation
  *[CMA-ES]: derandomised evolution strategy with covariance matrix adaptation