Embedding Bifurcations into Pneumatic Artificial Muscle

Abstract Harnessing complex body dynamics has long been a challenge in robotics, particularly when dealing with soft dynamics, which exhibit high complexity in interacting with the environment. Recent studies indicate that these dynamics can be used as a computational resource, exemplified by the McKibben pneumatic artificial muscle, a common soft actuator. This study demonstrates that bifurcations, including periodic and chaotic dynamics, can be embedded into the pneumatic artificial muscle, with the entire bifurcation structure using the framework of physical reservoir computing. These results suggest that dynamics not present in training data can be embedded through bifurcation embedment, implying the capability to incorporate various qualitatively different patterns into pneumatic artificial muscle without the need to design and learn all required patterns explicitly. Thus, this study introduces a novel approach to simplify robotic devices and control training by reducing reliance on external pattern generators and the amount and types of training data needed for control.


INTRODUCTION
Recent studies have revealed that mechanical devices can be designed to use their body dynamics for desired information processing, such as a mechanical random number generator (1) and mechanical neural networks (2).Furthermore, the natural dynamics of mechanical bodies not designed for computation can be used as an information processing resource.The complex dynamics in soft robotic arms, which are inspired by the octopus, can be used for real-time computation, embedding a timer, and controlling the arm by employing the approach of physical reservoir computing (PRC) (3)(4)(5)(6)(7).Reservoir computing (8)(9)(10) is a recurrent neural network framework characterized by use of a highdimensional neural network with nonlinearity and memory.In PRC (11), the network is replaced with physical dynamics.There are various types of robotic bodies, such as mechanical spring-mass-dampers (12), tensegrity (13), quadruped robots (14), and fish robots (15).This suggests that body dynamics can be directly exploited for information processing and control without external memory and nonlinearity.
A pneumatic artificial muscle (PAM) is a typical soft actuator that realizes expansion/contraction or bending dynamics through air pressurization.PAM been studied since the dawn of soft robotics (16,17).The McKibben PAM (18,19) is a central component of soft machines and devices, such as wearable devices (20) and robotic arms (21), and has the advantages of high durability against impact and vibration, a high forceto-weight ratio, and low manufacturing costs.In addition, PAM has been studied as a physical reservoir.PAM length sensors can be emulated by other sensory values in the PAM using the PRC framework (22).The air pressure of a rubber tube connected to a PAM that is attached to an assistive walking device could estimate the posture of the wearer using PRC (23).Moreover, control led to the assistant timing of the walking device based on the estimated information (23).Periodic dynamics have been embedded into a robotic arm composed of PAMs with PRC closed-loop control (24).
The bifurcation structure in dynamical systems is characterized by a qualitative change in dynamics, such as periodic and chaotic, through changes in a parameter.Bifurcation structures appear in the dynamics of robot bodies (25,26).Bifurcation structures in central pattern generators have been shown to contribute to providing the capabilities of exploration and self-organized adaptation for robot control (27)(28)(29).Recently, it has been reported that artificial neural networks with rich information processing capabilities can reconstruct the entire bifurcation structure by learning a subset of dynamics included in the bifurcation structure, which we call bifurcation embedment (30)(31)(32)(33)(34).This study is the first attempt to realize bifurcation embedment into the robotic body.Embedding dynamics implies internalizing the central pattern generator into the body, which is usually externally attached to the robot.In addition, embedding bifurcation structures suggests that it is possible to control various qualitatively different patterns by learning several patterns, without the need to learn all patterns required in robot control.Such a generalization is more powerful than interpolation or extrapolation in traditional machine learning, as the properties of unseen dynamics in the bifurcation structure are qualitatively different from trained ones.Concretely, we demonstrated that periodic dynamics could be embedded into a PAM by training only chaotic dynamics, and vice versa.This study provides insight into reducing the weight and keeping the softness of robotic devices, such as wearable devices, by internalizing external pattern generators using computational capabilities in these morphologies.Furthermore, bifurcation embedment has the potential to significantly reduce the amount and types of training data for robot control.

Pneumatic artificial muscle
This study used the McKibben PAM (Fig. 1A), which consists of a cylindrical rubber tube covered by a braided cord.This PAM is forced by a nearly constant external load.If the PAM is pressurized, it expands in the radial direction and shrinks lengthwise.We used the PAM as a physical reservoir by injecting an input value as a control pressure and sensing its physical quantities.The measurement system is illustrated in Fig. 1B.Inner pressure, length, load, and electric resistance of the rubber were measured.Although traditional natural rubber has low conductivity, and it is difficult to measure its electric resistance value, this study increased the rubber conductivity from 1.0 × 10 −3 S/m to 20 S/m by mixing carbons with the rubber (35).The length of the PAM was 108 mm, with an outer diameter of 11 mm, an inner diameter of 9 mm, and an angle of braid π/6 rad in the equilibrium state (35).

Dynamics
Fig. 1C shows the typical behaviors of each sensor value used in PRC.The input value () represents the following piecewise constant periodic signal: where  and  are the input magnitude and bias, respectively, which tune the input to a suitable input range for the device. is the input interval, and  is the period of input.This study used the following input parameters: A = 0.2 MPa, B, = 0.3 MPa,  = 0.1 sec, and  = 1.2 sec.Sensor responses are presented in Fig. 1C.Sensor value sampling timing was the timing of updating the next input signal.Measured pressure showed nearly the same behavior as control pressure.Furthermore, length and resistance exhibited different curves between the compression and extension phases and possessed the nonlinearity and hysteresis for the input.The load value barely responded to the control pressure and was likely a noise signal.

Bifurcations of electric resistance through external load
The behavior of the electric resistance of the rubber changed drastically through external load and explained the mechanics of these bifurcations.This was accomplished by focusing on changes in the rubber's thickness.The external load added to the PAM changed by 5 N from 100 N to 250 N under the periodic input signal, as depicted in Fig. 1C and represented by Eq. 1.The relationship between length and resistance is presented in Fig. 2A.Resistance response was opposed after 106 mm, which is approximately the equilibrium length of the PAM ( 0 = 108 mm).Thus, these regions can be considered to correspond to contraction and expansion phases.The length gap between 106 mm and 108 mm is considered to occur due to offset measurements.The electric resistance of a rubber tube has the same tendency as the rubber thickness (36).Rubber thickness was derived from measurements of diameter and conservation of volume.The thickness model is provided in the Materials and Methods section.The thickness of the rubber peaks at the equilibrium length is depicted in Fig. 2A.In the expansion phase, the thickness   thinned out from the equilibrium thickness  0 because the inner diameter did not change.In the contraction phase, the thickness   thinned out from  0 due to expansion in the radial direction.Therefore, resistance peaks at the equilibrium length as rubber thickness.
Fig. 2B illustrates the response of length for applied load and pressure in the experiment.Using this method, three load regions were identified: compression phase alone, compression and extension mixing phase, and extension phase alone.The typical resistance time series of each phase is presented in Fig. 2C.Resistance in the compression phase responded to anti-phase pressure value.Resistance in the mixing phase changed to a two-peak behavior.In contrast to the compression phase, resistance in the extension phase responded to the pressure value in-phase.Fig. 2D depicts the bifurcation diagram, where local minimum values of resistance are plotted for the applied load values.Bifurcations were confirmed, revealing that the local minimum value changed from one to two at 160 N and became one again at 220 N.These bifurcation points correspond to change points of the compressing, mixing, and extension phases.

Computing scheme
Fig. 3 illustrates the computational scheme of PAM PRC.The input signal was injected as a control pressure, which is a one-dimensional value.The PAM acted as a physical reservoir by providing a nonlinear historical response to the input.We obtained reservoir variables by sensing these responses and constructing output values from a weighted sum of reservoir variables and bias term.
The input signal is a piecewise constant one-dimensional signal, which is represented in Eq. 1.The nonlinearity and memory of the physical reservoir can be tuned by tuning the input magnitude  and input interval  in Eq.
where the number of samples is represented by .This multiplexing method (37), referred to as time-multiplexing, boosts the computational power of the reservoir from a small number of variables and has been widely used in PRC (6,38,39).In this study,  was fixed at five in all the experiments.The output values  �  = � � ,1 , … ,  � , � ∈ ℝ  were generated as follows: where  out ∈ ℝ (+1) is output weight (+1 in the index means a bias term).The output weight is obtained by ridge regression: where  = �  wash +1 ⋯   wash + train � ∈ ℝ (+1) train is the training data matrix,  = �  wash +1 ⋯   wash + train � ∈ ℝ  train is a target data matrix corresponding, and  is the ridge parameter.The numbers of washout and training data are represented as  wash and  train , respectively.
The open-and closed-loop settings are depicted in Fig. 3. Open-loop represents a case in which the input signal   is external to the reservoir.Closed-loop represents the case in which the input signal   is the output value of the reservoir one time step prior.In the closed-loop setting, input values, which are control pressures, are restricted to a certain range [ min ,  max ] to prevent the breakdown of the system.

Information processing capacity of a single pneumatic artificial muscle
Information-processing capabilities, that is, the nonlinearity and memory of each sensor value in the PAM, were revealed.We used the information processing capacity (IPC) (40) criteria, which describes the function that dynamical system serves for an input signal from independent and identically distributed (i.i.d.) random variables (41).Memory and linear/nonlinear transformation capability of the reservoir can be obtained by checking the bases of this function.IPC limited linear components are called memory capacity (42).Detailed definitions and formulation of the IPC are provided in the Supplementary Materials 1.The IPC restricted the delay to ≤  and degree of polynomial functions to ≤ , as IPC[, ].IPC[, ] can be decomposed by function reconstruction capacities [], where  is an orthogonal basis function in the focusing functional space.The IPC[1, K] is referred to as memory capacity.For a number  of the linear independent states of reservoir variable, the following theoretical equation holds (40): where equality is established if the reservoir has echo state property (ESP).Here, ESP, which is an important property for reservoir computing, guarantees the reproducibility of the computing results (43).
First, we showed the IPCs of each single sensor value and multiplexing sensor values in the PAM, which included pressure, length, resistance, load, and all sensors combined and time-multiplexed.The IPCs are presented in Table 1.The IPC of the pressure was nearly one, and the pressure could be completely described by the input sequence.Therefore, the pressure could be a computational node that has the ESP for input.The IPCs of length and resistance were slightly lower than one, and these sensor values could nearly be described by the input sequence; however, they had few irreproducible components.The IPC of the load was nearly zero, and load moved nearly independently of the input.The IPC could be improved to approximately 10 by combining all types of sensors and using timemultiplexing.The total number of reservoir variables was 20 = 4 × 5.An IPC lower than the number of variables implied the existence of input-independent or linearly dependent components.
Further details, such as nonlinear and memory capacities of the IPCs in the PAM, were examined.Fig. 4 depicts the dependency of the IPCs on the external load.The capacities of the pressure did not change through external load.The capacities of delay zero in the length monotonically increased as external load increased from 50 N to 250 N. Therefore, a PAM with a smaller external load could produce information processing that requires more memory.However, the degree components in the resistance non-monotonically changed through the external load.In addition, the linear components monotonically increased as external load increased to 150 N, degree two components were dominant when the external load was 175 N and 200 N, and the linear components were dominant when the external load was greater than 225 N. Increasing and decreasing switches corresponded to the bifurcation points of resistance; therefore, these critical behaviors were derived from the intrinsic resistance bifurcations.

Open-loop experiments
We evaluated the performance of PAM PRC in an open-loop setting.We investigated PAM length sensor emulation (35), which is a practical task.PAM length sensor emulation is a real-world task that emulates the PAM length sensor value from the input pressure value.A laser displacement sensor, which is a standard length sensor for the PAM, is made of a rigid component that takes away the softness of the PAM.Therefore, it is a practically important method to ensure the softness of the PAM, emulate the length sensor by other sensory values, eliminate it from the PAM.Although the length dynamics of the PAM respond nonlinearly to hysteresis for the input pressure, PAM length time series can be predicted by a recurrent neural network (35,(44)(45)(46)(47)(48).Here, the input sequence in this task arose from uniformly random values and was transformed to [0, 0.5] MPa control pressure value, as (, ) = (0.5, 0) in Eq. 1.In addition, we provided the performance of PAM PRC for a typical benchmark task in the Supplementary Materials 3.
Furthermore, we evaluated the performance of the task using the normalized mean squared error (NMSE), as follows: where  eval is the number of evaluation data.In the following experiments, the number of washout data was  wash = 1,000, the number of training data was  train = 40,000, and the number of evaluation data was  eval = 9,000.We compared the performances between PAM PRC and echo state network (ENS) (9), which is a typical recurrent neural network in reservoir computing.The ENS had the same number of computational nodes as the number of reservoir variables in PAM PRC.The best hyperparameters of the ENS for each task were determined by the grid search.In addition, we compared the performance of the physical model in the length sensor emulation task.The formulations of the ESN and physical model are presented in the Materials and Methods section.
We multiplexed sensor values in the PAM using time-multiplexing  = 5.In one instance, we used 15 reservoir variables (which are time multiplexed pressure, resistance, and load), whereas in another, we used 10 reservoir variables, which were time multiplexed pressure and resistance.The number of nodes in the ESN was the same as in the case of using the load, so this is 15.
The NMSEs of the physical model, ESN, PAM PRC without loads, and PAM PRC with loads were 0.0893, 0.0436, 0.0302, and 0.0294, respectively.PAM PRC outperformed the ESN and physical model.Fig. 4A illustrates the time series of the input, reservoir variable, target, and output signals.The responses of pressure and resistance time series were induced by the random input sequence; however, the load time series was independent of the inputs, as was the result of the IPCs for the load signals in Table 1.Although the load was not induced into the input signal information, cases using load signals had a higher performance than those without the load.In addition, the performance using the load was higher than the upper bound of the performance of the ENS with the sufficient number of nodes, as depicted in Fig. 4B.The length value was not only driven by the input signal but also reflected the noise or intrinsic time-varying dynamics, as the IPC length was 0.975.Only the part of the dynamics that was driven by input can be reconstructed by external machine learning, such as an ESN, which transforms only the inputs.Conversely, PAM PRC with load could treat not only explicit input signals but also implicit inputs and the internal state of the PAM, as the load was independent of the input due to the load IPC (Table 1).
There is another advantage of PAM PRC in the case of few training data.Fig. 4C present performance comparison between an ESN with 600 nodes and PAM PRC through training data.The performance of the ESN significantly decreased when the training data were less than 1,000; however, PAM PRC archives nearly the same performance when the number of training data ranged from 100 to 10,000.This could be considered a problem derived from the dimension of machine-learning networks.An ESN with 600 nodes produces overfitting when the training data are limited; however, PAM PRC can work with as few computational nodes such as 15.This does not easily produce overfitting, even with a small amount of data.Therefore, selecting a suitable reservoir for the task has advantages not only in performance and calculation costs but also with regard to learning.

Closed-loop experiments: Attractor embedding
This study analyzed closed-loop control by PRC in a single PAM.First, we analyzed the potential to embed attractors in a PAM.We focused on limit cycle, strange attractor of the logistic map, and the Rössler attractor.The limit cycle defined by Eq. 2 was a onedimensional periodic dynamic.Such rhythm dynamics are important as a central pattern generator in robot control (49).In addition, not only periodic but also chaotic oscillators play an important role in robot control.As a central pattern generator, a chaotic oscillator can derive adaptive and exploratory behaviors by its complex dynamics (28,29).The logistic map, which is a one-dimensional dynamical system with discrete time, was defined by the following equation: where  is a model parameter and is set as a chaotic parameter  = 3.7.The embedding of logistic dynamics does not require memory, as the next step of the logistic map can be determined only by the current step of the logistic map.Furthermore, the Rössler attractor (50) of three-dimensional chaos with continuous time is defined by the following equation: (10) where , , and  are model parameters and set as typical chaotic parameters: (, , ) = (0.2, 0.2, 5.7).As chaos with continuous time does not occur unless it is a dynamical system with at least three dimensions (51), and the nonlinear term, which is essential for chaos, is only  1  3 , this model can be considered one of the simplest models in chaos with a continuous time.
In the training phase, we injected   into the reservoir using an open-loop and trained the output weight using yn+1 as a teacher signal (this training scheme is referred to as teacher forcing (52)).The input range of all experiments was set to [0.1, 0.5] MPa by tuning  and  in Eq. 1.The input signal of the Rössler system, which is a three-dimensional system, was  1 , which was discretized by a sampling interval of 0.5 sec.Although the input is onedimensional, the attractor can be reconstructed with all dimensions due to Takens' embedding theorem (53) if the reservoir has sufficient memory and nonlinearity.Furthermore, the input interval of the PAM control pressure is  = 0.1 sec in the experiments with the limit cycle and Rössler system, and  = 0.2 sec in the experiments with the logistic map.The embedding of the limit cycle and Rössler system require memory; however, the embedding of the logistic map did not (the relationship between the input interval and IPC of PAM PRC is discussed in the Supplementary Materials 4).In all of the experiments, we fixed the number of washout data at  wash = 1,000 and the number of training data at  train = 4,000.In the prediction phase, we switched the openloop and closed-loop after 1,000 time steps.In addition, we evaluate the embedding result using the output time series, the attractor in the delayed coordinate system, and power spectra.Fig. 6 presents the results of the attractor embedding.In limit cycle embedding, the embedded attractor had similar properties as a time series attractor, and spectra as the target limit cycle.In logistic map embedding, the output time series deviated from the target signal as time passed after switching from an open-to a closed-loop due to the initial state sensitivity of chaos.However, there was an output signal close to the target attractor, and the embedding attractor captured the characteristics of the target attractor.In addition, the output signal had a broad range of spectra that was the same as the target signal, rather than a clear perk of spectra such as the limit cycle.Finally, in the Rössler attractor embedding, the output attractor could reconstruct the highest peak spectra of the Rössler attractor; however, it could not reconstruct multiple peaks, which are characteristic of chaos, that is, of the Rössler attractor.
Next, we evaluated the robustness of the attractor embedding.For this, we injected a random signal from the target signal and confirmed that the output signal could quickly return to the target attractor after switching to a closed loop.We focused on the limit cycle as a target attractor, and the input signals in the open-loop were zero and random signals.
The results are depicted in Figs.6D and 6E.The output signals quickly returned to the target attractor after switching.

Closed-loop experiments: Bifurcation embedding
We confirmed that the IPC of the resistance in the PAM could drastically change through the change in external load in an open-loop setting.Moreover, we found the change in the output signal of PAM PRC through the external load in a closed-loop setting.The following training data were used: A) A limit cycle with a period of 1.2 sec with an external load of 100 N (same as the limit cycle in Fig. 6); B) Limit cycles with periods of 1.2 and 2.4 sec with external loads of 100 N and 250 N, respectively; C) The chaotic trajectory of the logistic map, where  = 3.7, with an external load of 100 N (same as in the case of the logistic map in Fig. 6); D) The period 4 trajectory of the logistic map, where  = 3.55, with an external load of 100 N.
We confirmed the change in the output signal in closed-loop control when the external load changed by 5 N from 100 N to 250 N every at 2,000 time steps.Fig. 7 depicts the results.In experiment A, the amplitude and frequency of the limit cycle continuously changed in the range from load 100 N to 200 N.However, the limit cycle structure of the output signal suddenly collapsed at an external load of 200 N, and the output signal changed to nearly static dynamics.This switching point was around the second bifurcation point of the resistance, as shown in Fig. 2D.Thus, dynamics may switch, as the bifurcation of the resistance propagated to the entire dynamics of the PAM via closed-loop control.Conversely, the result of experiment B indicated that it is possible to suppress the closed-loop bifurcations.In experiment B, we trained limit cycles with different frequencies when external loads were 100 N and 250 N. The results revealed that the frequency of the closed-loop dynamics with intermediate external loads was linearly interpolated.
The results of experiments C and D revealed that periodic and chaotic dynamics could be embedded simultaneously and that one of the dynamics could be generated from learning another one of them.In experiment C, we trained chaotic dynamics in the logistic map when the external load was 100 N. The dynamics switched from chaotic dynamics to period 2 dynamics when the external load was 170 N, which corresponded to the first step of the bifurcation of resistance.As the bifurcation diagram indicates, period 2 dynamics appeared intermittently, acting as a window of the period-doubling bifurcation.In experiment D, we trained period 4 dynamics in the logistic map when the external load was 100 N. The dynamics switched to chaotic dynamics with a one-dimensional attractor in the delay coordinate and broad spectra at the external load 200 N, which corresponded to the second bifurcation of the resistance.The chaotic attractor in the delay coordinate had an alternative shape, similar to a cubic function, to the logistic map.In addition, dynamics had an unstable fixed point near  +1 =   , as there was a hole at the intersection of  +1 =   and the attractor.
These bifurcation embedding results could be useful for robotics applications.For instance, the automatic switching conducted in experiment A could be used for an emergency stop when the external load exceeds the threshold and an idling stop that transitions to a stationary state while the main power is on.This presents the possibility of internalizing adaptive behavioral control that depend on changes in the environment.In addition, the results of experiment B revealed that this switching can be turned off by explicit training on both sides of the bifurcation of the inherent dynamics.Furthermore, the results of experiments C and D demonstrated that multiple qualitatively different dynamics, including chaos, could be switched according to changes in the environment.

DISCUSSION
This study demonstrated that various dynamics and bifurcation structures can be embedded into a soft robotic actuator through systematic analyses of PAM PRC.These results reveal potentials, limitations, and future directions of computing using the robot body.
In the open-loop setting, PAM PRC can outperform the external ESN with a sufficient number of computational nodes by using resistance and load sensor values.This performance is believed to derive from the fact that the resistance and load sensor values reflect the PAM internal state, such as time-variant components and extra inputs, which cannot be represented by the pressure input.The evaluation of PAM IPCs should be extended to multi-inputs and time-variant form to test this hypothesis (41).External machine-learning networks can achieve the same level of performance as PAM PRC if resistance and load sensor values are injected into them.In addition, this study revealed that PAM PRC can obtain high performance of the sensor emulation from small-size training data.This aspect is an important advantage in the information processing of soft materials, as soft materials generally have lower durability than rigid materials, and their material properties can be easily changed over a long time period.
Furthermore, in the closed-loop control, we succeeded in embedding the attractor of the sinusoidal wave and logistic map but failed to embed the Rössler attractor.The dimension of the attractor is believed to have caused this failure.The sinusoidal wave and logistic map have one-dimensional attractor; however, the Rössler attractor that is embedded in a three-dimensional space has a fractal dimension ranging from two to three.The number of inputs in this study is one as the control pressure.To embed high-dimensional dynamics, such as Rössler, it is necessary to utilize the memory of the reservoir due to the Takens delay embedding theorem (53).Therefore, the failure could be because the memory of a single PAM is insufficient to reconstruct the target variables that are not injected as inputs.
In addition, we demonstrated that the PAM can be embedded qualitatively different attractors from the training attractor in the closed-loop experiments.These results suggest that bifurcations in the morphology may have a potential to be exploited to embed the bifurcation structure of the target dynamical system.However, the mechanism to embed the bifurcation structure into the reservoir has not been fully understood to date (33).Moreover, the necessity of the intrinsic bifurcations of the reservoir for the bifurcation embedding remains unknown.
This bifurcation embedding into the body suggests the strong potential of the control of robotics.For instance, if we can embed the period-doubling bifurcation in the morphology of the robot, it may be possible for the robot to generate all the arbitrary periodic dynamics and chaos underlying the Li-Yorke chaos (54) from learning only finite period patterns.
This study indicated that the body dynamics of the single PAM have high computational capability.We believe that these results can be expanded to practical situations.The structures that consist of multiple PAMs, such as a robot arm and wearable assistance suit, may have the potential to embed higher-dimensional and more complex dynamics than single PAMs.For example, if we can embed chaotic itineracy in the robot body, the robot can switch many primitive patterns autonomy and stochastically (55).Moreover, the embedded bifurcation structure could serve as an adaptive pattern switch for the environment, such as an anomaly detection and failure prevention of robot, as the bifurcation points correspond to the change points of the dynamic phase of the body dynamics, such as contraction and extension phases in the PAM.

MATERIALS AND METHODS Pneumatic artificial muscle thickness model
The thickness model of the PAM is presented in Fig. 2A.Thickness was calculated using the following equation:  =  −  (11) where we assume that the rubber tube in the PAM is a uniform cylinder and that  and  are the outer and inner radius of the cylinder, respectively.Furthermore, we assumed the below linear relationship between length and outer radius, as the coefficient of determination between the length and outer radius was 0.9934, which was obtained from the 9 values of length and thickness of the rubber tube with an external load of 50 N.  = −0.3382+ 47.525 (12) where the length of the rubber tube is represented by .We obtained the inner radius r using the following equation due to the constraint of the constant volume of rubber and restriction of both ends of the tube: where  and  are the volume and cross-sectional area of rubber, respectively, and  0 and  0 are the inner radius and length in equilibrium length, respectively.Fig. 2A presents the length and thickness that were calculated using the above equations.

Echo state network
The architectures of the ESN were compared with PAM PRC.
where the activation function is given by , which is the hyperbolic tangent; each node of the input weight  in = � , in � comprises a uniform distribution with [−1, 1], and each node of the internal weight  = � , � comprises a uniform distribution with [−1, 1] and is normalized to make the spectral radius one.The coupling magnitude  cp coincides with the spectral radius of  cp .The bias term   0 is set as   0 = 1.The output weight  out = (  ) is tuned by training.We fix  in = 1 and optimize  cp by grid search for the range [0, 1.2] in each task.

Dynamical model of the PAM
We estimated the length dynamics of PAM from the injected pressure and load for the control.The models of PAM dynamics have been widely investigated in previous studies (56,57).Based on these studies, we used the following length model of the PAM: ̈= − elas () −  fric () −  pre �, ()� +  ex () ( 16) where the displacement of the PAM length is represented as x, and the mass of the PAM is represented as ; the elastic force of the rubber, friction of the rubber, and tension of volume change by pressure are represented as  elas (),  fric (), and  pre �, ()�, respectively; and the input pressure and input load are represented as () and  ex ().Here, tension by pressure is derived from the following Schulze equation: where the strain of the PAM is represented as  = ( 0 − )/ 0 and the equilibrium length, inner radius, and angle of the braided code are represented as  0 ,  0 , and  0 , respectively.Note that the Schulze equation assumes that the PAM is a uniform cylinder with zero thickness.When the real PAM is compressed, it is not a cylinder but a bent shape, as both ends of the PAM are fixed.The model that considers the non-uniform and bent shape of the PAM has been previously proposed in extant literature (57).We ensured the linear elasticity  elas () ∝  is the elasticity of the PAM.We could accurately estimate the length of the PAM by solving the equation of the equilibrium of  pre (),  elas (), and  ex in the static state.However, in the dynamic state, in which the PAM continues to move, it is difficult to estimate the PAM dynamic, as the Schulze equation cannot consider the hysteresis depicted in Fig. 2 in the main text.The causes of the hysteresis can be considered as the effects of Coulomb and viscous frictions (58)(59)(60).Therefore, Eq. 25 can be rewritten using the following equation: ̈= − − ̇− sgn() +  �− pre �, ()� +  ex ()� ( 18) Here, , , , and  are the parameters of the model.We optimized these parameters using grid search, and the parameters used in the Section 2.6 were (, , , ) = (6353, 80.05, 10, 0.635).The measured length values were offset due to a measurement error; thus, we added a bias to the length-predicting value from the physical model to coincide with the average values of the measured and predicted lengths.

Fig. 1 .
Fig. 1.Pneumatic artificial muscles, measurement systems, and pneumatic artificial muscle dynamics.(A) No pressurized and pressurized pneumatic artificial muscles.(B) Pneumatic artificial muscle measurement systems.(C) Sensor responses for a sinusoidal wave input.

Fig. 2 .
Fig. 2. Intrinsic bifurcations of pneumatic artificial muscle.(A) The top graph plots length versus resistance; the bottom figures depict the schematic illustration of the thickness change.(B) Color map of the length for control load and pressure.(C) Resistance and pressure time series in four load conditions.(D) Bifurcation diagram of control load versus local minimum resistance.

Fig. 3 .
Fig. 3. Schematics of reservoir computing.(A) The right-hand side depicts a schematics of typical reservoir computing.The left-hand side depicts a schematics of Physical reservoir computing for pneumatic artificial muscle.

Fig. 4 .
Fig. 4. Information processing capacities of pneumatic artificial muscle sensory values.The bars indicate the decompositions of the IPCs through the degree components and memory components (The method for the decomposition is provided in the Supplementary Materials 1).

Fig. 5 .
Fig. 5. Open-loop length sensor emulation.(A) Time series of the input, reservoir variables, output, and prediction in the length sensor emulation task.The red line is the prediction signal of physical reservoir computing for pneumatic artificial muscle.(B) Normalized mean squared error through the number of echo state network nodes.(C) Normalized mean squared error using the number of training data.Echo state network nodes: n = 600.

Fig. 6 .
Fig. 6. Results of the closed-loop at the attractor embedding.(A) Time series of the target and physical reservoir computing output signals.(B) Attractors of the target and physical reservoir computing output.(C) Spectra of the target and physical reservoir computing output.(D) Time series and attractor in the limit cycle, embedding from zero inputs.(E) Time series and attractor in the limit cycle, embedding from a random input.
1.The dependency of nonlinearity and PAM physical reservoir memory on input magnitude  and input interval  is discussed in the Supplementary Materials 2. Sensor values at time  are represented in the form () = � 1 (), … ,   ()�, where M is the number of sensor values used as reservoir variables.We obtained reservoir variable   , which corresponds to input   , based on sensing  times from input injected time  to input updated time  + .Eq. 2 presents   ∈ ℝ +1 : The results of experiments C and D did not indicate the desired bifurcation structure but rather a bifurcation structure based on training data and reservoir dynamics.We present the embedding a bifurcation structure, which includes desired qualitative different signals, by training dynamics on both sides of bifurcation explicitly in the Supplementary Materials 4.
The th computational node at time  is represented as    , the th input node is represented as   , and the th output node at time  is represented as  �   .The computational nodes and outputs of the ESN are given by