Further author information: (Send correspondence to P.B.)
P.B.: E-mail<EMAIL_ADDRESS>
S.L.: E-mail<EMAIL_ADDRESS>
# Physical Reservoir Computing with Origami and its Application to Robotic
Crawling

Priyanka Bhovad  Department of Mechanical Engineering, Clemson University,
Clemson, SC, US  Suyi Li  Department of Mechanical Engineering, Clemson
University, Clemson, SC, US

###### Abstract

A new paradigm called physical reservoir computing has recently emerged, where
the nonlinear dynamics of high-dimensional and fixed physical systems are
harnessed as a computational resource to achieve complex tasks. Via extensive
simulations based on a dynamic truss-frame model, this study shows that an
origami structure can perform as a dynamic reservoir with sufficient computing
power to emulate high-order nonlinear systems, generate stable limit cycles,
and modulate outputs according to dynamic inputs. This study also uncovers the
linkages between the origami reservoir’s physical designs and its computing
power, offering a guideline to optimize the computing performance.
Comprehensive parametric studies show that selecting optimal feedback crease
distribution and fine-tuning the underlying origami folding designs are the
most effective approach to improve computing performance. Furthermore, this
study shows how origami’s physical reservoir computing power can apply to soft
robotic control problems by a case study of earthworm-like peristaltic
crawling without traditional controllers. These results can pave the way for
origami-based robots with embodied mechanical intelligence.

###### keywords:

Physical reservoir computing, Origami, Morphological computation, Soft
robotics, Peristaltic locomotion

##  1 INTRODUCTION

The animal kingdom is an endless source of inspiration for soft robotics [1,
2]. Researchers have constructed compliant robots that can mimic all kinds of
animal motions, like octopus locomotion [3], elephant trunk grasping [4],
insect flying [5], jellyfish and fish swimming [6, 7, 8], as well as snake and
insects crawling [9, 10, 11]. These robots share many similarities with
animals regarding their shape and motion kinematics; however, their underlying
sensing, actuation, and control architectures could be fundamentally
different. Our engineered soft robots typically rely on a centralized
controller (aka. an “electronic brain”) that takes up all computing work to
process sensor information, generate control commands, and make decisions.
This approach often struggles to achieve high actuation speed and control
effectiveness as soft robots exhibit virtually infinite degrees of freedom and
complicated dynamic characteristics. On the other hand, animals have highly
interconnected networks of nerves and muscles that can share the workload with
the brain [12, 13, 14]. The animal body’s morphology is an integral part of
its actuation, control, and ultimately its “brain’s” decision-making process,
leading to far superior efficiency than our engineered soft robots.

Motivated by this disparity, an increasing number of researchers have embraced
soft bodies’ nonlinear dynamics as a computational resource to create an
embodied intelligence and control [15, 16, 17, 18, 19, 20, 21]. As a result, a
new computational paradigm called morphological computation has emerged in
which the physical body of the robot itself takes part in performing low-level
control tasks, such as locomotion coordination and modulation, to simplify the
overall control architecture significantly [15, 16, 18, 17, 22]. The
contributions of body morphology to cognition and control involve three major
categories [20]: (1) Morphology facilitating control: wherein the physical
design enables certain behaviors such as motion sequencing (e.g., passive
dynamic walker [23]). (2) Morphology facilitating perception: wherein the
physical design enables sensing (e.g., the nonuniform distribution of cells in
the compound eyes of fly [24]). (3) Morphological computation, such as the
_physical reservoir computing_ (PRC), wherein a physical body performs genuine
computations. Among these, physical reservoir computing shows promising
potentials because of its balanced simplicity and versatility to perform
applicable computation with encoding and decoding [20].

Reservoir computing is a computational framework based on artificial recurrent
neural networks (RNNs), which have been used extensively for problems
involving time-series prediction like the stock market and weather
forecasting, robotic motion planning and control, text and speech recognition
[25, 26, 27, 28, 29, 30, 21, 31]. In RNNs, the output of the current time step
depends on the results from the previous time step in addition to the current
input. Since RNNs involve both forward and back-propagation of input data,
training them became a challenging task. To address this difficulty, Jaeger
introduced the concept of a _fixed_ recurrent neural network as Echo State
Networks (ESNs) [25], and Maass introduced Liquid State Machines (LSMs) [26].
Later, these two concepts merged under the umbrella of reservoir computing
(RC). In RC, the neural network (aka. the “reservoir”) has fixed
interconnections and input weights, and only the linear output readout weights
are trained by simple techniques like linear or ridge regression. These
reservoirs’ dynamics transform the input data stream into a high-dimensional
state space, capturing its non-linearities and time-dependent information for
computation tasks.

More importantly, the reservoir’s fixed nature opens up the possibility of
using physical bodies — such as a random network of nonlinear spring and mass
oscillators [18, 32, 33], tensegrity structures [15, 16, 34, 17], and soft
robotic arms [19, 35, 36] — to conduct computation, hence the paradigm of
Physical Reservoir Computing. These physical systems have shown to possess
sufficient computational power to achieve complex computing tasks like
emulating other non-linear dynamic systems, pattern generation [34, 17, 18,
32, 19, 21], speech recognition [37], and machine learning [36, 21, 31, 33].
More importantly, robotic bodies with sufficient nonlinear dynamics can also
perform like a physical reservoir and directly generate locomotion gait
without using the traditional controllers [17, 21, 38, 39, 40].

In this study, we investigate the use of origami as a physical reservoir and
harness its computing power for robotic locomotion generation. Origami is an
ancient art of folding paper into sophisticated and three-dimensional shapes.
Over the past decades, it has evolved into an engineering framework for
constructing deployable structures [41, 42, 43], advanced materials [44, 45,
46, 47, 48, 49], and robotics [50, 51, 52, 53, 54, 55, 56]. Origami has many
appealing advantages for use in robotics. It is compact, easy to fabricate,
and scale-independent (aka. Origami robots can be fabricated at different
scales but still follow similar folding principles [57, 50, 58, 59]).
Moreover, the nonlinear mechanics and dynamics induced by folding could
enhance robotic performance [60, 61].

We show that origami’s nonlinear folding dynamics also possess significant
computing power. A mechanical system must exhibit several essential properties
to perform as a reservoir [21]. The first one is high-dimensionality, which
allows the reservoir to gather as much information possible from the input
data stream, separating its spatio-temporal dependencies and projecting it
onto a high-dimensional state-space. The second one is non-linearity so that
the reservoir acts as a nonlinear filter to map the information from the input
stream. All the computation complexity is associated with this nonlinear
mapping, thus training the linear static readout becomes a straightforward
task. The third one is fading memory (or short-term memory), ensuring that
only the recent input history influences the current output. The fourth one is
separation property to classify and segregate different response signals
correctly, even with small disturbances or fluctuations. Moreover, if two
input time series differed in the past, the reservoir should produce different
states at subsequent time points [62]. Our physics-informed numerical
simulations prove that origami inherently satisfies these four requirements
and can complete computation tasks like emulation, pattern generation, and
output modulation.

Moreover, we conduct extensive numerical simulations to uncover the linkage
between origami design and its computing power, providing the guideline to
optimize computing performance. Finally, we demonstrate how to directly embed
reservoir computing in an origami robotic body to generate earthworm-like
peristalsis crawling without using any traditional controllers. This study’s
results could foster a new family of origami-based soft robots that operate
with simple mechatronics, interact with the environment through distributed
sensor and actuator networks, and respond to external disturbances by
modulating their activities.

In what follows: Section (2) details the construction of an origami reservoir,
including the lattice framework used to simulate its nonlinear dynamics.
Section 3 elucidates the origami reservoir’s computing power through various
numerical experiments. Section 4 discusses the parametric analysis that
uncovers the linkages between computing performance and physical design.
Section 5 applies the reservoir computing to an origami robot’s crawling
problem. Finally, Section 6 concludes this paper with a summary and
discussion.

##  2 Constructing The Origami Reservoir

In this study, we construct a physical reservoir using the classical Miura-ori
sheets. It is essentially a periodic tessellation of unit cells, each
consisting of four identical quadrilateral _facets_ with _crease_ lengths $a$
$b$ and an internal sector angle $\gamma$ (Figure 1 (a)) [63, 44]. The folded
geometry of Miura-ori can be fully defined with a dihedral _folding angle_
$\theta$ ($\in{[-\pi/2,\pi/2]}$) between the $x$-$y$ reference plane and its
facets. The reservoir size is defined as $n\times m$, where $n$ and $m$ are
the number of origami _nodes_ (aka. vertices where crease lines meet) in $x$
and $y$-directions, respectively. $N$ is the total number of creases in the
origami reservoir.

###  2.1 Dynamics Modeling of the Origami

To investigate this origami reservoir’s computing capacity, one must first
obtain its time responses under dynamic excitation. To this end, we adopt and
expand the lattice framework approach to simulate its nonlinear dynamics [63,
64, 65]. In this approach, origami creases are represented by pin-jointed
stretchable truss elements with prescribed spring coefficient $K_{s}$. Folding
(or bending) along the crease line is simulated by assigning torsional spring
coefficient $K_{b}$ (Figure 1 (b)). We further triangulate the quadrilateral
facets with additional truss elements to estimate the facet bending with
additional torsional stiffness (typically, $K_{b}$ across the facets is larger
than those along the creases). Therefore, this approach discretizes the
continuous origami sheet into a network of pin-jointed truss elements
connected at the nodes. A typical reservoir consists of an interconnected
network of units governed by nonlinear dynamics, and the origami reservoir, in
this case, consists of a network of nodes with their interconnections defined
by the underlying crease pattern. The corresponding governing equations of
motion, in terms of node #p’s displacement ($\mathbf{x}_{p}$) as an example,
are:

$m_{p}\ddot{\mathbf{x}}_{p}^{(j)}=\mathbf{F}_{p}^{(j)},$ (1)

where the superscript “$(j)$” represents the $j^{\text{th}}$ time step in
numerical simulation, and $m_{p}$ is the equivalent nodal mass. Unless noted
otherwise, the mass of the origami sheet is assumed to be equally distributed
to all its nodes. $\mathbf{F}_{p}^{(j)}$ is the summation of internal and
external forces acting on this node in that

$\mathbf{F}_{p}^{(j)}=\sum\mathbf{F}_{s,p}^{(j)}+\sum\mathbf{F}_{b,p}^{(j)}+\mathbf{F}_{d,p}^{j}+\mathbf{F}_{a,p}^{(j)}+m_{p}\textbf{g},$
(2)

where the five terms on the right hand side are the forces from truss
stretching, crease/facet bending, equivalent damping, external actuation, and
gravity, respectively. The formulation of these forces are detailed below.

Figure 1: The nonlinear Truss-frame approach for simulating the origami
dynamics. (a) The crease pattern of the classical Miura-ori, with a unit cell
highlighted. (b) The rigid-folding kinematics of the Miura-ori. (c) The truss-
frame approach discretizes the Miura-ori unit cell, showing the distribution
of truss elements along the creases and across the facets, as well as the
nodal masses. (d) Detailed kinematics and mechanics set up to analyze the
bending and stretching along the truss #$pq$. Notice that $\mathbf{m}^{(j)}$
and $\mathbf{n}^{(j)}$ are the current surface normal vectors defined by
triangles #$pqr$ and #$pqv$, respectively. (e) The bending of the Miura-ori
sheet under its weight. This simulation serves to validate appropriate
material property assignments.

Truss stretching forces: The truss elements are essentially elastic springs
with axial stretching stiffness ($K_{s}^{(j)}=EA/l^{(j)}$). Here, $EA$ is the
material constant, and $l^{(j)}$ is the truss element’s length at the current
$j^{\text{th}}$ time step. Thus, the axial stiffness is updated at each time-
step, accommodating the truss element’s increase in stiffness as it is
compressed and vice-a-versa. The stretching forces from a truss connecting
node #p and one of its neighbouring nodes #$q$ is,

$\mathbf{F}_{s,p}^{(j)}=-K_{s}^{(j)}\left(l_{pq}^{(j)}-l_{pq}^{(0)}\right)\frac{\mathbf{r}_{p}^{(j)}-\mathbf{r}_{q}^{(j)}}{|\mathbf{r}_{p}^{(j)}-\mathbf{r}_{q}^{(j)}|}$
(3)

where $l_{pq}^{(0)}$ is the truss length at its initial resting state.
$\mathbf{r}_{p}^{(j)}$ and $\mathbf{r}_{q}^{(j)}$ are the current position
vectors of these two nodes, respectively. To calculate the total truss
stretching forces acting on node #$p$, similar equations apply to all of its
neighbour nodes through trusses (e.g., node $q$, $r$, $s$, $t$, $u$, and $v$
in Figure 1(c)).

Crease/facet bending forces: The crease folding and facet bending are
simulated with torsional spring coefficient ($K_{b}^{(j)}=k_{b}l^{(j)}$),
where $k_{b}$ is torsional stiffness _per unit length_. Here, we adopt the
formulation developed by Liu and Paulino [64]. For example, the force acting
on nodes #$p$ due to the crease folding along the truss between #$p$ and #$q$
is:

$\mathbf{F}_{b,p}^{(j)}=-K_{b}^{(j)}(\varphi_{pq}^{(j)}-\varphi_{pq}^{(0)})\frac{\partial\varphi_{pq}^{(j)}}{\partial\mathbf{r}_{p}^{(j)}}$
(4)

where $\varphi_{pq}^{(j)}$ is the current dihedral angle along truss $pq$
(aka. the dihedral angle between the triangles #$pqr$ and #$pqv$ in 1(d)), and
$\varphi_{pq}^{(0)}$ is the corresponding initial value. $\varphi_{pq}^{(j)}$
can be calculated as

$\displaystyle\varphi_{pq}^{(j)}$
$\displaystyle=\eta\arccos\left(\frac{\mathbf{m}^{(j)}\cdot\mathbf{n}^{(j)}}{|\mathbf{m}^{(j)}||\mathbf{n}^{(j)}|}\right)\text{
modulo }2\pi$ (5) $\displaystyle\eta$
$\displaystyle=\begin{cases}\text{sign}\left(\mathbf{m}^{(j)}\cdot\mathbf{r}_{pv}^{(j)}\right),&\mathbf{m}^{(j)}\cdot\mathbf{r}_{pv}^{(j)}\neq
0\\\ 1.&\mathbf{m}^{(j)}\cdot\mathbf{r}_{pv}^{(j)}=0\end{cases}$ (6)

Here, $\mathbf{m}^{(j)}$ and $\mathbf{n}^{(j)}$ are current surface normal
vector of the triangles #$pqr$ and #$pqv$, respectively, in that
$\mathbf{m}^{(j)}=\mathbf{r}_{rq}^{(j)}\times\mathbf{r}_{pq}^{(j)}$ and
$\mathbf{n}^{(j)}=\mathbf{r}_{pq}^{(j)}\times\mathbf{r}_{pv}^{(j)}$. In
addition, $\mathbf{r}_{pq}^{(j)}=\mathbf{r}_{p}^{(j)}-\mathbf{r}_{q}^{(j)}$,
$\mathbf{r}_{rq}^{(j)}=\mathbf{r}_{r}^{(j)}-\mathbf{r}_{q}^{(j)}$, and
$\mathbf{r}_{pv}^{(j)}=\mathbf{r}_{p}^{(j)}-\mathbf{r}_{v}^{(j)}$. This
definition of $\varphi_{pq}^{(j)}$ ensures that the folding angle for valley
crease lies in $(0,\pi]$ and the folding angle for mountain crease lies in
$(\pi,2\pi]$. The derivative between folding angle $\varphi_{pq}^{(j)}$ and
the nodal #$p$’s current position vector is

$\frac{\partial\varphi_{pq}^{(j)}}{\partial\mathbf{r}_{p}^{(j)}}=\left(\frac{\mathbf{r}_{pv}^{(j)}\cdot\mathbf{r}_{pq}^{(j)}}{|\mathbf{r}_{pq}^{(j)}|^{2}}-1\right)\frac{\partial\varphi_{pq}^{(j)}}{\partial\mathbf{r}_{v}^{(j)}}-\left(\frac{\mathbf{r}_{rq}^{(j)}\cdot\mathbf{r}_{pq}^{(j)}}{|\mathbf{r}_{pq}^{(j)}|^{2}}\right)\frac{\partial\varphi_{pq}^{(j)}}{\partial\mathbf{r}_{r}^{(j)}}$
(7)

where

$\displaystyle\frac{\partial\varphi_{pq}^{(j)}}{\partial\mathbf{r}_{r}^{(j)}}$
$\displaystyle=\frac{|\mathbf{r}_{pq}^{(j)}|}{|\mathbf{m}^{(j)}|^{2}}\mathbf{m}^{(j)},$
(8)
$\displaystyle\frac{\partial\varphi_{pq}^{(j)}}{\partial\mathbf{r}_{v}^{(j)}}$
$\displaystyle=-\frac{|\mathbf{r}_{pq}^{(j)}|}{|\mathbf{n}^{(j)}|^{2}}\mathbf{n}^{(j)}.$
(9)

Again, to calculate the total crease folding and facet bending forces acting
on node #$q$, similar equations apply to trusses connected to this node (e.g.,
truss $pq$, $pr$, $ps$, $pt$, $pu$, and $pv$ in Figure 1(b)).

Damping forces: Estimating damping ratio and damping force is essential to
achieve realistic dynamic responses and reduce numerical simulation error
accumulation. In this study, we follow the formulation developed in [66, 65].
This formulation first calculates the average velocity of a node with respect
to its neighboring nodes ($\mathbf{v}_{\text{avg}}^{(j)}$) to effectively
remove the rigid body motion components from the relative velocities and
ensure that these components are not damped. Then damping force
$\mathbf{F}_{d,p}^{(j)}$ applied on node #$p$ is given by

$\displaystyle\mathbf{F}_{d,p}^{(j)}$
$\displaystyle=-c_{d}^{(j)}(\mathbf{v}_{p}^{(j)}-\mathbf{v}_{\text{avg}}^{(j)})$
(10) $\displaystyle c_{d}^{(j)}$ $\displaystyle=2\zeta\sqrt{K_{s}^{(j)}m_{p}}$
(11)

where $c_{d}^{(j)}$ is the equivalent damping coefficient, and $\zeta$ is the
damping ratio.

Actuation force: In the origami reservoir, two types of creases receive
actuation. The first type is “input creases,” and they receive input signal
$u(t)$ required for emulation and output modulation tasks. The second type is
“feedback creases,” and they receive reference or current output signal $z(t)$
required by all computing tasks in this study except for the emulation task
(more on the applications of input and feedback creases in Section 2.2). In
the case of multiple outputs, different groups of feedback creases are
present. Here, the selection of input and feedback creases are random. There
are many methods to implement actuation to deliver input $u(t)$ and
reference/feedback signal $z(t)$ to the reservoir. For example, the actuation
can take the form of nodal forces on a mass-spring-damper network [18, 32],
motor generated base rotation on octopus-inspired soft arm [19], or spring
resting length changes in a tensegrity structure [34]. In origami, the
actuation can take the form of moments that can fold or unfold the selected
creases. We assume that the resting angle $\varphi^{(0)}$ of the input and
feedback creases will change — in response to the actuation at every time step
— to a new equilibrium $\varphi_{a,0}^{(j)}$ in that [67, 34]

$\displaystyle\varphi_{a,0}^{(j)}$
$\displaystyle=W_{\text{in}}\tanh(u^{(j)})+\varphi^{(0)}\quad\text{for input
creases;}$ (12) $\displaystyle\varphi_{a,0}^{(j)}$
$\displaystyle=W_{\text{fb}}\tanh(z^{(j)})+\varphi^{(0)}\quad\text{for
feedback creases.}$ (13)

where $W_{\text{in}}$ and $W_{\text{fb}}$ are the input and feedback weight
associated with these actuated creases. They are assigned before the training
and remain unchanged after that. $u^{(j)}$ and $z^{(j)}$ are the input and
feedback signal at the $j^{\text{th}}$ time step. The magnitude of
$W_{\text{in}}$ and $W_{\text{fb}}$ are selected such that
$\varphi_{a,0}^{(j)}\in[0,2\pi]$ and consistent with the folding angle
assignment. This approach of assigning new equilibrium folding angles is
similar to traditional neural network studies that use $\tanh$ as a nonlinear
activation function to transform function $z(t)$ into a new one with
magnitudes between $[-1,1]$. Additionally, it prevents actuator saturation due
to spurious extreme values of $z(t)$. Denote the torsional stiffness of
actuated creases by $K_{b,a}^{(j)}$, and we can update Equation (4) for the
actuated creases (using node #$p$ as an example)

$\displaystyle\mathbf{F}_{a,p}^{(j)}=-K_{b,a}^{(j)}\left(\varphi_{pq}^{(j)}-\varphi_{a,0,pq}^{(j)}\right)\frac{\partial\varphi_{pq}^{(j)}}{\partial\mathbf{r}_{p}^{(j)}},$
(14)

The calculation of other terms in this equation are the same as those in the
force from crease folding and facet bending. Once the governing equations of
motion are formulated, they are solved using MATLAB’s ode45 solver with
$10^{-3}$ second time-steps. Although the governing equation of motions use
nodal displacement $\mathbf{x}^{(j)}$ as the independent variables, we use the
dihedral crease angles $\varphi^{(j)}$ as the _reservoir state_ variables to
characterize the origami’s time responses. This is because measuring crease
angles is easier to implement by embedded sensors, and $\varphi^{(j)}$ can be
directly calculated from $\mathbf{x}^{(j)}$ via the Equations 5 and 6.

###  2.2 Setting Up Reservoir Computing

Similar to the actuated creases (aka. input creases and feedback creases), we
designate “sensor creases” for measuring the reservoir states. We denote
$N_{a}$ as the number of actuated creases, and $N_{s}$ for sensor creases. It
is worth noting that, the actuated creases are typically small subset of all
origami creases (i.e., $N_{a}<N$). The sensor creases, on the other hand, can
be all of the origami creases ($N_{s}=N$) or a small subset as well
($N_{s}<N$).

Once the selections of input, feedback, and sensor creases are completed, one
can proceed to the computing. Physical reservoir computing for tasks that
require feedback (e.g., pattern generations in Section 3.2, and output
modulation in 3.3) consists of two phases: The “training phase” and “closed-
loop phase.” While the emulation tasks require the training phase only
(Section 3.1).

Training phase: In this phase, we use the teacher forcing to obtain the
readout weights $W_{i}$ corresponding to every reservoir state (aka. the
dihedral angles of the sensor creases). Suppose one wants to train the
reservoir to generate a nonlinear time series $z(t)$ (aka. the reference
output). The feedback creases receive the reference output and it dynamically
excites the origami reservoir under an open-loop condition without feedback
(Figure 2(a)). The reservoir states $\varphi^{(j)}$ at every time step are
measured and then compiled into a matrix $\mathbf{\Phi}$.

Once the numerical simulation is over, we segregate the reservoir state matrix
$\mathbf{\Phi}$ into the washout step, training step, and testing step. The
washout step data is discarded to eliminate the initial transient responses.
We then calculate the output readout weights $W_{i}$ using the training step
data via simple linear regression:

$\displaystyle\mathbf{W}_{\text{out}}=[\mathbf{1}\;\mathbf{\Phi}]^{+}\mathbf{Z}=\mathbf{\bar{\Phi}}^{+}\mathbf{Z}$
(15)

where, $[.]^{+}$ refers to the Moore-Penrose pseudo-inverse to accommodate
non-square matrix. $\mathbf{1}$ is a column of ones for calculating the bias
term $W_{\text{out},0}$ to shift the fitted function when necessary.
$\mathbf{Z}$ contains the reference signals at each time step, and it is a
matrix if more than one references present. Lastly, we use testing step data
to verify reservoir performance. It is worth noting that white noise of
amplitude $10^{-3}$ is superimposed on the reservoir state matrix during
training to ensure the robustness of the readout result against numerical
imperfections, external perturbations [32], and instrument noise in “real-
world” applications.

Closed-loop phase: Once the training phase is over and readout weights are
obtained, we run the reservoir in the closed-loop condition. That is, instead
of using the reference output $z(t)$, the current output $z^{*}(t)$ is sent to
the feedback creases (Figure 2(b)), and

$\displaystyle
z^{*}(t)=W_{\text{out},0}+\sum_{i=1}^{N_{s}}W_{\text{out},i}\varphi_{i}(t)=\mathbf{W}_{out}^{T}\bar{\mathbf{\Phi}}$
(16)

where, $N_{s}$ is the number of sensor creases, and
$\bar{\mathbf{\Phi}}=[\mathbf{1}\;\mathbf{\Phi}]$. Thus, the reservoir runs
autonomously in the closed-loop phase without any external interventions.

Figure 2: The setup of physical reservoir computing with origami. (a) The
training phase. The feedback creases receive the reference (or targeted)
output $z(t)$; while white noise is added to the reservoir state vector
$\mathbf{\Phi}(t)$ before calculating output weights
$\mathbf{W}_{\text{out}}$; (b) The closed-loop phase. The output weights
obtained in the training phase are used to calculate the current output, which
is fed back to the feedback creases.

We study the closed loop performance of reservoir by calculating the Mean
Squared Error (MSE) using M time-steps as follows:

$\text{MSE}=\frac{1}{M}\sum_{j=1}^{M}\left(z(j)-z^{*}(j)\right)^{2}$ (17)

To estimate performance when multiple reference outputs are present, we
combine the MSEs by taking a norm over the individual MSEs.

##  3 Computation Tasks By the Origami Reservoir

In this section, we use the origami reservoir to emulate multiple non-linear
filters simultaneously, perform pattern generation, and modulate outputs. The
baseline variables for the origami geometric design, material properties, and
reservoir parameters are given in Table 1.

Figure 3: Emulation tasks with the origami reservoir. (a) The Miura-ori
reservoir used for this task with input creases highlighted. Appropriate
boundary conditions are also necessary. (b) Examples of trajectories generated
in the emulation task including (from top to bottom) input signal $u(t)$, 2nd
order, 10th order system, and Volterra series. Dashed curves are the targeted
trajectories, and solid curves are the result of the reservoir. (c) Error
analysis of the emulation tasks. Circles are the standard deviation of MSE,
and horizontal bars are the corresponding extreme values. Table 1: Design of a
baseline origami reservoir in this study

Reservoir size and material properties
---
Parameter |  Value
Size |  9$\times$9
Nodal Mass |  7 g
$k_{s}$ |  100 N/m
$k_{c}^{a}$ |  1 N/(m-rad)
$k_{c}$ |  0.2525 N/(m-rad)
$K_{f}$ |  10 N/(m-rad)
$\zeta$ |  0.2
Geometric design of Miura-ori
Parameter |  Value
$a$ |  16 mm
$b$ |  10 mm
$\gamma$ |  48$\degree$
$\theta$ |  60$\degree$
Actuator and sensor creases
Parameter |  Value
No. of sensors ($N_{s}$) |  $N$
No. of actuators ($N_{a}$) |  0.45$N$
No. of Feedback creases |  0.3$N$
No. of Input creases |  0.15$N$

Table 2: Emulation task functions

Type | Functions in discretized form (at $j^{th}$ time step)
---|---
Input | $u(j)=0.2\sin(2\pi f_{1}j\Delta t)\sin(2\pi f_{2}j\Delta t)\sin(2\pi f_{3}j\Delta t)$
|  $f_{1}=2.11$ Hz, $f_{2}=3.73$ Hz, $f_{3}=4.33$ Hz
$2^{\text{nd}}$ order system | $z_{1}(j+1)=0.4z_{1}(j)+0.4z_{1}(j)z_{1}(j-1)+0.6(u(j\Delta t))^{3}+0.1$
$10^{\text{th}}$-order system | $\displaystyle z_{2}(j+1)=0.3z_{2}(j-1)+0.05z_{2}(j-1)\sum_{i=1}^{10}z_{2}(j-i)$
| $+1.5u((j-10)\Delta t)u((j-1)\Delta t)+0.1$
Discrete Volterra series | $\displaystyle z_{3}(j+1)=100\sum_{\tau_{1}=0}^{T}\sum_{\tau_{2}=0}^{T}h(\tau_{1},\tau_{2})u(j-\tau_{1})u(nj-\tau_{2})$
| $\displaystyle h(\tau_{1},\tau_{2})=\exp\left(\frac{(\tau_{1}\Delta
t-\mu_{1})^{2}}{2\sigma_{1}^{2}}+\frac{(\tau_{2}\Delta
t-\mu_{2})^{2}}{2\sigma_{2}^{2}}\right)$
| $\mu_{1}=\mu_{2}=0.1,\sigma_{1}=\sigma_{2}=0.05,\Delta t=10^{-3}$

###  3.1 Emulation Task

This sub-section shows that the origami reservoir can emulate multiple
nonlinear filters simultaneously using a single input. Such emulation is a
benchmark task for evaluating the performance in RNN training [68] and prove
the multi-tasking capability of physical reservoirs [18, 19]. Note that the
emulation task involves only the training phase, so there are no feedback
creases in this case. Consequently, we excite the reservoir by sending the
input function $u(t)$ to the input creases and train it to find three sets of
readout weights in parallel via linear regression. Here, $u(t)$ is a product
of three sinusoidal functions with different frequencies, and the three target
non-linear filters are $2^{\text{nd}}$-order non-linear dynamic system
$z_{1}(t)$, a $10^{\text{th}}$-order non-linear dynamic system $z_{2}(t)$, and
discrete Volterra series $z_{3}(n)$ (detailed in Table 2).

We use a $9\times 9$ Miura-ori reservoir in this task, exciting the reservoir
from complete rest and training it for 100 seconds. We discard the first 50
seconds of data as the washout step, use the data from the next 45 seconds to
calculate the optimum static readout weights, and then use the last 5 seconds
of data to calculate the MSE for performance assessments. Results in Figure 3
show that the origami reservoir can emulate these three nonlinear filters. As
the nonlinearity and complexity of the nonlinear filter increases, MSE also
increases (Figure 3(b)).

Moreover, we compare the emulation performance when all $N$ creases are used
as sensor creases versus when only actuated creases are used as sensors
($N_{s}=N_{a}=pN$). The increase in MSE is marginal in the latter case.
Therefore, the origami satisfies the previously mentioned nonlinearity and
fading memory requirements to be a physical reservoir, and one only needs to
use the input creases angles as the reservoir states to simplify the reservoir
setup.

###  3.2 Pattern Generation Task

Pattern generation tasks are essential for achieving periodic activities such
as robotic locomotion gait generation and manipulator control where persistent
memory is required. That is, by embedding these patterns (or limit cycles) in
the origami reservoir, one can generate periodic trajectories in the closed-
loop. We again use a $9\times 9$ Miura-ori reservoir and randomly select
$30\%$ of its creases as the feedback creases (this task does not require
input creases). These feedback creases are divided into two groups for the two
components of 2D trajectories. We run the training phase for 100 seconds for
each pattern, discard the initial 15 seconds of data as the washout step and
use the next 51 seconds’ data to calculate the optimum output readout weights.

Generating non-linear Limit cycles: In the following results, the origami
reservoir demonstrates its computation capability via generating quadratic
limit cycles (LC), Van der Pol limit cycles, and the Lissajous curve in
closed-loop. The quadratic limit cycle is defined by two differential
equations:

$\displaystyle\dot{x}_{1}$
$\displaystyle=x_{1}+x_{2}-\epsilon(t)x_{1}\left(x_{1}^{2}+x_{2}^{2}\right),$
(18) $\displaystyle\dot{x}_{2}$
$\displaystyle=-2x_{1}+x_{2}-x_{2}\left(x_{1}^{2}+x_{2}^{2}\right),$ (19)

where the parameter $\epsilon(t)$ determines the shape of the limit cycle
($\epsilon(t)=1$ in this case). The Van der Pol limit cycle is defined by:

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{2},$ (20)
$\displaystyle\dot{x}_{2}$
$\displaystyle=-x_{1}+\left(1-x_{1}^{2}\right)x_{2}.$ (21)

The Lissajous curve is a graph of two sinusoidal signals parameterized by
their frequency ratio ($f_{1}/f_{2}=0.5$) and phase difference
($\delta=\pi/2$):

$\displaystyle x_{1}$ $\displaystyle=\sin\left(f_{1}t+\delta\right)$ (22)
$\displaystyle x_{2}$ $\displaystyle=\sin\left(f_{2}t\right)$ (23)

Figure 4: Stable pattern generation under closed-loop using the Miura-ori
reservoir. (a) This task’s origami reservoir includes two groups of feedback
creases required to generate 2D limit cycles. (b-d) The closed-loop
trajectories of quadratic limit cycle, Van der Pol oscillator, and the
Lissajous curve, respectively. In these plots, the first row of time responses
shows the closed-loop output after 100s of training. The third row of time
responses shows how the trained reservoir can recover the targeted limit
cycles from an initial resting condition. The corresponding phase portraits
are as shown in the second row. Here, the dashed curves are targeted
trajectories, and the solid curves are the reservoir’s outputs. (e) Van der
Pol limit cycle recovery after the temporary failure of sensor and actuator
creases. The two simulations are the same except for the number of sensor
creases ($N_{s}=N$ for the first test, $N_{s}=0.3N$ for the second). The
insert figures show the corresponding phase-portraits.

As shown in Figure 4(b), the origami reservoir can generate all three periodic
trajectories just by changing the output readout weights. The MSE for
Quadratic LC, Van der Pol LC, and Lissajous curves, calculated using the data
for first 10 seconds’ closed-loop run (M = 10000), are $3.28\times 10^{-7}$,
$2.03\times 10^{-5}$, and $5.5\times 10^{-4}$, respectively. As expected, MSE
increases as the complexity of the curve increases.

Stability and robustness of the pattern generation: After finding the readout
weights, we test the stability of these three limit cycles by starting the
origami reservoir from total rest in the close-loop and running it for more
than 1000 seconds. The limit cycle is stable if and only it can recover the
pattern from zero initial conditions and stays on target for at least 1000
seconds of simulation [32, 19]. The results in Figure 4(c) indicates that the
torsional moments generated from the feedback signals on the feedback creases
are sufficient to recover and maintain the three limit cycles from total rest.
Small phase differences occur between generated trajectories and the targets
because the reservoir takes a slightly different path than the target, and the
Lissajous curve takes more than 15 seconds to recover fully. Nonetheless, the
origami reservoir successfully passes this test.

To further analyze the robustness of reservoir-generated limit cycles, we
simulate actuator and sensor failures. As the origami reservoir generates the
Van der Pol limit cycles in these tests, all feedback and sensor creases stop
working (aka. their signals set to zero) for 10 seconds. We conduct these
tests when all creases are used as sensor creases ($N_{s}=N$) and when only
feedback creases are sensor creases ($N_{s}=N_{a}=0.3N$). The simulation
results in Figure 4(e) show that, although the reservoir diverges to a
trajectory far away from the target during the actuator and sensor failure, it
can immediately recover the Van der Pol limit cycles after the end of these
failures.

Figure 5: Results of modulation task under closed-loop using The Miura-ori
reservoir. (a) This task’s origami reservoir includes two groups of feedback
creases and input creases. (b) Quadratic limit cycle trajectories under
closed-loop and the corresponding input signal $\epsilon(t)$. The results are
obtained after 500 seconds of training. (c) Closed-loop trajectory recovery
from the initial resting conditions. (d) The corresponding phase-portraits,
where the targeted trajectories are overlaid on top of the reservoir output.

###  3.3 Output Modulation Task

Output modulation capability allows the reservoir to adjust its output
according to a randomly varying input signal without changing the readout
weights. This ability is also essential for soft robotic control applications
because it allows the robot to switch behaviors according to external stimuli
or environmental changes. In this task, we randomly select input creases,
which account for $15\%$ of the total creases, in addition to the feedback
creases (Figure 5(a)). Moreover, all creases are used as sensor creases
($N_{s}=N$). The simulation results in Figure 5(b, c) shows the generated
quadratic limit cycles with modulated input (Equation (18, 19)). The origami
reservoir can react to the input and modulate the magnitude of the quadratic
limit cycles. The MSE is $3.8\times 10^{-4}$, which is remarkably small,
considering this task’s complexity.

##  4 Correlating Physical Design and Computing Performance

In this section, we use the mean squared error (MSE) as the metric to examine
the connections between the origami reservoir’s design and computing
performance. In particular, This analysis aims to investigate the sensitivity
of MSE to different parameter changes and identify the optimal origami
designs. To this end, in-depth parametric analyses are conducted to examine
the effect of (1) reservoir size and material properties, (2) crease pattern
geometry, and (3) feedback and sensor crease distribution. We use both Van der
Pol and quadratic limit cycle generation tasks to ensure the broad
applicability of parametric study results.

Table 3: Variables for reservoir size and material properties parametric study

Parameter | Base value | Distribution
---|---|---
Nodal mass (g) | 7 | [1,50]
Geometric | Standard | $\sigma=\chi\exp(\frac{-||(N_{i}-N_{j})||}{l})$
imperfections | Miura-ori |  $\mu=0$, $\chi=0.4a$, $l=4a$
Truss torsional |  $K_{b}^{a}=1$, |  $K_{b}^{a}=1$,
stiffness N/(m-rad) | $K_{b}=0.2525$ | $K_{b}\in[0.005,0.5]$

###  4.1 Reservoir Size, Material Properties, and Vertices Perturbation

We observe that feedback crease distribution affects reservoir computing
performance quite significantly. In particular, poorly distributed feedback
creases might result in failed pattern generating tasks. Therefore, we first
conduct numerical simulations by randomly changing the feedback crease
distributions (72 unique designs in total) and identifying the best performing
one (with the least MSE). We refer to this best performing feedback crease
distribution as the _base design_ (Figure 6(a, c)) for the following
parametric studies. Then, we conduct another parametric study regarding the
nodal mass, crease stiffness, and vertices perturbation. We vary these three
parameters, one at a time, for 72 randomly selected designs (six batches of
jobs in parallel on a computer with 12 cores). The baseline values and range
of the parameters are in Table 3.

The origami reservoir performance turns out to be highly sensitive to the
nodal mass variation. As opposed to the uniform nodal mass in base design, a
randomly distributed nodal mass can significantly increase or decrease the MSE
for both pattern generation tasks. However, randomly distributing mass in an
origami sheet is quite challenging in practical applications. So the use of
varying mass distribution should be judicially done based on the particular
application at hand. On the other hand, the origami performance is much less
sensitive to the crease torsional stiffness. By randomly changing the
stiffness, one can achieve performance at par with the base design.

Moreover, we investigate the effects of random geometric imperfection in the
base designs of origami reservoir. To this end, we adopt the formulation
introduced by Liu et al. [69], which introduce small perturbations to the
nodal positions in folded origami. Such imperfections are inevitable in
practice due to various manufacturing defects. It is found that these small
imperfections do not worsen the MSE significantly and in fact could reduce the
MSE by a moderate degree (Figure 6(a),(b)).

It is also worth noting that the larger $9\times 9$ Miura origami reservoir
performs better than the smaller one because larger origami contains more
folding angles to constitute the reservoir state matrix. Therefore, the high-
dimensionality of a reservoir is desirable to produce smaller MSE.

Figure 6: Effect of reservoir size and material properties on the reservoir
computing performance. (a) The distribution of MSE from the Quadratic limit
cycle simulations using random feedback crease distributions and different
design parameter distributions. Here “FB” stands for feedback crease
distribution, “M” stands for nodal mass distribution, “V” stands for origami
vertices geometry perturbation, and “$K_{f}$” stands for crease torsional
stiffness distribution. It is worth emphasizing that the “FB” results come
from one parametric study of 72 unique designs, and the “M,” “V,” and
“$K_{f}$” are results of the subsequent simulation. The bar charts depict the
average value, standard deviation (circles), and extreme values (horizontal
bars) of MSE. (b) A similar result from the Van der Pol limit cycle generation
task. (c) The feedback crease distributions of the four different baseline
designs used in this parametric study.

###  4.2 Origami Design

Figure 7: Effect of Miura-ori geometric design on the reservoir performance.
(a-c) The Miura-ori geometry and the corresponding landscape of MSE
distribution when $\theta=50\degree$, $60\degree$, and $70\degree$,
respectively. The lighter and darker regions correspond to larger and smaller
errors, respectively, while the white regions depict origami designs that
failed the computing task. (d) The unit cell geometry of four representative
designs with the same crease $a$ length but different sector angles $\gamma$
and crease length ratios $a/b$.

A unique advantage of origami based structures and materials is their
considerable freedom to tailor the geometric design. To this end, we start
from the Base Design I of $9\times 9$ Miura-ori reservoir, vary its crease
length ratio $(a/b)$ and internal sector angle $(\gamma)$, and then run the
quadratic limit cycle task with 100 crease length and sector angle
combinations at three folding angles $(\theta=50\degree,60\degree,70\degree)$.
The results of the parametric analysis are shown in Figure 7. We observe that,
at lower folding angles (flatter origami), the origami reservoir has a higher
possibility to fail the pattern generation tasks. The computing performance
improves significantly with a reduced MSE as the origami folds more (or as
$\theta$ increases). This trend is probably because highly folded origami
offers an increased range of folding motion.

Moreover, there are two design sets with the lowest MES: $a/b\approx 1.5$,
$\gamma\approx 45\degree$, and $a/b\approx 2.5$, $\gamma\approx 60\degree$.
Generally speaking, a moderate to high crease-length ratio and small sector
angles can create “skewed” origami patterns that appear to give better
computing performance across all values folding angles. The best designs here
have MSEs at the order of $10^{-7}$, which is of the same magnitude as we
found previously by tailoring the nodal mass and crease stiffness.

Figure 8: Effect of varying the number of actuator and sensor creases.

###  4.3 Actuator and Sensors Distribution

Finally, it is important, for practical applications, to find the minimum
amount of input/feedback and sensor creases required for achieving acceptable
computing performance. To this end, we start with the $9\times 9$ Miura-ori
reservoir and conduct two tests. In the first test, we vary the percentage of
feedback creases ($N_{a}=0.2N,0.3N,0.4N,0.5N$, each with 24 randomly generated
crease distributions) while using all crease dihedral angles to constitute the
reservoir state matrix (i.e., $N_{s}=N$). In the second test, we use the same
feedback crease design and only use these feedback creases’ dihedral angles to
formulate the reservoir state matrix (i.e., $N_{s}=N_{a}$).

We find that if only $20\%$ of crease are used for feedback, the origami
reservoir might fail the quadratic limit cycle task. On the other hand, the
MSE reduces only marginally as we increase the percentage of feedback creases
beyond $30\%$ (Figure 8). Therefore, we can conclude that using only
$30\%-40\%$ of total creases as the feedback and sensors crease will provide
us an adequate computing performance. This result is significant because it
shows that, even though a large size (high-dimensionality) of the reservoir is
essential for computing performance, one does not need to measure (readout)
every reservoir state. In this way, the practical implementation of the
origami reservoir can be significantly simplified.

In conclusion, the parametric analyses lay out the strategy to optimize the
origami reservoir performance by tailoring the underlying physical and
computational design. A larger origami with a higher-dimension can ensure low
computational error, but one only needs to use $30\%$ $40\%$ of its creases as
the feedback and sensor creases to tap into the origami’s computing capacity.
Meanwhile, the distribution of these feedback and sensor creases must be
carefully chosen with extensive simulations. To further improve computing
performance, one can tailor the origami’s mass distribution, crease stiffness,
and geometric design. Among these options, optimizing the folding geometry
should be the most effective because it is easy to implement in practical
applications.

##  5 Application to soft robotic crawling

This section demonstrates the application of origami reservoir computing to
generate an earthworm-inspired peristaltic crawling gait in a robotic system.
The earthworm uses peristalsis to navigate uneven terrain, burrow through
soil, and move in confined spaces. The lack of complex external appendages
(aka., legs or wings) makes earthworm-inspired robots ideal for field
exploration, disaster relief, or tunnel drilling [70, 71, 72]. The body of an
earthworm consists of segments (metamerism) grouped into multiple “driving
modules” [73, 60]. Each driving module includes contracting, anchoring, and
extending segments actuated by antagonistic muscles (Figure 9(a)). During
peristaltic locomotion, these segments alternately contract, anchor (to the
environment with the help of setae), and extend to create a propagating
peristalsis wave, thus moving the body forward.

We design an earthworm-inspired origami robot consisting of two $3\times 9$
Miura-ori reservoir connected via a stiff central bridge (9(b)). The left and
right half of the robots are symmetric in design, and the central bridge
design allows differential motion between the two halves to facilitate turning
in response to the external input. In each origami reservoir, we embed two
groups of feedback creases (Figure 9(b)) with feedback weights assigned such
that their values for the front and back-half are equal but opposite to each
other. This arrangement reduces the number of reference outputs needed to
generate a crawling gait. To create a peristalsis locomotion gait, we train
the origami reservoirs to generate multiple harmonic signals with a phase
difference of $\pi/2$ among them (aka. a pattern generation task shown Figure
9(b)). We train the robot for 100 seconds and discard the first 15 seconds of
data as the washout step.

Figure 9: Reservoir computing powered crawling origami robot. (a) The
kinematics of a peristaltic locomotion cycle in an earthworm. For clarity, the
earthworm body is simplified and consists of six identical segments organized
into two driving modules. The earthworm body moves forward while the
peristaltic wave of anchoring segments (or driving modules) propagates
backward. (b) The design of an earthworm inspired origami crawling robot that
features two stripes of Miura-ori connected by a zig-zag shaped “ridge.” This
robot has four groups of feedback creases. (c) The closed-loop trajectory
generated by the feedback creases after training. (d) Peristaltic locomotion
cycle in the origami robot as a result of the generated trajectory.

Also, we apply ideal anchors to the bottom origami creases that are in contact
with the surface below. These anchors are assumed to be kinematically attached
to the ground when the associated origami crease folds and relaxed as the
crease unfolds (or flattens). Such anchor design is feasible by leveraging the
origami facets’ folding motion, as shown in the author’s previous study [60].

Figure 9(d) illustrates the robotic locomotion generated by reservoir
computing, while Figure 9(c) depicts the closed-loop response and the limit
cycle recovery from total rest (MSE is $3.9\times 10^{-04})$. As the origami
reservoir generates the multiple harmonic signals with a phase difference, its
folding motion naturally “synchronizes” to these signals, generating a
peristaltic wave of folding and unfolding. As a result, the robot crawls
forward like an earthworm, without using any traditional controllers.

##  6 Summary and Conclusion

We demonstrate the physical reservoir computing capability of origami via
extensive benchmark simulations and parametric studies. First, we develop a
simulation environment to study the nonlinear origami dynamics and detail the
origami reservoir setup. This reservoir successfully achieves many computing
tasks such as emulation, pattern generation, and modulation, all of which are
relevant to robotic applications. We also conduct comprehensive parametric
analysis to uncover the linkage between origami reservoir design and its
computing performance. This new knowledge base offers us a guideline to
optimize computing performance. To the authors’ best knowledge, this is the
first study to rigorously examine the performance of physical reservoir
computer from the lens of the physical design. Finally, we demonstrate how to
embed reservoir computing into an origami robot for control without
traditional controllers through the example of peristaltic crawling.

We list four requirements for a mechanical system to be a reservoir in the
introduction, and origami satisfies all these requirements. The tessellated
origami structures are inherently high-dimensional. For example, a $7\times 7$
Miura-ori with 49 nodes contains $N=60$ crease dihedral angles, all or a small
portion of them can serve as the reservoir states. The nonlinearity of origami
partly originates from the nonlinear kinematic relationships between these
crease angles and external geometry. Also, since origami patterns are highly
structured (ordered), small perturbations in the material properties,
imperfections of crease geometry, and the introduction of local actuation are
sufficient to destroy the regularity and create disorder. These properties
make origami highly nonlinear dynamic reservoirs. The origami reservoir’s
performance in the emulation task proves that it can act as a nonlinear filter
and satisfies fading memory property. Nonlinear patterns can be embedded into
the origami reservoir, and the resulting pattern generation is robust against
external disturbances and recoverable under different initial conditions,
proving separation property. Finally, adding the feedback can create
persistent memory, which is conducive to learning new tasks.

For future robots to work autonomously in unstructured and dynamic
environments, the robot body and brain have to work together by continuously
exchanging information about the current condition, processing this
information, and taking appropriate actions. The physical reservoir computing
embodied robots shown in this study presents a step toward this vision. The
reservoir embedded in the robot body directly gathers information from the
distributed sensor-actuator network to perform low-level control tasks like
locomotion generation. The resulting soft robot can generate the global target
behavior autonomously without controlling every element individually.
Moreover, the generated trajectories could be robust against external
disturbances and modulated according to changing working conditions.

A challenge in implementing physical reservoir computing is the many sensors
and actuators required, even though these sensors and actuators can be simple
individually. Our results contribute in this regard by showing that only a
small portion of origami creases are required to be equipped with sensors and
actuators to tap into the reservoir computing power.

In summary, origami reservoir computing provides an attractive pathway for
facilitating synergistic collaboration between the soft robot’s body and the
brain. The reservoir computing, coupled with unique mechanical properties that
origami can offer — multi-stability [47, 49, 74], nonlinear stiffness [47, 44,
45, 49], and negative Poisson’s ratio [47, 44, 49] — opens up new avenues to
the next generation of soft robots with embedded mechanical intelligence.

## ACKNOWLEDGMENTS

The authors acknowledge the support from the National Science Foundation (CMMI
– 1933124), as well as the Clemson University for the generous allotment of
computing time on Palmetto cluster.

## References

  * [1] Miriyev, A. and Kovač, M., “Skills for physical artificial intelligence,” Nature Machine Intelligence 2(11), 658–660 (2020).
  * [2] Laschi, C., Mazzolai, B., and Cianchetti, M., “Soft robotics: Technologies and systems pushing the boundaries of robot abilities,” Science Robotics 1(1), eaah3690 (2016).
  * [3] Cianchetti, M., Calisti, M., Margheri, L., Kuba, M., and Laschi, C., “Bioinspired locomotion and grasping in water: The soft eight-arm OCTOPUS robot,” Bioinspiration and Biomimetics 10(3), 035003 (2015).
  * [4] Hannan, M. W. and Walker, I. D., “Kinematics and the Implementation of an Elephant’s Trunk Manipulator and Other Continuum Style Robots,” Journal of Robotic Systems 20(2), 45–63 (2003).
  * [5] Ma, K. Y., Chirarattananon, P., Fuller, S. B., and Wood, R. J., “Controlled Flight of a Biologically Inspired, Insect-Scale Robot,” Science 340, 603–607 (may 2013).
  * [6] Joshi, A., Kulkarni, A., and Tadesse, Y., “Fludojelly: experimental study on jellyfish-like soft robot enabled by soft pneumatic composite (spc),” Robotics 8(3), 56 (2019).
  * [7] Ren, Z., Hu, W., Dong, X., and Sitti, M., “Multi-functional soft-bodied jellyfish-like swimming,” Nature Communications 10(1) (2019).
  * [8] Katzschmann, R. K., Marchese, A. D., and Rus, D., “Hydraulic autonomous soft robotic fish for 3D swimming,” Springer Tracts in Advanced Robotics 109, 405–420 (2016).
  * [9] Rafsanjani, A., Zhang, Y., Liu, B., Rubinstein, S. M., and Bertoldi, K., “Kirigami skins make a simple soft actuator crawl,” Science Robotics 3(15), 1–8 (2018).
  * [10] Wu, Y., Yim, J. K., Liang, J., Shao, Z., Qi, M., Zhong, J., Luo, Z., Yan, X., Zhang, M., Wang, X., et al., “Insect-scale fast moving and ultrarobust soft robot,” Science Robotics 4(32), eaax1594 (2019).
  * [11] Seok, S., Onal, C. D., Cho, K. J., Wood, R. J., Rus, D., and Kim, S., “Meshworm: A peristaltic soft robot with antagonistic nickel titanium coil actuators,” IEEE/ASME Transactions on Mechatronics 18(5), 1485–1497 (2013).
  * [12] Hoffmann, M. and Müller, V. C., “Simple or complex bodies? Trade-offs in exploiting body morphology for control,” Studies in Applied Philosophy, Epistemology and Rational Ethics 28(January 2018), 335–345 (2017).
  * [13] Laschi, C. and Mazzolai, B., “Lessons from animals and plants: The symbiosis of morphological computation and soft robotics,” IEEE Robotics and Automation Magazine 23(3), 107–114 (2016).
  * [14] Trivedi, D., Rahn, C. D., Kier, W. M., and Walker, I. D., “Soft robotics: Biological inspiration, state of the art, and future research,” Applied Bionics and Biomechanics 5(3), 99–117 (2008).
  * [15] Paul, C., Investigation of Morphology and Control in Biped Locomotion, PhD thesis, University of Zurich, Switzerland (2004).
  * [16] Paul, C., “Morphological computation. A basis for the analysis of morphology and control requirements,” Robotics and Autonomous Systems 54(8), 619–630 (2006).
  * [17] Caluwaerts, K., D’Haene, M., Verstraeten, D., and Schrauwen, B., “Locomotion without a brain: Physical reservoir computing in Tensegrity structures,” Artificial Life 19(1), 35–66 (2013).
  * [18] Hauser, H., Ijspeert, A. J., Füchslin, R. M., Pfeifer, R., and Maass, W., “Towards a theoretical foundation for morphological computation with compliant bodies,” Biological Cybernetics 105(5-6), 355–370 (2011).
  * [19] Nakajima, K., Hauser, H., Kang, R., Guglielmino, E., Caldwell, D. G., and Pfeifer, R., “A soft body as a reservoir: case studies in a dynamic model of octopus-inspired soft robotic arm,” Frontiers in Computational Neuroscience 7(July), 1–19 (2013).
  * [20] Müller, V. C. and Hoffmann, M., “What Is Morphological Computation? On How the Body Contributes to Cognition and Control,” Artificial Life 23, 1–24 (feb 2017).
  * [21] Tanaka, G., Yamane, T., Héroux, J. B., Nakane, R., Kanazawa, N., Takeda, S., Numata, H., Nakano, D., and Hirose, A., “Recent advances in physical reservoir computing: A review,” Neural Networks 115, 100–123 (2019).
  * [22] Füchslin, R. M., Dzyakanchuk, A., Flumini, D., Hauser, H., Hunt, K. J., Luchsinger, R. H., Reller, B., Scheidegger, S., and Walker, R., “Morphological computation and morphological control: Steps toward a formal theory and applications,” Artificial Life 19(1), 9–34 (2013).
  * [23] Collins, S. H., Wisse, M., and Ruina, A., “A Three-Dimensional Walking Robot with Two Legs and Knees,” The International Journal of Robotics Research 20(7), 607–615 (2001).
  * [24] Floreano, D., Pericet-Camara, R., Viollet, S., Ruffier, F., Brückner, A., Leitel, R., Buss, W., Menouni, M., Expert, F., Juston, R., Dobrzynski, M. K., L’Eplattenier, G., Recktenwald, F., Mallot, H. A., and Franceschini, N., “Miniature curved artificial compound eyes,” Proceedings of the National Academy of Sciences of the United States of America 110(23), 9267–9272 (2013).
  * [25] Jaeger, H., “The “echo state” approach to analysing and training recurrent neural networks-with an erratum note,” Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148(34), 13 (2001).
  * [26] Maass, W. W., Markram, H., and Natschläger, T., “The “ liquid computer”: A novel strategy for real-time computing on time series,” Special Issue on Foundations of Information Processing of TELEMATIK 8(1), 39–43 (2002).
  * [27] Maass, W., Joshi, P., and Sontag, E. D., “Computational aspects of feedback in neural circuits,” PLoS Computational Biology 3(1), 0015–0034 (2007).
  * [28] Maass, W., “Liquid state machines: motivation, theory, and applications,” in [Computability in context: computation and logic in the real world ], 275–296, World Scientific (2011).
  * [29] Lukoševičius, M. and Jaeger, H., “Reservoir computing approaches to recurrent neural network training,” Computer Science Review 3(3), 127–149 (2009).
  * [30] Schrauwen, B., Verstraeten, D., and Van Campenhout, J., “An overview of reservoir computing: theory, applications and implementations,” in [Proceedings of the 15th european symposium on artificial neural networks. p. 471-482 2007 ], 471–482 (2007).
  * [31] Nakajima, K., “Physical reservoir computing-an introductory perspective,” Japanese Journal of Applied Physics 59(6) (2020).
  * [32] Hauser, H., Ijspeert, A. J., Füchslin, R. M., Pfeifer, R., and Maass, W., “The role of feedback in morphological computation with compliant bodies,” Biological Cybernetics 106(10), 595–613 (2012).
  * [33] Morales, G., Huam, S. G., Telles, J., Yamanaka, Y., Nakajima, K., Yaguchi, T., and Hauser, H., “Mass-Spring Damper Array as a Mechanical Medium for Computation,” International Conference on Artificial Neural Networks 1, 208–217 (2018).
  * [34] Caluwaerts, K. and Schrauwen, B., “The body as a reservoir: locomotion and sensing with linear feedback,” in [2nd International conference on Morphological Computation (ICMC 2011) ], (2011).
  * [35] Li, T., Nakajima, K., Cianchetti, M., Laschi, C., and Pfeifer, R., “Behavior switching using reservoir computing for a soft robotic arm,” Proceedings - IEEE International Conference on Robotics and Automation 1, 4918–4924 (2012).
  * [36] Nakajima, K., Hauser, H., Li, T., and Pfeifer, R., “Exploiting the dynamics of soft materials for machine learning,” Soft Robotics 5(3), 339–347 (2018).
  * [37] Fernando, C. and Sojakka, S., “Pattern recognition in a bucket,” in [Advances in Artificial Life ], Banzhaf, W., Ziegler, J., Christaller, T., Dittrich, P., and Kim, J. T., eds., 588–597, Springer Berlin Heidelberg, Berlin, Heidelberg (2003).
  * [38] Degrave, J., Caluwaerts, K., Dambre, J., and Wyffels, F., “Developing an embodied gait on a compliant quadrupedal robot,” IEEE International Conference on Intelligent Robots and Systems 2015-Decem, 4486–4491 (2015).
  * [39] Agogino, A. K., SunSpiral, V., and Atkinson, D., “Super ball bot-structures for planetary landing and exploration,” NASA Technical Report (2018).
  * [40] Urbain, G., Degrave, J., Carette, B., Dambre, J., and Wyffels, F., “Morphological properties of mass-spring networks for optimal locomotion learning,” Frontiers in Neurorobotics 11(MAR), 1–13 (2017).
  * [41] Filipov, E. T., Paulino, G. H., and Tachi, T., “Origami tubes with reconfigurable polygonal cross-sections,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 472(2185), 20150607 (2016).
  * [42] Morgan, J., Magleby, S. P., and Howell, L. L., “An Approach to Designing Origami-Adapted Aerospace Mechanisms,” Journal of Mechanical Design 138(5), 052301 (2016).
  * [43] Dang, X., Feng, F., Plucinsky, P., James, R. D., Duan, H., and Wang, J., “Deployable origami structures that approximate a general surface,” (2020).
  * [44] Schenk, M. and Guest, S. D., “Geometry of Miura-folded metamaterials,” Proceedings of the National Academy of Sciences 110, 3276–3281 (feb 2013).
  * [45] Silverberg, J. L., Evans, A. A., McLeod, L., Hayward, R. C., Hull, T. C., Santangelo, C. D., and Cohen, I., “Using origami design principles to fold reprogrammable mechanical metamaterials,” Science 345, 647–650 (aug 2014).
  * [46] Yasuda, H., Gopalarethinam, B., Kunimine, T., Tachi, T., and Yang, J., “Origami‐Based Cellular Structures with In Situ Transition between Collapsible and Load‐Bearing Configurations,” Advanced Engineering Materials 1900562, 1900562 (2019).
  * [47] Li, S., Fang, H., Sadeghi, S., Bhovad, P., and Wang, K. W., “Architected Origami Materials: How Folding Creates Sophisticated Mechanical Properties,” Advanced Materials 31(5), 1–18 (2019).
  * [48] Yan, Z., Zhang, F., Wang, J., Liu, F., Guo, X., Nan, K., Lin, Q., Gao, M., Xiao, D., Shi, Y., Qiu, Y., Luan, H., Kim, J. H., Wang, Y., Luo, H., Han, M., Huang, Y., Zhang, Y., and Rogers, J. A., “Controlled Mechanical Buckling for Origami-Inspired Construction of 3D Microstructures in Advanced Materials,” Advanced Functional Materials 26, 2629–2639 (apr 2016).
  * [49] Kamrava, S., Mousanezhad, D., Ebrahimi, H., Ghosh, R., and Vaziri, A., “Origami-based cellular metamaterial with auxetic, bistable, and self-locking properties,” Scientific Reports 7, 46046 (2017).
  * [50] Rus, D. and Tolley, M. T., “Design, fabrication and control of origami robots,” Nature Reviews Materials 3, 1 (jun 2018).
  * [51] Miyashita, S., Guitron, S., Li, S., and Rus, D., “Robotic metamorphosis by origami exoskeletons,” Science Robotics 2(10), eaao4369 (2017).
  * [52] Belke, C. H. and Paik, J., “Mori: A Modular Origami Robot,” IEEE/ASME Transactions on Mechatronics 22(5), 2153–2164 (2017).
  * [53] Onal, C. D., Tolley, M. T., Wood, R. J., and Rus, D., “Origami-Inspired Printed Robots,” IEEE/ASME Transactions on Mechatronics 20(5), 2214–2221 (2015).
  * [54] Onal, C. D., Wood, R. J., and Rus, D., “An origami-inspired approach to worm robots,” IEEE/ASME Transactions on Mechatronics 18(2), 430–438 (2013).
  * [55] Yan, R., Luo, M., Wan, Z., Qin, Y., Santoso, J., Skorina, E., and Onal, C., “OriSnake: Design, Fabrication and Experimental Analysis of a 3-D Origami Snake Robot,” IEEE Robotics and Automation Letters 3(3), 1–1 (2018).
  * [56] Novelino, L. S., Ze, Q., Wu, S., Paulino, G. H., and Zhao, R., “Untethered control of functional origami microrobots with distributed actuation,” Proceedings of the National Academy of Sciences 117(39), 24096–24101 (2020).
  * [57] Peraza-Hernandez, E. A., Hartl, D. J., Malak Jr, R. J., and Lagoudas, D. C., “Origami-inspired active structures: a synthesis and review,” Smart Materials and Structures 23(9), 094001 (2014).
  * [58] Ning, X., Wang, X., Zhang, Y., Yu, X., Choi, D., Zheng, N., Kim, D. S., Huang, Y., Zhang, Y., and Rogers, J. A., “Assembly of advanced materials into 3d functional structures by methods inspired by origami and kirigami: a review,” Advanced Materials Interfaces 5(13), 1800284 (2018).
  * [59] Morris, E., McAdams, D. A., and Malak, R., “The state of the art of origami-inspired products: A review,” in [International Design Engineering Technical Conferences and Computers and Information in Engineering Conference ], 50169, V05BT07A014, American Society of Mechanical Engineers (2016).
  * [60] Bhovad, P., Kaufmann, J., and Li, S., “Peristaltic locomotion without digital controllers: Exploiting the origami multi-stability to coordinate robotic motions,” Extreme Mechanics Letters 32, 100552 (2019).
  * [61] Zhakypov, Z., Falahi, M., Shah, M., and Paik, J., “The design and control of the multi-modal locomotion origami robot, tribot,” in [2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) ], 4349–4355 (2015).
  * [62] Legenstein, R. and Maass, W., “What makes a dynamical system computationally powerful,” New directions in statistical signal processing: From systems to brain , 127–154 (2007).
  * [63] Schenk, M. and Guest, S. D., “Origami folding: A structural engineering approach,” in [Proceedings of The Fifth International Meeting of Origami Science Mathematics and Education (5OSME) ], 291–303 (2011).
  * [64] Liu, K. and Paulino, G. H., “Nonlinear mechanics of non-rigid origami: an efficient computational approach,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 473(2206), 20170348 (2017).
  * [65] Ghassaei, A., Demaine, E. D., and Gershenfeld, N., “Fast, interactive origami simulation using gpu computation,” Origami 7, 1151–1166 (2018).
  * [66] Hiller, J. and Lipson, H., “Dynamic Simulation of Soft Multimaterial 3D-Printed Objects,” Soft Robotics 1(1), 88–101 (2014).
  * [67] Paul, C., Valero-Cuevas, F. J., and Lipson, H., “Design and control of tensegrity robots for locomotion,” IEEE Transactions on Robotics 22(5), 944–957 (2006).
  * [68] Atiya, A. F. and Parlos, A. G., “New results on recurrent network training: unifying the algorithms and accelerating convergence,” IEEE Transactions on Neural Networks 11(3), 697–709 (2000).
  * [69] Liu, K., Novelino, L. S., Gardoni, P., and Paulino, G. H., “Big influence of small random imperfections in origami-based metamaterials,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476(2241), 20200236 (2020).
  * [70] Calderón, A. A., Ugalde, J. C., Chang, L., Cristóbal Zagal, J., and Pérez-Arancibia, N. O., “An earthworm-inspired soft robot with perceptive artificial skin,” Bioinspiration and Biomimetics 14, 0–12 (mar 2019).
  * [71] Kamata, M., Yamazaki, S., Tanise, Y., Yamada, Y., and Nakamura, T., “Morphological change in peristaltic crawling motion of a narrow pipe inspection robot inspired by earthworm’s locomotion,” Advanced Robotics 32(7), 386–397 (2018).
  * [72] Fang, H., Zhang, Y., and Wang, K. W., “Origami-based earthworm-like locomotion robots,” Bioinspiration and Biomimetics 12(6), 065003 (2017).
  * [73] Quillin, K. and Quillin, “Kinematic scaling of locomotion by hydrostatic animals: ontogeny of peristaltic crawling by the earthworm lumbricus terrestris,” The Journal of experimental biology 202 (Pt 6), 661–74 (mar 1999).
  * [74] Waitukaitis, S., Menaut, R., Chen, B. G., and van Hecke, M., “Origami multistability: from single vertices to metasheets,” Physical Review Letters 114, 055503 (feb 2015).