# Size-dependent ferroelectric-to-paraelectric sliding transformations and
antipolar-to-ferroelectric topological phase transitions in binary
homobilayers

Alejandro Pacheco-Sanjuan,1 Pradeep Kumar,2 and Salvador Barraza-Lopez2,3
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1\. Department
of Mechanical Engineering, Universidad Técnica Federico Santa María,
Valparaíso, Chile
2\. Department of Physics, University of Arkansas, Fayetteville, AR 72701, USA
3\. MonArk NSF Quantum Foundry, University of Arkansas, Fayetteville, AR
72701, USA

###### Abstract

The recent discovery of ferroelectric behavior in few-layer materials,
accompanied by the observation of antipolar domains in hexagonal boron nitride
and transition metal dichalcogenide moiré bilayers, is paving the way for
revolutionary advancements in the generation and manipulation of intrinsic
electric dipoles through stacking. In addition, these cutting-edge quantum
materials are reshaping our comprehension of phase transitions. Within the
present study, we unveil a hitherto unreported size-dependent sliding behavior
that marks a significant departure from conventional ferroelectrics. We also
shed light on thermally induced spontaneous hyperlubric sliding within moiré
bilayers, which can be used as a signal to distinguish topological phase
transitions from an antipolar onto a ferroelectric bilayer. Our findings also
suggest that the (topological) pinning of AA nodes in antipolar moiré
homobilayers prevents the occurrence of an antipolar-to-paraelectric
transformation.

##  I Introduction

Two-dimensional (2D) ferroelectrics are fundamentally different from ternary
bulk ferroelectric oxides. Those differences stem from (i) a chemistry that
now includes binary [1, 2, 3, 4] or even elemental [5, 6, 7, 8, 9] materials,
(ii) the presence of new mechanisms for removing centers of inversion such as
relative rotations within a homobilayer [2] or by the simple creation of a
heterobilayer, and (iii) by the presence of ferroelectric-to-paraelectric
phase transitions [10] by relative sliding [11, 12] (whereby the putative
paraelectric structure is a time average over ferroelectric configurations
that change the direction of their intrinsic electric dipole at discrete
sliding events). Bilayer ferroelectrics are created by a $60^{\circ}$ rotation
of one binary monolayer on a bilayer stack, away from a ground state
configuration [2], and they feature a macroscopically large number of
degenerate minima [12].

Despite the existence of works displaying double-well energy potentials [12,
11] and of analytical calculations describing ferroelectric-to-paraelectric
phase transitions in those homobilayers, the size-dependence of the propensity
for sliding has not been addressed to date. Furthermore, it is difficult to
realize a precise $60^{\circ}$ rotation among bilayers [13], and a more likely
outcome is to achieve rotations by $60^{\circ}-\delta$, where $\delta$ is a
small angle. As it is well-known by now, such mismatch gives rise to moirés,
which have originally been studied for their unique electronic properties [14,
15, 16, 17, 18, 19, 20, 21, 22, 23] but also display antipolar domains on
bilayers made from binary compounds [24, 25, 26, 10, 27, 28, 29]. According to
Bennett and coworkers, binary moiré homobilayers are antipolar, not
antiferroelectric [30, 31]. As superbly described by Cazeaux and coworkers,
such antipolar structures feature AB and BA triangular domains separated by
dislocation lines meeting at AA-nodes [29]. Bennett and coworkers argue that
the swap of polarization among AB and BA domains is tied to a topological
winding of the local intrinsic dipole moment (which thus also includes in-
plane components), with a meron and antimeron texture at AB and BA domains,
respectively [32]. Such topological behavior may be readable through
superlubricity, but a direct comparison of the sliding behavior of
ferroelectric and antipolar bilayers has not been presented yet.

In addition, there should be a coupling of such intrinsic electric dipole with
non-trivial topology and the inherent topology induced by strain in moiré hBN
bilayers [33, 34, 35, 26], which have one AA node per moiré supercell. The
density of domains cannot be changed spontaneously but only by a change of the
angle of rotation, strain, or a combination of both [28, 29]. Atomistic
reconstructions take place such that AB and BA domains occupy most of the
bilayer’s area, while AA nodes [29] are topological in nature [33, 34, 35, 36,
37, 38, 39, 40, 41], and their number cannot be modified once a relative angle
of rotation has been set up. Molecular dynamics calculations were carried out
at the highest allowed temperature in a molecular dynamics code to shed light
on this question.

This study is based on molecular dynamics calculations, whose methodology is
described in Sec. II. The results appear in Sec. III and conclusions in Sec.
IV.

Figure 1: Initial hBN bilayer structures after atomistic optimization at 1 K;
false color indicates the relative height among layers, and the dimmest yellow
locations indicate AA-nodes (which are topological defects; see Ref. [28]).
(a) A ferroelectric ($\delta=0^{\circ}$) sample initially set into an AA
relative configuration. Moiré hBN bilayers with (b) $\delta=1.02^{\circ}$ and
50444 atoms, (c) $\delta=2.01^{\circ}$ and 12944 atoms, (d)
$\delta=3.15^{\circ}$ and 5216 atoms, and (e) $\delta=3.89^{\circ}$ and 3408
atoms. All structures were drawn using the same length scale for a direct
comparison.

##  II Methods

Atomistic models for hBN bilayers were constructed with an initial interatomic
distance of 1.47 Å and an interlayer distance of 3.33 Å. The models consisted
of commensurate periodic hBN bilayers with a relative rotation angle of
$60^{\circ}-\delta$ degrees, with $\delta=0^{\circ}$ for a ferroelectric
configuration, or $\delta=1.02^{\circ}$, $2.01^{\circ}$, $3.15^{\circ}$, and
$3.89^{\circ}$, for moiré antipolar configurations. The studied structures
contain 50,444, 12,944, 5,216, and 3,408 atoms, respectively. Two additional
ferroelectric bilayers with either 200 or 400 atoms were also considered to
better understand the drastically different sliding behaviors of ferroelectric
and antipolar bilayers.

Structures were initially set into an AA configuration, and forces were
allowed to equilibrate in this position of (unstable) equilibrium in LAMMPS
[42]. Molecular dynamics simulations were performed with the same numerical
tool afterwards. Intr-alayer energetics were described by the Tersoff
potential with parameters fitted for BN structures [43, 44]. This interatomic
potential accurately replicates the strain energy response, equilibrium
lattice constant, and phonon dispersion relations of boron nitride
nanostructures with data obtained from $x-$ray scattering experiments and
density functional theory (DFT) calculations. The anisotropic interlayer
potential (ILP) for h-BN was utilized [45] to represent the complex interlayer
interactions. The velocity Verlet scheme was used for the time-integration of
the resulting equations of motion with a constant time step of 0.001 ps. The
interlayer energy $E_{int}$ was minimized until reaching a force threshold of
$10^{-6}$ eV/Å. This permits elucidating AA nodes, which occur at the longest
interlayer separation. An NPT equilibration stage (in which the number of
atoms is fixed, pressure is set to 0 Pa, and temperature oscillate around 700
K) of $5\times 10^{4}$ time steps was performed using the Nosé-Hoover
thermostat [46, 47]. Last, a production stage was set in an NVE ensemble
(where the number of atoms, the volume, and total energy remain constant) for
$5\times 10^{6}$ time steps (5 ns) and a sampling rate of 0.5 ps. Each
monolayer’s center of mass’ coordinates were recorded to detect their relative
motion. Additional molecular dynamics calculations were carried out at 300 K
on a moiré bilayer with a $3.15^{\circ}$ relative rotation to determine the
evolution of diffusion coefficients with temperature.

Figure 2: Sliding of ferroelectric ($\delta=0^{\circ}$) hBN bilayers at 700 K
starting from an AA relative configuration. Bilayers containing 200 and 400
atoms (subplots (a) and (b)) undergo telegraph-noise-like discrete
displacements $\Delta x$ and $\Delta y$ that sample the honeycomb lattice,
similar to those described in Ref. [12]. The right panels contain the squared
relative displacement of the monolayer’s centers of mass as a function of time
$\langle\Delta r(\tau)\rangle^{2}$. As bilayers increase their size (subplots
(c) through (f)), the single shear mode underpinning a whole concerted
relative motion by one interatomic distance becomes suppressed, as seen by the
straight lines and the comparatively smaller in-plane sliding scale. This
means that a critical ferroelectric-to-paraelectric critical temperature
cannot be uniquely assigned, as it depends on the bilayer’s size. In subplots
(c) through (f), the bilayer quickly moves from the unstable AA configuration
into a stable local minima (either AB or BA), and it does not move from that
local minima for the extent of the simulations; inserts show the relative
displacement of one monolayer as the bilayer falls into the energy minima.

##  III Results

###  III.1 Size-dependent sliding propensity

Structures optimized in forces at zero temperature are depicted in Fig. 1.
Fig. 1(a) is a ferroelectric hBN bilayer of size similar to a moiré bilayer
with $\delta=3.15^{\circ}$. The green color indicates a constant relative
height among monolayers for the ferroelectric configuration. Figures 1(b)
through 1(e) depict moiré bilayers with increasing $\delta$ (and a concomitant
decrease in the number of atoms $n$). The regions in which individual
constituent monolayers are farthest apart are colored in yellow. Those are
energetically unfavorable AA-stacked nodes [48, 29].

A precision concerning size effects and periodic boundary conditions is
necessary: Even when under periodic boundary conditions, a system with $N$
atoms can only support 3N vibrational modes. Within the present context, the
periodic boundary serves the goals of avoiding (i) dangling chemical bonds and
(ii) large out-of-plane excursions commonly seen at the edges of (non-
periodic) finite-size flakes.

Figure 2 contains the first original result from this work: Refs. [11] and us
[12] stated specific values for a critical temperature ($T_{C}$) on
ferroelectric sliding bilayers by the onset of sliding events and–starting
from an unstable AA configuration–Figs. 2(a) and 2(b) do show discrete,
telegraph-noise in-plane sliding ($\Delta x$,$\Delta y$) events (see left
panels in Figures 2(a) and 2(b)) consistent with discrete displacements
resembling the honeycomb lattice on molecular dynamics calculations at 700 K
when the number of atoms is either 200 or 400 (see the sampling of the
honeycomb lattice on the right subplots in Figs. 2(a) and 2(b)). Nevertheless,
the frequency of those steps decreases with the number of atoms (as do the
relative displacements of the center of mass) until the frequency of those
displacements turns into a complete halt (see the straight lines that imply no
discrete jumps on Figs. 2(c) through 2(f), where the number of atoms ranges
from 3,408 to 50,444). This is to say that ferroelectric bilayers by sliding
should have a propensity to slide that depends not just on an inherent energy
barrier, but also on the sample’s size. Similar size-dependent sliding events
have been reported on graphene bilayers [49].

The observed phenomena can be explained as follows: assuming that barrier
crossing can be described by one dimensional Smoluchowski equation [50], the
rate of sliding events of ferroelectric bilayers $r_{s}$ could be written as:

$r_{s}=\frac{D\beta\omega_{\rm min}\omega_{\rm max}}{2\pi}e^{-\beta\Delta U},$
(1)

where $D$ is the diffusion coefficient (to be determined momentarily),
$\beta=1/k_{B}T$ is the inverse temperature, $\omega_{\rm min}$ and
$\omega_{\rm max}$ are the frequencies at the minimum and the maximum of the
potential well, and $\Delta U$ is the barrier height. The dependency of $D$ on
the size of the bilayer, and on its ferroelectric or moiré antipolar phase,
will be studied soon.

To quantify the relative sliding of monolayers seen in Figs. 2(a) and 2(b), we
calculated the square of the displacement of the atom nearest to the center of
mass on one monolayer as a function of time $\langle\Delta r(\tau)\rangle^{2}$
defined as:

$\langle\Delta
r(\tau)\rangle^{2}=\langle|\mathbf{r}_{0}(\tau+t)-\mathbf{r}_{0}(t)|\rangle^{2}$
(2)

(where $\mathbf{r}_{0}(t)$ is the relative sliding of the center of mass among
monolayers at time $t$), and shown within insets in the right subpanels of
Figs. 2(a) and 2(b).

###  III.2 Superlubric moiré bilayers

The main difficulty in assembling moirés is the sudden, spontaneous unwanted
relative motion of monolayers [13] due to a loss of commensuration [51] upon
rotation [52]. A direct comparison among Figs. 2 and 3 illustrates the
different sliding behavior of ferroelectric and antipolar binary bilayers,
which provides a mechanism to tell them apart.

To begin with, ferroelectric structures with 3408, 5216, 12944, and 50444
atoms have all commensurate unit cells containing four atoms which turn them
“rougher” against relative thermally-induced displacements. This is the
meaning of the straight lines on Fig. 2(c)-(f). The left subplots in Fig. 3
display the in-plane sliding of the two monolayers in a moiré configuration.
Even though they have the same number of atoms than the structures in Figs.
2(c)-(f), a relative sliding without discrete telegraph-like jumps can be
clearly observed. Even more, the displacements occur over a range of 100 Å,
much larger than that seen in ferroelectric bilayers. The larger sliding of
moiré bilayers with a lack of registry of the honeycomb lattice is evident in
the right panels of Fig. 3, in which a chosen atom moves in what appears to be
a continuous fashion.

Figure 3: Sliding of moiré hBN bilayers at 700 K. Similar to Figs. 2(a) and
2(b), the left sub-panels show in-plane sliding displacements. The relative
motion of the centers of mass of the individual monolayers does not track the
honeycomb lattice anymore (compare right panels to those in Fig. 2(a) and
2(b)), and the diffusion is orders of magnitude larger due to a lack of atomic
registry upon relative rotation. Panels to the right contain $\langle\Delta
r(\tau)\rangle^{2}$ as well and permit observing a suppressed diffusion as
size increases. Figure 4: Sliding of the moiré hBN bilayers with
$\delta=3.15^{\circ}$ at 300 K. The left sub-panels show in-plane sliding
displacements. The right panel contains $\langle\Delta r(\tau)\rangle^{2}$
and, when contrasted against Fig. 3(b), permit observing a suppressed
diffusion as temperature decreases.

We produced additional molecular dynamics calculations at 300 K for the
bilayer at $3.15^{\circ}$ (containing 5216 atoms) and present those results in
Fig. 5. Contrasted against the results shown in Fig. 3(b), the motion is
restricted at the lower temperatures. Indeed, the diffusion coefficient $D$ is
reduced from 0.0879 Å2/ps at 700 K, down to 0.0727 Å2/ps at 300 K.

Figure 5: Diffusion coefficient $D(n)$ obtained from the slope of
$\langle\Delta r(\tau)\rangle^{2}=4D(n)\tau$. A power-law fit is included. The
diffusion coefficients for two smaller ferroelectric bilayers are shown as
well; note how they are orders of magnitude comparatively smaller for a given
number of atoms.

The differences in sliding can be further ascertained by a comparison of the
slopes of the $\langle\Delta r(\tau)\rangle^{2}$ versus $\tau$ plots, which
depend on the number of atoms $n$. The idea is that, for diffusive dynamics,
$\langle\Delta r(\tau)\rangle^{2}$ is related to the 2D diffusion constant
$D(n)$ by:

$\langle\Delta r(\tau)\rangle^{2}=4D(n)\tau.$ (3)

$D(n)$, plotted in Fig. 5, further confirms the striking difference in
relative sliding among ferroelectric and antipolar bilayers. In connection to
the ferroelectric-to-paraelectric transition in ferroelectric bilayers by
sliding, one observes a quick decay of $D$ to zero as the number of atoms
increases on the ferroelectric bilayer (inset on Fig. 5). The frequency of
sliding effects $r_{s}$ is then shown to go to zero in Eqn. (1) along with
$D(n)$, such that the critical temperature for a ferroelectric-to-paraelectric
transition by sliding is size-dependent.

An analysis of friction of moiré bilayers is provided next.

###  III.3 Friction within moiré bilayers

We first calculate the coefficient of friction $\mu$ between the monolayers
from our simulations (i.e., at $T=700$ K). $\mu$ is defined here as:

$\mu=\left<\frac{F_{f}}{F_{z}}\right>,$ (4)

where $F_{z}\geq 0$ and $F_{f}$ are the instantaneous normal and friction
forces, respectively, and $\left<\cdot\right>$ denotes a time average.

Figure 6: Dependence of (a) average friction force $\left<F_{f}\right>$, (b)
average normal force $\left<F_{N}\right>$, (c) average interlayer energy
$\left<E_{\rm int}\right>$ and (d) coefficient of friction $\mu$ on the number
of atoms (i.e., on the angle $\delta$) for moiré hBN bilayers. The recorded
values of $\mu$ imply that these are superlubric.

The variation of $\left<F_{z}\right>$, $\left<F_{f}\right>$, average
interlayer energy $\left<E_{\rm int}\right>$, and the coefficient of friction
$\mu$ as a function of the number of atoms $n$ (and, implicitly, on the
relative rotation angle $\delta$) for moiré hBN bilayers are depicted in Fig.
6. The friction and normal forces increase with $n$, but they follow different
trends. The coefficient of friction remains in the superlubricity regime
($\mu<0.01$) [53] for all values of $n$ ($\delta$). Rue et al. determined
$\mu=0.005$ for MoS2/MoSe2 heterostructures [54] at 300 K. for
$\delta=5^{\circ}$ and $\langle F_{z}\rangle=0.2$ nN/atom. In our calculations
at 700 K and for the largest angle we studied ($\delta=3.89^{\circ}$),
$\langle F_{z}\rangle=0.0183$ nN/atom, and the coefficient of friction was
$0.0016$, which is comparable with that of the MoS2/MoSe2 heterostructure.

At this point, we use the Green-Kubo formalism to calculate the inter-layer
friction constant $\lambda$. There, the interlayer friction constant $\lambda$
is defined as [55, 56, 57, 58, 59]:

$\lambda=\lim_{t\to\infty}\frac{1}{2k_{B}TA}\int_{0}^{t}c(t^{\prime})dt^{\prime},$
(5)

where $k_{B}$ is the Boltzmann constant, $T$ is the temperature, $A$ is the
contact area, and $c(t)$ is the lateral (in-plane) force correlation function
defined as:

$c(t)=\left<\mathbf{F}_{||}(t)\mathbf{F}_{||}(0)\right>.$ (6)

The force autocorrelation function $c(t)$ is shown in Fig. 7(a) for moiré hBN
bilayers of different sizes ($\delta$) at $T=700$ K. For all values of $n$,
$c(t)$ decreases rapidly and decays to zero at the time-scale of order of
$200$ ps. In order to investigate the frequency response of the friction
constant, we next looked the Fourier transform of $c(t)$:

$c(\omega)=\int_{0}^{\infty}c(t)e^{-i\omega t}dt.$ (7)

$c(\omega)$ was calculated for different values of $\delta$ in Fig. 7(b). We
find that $c(\omega)$ exhibits two peaks–one at lower frequency ($\approx 0.1$
THz), and another at slightly higher frequency ($\approx 2.0$ THz). The
amplitude of the lower frequency peak becomes smaller and moves to slightly
higher frequency, while the amplitude of the high-frequency peak becomes
weaker upon decreasing $n$. In Fig. 7(c), we show the behavior of $\lambda(t)$
for four different values of $n$. The values of the friction constant obtained
from the long-time behavior of $\lambda(t)$ are shown as a function of $n$ in
Fig. 7(d); $\lambda$ increases monotonically with increasing $n$. The behavior
of the friction constant is consistent with the increasing value of friction
force with $n$ discussed in Fig. 6(a).

Figure 7: (a) Time-autocorrelation of the lateral force between the
monolayers. (b) Fourier transform of the force autocorrelation function. (c)
Time dependence of $\lambda(t)$. (d) Friction constant $\lambda$ as a function
of the number of atoms $n$. Figure 8: The number of topological defects (dark-
colored patches corresponding to AA-stacked sections [28, 29]) for the moiré
hBN-bilayer with $\delta=1.02^{\circ}$ at (a) 2.05 ns, (b) 2.15 ns, and (c)
2.30 ns and at 700 K. Subplots (d) through (f) display the moiré bilayer at
700 K at $2.01^{\circ}$ at 4.295 ns, 4.345 ns, and 4.395 ns, respectively. In
both cases, a diameter of about 50 Å can be observed for the AA-nodes.

###  III.4 Observations on an antipolar to paraelectric phase transitions

The network of AB and BA domains forms a network separated by dislocation
lines meeting at AA nodes [28, 29]. Engelke and coauthors posited that one can
map the local relative sliding configurations in individual unit cells though
the moiré into a closed space akin to a punctured torus, where the punctured
section is precisely the AA-node. In this map, the torus is the closed surface
in which the winding of a local displacement field turns out to be topological
[28]. And so, the last item for discussion is the possibility of a mechanism
by which the density of topologically-protected AA-nodes alluded on the
previous paragraph [28, 29] could be thermally altered. That possibility was
the reason to perform molecular dynamics calculations at 700 K in the first
place (which was the highest we could achieve without facing numerical
instabilities in the code).

Figure 8 represents a qualitative answer to the question. Figures 8(a-c) show
the relative height among the two monolayers in the $\delta=1.02^{\circ}$
configuration at three different times (2.05 ns, 2.15 ns, and 2.30 ns) along
the molecular dynamics evolution. The darkest hue is a 3.5 Å cutoff, above
which the local bilayer has a local AA-stacking. The point is that, though
those move around, four AA-stacked sections always remain visible on the
snapshots.

We expect the antipolar behavior described by Bennett et al. [32] to still be
present at high temperature. Then, we showed by molecular dynamics
calculations that the AA-nodes cannot be thermally removed. Therefore, to
produce a transition onto a ferroelectric or paraelectric phase, the only
options are either (i) a quantum phase transition, or (ii) a topological one.
As the temperature is high, a quantum phase transition is discarded as well,
and a topological transition signified by the removal of AA-nodes is the only
possibility left. Indeed, the rotation by $-\delta$ renders the antipolar
moiré onto a ferroelectric configuration which may be unstable to sliding
events depending on the samples’ size and temperature.

##  IV Conclusion

To conclude, we undertook a study of hBN bilayers in the ferroelectric
($\delta=0^{\circ}$) and moiré antipolar ($\delta>0^{\circ}$) configurations.
We first established the fact that ferroelectric-to-paraelectric transitions
by sliding occur at temperatures dependent on system size. Then, using
configurations with identical number of atoms, we establish an easier
propensity to slide of antipolar moiré bilayers. This study included the
calculation of the diffusion coefficient, and of parameters $\mu$ and
$\lambda$ that indicate superlubric behavior. The work ends by the observation
that the area density of topological defects (given by AA-nodes) remains
protected at the highest temperature of 700 K, in which molecular dynamics
simulations could run without developing numerical instabilities. Therefore,
the only possibility for the removal of AA-nodes is a rotation out of the
topologically-nontrivial moiré conformation onto a (trivial) configuration
with a single realization of relative sliding among monolayers across all unit
cells.

##  V Acknowledgments

S.B.-L. acknowledges funding from the U.S. Department of Energy (Award DE-
SC0022120).

## References

  * Chang et al. [2016] K. Chang, J. Liu, H. Lin, N. Wang, K. Zhao, A. Zhang, F. Jin, Y. Zhong, X. Hu, W. Duan, et al., Science 353, 274 (2016).
  * Li and Wu [2017] L. Li and M. Wu, ACS Nano 11, 6382 (2017).
  * Yang et al. [2018a] Q. Yang, M. Wu, and J. Li, J. Phys. Chem. Lett. 9, 7160 (2018a).
  * Barraza-Lopez et al. [2021] S. Barraza-Lopez, B. M. Fregoso, J. W. Villanova, S. S. P. Parkin, and K. Chang, Rev. Mod. Phys. 93, 011001 (2021).
  * Xiao et al. [2018] C. Xiao, F. Wang, S. A. Yang, Y. Lu, Y. Feng, and S. Zhang, Adv. Funct. Mater. 28, 1707383 (2018).
  * Gou et al. [2023] J. Gou, H. Bai, X. Zhang, Y. L. Huang, S. Duan, A. Ariando, S. A. Yang, L. Chen, Y. Lu, and A. T. S. Wee, Nature 617, 62 (2023).
  * Garcia-Ruiz et al. [2023] A. Garcia-Ruiz, V. Enaldiev, A. McEllistrim, and V. I. Fal’ko, Nano Lett. 23, 4120 (2023).
  * Yang et al. [2023] L. Yang, S. Ding, J. Gao, and M. Wu, Phys. Rev. Lett. 131, 096801 (2023).
  * Salleh Atri et al. [2023] S. Salleh Atri, W. Cao, B. Alon, N. Roy, M. Vizner Stern, V. Fal’ko, M. Goldstein, L. Kronik, M. Urbakh, O. Hod, et al. (2023), eprint 22305.10890.
  * Liu et al. [2022] Y. Liu, S. Liu, B. Li, W. J. Yoo, and J. Hone, Nano Lett. 22, 1265 (2022).
  * Tang and Bauer [2023] P. Tang and G. E. W. Bauer, Phys. Rev. Lett. 130, 176801 (2023).
  * Marmolejo-Tejada et al. [2022] J. M. Marmolejo-Tejada, J. E. Roll, S. P. Poudel, S. Barraza-Lopez, and M. A. Mosquera, Nano Lett. 22, 7984 (2022).
  * Lau et al. [2022] C. N. Lau, M. W. Bockrath, K. F. Mak, and F. Zhang, Nature 602, 41 (2022).
  * Lopes dos Santos et al. [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Li et al. [2010] G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong, and E. Y. Andrei, Nat. Phys. 6, 109 (2010).
  * Trambly de Laissardiére et al. [2010] G. Trambly de Laissardiére, D. Mayou, and L. Magaud, Nano Lett. 10, 804 (2010).
  * Suárez Morell et al. [2010] E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
  * Naik and Jain [2018] M. H. Naik and M. Jain, Phys. Rev. Lett. 121, 266401 (2018).
  * Wu et al. [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. 121, 026402 (2018).
  * He et al. [2014] J. He, K. Hummer, and C. Franchini, Phys. Rev. B 89, 075409 (2014).
  * Zhang et al. [2020] Y. Zhang, N. F. Q. Yuan, and L. Fu, Phys. Rev. B 102, 201115 (2020).
  * Wang et al. [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, et al., Nat. Mater. 19, 861 (2020).
  * Xian et al. [2019] L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, Nano Lett. 19, 4934 (2019).
  * Yang et al. [2018b] Q. Yang, M. Wu, and J. Li, J. Phys. Chem. Lett. 9, 7160 (2018b).
  * Park et al. [2019] J. Park, I. W. Yeu, G. Han, C. S. Hwang, and J.-H. Choi, Sci. Rep. 9, 14919 (2019).
  * Wang et al. [2022] X. Wang, K. Yasuda, Y. Zhang, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, L. Fu, and P. Jarillo-Herrero, Nat. Nanotechnol. 17, 367 (2022).
  * Wu and Li [2021] M. Wu and J. Li, Proc. Natl. Acad. Sci. (USA) 118, e2115703118 (2021).
  * Engelke et al. [2023] R. Engelke, H. Yoo, S. Carr, K. Xu, P. Cazeaux, R. Allen, A. M. Valdivia, M. Luskin, E. Kaxiras, M. Kim, et al., Phys. Rev. B 107, 125413 (2023).
  * Cazeaux et al. [2023] P. Cazeaux, D. Clark, R. Engelke, P. Kim, and M. Luskin, J. Elast. 154, 443 (2023).
  * Bennett and Remez [2022] D. Bennett and B. Remez, npj 2D Mater. Appl. 6, 1 (2022).
  * Bennett [2022] D. Bennett, Phys. Rev. B 105, 235445 (2022).
  * Bennett et al. [2023] D. Bennett, G. Chaudhary, R.-J. Slager, E. Bousquet, and P. Ghosez, Nat. Commun. 14, 1629 (2023).
  * Woods et al. [2021] C. R. Woods, P. Ares, H. Nevison-Andrews, M. J. Holwill, R. Fabregas, F. Guinea, A. K. Geim, K. S. Novoselov, N. R. Walet, and L. Fumagalli, Nat. Commun. 12, 347 (2021).
  * Yasuda et al. [2021] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Science 372, 1458 (2021).
  * Vizner Stern et al. [2021] M. Vizner Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. Ben Shalom, Science 372, 1462 (2021).
  * Yoo et al. [2019] H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, et al., Nat. Mater. 18, 448 (2019).
  * Xu et al. [2019] S. G. Xu, A. I. Berdyugin, P. Kumaravadivel, F. Guinea, R. Krishna Kumar, D. A. Bandurin, S. V. Morozov, W. Kuang, B. Tsim, S. Liu, et al., Nat. Commun. 10, 4008 (2019).
  * Moore et al. [2021] S. L. Moore, C. J. Ciccarino, D. Halbertal, L. J. McGilly, N. R. Finney, K. Yao, Y. Shao, G. Ni, A. Sternbach, E. J. Telford, et al., Nat. Commun. 12, 5741 (2021).
  * Weston et al. [2020a] A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, et al., Nat. Nanotech. 15, 592 (2020a).
  * Weston et al. [2020b] A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, et al., Nat. Nanotechnol. 15, 592 (2020b).
  * Van Winkle et al. [2023] M. Van Winkle, S. Craig, Isaac M.and Carr, M. Dandu, K. C. Bustillo, J. Ciston, C. Ophus, T. Taniguchi, K. Watanabe, A. Raja, S. M. Griffin, et al., Nat. Commun. 14, 2989 (2023).
  * Plimpton [1995] S. Plimpton, J. Comp. Phys. 117, 1 (1995).
  * Kınac ı et al. [2012] A. Kınac ı, J. B. Haskins, C. Sevik, and T. Çağ ın, Phys. Rev. B 86, 115410 (2012).
  * Sevik et al. [2011] C. Sevik, A. Kinaci, J. B. Haskins, and T. Çağ ın, Phys. Rev. B 84, 085409 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.085409.
  * Ouyang et al. [2018] W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod, Nano Lett. 18, 6009 (2018).
  * Hoover [1985] W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
  * Nosé [1984] S. Nosé, J. Chem. Phys. 81, 511 (1984).
  * Zhou et al. [2015] S. Zhou, J. Han, S. Dai, J. Sun, and D. J. Srolovitz, Phys. Rev. B 92, 155438 (2015).
  * Bagchi et al. [2020] S. Bagchi, H. T. Johnson, and H. B. Chew, Phys. Rev. B 101, 054109 (2020).
  * Zwanzig [2001] R. Zwanzig, _Nonequilibrium statistical mechanics_ (Oxford university press, 2001).
  * Hirano and Shinjo [1990] M. Hirano and K. Shinjo, Phys. Rev. B 41, 11837 (1990).
  * Dienwiebel et al. [2004] M. Dienwiebel, G. S. Verhoeven, N. Pradeep, J. W. Frenken, J. A. Heimberg, and H. W. Zandbergen, Phys. Rev. Lett. 92, 126101 (2004).
  * Müser [2014] M. H. Müser, in _Fundamentals of Friction and Wear on the Nanoscale_ (Springer, 2014), pp. 209–232.
  * Ru et al. [2020] G. Ru, W. Qi, K. Tang, Y. Wei, and T. Xue, Tribol. Int. 151, 106483 (2020).
  * Kubo [1957] R. Kubo, J. Phys. Soc. Japan 12, 570 (1957).
  * Kirkwood [1946] J. G. Kirkwood, J. Chem. Phys. 14, 180 (1946).
  * Bocquet and Barrat [2013] L. Bocquet and J.-L. Barrat, J. Chem. Phys. 139 (2013).
  * Panoukidou et al. [2021] M. Panoukidou, C. R. Wand, and P. Carbone, Soft Matter 17, 8343 (2021).
  * Shkulipa et al. [2005] S. Shkulipa, W. K. den Otter, and W. J. Briels, Biophys. J. 89, 823 (2005).


# Online Modeling and Control of Soft Multi-fingered Grippers
via Koopman Operator Theory

Lu Shi, Caio Mucchiani, and Konstantinos Karydis The authors are with the
Dept. of Electrical and Computer Engineering, University of California,
Riverside. Email: {lshi024, caiocesr<EMAIL_ADDRESS>We gratefully
acknowledge the support of NSF under grants # IIS-1910087 and # CMMI-2133084.
Any opinions, findings, and conclusions or recommendations expressed in this
material are those of the authors and do not necessarily reflect the views of
the funding agencies.

###### Abstract

Soft grippers are gaining momentum across applications due to their
flexibility and dexterity. However, the infinite-dimensionality and non-
linearity associated with soft robots challenge modeling and closed-loop
control of soft grippers to perform grasping tasks. To solve this problem,
data-driven methods have been proposed. Most data-driven methods rely on
intensive model learning in simulation or offline, and as such it may be hard
to generalize across different settings not explicitly trained upon and in
physical robot testing where online control is required. In this paper, we
propose an online modeling and control algorithm that utilizes Koopman
operator theory to update an estimated model of the underlying dynamics at
each time step in real-time. The learned and continuously updated models are
then embedded into an online Model Predictive Control (MPC) structure and
deployed onto soft multi-fingered robotic grippers. To evaluate the
performance, the prediction accuracy of our approach is first compared against
other model-extraction methods among different datasets. Next, the online
modeling and control algorithm is tested experimentally with a soft 3-fingered
gripper grasping objects of various shapes and weights unknown to the
controller initially. Results indicate high success ratio in grasping
different objects using the proposed method. Sample trials can be viewed at
https://youtu.be/i2hCMX7zSKQ.

##  I Introduction

Grasping and manipulation are fundamental functions of robotics. Traditional
robotic grippers consist of a set of mostly rigid joints and links for a
variety of tasks such as perching [1, 2], digging [3, 4], sorting [5, 6], and
others [7, 8]. Soft grippers offer flexibility and dexterity, and have gained
momentum recently (e.g., [9, 10]). Various applications have been increasingly
using soft grippers, such as minimally invasive surgery [11], agriculture
[12], human rehabilitation [13], as well as manufacturing techniques (3d
printing [14, 15], silicone molding [16, 17], surface and shape deposition
[18]), and others [19, 20, 21].

Soft robots (including soft grippers) have garnered significant attention and
grown into an important research topic, one of the prevailing challenges
relates to their modeling and control [22]. Due to the infinite degrees of
freedom of soft robots, accurate first-principles-based models are hard to
derive. Moreover, in most cases, system non-linearities pose additional
difficulties to control design. One way to resolve such constraints is by
considering data-driven methods.

Various methods based on neural networks have been proposed for modeling and
control of soft grippers. The authors in [23] show the efficiency of recurrent
neural networks in terms of grasp success and stability using a shape-memory
soft gripper with integrated hybrid microelectronic sensing skin. Uniformly
sampled x-y positions and z axis orientations are considered as candidates in
the dataset while no closed-loop control is implemented. In [24], a highly
compliant multi-fingered hand is controlled by inflating or deflating chambers
to demonstrate tasks that human’s hand could do, e.g., turning a valve and
grasping a bottle from a table. The control algorithm is designed based on
imitation learning and reinforcement learning that utilize measurements and
experience to update the policy. Another work [25] tries to estimate objects’
instability and slippage due to acceleration of the robot or insufficient
grasp strength using neural networks to improve grasp performance.

Figure 1: Instance of the soft gripper considered herein grasping a soft ball
using the proposed online learning-based method. The method is designed to be
agnostic to the shape and/or weight of the object to be grasped.

An alternative data-driven direction approach, which is also pursued in this
paper, is based on the Koopman operator theory. The Koopman operator theory
can be used to map a complicated nonlinear system into a higher-dimensional
linear system approximating the unknown nonlinear dynamics directly from data.
Thus, it is appropriate for modeling and control of soft robots [26]. For
example, the highly nonlinear dynamics of a soft robot consisting of three
pneumatically-driven actuators is identified by Koopman operator [27]. The
method constructs the relationship between PWM signals corresponding to
actuator pressure and position/velocity of the end-effector; learned model of
this form can then be used directly for control design. In [28], Hankel
Dynamic Mode Decomposition is used to estimate the Koopman operator. The
learned model is then used in a Linear Quadratic Regulator (LQR) to control a
soft robot arm to reach a target position. The Koopman operator theory is also
implemented in the modeling of a soft robotic swimmer in simulation [29].

In this paper, we propose an online modeling and control approach for soft
multi-fingered robotic grippers based on Koopman operator theory. First,
dictionaries of lifting functions are designed based on the workspace of the
soft gripper. Then, instead of implementing the meta-learning based Koopman
identifier [30], ACD-EDMD [31] is used to learn and update the model at every
time step in real-time without any complicated dictionary learning. Updated
models, in turn, provide the dynamics constraints within an MPC structure to
help achieve desired control objectives.

The paper’s contribution are twofold:

  * •

We propose an online modeling approach that utilizes the Koopman operator
theory. An online control algorithm is designed that uses the learned and
continuously-refined models in real-time and adapts control inputs for an MPC
structure accordingly.

  * •

We evaluate and compare via physical experiments the prediction accuracy of
different modeling methods among different datasets using two soft grippers.
The proposed online modeling and control algorithm is implemented in the
grippers to grasp different objects of varying shape and weight, without any
explicit information provided to the algorithm a-priori.

##  II Preliminary Material

###  II-A Koopman Operator Theory and Extended Dynamic Mode Decomposition
(EDMD)

The Koopman operator theory [32] details how to construct an infinite-
dimensional linear operator that governs the evolution of observables
$g(\xi_{t})$; the latter is also called the “lifted” variables of original
states $\xi_{t}$. Then, the original nonlinear system $f$ can be represented
by Koopman modes, eigenvalues and eigenfunctions. Consider the nonlinear
system $\xi_{t+1}=f(\xi_{t})$, with $\xi\in\mathbb{R}^{n_{x}}$. The evolution
of observables $g$ using Koopman operator
$\mathcal{K}:\mathcal{F}\to\mathcal{F}$ can be written as
$\mathcal{K}g(\xi_{t})=g(f(\xi_{t})).$ Decomposition of the full state
observable $g(\xi)=\xi$ using $N$ Koopman modes $v_{n}$, eigenvalues
$\lambda_{n}$ and eigenfunctions $\varphi_{n}$, yields [33]

$\xi_{t+1}=g(f(\xi_{t}))=\mathcal{K}g(\xi_{t})\to\xi_{t+1}=\sum_{n=1}^{N}v_{n}\lambda_{n}\varphi_{n}(\xi_{t})\enspace.$
(1)

Given state history $\Xi=[\xi_{1},\xi_{2},\dots,\xi_{M},\xi_{M+1}]$ (commonly
referred to as “snapshots”), $\mathcal{K}$ can be expressed as a finite-
dimensional approximation $K:\mathcal{F}_{N}\to\mathcal{F}_{N}$ of the Koopman
operator $\mathcal{K}:\mathcal{F}\to\mathcal{F}$ via EDMD. To do so, we need a
dictionary of functions that lift state variables to a higher-dimensional
space where observable dynamics is approximately linear.

Given a dictionary
$\mathbf{\Psi}(\xi_{m})=[\psi_{1}(\xi_{m}),\ldots,\psi_{N}(\xi_{m})]$, the
Koopman operator can be approximated by minimizing the total residual between
snapshots,
$J=\frac{1}{2}\sum_{m=1}^{M}\left(\mathbf{\Psi}\left(\xi_{m+1}\right)-\mathbf{\Psi}\left(\xi_{m}\right)K\right)^{2}$.
This can be solved by truncated singular value decomposition, yielding

$K\triangleq\boldsymbol{G}^{\dagger}\boldsymbol{A},\enspace\text{where
}\begin{cases}\boldsymbol{G}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{\Psi}_{m}^{*}\mathbf{\Psi}_{m}\enspace,\\\
\boldsymbol{A}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{\Psi}_{m}^{*}\mathbf{\Psi}_{m+1}\enspace,\end{cases}$
(2)

with †, $T$ and ∗ denoting pseudoinverse, transpose and conjugate transpose,
respectively. With $K$ via (2), we obtain

$\begin{cases}v_{n}=(w_{n}^{*}B)^{T}\enspace,\\\
\lambda_{n}\eta_{n}=K\eta_{n}\enspace,\\\
\varphi_{n}=\mathbf{\Psi}_{t}\eta_{n}\enspace,\end{cases}$ (3)

where $\eta_{n}$ is the $n$-th eigenvector, $w_{n}$ is the $n$-th left
eigenvector of $K$ scaled so $w_{n}^{T}\eta_{n}=1$, and $B$ is a weight matrix
such that $\xi=(\mathbf{\Psi}B)^{T}$ [33].

The evolution of the original nonlinear system using the estimated Koopman
operator is described by replacing expressions (3) into (1). Control inputs
can be readily incorporated to the definition of $\mathbf{\Psi}$ as an
augmented state [34].

###  II-B ACD-EDMD

One of the most critical steps that affects the performance of Koopman
operator-based methods is how to choose a proper dictionary for lifting the
original states,
$\mathcal{D}=span\left\\{\psi_{1},\psi_{2},\ldots,\psi_{N}\right\\}$. In our
previous work [31], we have proposed ACD-EDMD as a general and analytical
approach to formalize the construction of the dictionary of lifting functions
based on fundamental dynamical system properties, namely the configuration
space of a rigid robot and the workspace for a soft robot. Hermite polynomials
are the basis functions for Euclidean spaces, whereas for non-Euclidean
spaces, state variables are mapped into a higher Euclidean space. Then, the
dictionary of states $D(\xi)$ is expressed as the Kronecker product of the
lifting functions of all these lower-dimensional spaces. For the control
input, $D\mathbf{u})$ is computed by the zero- and first-order Hermite
polynomials and Kronecker products thereof similarly to the
$\mathbb{R}^{n_{x}}$ cases for states. Finally, the complete dictionary is
formed as the Kronecker product between lifting functions for states and
inputs, i.e. $\mathcal{D}=kron(D(\xi),D(\mathbf{u}))$. With the proper design
of the lifting functions, we can then apply the algorithm introduced before in
Section II-A to extract a model of the robotic system at hand using the
Koopman operator theory.

##  III Soft Gripper System Design and Properties

In this work we consider soft multi-fingered robotic grippers that are
pneumatically actuated. A rigid sub-system is used for gross transport and
positioning of the wrist plate, which contains the soft fingers (actuators)
that are responsible for object grasping (Fig. 1). Note that motion of the
rigid sub-system is not included in the model we aim to learn. Instead, we
only consider the grasping behavior of the soft fingers in this paper.

###  III-A Soft Gripper Design and Integration

The soft gripper comprises three identical bending actuators [35] (Fig. 2)
which are retrofitted PneuNets [36]. These actuators use pneumatic
pressurization and depressurization applied to in-series connected deformable
chambers to create rotational motion. The soft actuators are made via silicone
(Smooth-On Dragon Skin 30) casting in 3D-printed molds (with carbon-fiber
reinforced Onyx filament), and are attached to 3D-printed plates. We consider
two distinct plates for arranging the soft fingers, one where the three
fingers are arranged in a symmetric manner and one that is asymmetric (Fig.
3). The soft gripper assembly is attached to the ReactorX 150, a rigid
manipulator with four links and four joints.111 The original fifth joint of
the robot is a rigid pincher which is removed for the purposes of this work.
The plate acts as a wrist affixed to the manipulators endpoint.

The rigid sub-system serves a dual purpose: first to move the soft gripper to
a desired position and then to lift the soft gripper after it has grasped an
object. Note that the focus of this work is on the soft gripper, assessing how
our proposed Koopman operator-based algorithm can enable the soft gripper to
grasp objects of different shapes and weight without explicit knowledge of. To
evaluate grasping success, we quantify if the gripper can hold the object as
it is being lifted up by the rigid manipulator.

Figure 2: Soft finger depressurized (left), idle (middle) and pressurized
(right).

Figure 3: CAD renderings of the rigid attachment plates considered herein that
afford symmetric (left) and asymmetric (right) soft finger arrangements.

###  III-B Pneumatic Actuation

The soft fingers are pneumatically actuated using Programmable-Air boards
(Fig. 4). The board is based on the Arduino nano (ATMega328P), and has a
maximum pressure of $\pm\;50\;kPa$ and flow rate of $2\;liters/min$.
Pressurization is controlled by changing the PWM signal sent to pumps. Two
such boards were considered for the experiments; the detailed configuration is
discussed in Section V.

Figure 4: The Programmable-Air Pneumatic Control Board.

###  III-C Workspace Construction

As discussed in [31], the workspace of soft robots needs to be defined to
obtain the dictionary of lifting functions, which is used in the Koopman
operator-based modeling approach. The workspace of a robot is the set of
positions that it can reach. For soft multi-fingered grippers, the position of
the tips of the fingers can be considered to construct the workspace. The
state of the soft gripper is thus defined by the concatenation of the
positions of each tip as $\xi_{1}=[x_{1},y_{1},z_{1}]$,
$\xi_{2}=[x_{2},y_{2},z_{2}]$, and $\xi_{3}=[x_{3},y_{3},z_{3}]$, expressed in
terms of an inertial frame set at the center of the plate (Fig. 5).222
Although this work uses a $3$-fingered soft gripper, the proposed approach can
be directly extended to account for general $n$-fingered soft grippers.

We can further simplify the workspace used for operator learning by exploiting
the physical constraints underlying the $\\{x_{i},y_{i},z_{i}\\},i=1,2,3$
state vectors. The types of actuators considered herein are constrained to
create planar motion and can only bend about the normal vector to the plane
they are constrained to create motion at (which is determined herein based on
the arrangement of the rigid attachment plate). This way, only one of
$\\{x_{i},y_{i},z_{i}\\}$ for each finger suffices, thus giving rise to a
workspace for the whole soft gripper of
$\mathbb{R}^{1}\times\mathbb{R}^{1}\times\mathbb{R}^{1}$. Without loss of
generality, here we choose $z_{i}$ as the representations of sub-workspace for
each of the fingers, but any other state component would work as well. Besides
the Kronecker product of sub-dictionaries, we also include in the library
observations of other states and inputs for full-observability. We define
$D_{\xi}=[1,z_{1},z_{2},z_{3},z_{1}z_{2},z_{1}z_{3},z_{2}z_{3},z_{1}z_{2}z_{3}]$.
The dictionary is computed as
$D=[\xi_{1},\xi_{2},\xi_{3},kron(\mathbf{u},D_{\xi})]$. After eliminating
redundant terms, the final dictionary contains $29$ terms.333 It is beneficial
to seek to construct the smallest possible dictionary that best describes the
system at hand for computational expediency. At times, doing so relies on
exploiting engineering intuition and/or any system symmetries. In this work we
take advantage of the bending-only motion afforded by the specific fingers we
consider. However, our method is general and can directly apply to any type of
soft fingers that afford unconstrained position (and orientation) control of
their tips in full 3D space, albeit at the expense of increasing the size of
the constructed dictionary.

##  IV Online Modeling and Control

Unlike rigid robotic systems whereby it is known how to construct their
dynamics, soft robotic systems possess complicated and, most often, infinite-
dimensional dynamics. Thus, data-driven approaches may be a better fit for
extracting models for soft robotic systems directly from data. However, models
learned offline or in simulation (a technique often employed in relevant
works) may be invalidated at runtime due to the presence of uncertainty [37],
unless a very delicately designed experimental environment and operation
process are imposed for the dynamics learned offline to remain valid for
online model-based control. This process can get very tedious and is not
scalable. Instead, online control should require limited data and limited time
to enable online computation, which are key properties enabled via Koopman
operator-based methods [30, 38]. Thus, we propose an online modeling approach
using Koopman operator theory to update dynamics constraints in real time and
adjust inputs accordingly within an MPC structure.444 While we consider MPC in
this work, the continuous refined learned model can be used in any other
online model-based control algorithm.

###  IV-A Online Model Predictive Control (MPC) Formulation

MPC determines control inputs at each time step by solving an optimization
problem among a range of prediction horizon $H_{p}$ constrained by robot
dynamics. The general optimization algorithm can be formulated as

$\min\limits_{u_{i}}\enspace
J((u_{i})_{i=0}^{H_{p}-1},(\xi_{i})_{i=0}^{H_{p}},),$
$\begin{cases}\xi_{i+1}=f(\xi_{i},u_{i}),\enspace i=0,\dots,H_{p}-1,\\\
\xi_{0}=\xi_{initial}.\end{cases}$

In this work, the soft grippers’ dynamic constraints are obtained via ACD-EDMD
in real time. In other words, at each time step, we solve an optimization
problem using an updated learned model. The cost function $J$ is usually user-
designed. Herein, we set a desired trajectory $r$ and minimize the difference
between predicted states and the desired points. With these considerations in
place, the MPC problem solved here at each time step $t$ is expressed as

$\min\limits_{u_{i}}\enspace
J=((u_{i})_{i=0}^{H_{p}-1},(\xi_{i}-r_{i})_{i=0}^{H_{p}},),$ (4)
$\begin{cases}\xi_{i+1}=\text{ACD-EDMD}(\xi_{i},u_{i}),\enspace
i=0,\dots,H_{p}-1,\\\ \xi_{0}=\xi_{t}.\end{cases}$

###  IV-B Online Modeling and Control Algorithm

We are now ready to present the proposed online modeling and control algorithm
(Alg. 1). After initialization, at each time step $t$, control signals are
obtained from the MPC controller with the updated model constraints learned
based on data from the previous time step. Then, the physical robot is driven
by the control inputs and observations are collected to the online training
dataset. A refined model estimated from the updated dataset via ACD-EDMD is
then used to create updated dynamic constraints for the following step. The
process is repeated until the control task (here grasping different objects)
is achieved.

##  V Experimental Testing and Results

To test the efficiency of our online modeling and control algorithm, we
evaluate separately modeling accuracy and control performance. First, we
compare the performance of different model-estimation approaches in online and
offline datasets to show the benefits of training with online observations. We
further compare the ACD-EDMD based online modeling estimation accuracy with
SINDy (Sparse Identification of Nonlinear Dynamical Systems (SINDy) [39]) and
LSTM (Long Short-Term Memory [40]) using online training datasets. Then, we
implement our online control algorithm in the soft gripper and evaluate it in
grasping different objects.

Figure 5: Coordinate systems and motion capture marker configuration.

1 initialize: Evolve the system with random valid inputs for the first $M$
time steps. Then, compute the estimated Koopman operator from the
observations: $K=\text{ACD-EDMD}(\\{\xi\\}_{1}^{M+1},\\{u\\}_{1}^{M})$

2 for  _$t\geq M$_ do

3 repeat

4 MPC Controller: Calculate the control inputs with the model constraints
learned in the previous step by solving optimization (4).

5 Plant: Propagate the soft robotic gripper with $u=[u_{1},u_{2}]$:
$\xi_{t+1}=f(\xi_{t},u)$ , where $f$ is the evolving law of real plant.

6 Model: Collect new data from observations and update the ACD-EDMD model.

7 $t\leftarrow t+1$

8 until  _control task is finished_ ;

Algorithm 1 Online Modeling and Control TABLE I: Prediction accuracy of different approaches _Method_ |  _Dataset_ |  MSE($|\xi_{1}|$) |  MSE($|\xi_{2}|$) |  MSE($|\xi_{3}|$) |  (magnitude)
---|---|---|---|---|---
SINDy |  Offline (Unloaded) |  $[0.00,556.66,0.062]$ |  $[340.17,0.44,188.47]$ |  $[12.84,327.78,0.014]$ |  $\times 10^{3}$
ACD-EDMD |  Offline (Unloaded) |  $[0.26,256.62,5.76]$ |  $[180.83,7.12,133.35]$ |  $[7.91,141.33,0.69]$ |  $\times 10^{2}$
SINDy |  Offline (Loaded) |  $[364.17,38.57,0.74]$ |  $[0.95,0.99,26.59]$ |  $[122.65,0.13,3.25]$ |  $\times 10^{1}$
ACD-EDMD |  Offline (Loaded) |  $[211.90,75.52,0.0026]$ |  $[0.43,0.13,0.47]$ |  $[62.36,15.21,0.05]$ |  $\times 10^{0}$
SINDy |  Online |  $[0.0047,0.0113,0.0006]$ |  $[1.36,1.24,0.09]$ |  $[0.44,0.0022,0.03]$ |  $\times 10^{0}$
LSTM |  Online |  $[0.58,1.40,0.0136]$ |  $[0.38,0.76,0.60]$ |  $[1.45,1.81,0.87]$ |  $\times 10^{-2}$
ACD-EDMD |  Online |  $[0.32,2.16,1.58]$ |  $[0.17,0.93,0.20]$ |  $[0.40,1.35,0.15]$ |  $\times{\bf 10^{-4}}$

Motion capture system (12-camera Prime13 Optitrack) is used to collect
position data ($x$, $y$, $z$) of the tip of each finger. To ensure consistency
among all training data, we set the origin of the inertial frame at the
geometric center of the attachment plate (Fig. 5). In both of the grippers,
one control board controls soft fingers $1$ and $3$, and the other control
board controls finger $2$. Thus, together with the PWM signals of two control
boards, $9$ states and $2$ inputs of the robotic soft gripper are observed.

###  V-A Settings of Different Approaches to Compare Against

####  V-A1 SINDy

SINDy is a data-driven approach to extract governing equations of nonlinear
systems. It builds on the assumption that only a few key nonlinear candidates
governing the evolution of original system. Sparse regression determines
weights of those nonlinear terms from observation data. Here we adopt the
example in [39] as the library $\Theta(X)$ of candidate nonlinear functions,
consisting of constant, polynomial and trigonometric terms, that is

$\Theta(X)=\left[\begin{matrix}1&X&X^{2}&\sin(X)&\cos(X)\end{matrix}\right].$

The sparse regression problem is solved by least absolute shrinkage and
selection operator (LASSO) [41].

####  V-A2 LSTM

LSTM is an artificial recurrent neural network that offers feedback
connections and the ability to deal with entire datasets rather than single
data points. LSTM networks have been used in robotics applications (e.g.,
[42]), notably in-hand robotic manipulation [43]. In this paper, we use the
LSTM layer in MATLAB R2020a Deep Learning Toolbox to construct the network. A
single layer of 20 hidden units is applied as designed in [23].

###  V-B Model Estimation Accuracy

Different datasets for training are collected as follows.

  1. 1.

Unloaded Offline training data: In offline training, constant PWM signals are
given by the control board. Each signal lasts for $16$ sec with a sampling
rate of $2$ Hz ($0.5$ sec). Preliminary testing revealed that the finger can
function effectively when the PWM is within the range of $20\%$ to $35\%$.
Thus, we collect the position and voltage signals data by either pressurizing
or depressurizing all the fingers when the gripper is unloaded within that
range. A total number of $32$ sets are collected as the ‘Unloaded Offline’
training data. A moving-average filter with window size of $5$ is applied to
the offline training data.

  2. 2.

Loaded Offline training data: Data from four trials of grabbing plush oranges
(object $1$ in Fig. 6), two trials of grabbing the smaller rectangular box
(object $2$ in Fig. 6) and two trials of grabbing the larger square box
(object $3$ in Fig. 6) are collected with length of $26$ points for each
trial. The ‘Loaded Offline data’ case is constructed as a combination of these
$8$ data sets.

  3. 3.

Online training data: Online observations are collected when the gripper
attempts to grasp the plush orange. Let $N_{T}$ be the training data length;
every $\xi_{t-N_{T}-1}^{t}$ observations are used to predict $\xi_{t+1}$.

The three aforementioned datasets are used to predict $\xi_{t+1}$ for the
whole online data trial of length $L$. We calculate the Mean Squared Error
(MSE) between predicted (superscript $p$) and true (superscript $t$) states
per
$MSE=\frac{1}{L-N_{T}}\sum_{N_{T}}^{L}([\xi_{1}^{p},\xi_{2}^{p},\xi_{3}^{p}]-[\xi_{1}^{t},\xi_{2}^{t},\xi_{3}^{t}])^{2}$
for the whole trajectory. Results are contained in Table I. For all
algorithms, predictions computed from the offline loaded data are more
accurate than those made based on unloaded data. Further, online approaches
lead to significantly smaller error than all offline approaches.

We also compare the efficiency of all methods with online training data by
evaluating repeated grasping of object $1$ for $8$ times to avoid bias. We
tested all three approaches with different training data size,
$N_{T}=\\{3,5,10\\}$.

Results are depicted in Table II (training time) and Fig. 7 (prediction
accuracy). Our algorithm can run $2\times$ and $3\times$ orders of magnitude
faster than SINDy and LSTM, respectively, per each time step. This performance
is attainable because ACD-EDMD leverages fundamental characteristics of the
dynamical system at hand (herein the workspace of the soft gripper) which
leads to a simultaneously smaller and more descriptive dictionary to use.
Further, with reference to Fig. 7, our method can achieve similar accuracy
with LSTM for smaller training data size, and comparable accuracy for larger
training data size. (More specifically, our method can lead to more variable
accuracy compared to LSTM for $N_{T}=10$, but mean values are close.) The
combination of comparable accuracy but substantially faster training time
makes our proposed approach appropriate for online learning-based modeling and
control tasks.

Figure 6: Objects of different shape and/or weight considered in this work. Sample trials can be viewed at https://youtu.be/i2hCMX7zSKQ. TABLE II: Training times per each time step for different approaches Training Size \ Approach | ACD-EDMD | SINDy | LSTM
---|---|---|---
$N_{T}=3$ |  $0.00028$ s |  $0.0260$ s |  $0.7240$ s
$N_{T}=5$ |  $0.00030$ s |  $0.0583$ s |  $0.7837$ s
$N_{T}=10$ |  $0.00044$ s |  $0.1542$ s |  $0.8468$ s

The bottom row panels of Fig. 7 show that the prediction error obtained by
SINDy is much larger than the other two approaches.555 For error bars to be
visible, we omitted those of SINDy in $N_{T}=\\{3,5\\}$. SINDy’s worse than
expected performance may be because of the small training data size used in
the online procedure, or perhaps the choice of candidate nonlinear functions
selected based on the original publication [39]. This is in fact one key point
of departure from SINDy in that it requires selection of candidate nonlinear
functions and that selection can affect the overall performance.

Figure 7: $1$-$\sigma$ error prediction accuracy of all online methods. SINDy
(blue) has large errors; to keep the bars visible, one case ($N_{T}=10$) is
shown.

###  V-C Online Modeling and Control: Grasping

For the final evaluation of our method, we apply Alg. 1 to grasp objects using
the two soft grippers with different finger configurations. The online
training size $N_{T}$ is $5$ and the prediction horizon $H_{p}$ is set at $3$.
In this experiment, we manually designed desired trajectories for each
finger.666 Note that this experiment is only one example to show how our
online modeling and control approach can be implemented in practice. Other
model-based controllers and desired trajectories can be used as well, and one
can embed our online modeling method into their design directly, which is one
key benefit of our proposed method. The same desired trajectory is used for
grasping objects of different shape and same shape but different weight (Fig.
6), without any information about the objects. Each grasping task is repeated
$10$ times to avoid bias. Successful grasps are considered those that the
gripper grasps the object, holds it until fully lifted and keeps it for $3$
sec after lifting.

Results show that the asymmetric gripper functions better than the symmetric
one in terms of grasping heavier objects, whereas the symmetric gripper
performs better on grasping objects of varying shape (Table III). However, we
can still deduce that our algorithm can be easily implemented into different
robots to grasp different objects without requiring any knowledge about the
robot and the objects.

##  VI Conclusions

Soft multi-fingered grippers have been gaining momentum across applications.
However, because of their infinite-dimensionality and non-linearities,
necessary modeling and control algorithms are difficult to derive. The online
modeling and control method proposed in this paper relies on Koopman operator
theory for data-driven modeling and online MPC for control. Our proposed
method can achieve solid prediction accuracy while training significantly
faster compared to related data-driven methods currently in-use.

Our method can be directly embedded into model-based online control of
robotics systems, without requiring explicit knowledge of system dynamics
except for their configuration space or workspace. Results by testing with two
soft 3-fingered robotic grippers attest that our method can grasp a variety of
objects of different shape and/or weight without any a-priori knowledge. This
work can be beneficial in applications like smart manufacturing and precision
agriculture.

TABLE III: Success rates of grasping different objects  Object | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
Gripper | Symmetric Gripper
Success Rate | $100\%$ | $100\%$ | $90\%$ | $10\%$ | $0\%$ | $0\%$
Gripper | Asymmetric Gripper
Success Rate | $90\%$ | $80\%$ | $100\%$ | $90\%$ | $70\%$ | $0\%$

Future work aims to combine the online modeling approach with other model-
based controllers (e.g., [44]), and develop an approach for generating and
adapting online the desired trajectories to grasp a wider range of objects.

## References

  * [1] S. B. Backus, L. U. Odhner, and A. M. Dollar, “Design of hands for aerial manipulation: Actuator number and routing for grasping and perching,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2014, pp. 34–40.
  * [2] S. Liu, W. Dong, Z. Ma, and X. Sheng, “Adaptive aerial grasping and perching with dual elasticity combined suction cup,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 5, no. 3, pp. 4766–4773, 2020.
  * [3] Z. Tong, Y. H. Ng, C. H. Kim, T. He, and J. Seo, “Dig-grasping via direct quasistatic interaction using asymmetric fingers: An approach to effective bin picking,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 6, no. 2, pp. 3033–3040, 2021.
  * [4] D. Danfeng, M. Yan, G. Xiurong, and L. Huaimin, “Research on a forestation hole digging robot,” in _IEEE International Conference on Intelligent Computation Technology and Automation (ICICTA)_ , vol. 2, 2010, pp. 1073–1076.
  * [5] M. Gupta and G. S. Sukhatme, “Using manipulation primitives for brick sorting in clutter,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2012, pp. 3883–3889.
  * [6] Q. Wu, M. Li, X. Qi, Y. Hu, B. Li, and J. Zhang, “Coordinated control of a dual-arm robot for surgical instrument sorting tasks,” _Robotics and Autonomous Systems_ , vol. 112, pp. 1–12, 2019.
  * [7] C. Mucchiani and M. Yim, “Dynamic grasping for object picking using passive zero-dof end-effectors,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 6, no. 2, pp. 3089–3096, 2021.
  * [8] C. H. Kim and J. Seo, “Shallow-depth insertion: Peg in shallow hole through robotic in-hand manipulation,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 4, no. 2, pp. 383–390, 2019.
  * [9] J. Shintake, V. Cacucciolo, D. Floreano, and H. Shea, “Soft robotic grippers,” _Advanced materials_ , vol. 30, no. 29, p. 1707035, 2018.
  * [10] M. E. Giannaccini, I. Georgilas, I. Horsfield, B. Peiris, A. Lenz, A. G. Pipe, and S. Dogramadzi, “A variable compliance, soft gripper,” _Autonomous Robots_ , vol. 36, no. 1, pp. 93–107, 2014.
  * [11] G. Rateni, M. Cianchetti, G. Ciuti, A. Menciassi, and C. Laschi, “Design and development of a soft robotic gripper for manipulation in minimally invasive surgery: a proof of concept,” _Meccanica_ , vol. 50, no. 11, pp. 2855–2863, 2015.
  * [12] E. Navas, R. Fernández, D. Sepúlveda, M. Armada, and P. Gonzalez-de Santos, “Soft gripper for robotic harvesting in precision agriculture applications,” in _IEEE International Conference on Autonomous Robot Systems and Competitions (ICARSC)_ , 2021, pp. 167–172.
  * [13] H. Zhang, Y. Wang, M. Y. Wang, J. Y. H. Fuh, and A. S. Kumar, “Design and analysis of soft grippers for hand rehabilitation,” in _ASME International Manufacturing Science and Engineering Conference (MSEC)_ , vol. 50756, 2017, p. V004T05A003.
  * [14] C. Tawk, A. Gillett, M. in het Panhuis, G. M. Spinks, and G. Alici, “A 3d-printed omni-purpose soft gripper,” _IEEE Transactions on Robotics (T-RO)_ , vol. 35, no. 5, pp. 1268–1275, 2019.
  * [15] Z. Wang, D. S. Chathuranga, and S. Hirai, “3d printed soft gripper for automatic lunch box packing,” in _IEEE International Conference on Robotics and Biomimetics (ROBIO)_ , 2016, pp. 503–508.
  * [16] S. Liu, F. Wang, Z. Liu, W. Zhang, Y. Tian, and D. Zhang, “A two-finger soft-robotic gripper with enveloping and pinching grasping modes,” _IEEE/ASME Transactions on Mechatronics (TMECH)_ , vol. 26, no. 1, pp. 146–155, 2020.
  * [17] J. Zhou, S. Chen, and Z. Wang, “A soft-robotic gripper with enhanced object adaptation and grasping reliability,” _IEEE Robotics and automation letters (RAL)_ , vol. 2, no. 4, pp. 2287–2293, 2017.
  * [18] S. A. Suresh, D. L. Christensen, E. W. Hawkes, and M. Cutkosky, “Surface and shape deposition manufacturing for the fabrication of a curved surface gripper,” _Journal of Mechanisms and Robotics_ , vol. 7, no. 2, p. 021005, 2015.
  * [19] G. Gao, C.-M. Chang, L. Gerez, and M. Liarokapis, “A pneumatically driven, disposable, soft robotic gripper equipped with multi-stage, retractable, telescopic fingers,” _IEEE Transactions on Medical Robotics and Bionics (T-MRB)_ , vol. 3, no. 3, pp. 573–582, 2021.
  * [20] E. Roels, S. Terryn, J. Brancart, G. Van Assche, and B. Vanderborght, “A multi-material self-healing soft gripper,” in _2nd IEEE International Conference on Soft Robotics (RoboSoft)_ , 2019, pp. 316–321.
  * [21] P. Glick, S. A. Suresh, D. Ruffatto, M. Cutkosky, M. T. Tolley, and A. Parness, “A soft robotic gripper with gecko-inspired adhesive,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 3, no. 2, pp. 903–910, 2018.
  * [22] J. Hughes, U. Culha, F. Giardina, F. Guenther, A. Rosendo, and F. Iida, “Soft manipulators and grippers: a review,” _Frontiers in Robotics and AI_ , vol. 3, p. 69, 2016.
  * [23] J. Zimmer, T. Hellebrekers, T. Asfour, C. Majidi, and O. Kroemer, “Predicting grasp success with a soft sensing skin and shape-memory actuated gripper,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 7120–7127.
  * [24] A. Gupta, C. Eppner, S. Levine, and P. Abbeel, “Learning dexterous manipulation for a soft robotic hand from human demonstrations,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 3786–3793.
  * [25] P. M. Khin, J. H. Low, M. H. Ang Jr, and C. H. Yeow, “Development and grasp stability estimation of sensorized soft robotic hand,” _Frontiers in Robotics and AI_ , vol. 8, 2021.
  * [26] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Data-driven approximations of dynamical systems operators for control,” in _The Koopman Operator in Systems and Control_. Springer, 2020, pp. 197–234.
  * [27] D. Bruder, C. D. Remy, and R. Vasudevan, “Nonlinear system identification of soft robot dynamics using koopman operator theory,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 6244–6250.
  * [28] D. A. Haggerty, M. J. Banks, P. C. Curtis, I. Mezić, and E. W. Hawkes, “Modeling, reduction, and control of a helically actuated inertial soft robotic arm via the koopman operator,” _arXiv preprint arXiv:2011.07939_ , 2020.
  * [29] M. L. Castaño, A. Hess, G. Mamakoukas, T. Gao, T. Murphey, and X. Tan, “Control-oriented modeling of soft robotic swimmer with koopman operators,” in _IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM)_ , 2020, pp. 1679–1685.
  * [30] M. Mazouchi, S. Nageshrao, and H. Modares, “Finite-time koopman identifier: A unified batch-online learning framework for joint learning of koopman structure and parameters,” _arXiv preprint arXiv:2105.05903_ , 2021.
  * [31] L. Shi and K. Karydis, “Acd-edmd: Analytical construction for dictionaries of lifting functions in koopman operator-based nonlinear robotic systems,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 7, no. 2, pp. 906–913, 2021.
  * [32] B. O. Koopman, “Hamiltonian systems and transformation in hilbert space,” _Proceedings of the national academy of sciences of the united states of america (PNAS)_ , vol. 17, no. 5, p. 315, 1931.
  * [33] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” _Journal of Nonlinear Science_ , vol. 25, no. 6, pp. 1307–1346, 2015.
  * [34] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Generalizing koopman theory to allow for inputs and control,” _SIAM Journal on Applied Dynamical Systems_ , vol. 17, no. 1, pp. 909–930, 2018.
  * [35] E. Kokkoni, Z. Liu, and K. Karydis, “Development of a soft robotic wearable device to assist infant reaching,” _Journal of Engineering and Science in Medical Diagnostics and Therapy (JESMDT)_ , vol. 3, no. 2, 2020.
  * [36] R. F. Shepherd, F. Ilievski, W. Choi, S. A. Morin, A. A. Stokes, A. D. Mazzeo, X. Chen, M. Wang, and G. M. Whitesides, “Multigait soft robot,” _Proceedings of the national academy of sciences_ , vol. 108, no. 51, pp. 20 400–20 403, 2011.
  * [37] K. Karydis, I. Poulakakis, J. Sun, and H. G. Tanner, “Probabilistically valid stochastic extensions of deterministic models for systems with uncertainty,” _The International Journal of Robotics Research_ , vol. 34, no. 10, pp. 1278–1295, 2015.
  * [38] L. Shi, H. Teng, X. Kan, and K. Karydis, “A data-driven hierarchical control structure for systems with uncertainty,” in _IEEE Conference on Control Technology and Applications (CCTA)_ , 2020, pp. 57–63.
  * [39] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the national academy of sciences (PNAS)_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [40] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [41] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [42] D. Park, Y. Hoshi, and C. C. Kemp, “A multimodal anomaly detector for robot-assisted feeding using an lstm-based variational autoencoder,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 3, no. 3, pp. 1544–1551, 2018.
  * [43] O. M. Andrychowicz, B. Baker, M. Chociej, R. Jozefowicz, B. McGrew, J. Pachocki, A. Petron, M. Plappert, G. Powell, A. Ray, _et al._ , “Learning dexterous in-hand manipulation,” _The International Journal of Robotics Research (IJRR)_ , vol. 39, no. 1, pp. 3–20, 2020.
  * [44] C. Mucchiani, Z. Liu, I. Sahin, J. Dube, L. Vu, E. Kokkoni, and K. Karydis, “Closed-loop position control of a pediatric soft robotic wearable device for upper extremity assistance,” in _IEEE International Conference on Robot & Human Interactive Communication (RO-MAN)_, 2022.