Physical reservoir computing with origami and its application to robotic crawling

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.


INTRODUCTION
The animal kingdom is an endless source of inspiration for soft robotics. 1,2 0][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.3][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.
7][18]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][27][28][29][30][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). 26Later, 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][17]34 and soft robotic arms 19,35,36 -to conduct computation, hence the paradigm of Physical Reservoir Computing. Thse physical systems have shown to possess sufficient computational power to achieve complex computing tasks like emulating other non-linear dynamic systems, pattern generation, 17-19, 21, 32, 34 speech recognition, 37 and machine learning.21,31,33,36 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 isan ancient art of folding paper into sophisticated and threedimensional 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, ad scale-independent (aka.Origami robots can be fabricated at different scales but still follow similar folding principles 50,[57][58][59] ). Moreover, thenlinear 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 nonlinearity 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 Ou physics-informed numerical simulations prove that origami inherently satisfies these four requirements and can complete computation tasks like emulation, pattern generation, and output modulation. Moreoer, 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 tudy'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.

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 γ (Figure 1 (a)). 44,63 he folded geometry of Miura-ori can be fully defined with a dihedral folding angle θ (∈ [−π/2, π/2]) between the x-y reference plane and its facets.The reservoir size is defined as n × 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.

Dynamics Modeling of the Origami
To investigate this origami reservoir's computing capacity, one must first obtain its time responses under dynamic excitation.4][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 (x p ) as an example, are: where the superscript "(j)" represents the j 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.
p is the summation of internal and external forces acting on this node in that 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.
Truss stretching forces: The truss elements are essentially elastic springs with axial stretching stiffness (K Here, EA is the material constant, and l (j) is the truss element's length at the current j 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, where l (0) pq is the truss length at its initial resting state.r (j) p and r (j) q 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 where k b is torsional stiffness per unit length.Here, we adopt the formulation developed by Liu and Paulino. 64For example, the force acting on nodes #p due to the crease folding along the truss between #p and #q is: where ϕ pq is the current dihedral angle along truss pq (aka.the dihedral angle between the triangles #pqr and #pqv in 1(d)), and ϕ pq is the corresponding initial value.ϕ pq can be calculated as Here, m (j) and n (j) are current surface normal vector of the triangles #pqr and #pqv, respectively, in that pq and n (j) = r pv .In addition, r q , and r pq ensures that the folding angle for valley crease lies in (0, π] and the folding angle for mountain crease lies in (π, 2π].The derivative between folding angle ϕ (j) pq and the nodal #p's current position vector is where 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. 65,66 his formulation first calculates the average velocity of a node with respect to its neighboring nodes (v avg ) to effectively remove the rigid body motion components from the relative velocities and ensure that these components are not damped.Then damping force F (j) d,p applied on node #p is given by where c (j) d is the equivalent damping coefficient, and ζ 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. 34In origami, the actuation can take the form of moments that can fold or unfold the selected creases.We assume that the resting angle ϕ (0) of the input and feedback creases will change -in response to the actuation at every time step -to a new equilibrium ϕ a,0 in that 34,67 ϕ (j) a,0 = W fb tanh(z (j) ) + ϕ (0) for feedback creases.(13)   where W in and W 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 th time step.The magnitude of W in and W fb are selected such that ϕ a,0 ∈ [0, 2π] 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 (j) b,a , and we can update Equation (4) for the actuated creases (using node #p as an example) 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 x (j) as the independent variables, we use the dihedral crease angles ϕ (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 ϕ (j) can be directly calculated from x (j) via the Equations 5 and 6.

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. 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 ϕ (j) at every time step are measured and then compiled into a matrix Φ.
Once the numerical simulation is over, we segregate the reservoir state matrix Φ 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: where, [.] + refers to the Moore-Penrose pseudo-inverse to accommodate non-square matrix.We study the closed loop performance of reservoir by calculating the Mean Squared Error (MSE) using M time-steps as follows: To estimate performance when multiple reference outputs are present, we combine the MSEs by taking a norm over the individual MSEs.

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.

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 ote 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 nd -order non-linear dynamic system z 1 (t), a 10 th -order non-linear dynamic system z 2 (t), and discrete Volterra series z 3 (n) (detailed in Table 2).
We use a 9 × 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.

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 × 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 compu- tation 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: where the parameter ǫ(t) determines the shape of the limit cycle (ǫ(t) = 1 in this case).The Van der Pol limit 2 nd order system z 1 (j + 1) = 0.4z 1 (j) + 0.4z 1 (j)z 1 (j − 1) + 0.6(u(j∆t)) 3 + 0.1 10 th -order system z 2 (j + 1) = 0.3z 2 (j − 1) + 0.05z 2 (j − 1) +1.5u((j − 10)∆t)u((j − 1)∆t) + 0.1 Discrete Volterra series z 3 (j + 1) = 100 cycle is defined by: The Lissajous curve is a graph of two sinusoidal signals parameterized by their frequency ratio (f 1 /f 2 = 0.5) and phase difference (δ = π/2): x 2 = sin (f 2 t) As shown in Figure 4 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. 19,32 he 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.

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 × 10 −4 , which is remarkably small, considering this task's complexity.

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, indepth 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.

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 Fixed in z Phase Portrait of (b) Phase Portrait of (c) 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.
Truss torsional 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 × 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.

Origami Design
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 × 9 Miura-ori reservoir, vary its crease length ratio (a/b) and internal sector angle (γ), and then run the quadratic limit cycle task with 100 crease length and sector angle combinations at three folding angles (θ = 50 • , 60 • , 70 • ).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 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 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.

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 × 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., 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.

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.1][72] The body of an earthworm consists of segments (metamerism) grouped into multiple "driving modules". 60,73 ach 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 × 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 π/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.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. 60gure 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 × 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.

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 × 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, 44,45,47,49 and negative Poisson's ratio 44,47,49 -opens up new avenues to the next generation of soft robots with embedded mechanical intelligence.

Figure 1 :
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 m (j) and 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.
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).

Figure 2 :
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 Φ(t) before calculating output weights W 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.

10 - 4 NSFigure 3 :
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), 2 nd order, 10 th 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.
(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 × 10 −7 , 2.03 × 10 −5 , and 5.5 × 10 −4 , respectively.As expected, MSE increases as the complexity of the curve increases.

Figure 4 :
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.

Figure 5 :
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 ǫ(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.

Figure 6 :
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,

(
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.as the origami folds more (or as θ 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 ≈ 1.5, γ ≈ 45 • , and a/b ≈ 2.5, γ ≈ 60 • .

Figure 7 :
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 θ = 50 • , 60 • , and 70 • , 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 γ and crease length ratios a/b.

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

Figure 9 :
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.
1 is a column of ones for calculating the bias term W out,0 to shift the fitted function when necessary.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 32 the reservoir state matrix during training to ensure the robustness of the readout result against numerical imperfections, external perturbations,32and instrument noise in "real-world" applications.......

Table 1 :
Design of a baseline origami reservoir in this study

Table 2 :
Emulation task functions

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