Design and modeling of a novel soft parallel robot driven by endoskeleton pneumatic artificial muscles

Abstract

Owing to their inherent great flexibility, good compliance, excellent adaptability, and safe interactivity, soft robots have shown great application potential. The advantages of light weight, high efficiency, non-polluting characteristic, and environmental adaptability provide pneumatic soft robots an important position in the field of soft robots. In this paper, a soft robot with 10 soft modules, comprising three uniformly distributed endoskeleton pneumatic artificial muscles, was developed. The robot can achieve flexible motion in 3D space. A novel kinematic modeling method for variable-curvature soft robots based on the minimum energy method was investigated, which can accurately and efficiently analyze forward and inverse kinematics. Experiments show that the robot can be controlled to move to the desired position based on the proposed model. The prototype and modeling method can provide a new perspective for soft robot design, modeling, and control.

Key words： pneumatic artificial muscles; soft robot; modeling approach; principle of virtual work; external load

1 Introduction

Soft robots provide new ideas for solutions to problems regarding safety and flexibility, which are difficult to overcome by traditional robots in human–robot interaction and complex special environments [1]. Soft robots can be mainly divided into rope [2], pneumatic [3], hydraulic [4], magnetic field [5], smart material [6], chemical reaction [7], and bio-hybrid drives [8]. Owing to their light weight, high efficiency, non-polluting characteristic, good flexibility, and environmental adaptability, pneumatic soft robots have been widely studied by researchers.

Pneumatic actuators mainly include fiber-constrained pneumatic soft actuators [9,10], elastic gas chamber actuators [11,12], corrugated structure soft actuators [13,14], and origami-based pneumatic soft actuators [15,16]. Although the structural design forms of pneumatic soft actuators are diverse and different, the basic working principle is similar, that is, gas is used as the working medium and the elastic chamber produces directional expansion or contraction in a spatial dimension (e.g., axial, bending, and torsion) under the action of air (positive pressure or negative pressure) and structural constraints. The energy of the pneumatic soft actuator originates from external air pressure. Therefore, the volume-change capability of the soft actuator has a remarkable influence. The advantage of a high deformation rate when expanding, contracting, and folding has directed increased attention from researchers toward origami structures, which have been gradually applied to the design of soft actuators [17].

Kim et al. [18] proposed a dual-mode deformation origami actuator. By combining two deformation modes, origami folding deformation and an elastic air chamber, the actuator can achieve a very high deformation rate. Martinez et al. [19] designed a pneumatic soft origami actuator through different origami structures; when the actuator is inflated, the origami structure could achieve a variety of movements, such as elongation, bending, and oscillation. Yang et al. [20] proposed a new linear pneumatic soft actuator based on the principle that when the elastic chamber is a vacuum, the specially designed outer wall structure contracts under the pressure difference and generates a contraction force. Jiao et al. [21] designed a soft actuator driven by a vacuum inspired by an origami structure. Under the action of negative pressure, the actuator flexes in the form of folding motion similar to origami, realizing the compound motion of expansion, torsion, bending, and torsion. The above study required the use of the deformation of the chamber structure, but the deformation was uncertain and required a special design to achieve the intended motion. The introduction of an “endoskeleton” enabled a more deterministic form of actuator motion and improved the flexural stiffness of the actuator in other directions. Therefore, researchers have investigated pneumatic actuators with endoskeletons. Li et al. [22] proposed a soft actuator inspired by origami. The actuator only comprises a compressible origami skeleton and soft skin. The actuator can be programmed with a skeleton structure to achieve multiaxis motion, including contraction, bending, and twisting. Deshpande et al. [23] developed a new soft bidirectional pneumatic actuator based on the Yoshimura origami structure using a double-layer chamber. Under vacuum, the actuator axially contracts and can output a contraction force. Lee and Rodrigue [24] proposed a novel linear origami pneumatic muscle, using a metal frame as an “endoskeleton.” The chamber is internally supported by a uniformly spaced transverse metal frame to limit its lateral contraction, and the actuator linearly contracts under vacuum. Currently, the “endoskeleton” pneumatic actuator can realize axial motion by changing the length, and bending motion is mainly realized by asymmetric constraint structures, such as asymmetric wall thickness and asymmetric paper folding structure. However, this asymmetric structure can only achieve specific direction and angle bending, and cannot achieve a wide range of free, flexible bending motions. In this study, a soft robot based on an endoskeleton pneumatic artificial muscle (PAM) was developed, comprising 10 soft modules, 11 constraint disks, and a backbone. Each module is coupled with three PAMs in parallel, and when one of the PAMs is driven alone, the soft module bends on one side. The coupled motion of the three PAMs enables 3D flexible bending motion.

Theoretical modeling is necessary to perform motion control of soft robots. Unlike traditional rigid rod robots, soft robots have a theoretically infinite number of degrees of freedom, which poses a substantial challenge for theoretical modeling [25]. The kinematic analysis of soft robots usually discretizes the robot into a series of interest points, describes the spatial curves defined by these points using mathematical methods, and then fits the robot to the curves [26]. Modeling methods for soft robots can be divided into constant curvature and variable curvature approaches. The constant curvature approach discretizes the robot into a series of mutually tangential circular arcs, which is a simplified method. Owing to the advantages of simple models and easy solutions, the constant curvature approach has been widely applied in the field of soft robot modeling [27]. Sofla et al. [28] proposed a model based on the constant curvature approach and the concentrated mass method, which can describe the asymmetric hysteresis nonlinearity of pneumatic muscles. Tan et al. [29] derived a mechanical model of a soft actuator with three degrees of freedom based on a constant curvature approach and experimentally validated the workspace. Garriga-Casanovas et al. [30] combined the constant curvature approach and spiral curves to derive a closed-form solution for a concentric tube soft robot with spiral precurvature. Schiller et al. [31] extended the segmental constant curvature approach using the minimum energy method, and the proposed method can effectively estimate the forward kinematics of soft robots. Yang et al. [32] used the constant curvature approach and the principle of virtual work to establish a forward and inverse kinematic model for a rope-driven soft robot and analyzed the configuration properties in different cases. The constant curvature approach has some application limitations and does not exactly match the characteristics of soft robots. Robots with more sections and longer lengths have larger calculation errors. Therefore, various variable curvature methods have been developed. Zeng et al. [33] developed a model of a soft robot with two degrees of freedom based on the Euler–Bernoulli equation, which can predict the bending deformation of a robot in 2D workspace. Singh et al. [34] developed an inverse kinematic model of a soft robot with variable curvature based on Pythagorean hodograph curves, which can calculate deformation in 3D space. Godage et al. [35] developed a kinematic model of a multisegment soft robot based on mode shape functions to approximate complex, nonlinear parameter variations via simple mathematical functions. Yang et al. [36] developed a kinematic model of a cable-driven multisegment continuum robot based on the Lagrangian formulation that analyzed the elongation, bending, and torsional deformation of the robot. These modeling methods mainly considered the geometry of the soft robot and did not consider the external forces that the robot might be subjected to.

Some mechanics-based modeling approaches have considered actuation and external forces. Yuan et al. [37] established a static model of a soft robot based on Cosserat rod theory, which considers external forces and can describe the configuration of the robot in two dimensions. Huang et al. [38] established a kinematic model of a soft arm with two degrees of freedom based on the absolute nodal coordinate method. End loads and gravity were considered, and robot configuration was analyzed when the end points were moving in the plane. Bieze et al. [39] developed a deformation–space formulation of soft robot dynamics based on the finite element method, which considered bending, torsion, shear, and axial deformation caused by external loads. Sadati et al. [40] proposed a new discrete mechanic model based on Euler–Bernoulli beam segment and relative states, and verified it on a continuum robot prototype. Godage et al. [41] established a mechanic model of a pneumatic soft continuum robot based on the improved lumped parameter method and proved the superiority of the proposed method. The current kinematic analysis of robots under external forces mainly focuses on 2D space. The finite element method, which can be applied to 3D space, is usually computationally intensive. Therefore, this paper proposes a modeling method based on the minimum energy method, which can efficiently and accurately calculate the forward and inverse kinematics of a soft robot with external loads.

The main contributions of this paper are as follows. (1) A soft robot with 10 soft modules is developed based on endoskeleton PAMs. Each module comprises three axes of symmetrically distributed muscles. The robot can flexibly bend in different directions by controlling muscle stretching. (2) A modeling method based on the minimum energy method is investigated. First, the variational form of the total energy of the robot system is established, and the equilibrium equation is derived. Then, the equation is discretized by the finite difference method and solved via the nonlinear least-squares method and the Levenberg–Marquardt algorithm. Finally, the pose parameters of the point of interest are obtained by interpolation to control robot motion.

The remainder of this paper is organized as follows: The developed soft robot is introduced in Section 2. The forward and inverse kinematic model is investigated in Section 3. Experimental validation is conducted in Section 4, and the conclusion is provided in Section 5.

2 Design of pneumatic soft continuum robot

The endoskeleton PAM consists of an airbag and an endoskeleton encapsulated by an airbag. The airbag shortens when deflated and extends when inflated. Three axially symmetric PAMs form a parallel soft module. When driving a single muscle, the module bends on one side, and the module can be bent in different directions owing to the synergy of the three muscles. Ten soft modules are connected in series to form a soft robot, and its structure is shown in Fig.1. Two adjacent soft modules have a 60° relative turning angle, as shown in Fig.2. The fastening coordinate system of the tth constrained disk is N_t-xyz, and the position coordinates {\mathit{a}}_{t,i} of the ith muscle endpoint in the tth soft module are defined as:

Fig.1 Developed soft robot.

Full size|PPT slide

Fig.2 Relative deflection of adjacent modules.

Full size|PPT slide

(1){\mathit{a}}_{t,i}=\{\begin{array}{c}{\left[\begin{array}{ccc}{r}_{\mathrm{c}}\cos \left(i\cdot {\displaystyle \frac{2\text{π}}{3}}-{\displaystyle \frac{\text{π}}{6}}\right)& r_{\mathrm{c}}\sin \left(i\cdot {\displaystyle \frac{2\text{π}}{3}}-{\displaystyle \frac{\text{π}}{6}}\right)& 0\end{array}\right]}^{\mathrm{T}}\text{if}mod(t,2)=0,\\ {\left[\begin{array}{ccc}{r}_{\mathrm{c}}\cos \left(i\cdot {\displaystyle \frac{2\text{π}}{3}}+{\displaystyle \frac{\text{π}}{6}}\right)& r_{\mathrm{c}}\sin \left(i\cdot {\displaystyle \frac{2\text{π}}{3}}+{\displaystyle \frac{\text{π}}{6}}\right)& 0\end{array}\right]}^{\mathrm{T}}\text{if}mod(t,2)=1,\end{array}

where mod(t,2 ) denotes the remainder of t divided by 2, and r_{\mathrm{c}} denotes the constraint disk radius.

The structure of the endoskeleton has an important effect on the contraction length of the PAM, as shown in Fig.3. The number of creases is the main factor that affects the performance of the PAMs because the cross-section of the endoskeleton is a square. For few numbers of creases, when the PAMs contract, stress concentration occurs at the creases, and the airbag is easily fatigued and broken. The contraction length of the PAMs also changes nonlinearly, which is not easy to control. For a larger number of creases, the minimum contraction length of the PAM increases, and the working space of the robot decreases. Therefore, selecting the appropriate number of creases is necessary. By comparing the PAM contraction effects with different number of creases, 4 is selected as the most appropriate number.

Fig.3 Pneumatic artificial muscles with different numbers of creases. (a) Few numbers of creases, (b) appropriate number of creases, and (c) large number of creases.

Full size|PPT slide

The load carrying capacity of the PAM is related to air pressure inside ( P_{\mathrm{i}\mathrm{n}}) the airbag and the airbag cross-sectional area. The model assumes that the cross-section of the airbag is a square and the side length is d_{\mathrm{a}}. When the PAM is negative pressure, the load carrying capacity is F_{\mathrm{c}}=(P_{\mathrm{o}\mathrm{u}\mathrm{t}}-P_{\mathrm{i}\mathrm{n}})d_{\mathrm{a}}^2, where P_{\mathrm{o}\mathrm{u}\mathrm{t}} is the atmospheric pressure. The maximum load could be achieved when the PAM is a vacuum. The relationship between P_{\mathrm{i}\mathrm{n}}, d_{\mathrm{a}}, and F_{\mathrm{c}} is shown in Fig.4(a). Increasing d_{\mathrm{a}} can increase the load carrying capacity of PAM, but it also increases the diameter and volume of the robot prototype. Therefore, d_{\mathrm{a}}=15\mathrm{m}\mathrm{m} is selected under comprehensive consideration, and the maximum load of single PAM is 22.7 N.

Fig.4 (a) Relationship between P_{\mathrm{i}\mathrm{n}}, d_{\mathrm{a}}, and F_{\mathrm{c}}; (b) relationship between muscle length, inflation time, and deflation time.

Full size|PPT slide

The length of the PAM is only related to P_{\mathrm{i}\mathrm{n}}, which is difficult to adjust accurately. When the inflating/deflating speed is constant, P_{\mathrm{i}\mathrm{n}} is proportional to inflation/deflation time. Therefore, the length of the airbag could be accurately controlled by adjusting inflation/deflation time. The relationship between airbag length and inflation/deflation time is shown in Fig.4(b).

When the PAM is inflated/deflated, the PAM length and internal air pressure change, which results in the variation of stiffness. One end of the PAM is fixed, and the other end was with 0.2 N load, as shown in Fig.5(a). Then, the configuration changes of the PAM with different lengths were measured, and the results are shown in Fig.5(b). When the internal air pressure and length of the PAM decrease, the endpoint displacement of the PAM gradually increases, which indicates that PAM stiffness gradually decreases.

Fig.5 (a) PAM stiffness measured method; (b) experimental results of PAM configuration with different lengths.

Full size|PPT slide

According to selected design parameters, the production steps of the soft robot were as follows:

Step 1. Use polyethylene film to make an airbag with a length and width of 15 mm and height of 100 mm.

Step 2. Use a PVC sheet to make an origami structure with a cross-sectional side length of 14.5 mm as an endoskeleton.

Step 3. Place the endoskeleton of the origami structure into the airbag, and glue it to the airbag at the crease. Seal one end of the airbag, and seal the other end after inserting the air inlet tube (to ensure air tightness, the airbag can be further sealed with silicone at the airtight seal).

Repeat Steps 1–3 to produce multiple PAMs.

Step 4. Fix the three PAM axes symmetrically distributed on the constraint disk to form a soft module. Connect the PAM and the constraint disk by pneumatic connectors and fixed terminals.

Repeat Step 4 to produce 10 soft modules.

Step 5. Connect 10 flexible modules in series to form a soft robot.

Based on the previous structural design, a soft robot prototype was developed, as shown in Fig.6(a). The system consists of a power supply, control module, and a soft robot. The soft robot is covered with a woven mesh, and the structural details are shown. The control module consists of an upper computer (PC), a lower computer (STM32), 60 relays, and 60 two-position two-way solenoid valves, as shown in Fig.6(b). The air compressor and vacuum pump power the system. The upper computer communicates with the lower computer via Bluetooth. The upper computer sends commands to the lower computer to manage the switching status of the relays and solenoid valves, thus controlling the state of air pressure inside the PAM. Then, the configuration of the robot could be controlled bending and deforming in 3D space, as shown in Fig.7. The parameters of the robot are listed in Tab.1.

Tab.1 Parameters of the soft parallel robot

Parameter	Value
Robot length, L	1000 mm
Robot diameter, d_{\mathrm{r}}	40 mm
Elastic rod diameter, d	2.5 mm
Young’s modulus, E	6.2 MPa
Shear modulus, G	2.4 MPa
Robot gravity, f_{\mathrm{g}}	4 N

Fig.6 (a) Prototype of the designed soft parallel robot; (b) control module of the robot.

Full size|PPT slide

Fig.7 Soft robot prototype bending in 3D space (with and without load).

Full size|PPT slide

3 Variable curvature modeling approach

The kinematic analysis of the soft robot is equivalent to the analysis of the backbone. Therefore, this section serves as the kinematics of the robot by establishing the equilibrium equations of the backbone.

The inertial coordinate system is O- \xi \eta \zeta, and the spindle coordinate system at any point N on the skeleton is N-xyz. The unit vectors of the x, y, and z axes are {\mathit{e}}_1, {\mathit{e}}_2, and {\mathit{e}}_3, respectively. The position of O-\xi \eta \zeta after turning around the \zeta axis by angle \alpha is O-{\xi }_1{\eta }_1{\zeta }_1, the position after turning around the {\xi }_1 axis by angle \beta is O-{\xi }_2{\eta }_2{\zeta }_2, and the position after turning around the {\zeta }_2 axis by angle \theta is O-{\xi }_3{\eta }_3{\zeta }_3, as shown in Fig.8. O-{\xi }_3{\eta }_3{\zeta }_3 and N- xyz have the same pose. The direction cosine matrix R between the spindle coordinate system N- xyz and inertial coordinate system O-\xi \eta \zeta is

Fig.8 Force on the micro element arc segment of the backbone.

Full size|PPT slide

(2)\mathit{R}=\left[\begin{array}{ccc}\mathrm{c}\alpha \mathrm{c}\theta -\mathrm{c}\beta \mathrm{s}\alpha \mathrm{s}\theta & -\mathrm{c}\alpha \mathrm{s}\theta -\mathrm{c}\beta \mathrm{s}\alpha \mathrm{c}\theta & \mathrm{s}\beta \mathrm{s}\alpha \\ \mathrm{s}\alpha \mathrm{c}\theta +\mathrm{c}\beta \mathrm{c}\alpha \mathrm{s}\theta & -\mathrm{s}\alpha \mathrm{s}\theta +\mathrm{c}\beta \mathrm{c}\alpha \mathrm{c}\theta & -\mathrm{s}\beta \mathrm{c}\alpha \\ \mathrm{s}\beta \mathrm{s}\theta & \mathrm{s}\beta \mathrm{c}\theta & \mathrm{c}\beta \end{array}\right],

where s and c are abbreviations for the sin and cos functions, respectively.

The backbone is discretized into n units, with a total of (n+1) nodes ( N_0,N_1,... ,N_n). N_i and N_{i+1} are two adjacent points. The position vectors of N_i and N_{i+1} with respect to the origin of the inertial coordinate system are r and (\mathit{r}+\mathrm{\Delta }\mathit{r}), respectively, and the arc coordinates are s and (s+\mathrm{\Delta }s), respectively. The internal force at point N_i is F, the internal force at point N_{i+1} is ( \mathit{F}+\mathrm{\Delta }\mathit{F}), and the external force to N_i is {\mathit{f}}_i. When the micro element arc segment N_i{N}_{i+1} is in equilibrium, the force balance condition is

(3)\mathrm{\Delta }{\mathit{F}}_i+{\mathit{f}}_i\mathrm{\Delta }s=0.

Dividing Eq. (3) by \mathrm{\Delta }s and letting \mathrm{\Delta }s \to 0 obtains

(4){\mathit{F}}_i^{\prime }+{\mathit{f}}_i=0 ,

where {\mathit{F}}_i^{\prime } indicates the derivative of the internal force F with respect to the arc coordinate s.

Integrating Eq. (4),

(5)\mathit{F}(s)={\int }_0^s\left(-\mathit{f}\right)\mathrm{d}s.

The total elastic potential energy E_{\mathrm{t}} of the backbone is equal to the sum of elastic strain energy E_{\mathrm{e}} and external force potential energy E_{\mathrm{p}}:

(6)E_{\mathrm{t}}=E_{\mathrm{e}}+E_{\mathrm{p}}.

The variation of Eq. (6) is

(7)\delta E_{\mathrm{t}} =\delta E_{\mathrm{e}}+\delta E_{\mathrm{p}}.

When the robot is in equilibrium, the variable fraction of the total potential energy is 0.

(8)\delta E_{\mathrm{t}} =0.

The relative rotation angles of N_i- xyz and N_{i+1}-xyz is \mathrm{\Delta } \text{ϕ}. \text{ω} is the rate of change of \text{ϕ} with respect to s.

(9)\text{ω} =\underset{\mathrm{\Delta }s\to 0}{lim}\frac{\mathrm{\Delta }\text{ϕ}}{\mathrm{\Delta }s}.

The Young’s modulus of the backbone is E, the shear modulus is G, and the elastic rod diameter is d. The flexural stiffness of the x and y axes and torsional stiffness of the z axis are k_1, k_2, and k_3, respectively.

(10)\{\begin{array}{c}{k}_1=k_2=\frac{E\text{π}{d}^4}{64}\text{,}\\ k_3=\frac{G\text{π}{d}^4}{32}.\end{array}

The elastic strain energy is

(11)E_{\mathrm{e}}=\frac{1}{2}{\int }_0^L\left(k_1{\omega }_1^2+k_2{\omega }_2^2+k_3{\omega }_3^2\right)\mathrm{d}s ,

where L is the robot length. The variation of Eq. (11) is

(12)\delta E_{\mathrm{e}} =\delta \frac{1}{2}{\int }_0^L\left(k_1{\omega }_1^2+k_2{\omega }_2^2+k_3{\omega }_3^2\right)\mathrm{d}s .

External potential energy can be expressed in terms of the internal force as

(13)E_{\mathrm{p}}=-{\int }_0^L\left(\mathit{F}\cdot {\mathit{e}}_3\right)\mathrm{d}s.

The variable of Eq. (13) is

(14)\delta E_{\mathrm{p}} =-\delta {\int }_0^L\left(\mathit{F}\cdot {\mathit{e}}_3\right)\mathrm{d}s.

Substituting Eqs. (7), (12), and (14) into Eq. (8) obtains

(15)\delta {\int }_0^L\mathrm{\Gamma }\mathrm{d}s=0,

where

(16)\mathrm{\Gamma }=\frac{1}{2}\left(k_1{\omega }_1^2+k_2{\omega }_2^2+k_3{\omega }_3^2\right) -\mathit{F}\cdot {\mathit{e}}_3.

Taking \alpha, \beta, and \theta as the three generalized coordinates for determining the section pose, the generalized coordinate symbols, q_1=\alpha, q_2=\beta, and q_3=\theta, are introduced. \mathrm{\Delta }\text{ϕ}, \alpha, \beta, and \theta satisfy

(17)\{\begin{array}{c}\mathrm{\Delta }{\varphi }_i=\sum _{j=1}^3{\mathrm{\Psi }}_{ij}{q}_j,\\ q_j=\sum _{i=1}^3{\mathrm{\Phi }}_{ij}{\varphi }_i\text{,}\end{array}

(18)\{\begin{array}{c}{\omega }_i=\sum _{j=1}^3{\mathrm{\Psi }}_{ij}{q^{\prime }}_j,\\ {q^{\prime }}_j=\sum _{i=1}^3{\mathrm{\Phi }}_{ij}{\omega }_i\text{,}\end{array}

where

(19)\begin{array}{rl}\mathbf{\Psi }& =\left[\begin{array}{ccc}\mathrm{s}{q}_2\mathrm{s}{q}_3& \mathrm{c}{q}_3& 0\\ \mathrm{s}{q}_2\mathrm{c}{q}_3& -\mathrm{s}{q}_3& 0\\ \mathrm{c}{q}_2& 0& 1\end{array}\right]\mathbf{\Phi }={\mathbf{\Psi }}^{-1}\\ & =\left[\begin{array}{ccc}\csc q_2& \csc q_2\mathrm{c}{q}_3& 0\\ \mathrm{c}{q}_3& -\mathrm{s}{q}_3& 0\\ -\cot q_2\mathrm{s}{q}_3& \cot q_2\mathrm{c}{q}_3& 1\end{array}\right],i,j=1,2,3.\end{array}

{\mathrm{\Psi }}_{ij} and {\mathrm{\Phi }}_{ij} can be written as

(20)\{\begin{array}{c}{\mathrm{\Psi }}_{ij}=\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q^{\prime }}_j}\text{,}\\ {\mathrm{\Phi }}_{ij}=\frac{\mathrm{\partial }{q^{\prime }}_j}{\mathrm{\partial }{\omega }_i}.\end{array}

Substituting Eq. (20) into Eqs. (17) and (18), and then considering the variation in Eqs. (17) and (18) obtains

(21)\delta {\varphi }_i= \sum _{j=1}^3\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q^{\prime }}_j}\delta q_j,

(22)\delta {\omega }_i= \sum _{j=1}^3\left(\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q^{\prime }}_j}\delta {q^{\prime }}_j+\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q}_j}\delta q_j\right).

Calculating the derivative of \delta {\varphi }_i with respect to s yields

(23)\frac{\mathrm{d}}{\mathrm{d}s}( \delta {\varphi }_i)=\sum _{j=1}^3\left[\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q^{\prime }}_j}\right)\delta q_j+\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q^{\prime }}_j}\frac{\mathrm{d}}{\mathrm{d}s}\left(\delta q_j\right)\right].

Subtracting Eq. (22) from Eq. (23) and considering the exchangeability of the generalized coordinate differentiation and variational differentiation of the complete system obtain

(24)\frac{\mathrm{d}}{\mathrm{d}s}( \delta {\varphi }_i)-\delta {\omega }_i=\sum _{k=1}^3\sum _{j=1}^3\left[\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q^{\prime }}_j}\right)-\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q}_j}\right]\frac{\mathrm{\partial }{q^{\prime }}_j}{\mathrm{\partial }{\omega }_k}\delta {\varphi }_k.

From Eq. (20), the following is obtained

(25)\{\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q^{\prime }}_j}\right)\text{=}\sum _{k=1}^3\frac{\mathrm{\partial }{\mathrm{\Psi }}_{ij}}{\mathrm{\partial }{q}_k}{q^{\prime }}_k,\\ \frac{\mathrm{\partial }{\omega }_i}{\mathrm{\partial }{q}_j}\text{=}\sum _{k=1}^3\frac{\mathrm{\partial }{\mathrm{\Psi }}_{ik}}{\mathrm{\partial }{q}_j}{q^{\prime }}_k.\end{array}

Let

(26){\gamma }_{km}^i= \sum _{j=1}^3\sum _{l=1}^3\left(\frac{\mathrm{\partial }{\mathrm{\Psi }}_{ij}}{\mathrm{\partial }{q}_l}-\frac{\mathrm{\partial }{\mathrm{\Psi }}_{il}}{\mathrm{\partial }{q}_j}\right){\mathrm{\Phi }}_{jk}{\mathrm{\Phi }}_{lm} , i,k,m =1,2, 3.

Substituting Eqs. (25) and (26) into Eq. (24) obtains

(27)\delta {\omega }_i=(\delta {\varphi }_i{)}^{\prime }-\sum _{k=1}^3\sum _{m=1}^3{\omega }_m{\gamma }_{km}^i\delta {\varphi }_k.

By reducing Eq. (15) to a function of {\varphi }_i and {\omega }_i,

(28)\delta {\int }_0^L\mathrm{\Gamma }\mathrm{d}s=\sum _{i=1}^3\left[{\int }_0^L\left(\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_i}\delta {\omega }_i+\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }{\varphi }_i}\delta {\varphi }_i\right)\mathrm{d}s\right],

where

(29)\frac{\mathrm{\partial }}{\mathrm{\partial }{\varphi }_i}=\sum _{j=1}^3\frac{\mathrm{\partial }{q^{\prime }}_j}{\mathrm{\partial }{\omega }_i}\frac{\mathrm{\partial }}{\mathrm{\partial }{q}_j}.

Substituting Eq. (29) into Eq. (28) yields

(30)\begin{array}{rl}& \delta {\int }_0^L\mathrm{\Gamma }\mathrm{d}s=\\ \sum _{k=1}^m& \{{\int }_0^L\left[\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\varphi }_k}-\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_k}\right)-\sum _{i=1}^3\sum _{m=1}^3\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_i}{\omega }_m{\gamma }_{km}^i\right]\delta {\varphi }_k\mathrm{d}s\\ & +{\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_k}\delta {\varphi }_k|}_0^L\},\end{array}

where the edge value term can be reduced to a boundary condition,

(31){\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_k}\delta {\varphi }_k|}_0^L={\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_k}{\mathit{e}}_k^{\mathrm{T}}\delta \text{ϕ}|}_0^L=0.

Equation (30) can be transformed into

(32)\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_k}\right)-\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\varphi }_k}+\sum _{i=1}^3\sum _{m=1}^3\frac{\mathrm{\partial }\mathrm{\Gamma }}{\mathrm{\partial }{\omega }_i}{\omega }_m{\gamma }_{km}^i=0.

Substituting Eq. (16) into Eq. (32) yields

(33)\{\begin{array}{c}{k}_1{{\omega }^{\prime }}_{1,i}+(k_3-k_2){\omega }_{2,i}{\omega }_{3,i}-F_{2,i}=0,\\ k_2{{\omega }^{\prime }}_{2,i}+(k_1-k_3){\omega }_{1,i}{\omega }_{3,i}+F_{1,i}=0,\\ k_3{{\omega }^{\prime }}_{3,i}=0,\end{array}

where F_1 and F_2 are the projections of the internal force {\mathit{F}}_i on the x and y axes of the spindle coordinate system N_i-xyz, respectively.

To avoid singularities in the calculation of Euler angles, quaternions (x_1,x_2,x_3,x_4) are used to describe the cross-sectional poses. \text{ω} with quaternions satisfies

(34)\{\begin{array}{c}{\omega }_{1,i}=2\left(-x_{2,i}{x^{\prime }}_{1,i}+x_{1,i}{x^{\prime }}_{2,i}+x_{4,i}{x^{\prime }}_{3,i}-x_{3,i}{x^{\prime }}_{4,i}\right),\\ {\omega }_{2,i}=2\left(-x_{3,i}{x^{\prime }}_{1,i}-x_{4,i}{x^{\prime }}_{2,i}+x_{1,i}{x^{\prime }}_{3,i}+x_{2,i}{x^{\prime }}_{4,i}\right),\\ {\omega }_{3,i}=2\left(-x_{4,i}{x^{\prime }}_{1,i}+x_{3,i}{x^{\prime }}_{2,i}-x_{2,i}{x^{\prime }}_{3,i}+x_{1,i}{x^{\prime }}_{4,i}\right),\end{array}

(35)\{\begin{array}{c}{{\omega }^{\prime }}_{1,i}=2\left(-x_{2,i}{x}_{1,i}^{″}+x_{1,i}{x}_{2,i}^{″}+x_{4,i}{x}_{3,i}^{″}-x_{3,i}{x}_{4,i}^{″}\right),\\ {{\omega }^{\prime }}_{2,i}=2\left(-x_{3,i}{x}_{1,i}^{″}-x_{4,i}{x}_{2,i}^{″}+x_{1,i}{x}_{3,i}^{″}+x_{2,i}{x}_{4,i}^{″}\right),\\ {{\omega }^{\prime }}_{3,i}=2\left(-x_{4,i}{x}_{1,i}^{″}+x_{3,i}{x}_{2,i}^{″}-x_{2,i}{x}_{3,i}^{″}+x_{1,i}{x}_{4,i}^{″}\right).\end{array}

Substituting Eqs. (34) and (35) into Eqs. (33) yields

(36)\begin{array}{rl}{h}_{1,i}=& 2k_1(-x_{2,i}{x}_{1,i}^{″}+x_{1,i}{x}_{2,i}^{″}+x_{4,i}{x}_{3,i}^{″}-x_{3,i}{x}_{4,i}^{″})\\ & +4\left(k_3-k_2\right)\left(-x_{3,i}{x^{\prime }}_{1,i}-x_{4,i}{x^{\prime }}_{2,i}+x_{1,i}{x^{\prime }}_{3,i}+x_{2,i}{x^{\prime }}_{4,i}\right)\\ & \cdot \left(-x_{4,i}{x^{\prime }}_{1,i}+x_{3,i}{x^{\prime }}_{2,i}-x_{2,i}{x^{\prime }}_{3,i}+x_{1,i}{x^{\prime }}_{4,i}\right)-F_{2,i},\end{array}

(37)\begin{array}{rl}{h}_{2,i}=& 2k_2(-x_{3,i}{x}_{1,i}^{″}-x_{4,i}{x}_{2,i}^{″}+x_{1,i}{x}_{3,i}^{″}+x_{2,i}{x}_{4,i}^{″})\\ & +4(k_1-k_3)(-x_{2,i}{x^{\prime }}_{1,i}+x_{1,i}{x^{\prime }}_{2,i}+x_{4,i}{x^{\prime }}_{3,i}-x_{3,i}{x^{\prime }}_{4,i})\\ & \cdot (-x_{4,i}{x^{\prime }}_{1,i}+x_{3,i}{x^{\prime }}_{2,i}-x_{2,i}{x^{\prime }}_{3,i}+x_{1,i}{x^{\prime }}_{4,i})+F_{1,i},\end{array}

(38)h_{3,i}=2k_3(-x_{4,i}{x}_{1,i}^{″}+x_{3,i}{x}_{2,i}^{″}-x_{2,i}{x}_{3,i}^{″}+x_{1,i}{x}_{4,i}^{″}).

The sum of squares of quaternions is 1:

(39)h_{4,i}=x_{1,i}^2+x_{2,i}^2+ x_{3,i}^2+x_{4,i}^2 -1.

For nodes N_0, a three-point differential format is used. For the other nodes, a forward differential format is used, as shown in Eqs. (40) and (41):

(40)\{\begin{array}{c}{x^{\prime }}_{k,0}=\frac{-3x_{k,0}+4x_{k,1}-x_{k,2}}{\mathrm{\Delta }s},\\ x_{k,0}^{″}=\frac{-3{x^{\prime }}_{k,0}+4{x^{\prime }}_{k,1}-{x^{\prime }}_{k,2}}{\mathrm{\Delta }s}\text{,}\end{array}

(41)\{\begin{array}{c}{x^{\prime }}_{k,i}=\frac{x_{k,i}-x_{k,i-1}}{\mathrm{\Delta }s}\text{,}\\ x_{k,i}^{″}=\frac{{x^{\prime }}_{k,i}-{x^{\prime }}_{k,i-1}}{\mathrm{\Delta }s}\text{,}\end{array}

where k=1,2,3,4, and i=1,2,...,n.

The R in Eq. (2) can be expressed in quaternions as

(42)\mathit{R}\left(s_i\right)=\left[\begin{array}{ccc}{x}_{1,i}^2+x_{2,i}^2-x_{3,i}^2-x_{4,i}^2& 2(x_{2,i}{x}_{3,i}-x_{1,i}{x}_{4,i})& 2(x_{2,i}{x}_{4,i}+x_{1,i}{x}_{3,i})\\ 2(x_{2,i}{x}_{3,i}+x_{1,i}{x}_{4,i})& x_{1,i}^2-x_{2,i}^2+x_{3,i}^2-x_{4,i}^2& 2(x_{3,i}{x}_{4,i}-x_{1,i}{x}_{2,i})\\ 2(x_{2,i}{x}_{4,i}-x_{1,i}{x}_{3,i})& 2(x_{3,i}{x}_{4,i}+x_{1,i}{x}_{2,i})& x_{1,i}^2-x_{2,i}^2-x_{3,i}^2+x_{4,i}^2\end{array}\right].

Then,

(43)\left[\begin{array}{c}{F}_{1,i}\\ F_{2,i}\\ F_{3,i}\end{array}\right]\text{=}\mathit{R}{\left(s_i\right)}^{\mathrm{T}} \left[\begin{array}{c}{F}_{\xi ,i}\\ F_{\eta ,i}\\ F_{\zeta ,i}\end{array}\right].

Substituting Eqs. (40), (41), and (43) into Eqs. (36)–(39), when the robot is in equilibrium,

(44)\mathit{h}={\left[\begin{array}{cccc}{\mathit{h}}_1& {\mathit{h}}_2& {\mathit{h}}_3& {\mathit{h}}_4\end{array}\right]}^{\mathrm{T}}=\mathbf{0}.

The robot usually needs to move to a desired point (P_1,P_2, P_3) in space to perform tasks such as manipulation, so adding end position constraints to the robot is necessary. The position of each discrete node of the robot can be obtained by integrating the z axis.

(45)\mathit{r}(s_i)={\int }_0^{s_i}{\mathit{e}}_3\left(\sigma \right)\mathrm{d}\sigma ,

where \sigma is the integration variable.

From Eq. (42), the following is obtained:

(46)\mathit{r}(s)=\left[\begin{array}{c}\xi (s)\\ \eta (s)\\ \zeta (s)\end{array}\right]= \left[\begin{array}{c}2{\int }_0^s(x_2(\sigma )x_4(\sigma )+x_1(\sigma )x_3(\sigma ))\mathrm{d}\sigma \\ 2{\int }_0^s(x_3(\sigma )x_4(\sigma )-x_1(\sigma )x_2(\sigma ))\mathrm{d}\sigma \\ {\int }_0^s(x_1^2(\sigma )-x_2^2(\sigma )-x_3^2(\sigma )+x_4^2(\sigma ))\mathrm{d}\sigma \end{array}\right].

Interpolation of Eq. (46) by Simpson’s formula yields

(47)\{\begin{array}{rl}{\xi }_i=& \frac{\mathrm{\Delta }s}{3}((2x_{2,i-1}+x_{2,i})x_{4,i-1}+(2x_{2,i}+x_{2,i-1})x_{4,i}\\ & +(2x_{1,i-1}+x_{1,i})x_{3,i-1}+(2x_{1,i}+x_{1,i-1})x_{3,i})+{\xi }_{i-1},\\ {\eta }_i=& \frac{\mathrm{\Delta }s}{3}((2x_{3,i-1}+x_{3,i})x_{4,i-1}+(2x_{3,i}+x_{3,i-1})x_{4,i}\\ & -(2x_{1,i-1}+x_{1,i})x_{2,i-1}-(2x_{1,i}+x_{1,i-1})x_{2,i})+{\eta }_{i-1},\\ {\zeta }_i=& \frac{\mathrm{\Delta }s}{3}(x_{1,i-1}^2+x_{1,i-1}{x}_{1,i}+x_{1,i}^2-x_{2,i-1}^2-x_{2,i-1}{x}_{2,i}-x_{2,i}^2\\ & -x_{3,i-1}^2-x_{3,i-1}{x}_{3,i}-x_{3,i}^2+x_{4,i-1}^2+x_{4,i-1}{x}_{4,i}+x_{4,i}^2)+{\zeta }_{i-1}.\end{array}

The position constraint at the endpoint of the robot can be expressed as

(48){\mathit{h}}_5=\left[\begin{array}{c}{\xi }_n-P_1\\ {\eta }_n-P_2\\ {\zeta }_n-P_3\end{array}\right]= \mathbf{0}.

Sometimes the attitude of the endpoint of the robot also needs to be constrained:

(49){\mathit{h}}_6=\left[\begin{array}{c}{x}_{1,n}-p_1\\ x_{2,n}-p_2\\ x_{3,n}-p_3\\ x_{4,n}-p_4\end{array}\right]=\mathbf{0},

where p_1,p_2,p_3, and p_4 are the specified values.

Equation (44) can be expanded as

(50)\mathit{h}={\left[\begin{array}{cccc}\begin{array}{ccc}{\mathit{h}}_1& {\mathit{h}}_2& {\mathit{h}}_3\end{array}& {\mathit{h}}_4& {\mathit{h}}_5& {\mathit{h}}_6\end{array}\right]}^{\mathrm{T}}=\mathbf{0}.

Equation (50) is a soft robot model. Based on the least-squares method, the solution problem in Eq. (50) can be transformed into the following optimization problem. The objective function is as follows:

(51)\{\begin{array}{c}minH,\\ H=\frac{1}{2}\sum _{i=1}^{4n+7}{h}_i^2(x).\end{array}

The Levenberg–Marquardt algorithm can be used to solve Eq. (51), and the results of the quaternion solution are expressed in Eqs. (42) and (47) to obtain the pose and position of each discrete node. The spatial poses of the constrained disk scan are obtained by interpolation. The pose matrix of the tth constrained disk is

(52){\mathit{T}}_t=\left[\begin{array}{cc}{\mathit{R}}_t& {\mathit{r}}_t\\ 0& 1\end{array}\right].

Then, the position and pose matrix of the tth constrained disk with respect to the (t-1)\mathrm{t}\mathrm{h} constrained disk is

(53){\mathit{T}}_{t-1}^t={\mathit{T}}_{t-1}^{-1}{\mathit{T}}_t=\left[\begin{array}{cc}{\mathit{R}}_{t-1}^t& {\mathit{r}}_{t-1}^t\\ 0& 1\end{array}\right].

The length of the ith PAM in the tth soft module can be solved using the inverse kinematic equation of the parallel robot.

(54)l\mathrm{\_}{p}_{t,i}=\Vert \mathit{l}\text{\_}{\mathit{p}}_{t,i}\Vert =\sqrt{{(({\mathit{R}}_t^{\text{T}}{\mathit{R}}_{t+1}{\mathit{a}}_{t,i}+{\mathit{r}}_t^{t+1})-{\mathit{a}}_{t-1,i})}^{\mathrm{T}}(({\mathit{R}}_t^{\text{T}}{\mathit{R}}_{t+1}{\mathit{a}}_{t,i}+{\mathit{r}}_t^{t+1})-{\mathit{a}}_{t-1,i})},

where {\mathit{a}}_t is defined in Eq. (1).

By leveraging the relationship between muscle length and inflation and deflation time shown in Fig.4(b), the soft robot can be controlled to bend into the desired form by controlling the relay and solenoid valve on/off time.

4 Experimental validation of the proposed model

Soft robots usually need to move along a specific trajectory when performing a task. The trajectory can be discretized as a set of target points, and the motion along the predefined trajectory can be achieved by controlling the robot end points to reach these target points sequentially. Three sets of experiments, with no load, 0.2 N load, and 0.5 N load at the end of the robot were conducted to control the robot endpoints along a set of circulars, square, and cardioid curves. In each set of experiments, the trajectory was discretized into 20 points, and the robot was controlled to pass through these 20 points sequentially. A position sensor (3D Guidance trakSTAR) was used to measure the experimental position {\mathit{P}}_{\mathrm{a}} of the robot endpoint. The simulation position was denoted as P, and normalized error was defined as

(55)e=\frac{\Vert {\mathit{P}}_{\mathrm{a}}-\mathit{P}\Vert }{L}\times 100\text{\%}.

The configuration simulation results and endpoint position errors of the robot in the three sets of experiments, Experiment 1, Experiment 2, and Experiment 3, are shown in Fig.9 and Fig.10, Fig.11 and Fig.12, Fig.13 and Fig.14, respectively. In the trajectory tracking experiment, three groups of trajectory tracking results for the same path are used to reflect the effectiveness of the average error. The simulation and experimental results are very close. The errors of the robot endpoint were calculated, and the average errors of the three sets of experiments are 2.457%, 2.477%, and 2.514%. The robot endpoint load has a certain influence on the robot model solving and control accuracy, and calculation error gradually increases with increasing load.

Fig.9 Experiment 1: The endpoint trajectory is circular. (a) Spatial configuration of soft robot endpoints with different constraints, (b) \eta O\zeta projection view, (c) \eta O\xi projection view, and (d) comparison of endpoint simulation and experimental position.

Full size|PPT slide

Fig.10 Position errors of robot endpoints in Experiment 1.

Full size|PPT slide

Fig.11 Experiment 2: The endpoint trajectory is square. (a) Spatial configuration of soft robot endpoints with different constraints, (b) \eta O\zeta projection view, (c) \eta O\xi projection view, and (d) comparison of endpoint simulation and experimental position.

Full size|PPT slide

Fig.12 Position errors of robot endpoints in Experiment 2.

Full size|PPT slide

Fig.13 Experiment 3: The endpoint trajectory is a heart-shaped curve. (a) Spatial configuration of soft robot endpoints with different constraints, (b) \xi O\zeta projection view, (c) \xi O\eta projection view, and (d) comparison of endpoint simulation and experimental position.

Full size|PPT slide

Fig.14 Position errors of robot endpoints in Experiment 3.

Full size|PPT slide

5 Conclusions

In this paper, a soft robot containing 10 soft modules was developed based on endoskeleton PAMs. Compared with existing pneumatic soft robots with unstable motion and only specific directional and angle bending, the developed prototype can realize flexible bending motion in 3D space and can be controlled effectively. A novel kinematic modeling method for variable-curvature soft robots based on the minimum energy method was investigated, considering external loads. Based on this model, forward control of the soft robot was realized. The investigated model and the developed prototype were validated by three sets of experiments, and the results showed that the proposed modeling method could control the soft robot along the desired trajectory. This paper is expected to provide new research ideas for the design, modeling, and control of pneumatic soft robots.

6 Nomenclature

{\mathit{a}}_{t,i}	Position coordinates of the ith muscle endpoint in the tth soft module
d	Elastic rod diameter
d_{\mathrm{a}}	Side length of the airbag
d_{\mathrm{r}}	Robot diameter
e	Robot endpoint position error
E	Young’s modulus
E_{\mathrm{e}}	Elastic strain energy
E_{\mathrm{p}}	External force potential energy
E_{\mathrm{t}}	Total elastic potential energy
{\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3	Unit vectors of the x, y, and z axes, respectively
f_{\mathrm{g}}	Robot gravity
Fc	Load carrying capacity of PAM
{\mathit{f}}_i	External force of the ith discretized nodes N_i
F	Internal force
k_1,k_2	Flexural stiffnesses of the x and y axes, respectively
k_3	Torsional stiffness of the z axis
l_p	Pneumatic artificial muscle length
L	Robot length
N_i	Theith discretized nodes
P	Simulation position of the robot endpoint
P_{\mathrm{a}}	Experimental position of the robot endpoint
P_{\mathrm{i}\mathrm{n}}	Air pressure inside the PAM
P_{\mathrm{o}\mathrm{u}\mathrm{t}}	Atmospheric pressure
r_{\mathrm{c}}	Constraint disk radius
r	Position vector of the robot discrete nodes
(P_1,P_2,P_3)	Position of the desired point
(α,β,θ)	Euler angles
\text{ϕ}	Relative rotation angle of adjacent coordinate systems
\mathit{\omega }	Derivative of \text{ϕ} with respect to arc coordinate s

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant Nos. 51975566, 61821005, and U1908214), and the Key Research Program of Frontier Sciences, CAS, China (Grant No. ZDBS-LY-JSC011).

References

1 RusD, TolleyM T. Design, fabrication and control of soft robots. Nature, 2015, 521( 7553): 467– 475

2 LiuT L, XuW F, YangT W, LiY M. A hybrid active and passive cable-driven segmented redundant manipulator: design, kinematics, and planning. IEEE/ASME Transactions on Mechatronics, 2021, 26( 2): 930– 942

3 GuanQ H, SunJ, LiuY J, WereleyN M, LengJ S. Novel bending and helical extensile/contractile pneumatic artificial muscles inspired by elephant trunk. Soft Robotics, 2020, 7( 5): 597– 614

4 AubinC A, ChoudhuryS, JerchR, ArcherL A, PikulJ H, ShepherdR F. Electrolytic vascular systems for energy-dense robots. Nature, 2019, 571( 7763): 51– 57

5 KimY, YukH, ZhaoR K, ChesterS A, ZhaoX H. Printing ferromagnetic domains for untethered fast-transforming soft materials. Nature, 2018, 558( 7709): 274– 279

6 LiG R, ChenX P, ZhouF H, LiangY M, XiaoY H, CaoX N, ZhangZ, ZhangM Q, WuB S, YinS Y, XuY, FanH B, ChenZ, SongW, YangW J, PanB B, HouJ Y, ZouW F, HeS P, YangX X, MaoG Y, JiaZ, ZhouH F, LiT F, QuS X, XuZ B, HuangZ L, LuoY W, XieT, GuJ, ZhuS Q, YangW. Self-powered soft robot in the Mariana Trench. Nature, 2021, 591( 7848): 66– 71

7 ParkS J, GazzolaM, ParkK S, ParkS, Di SantoV, BlevinsE L, LindJ U, CampbellP H, DauthS, CapulliA K, PasqualiniF S, AhnS, ChoA, YuanH Y, MaozB M, VijaykumarR, ChoiJ W, DeisserothK, LauderG V, MahadevanL, ParkerK K. Phototactic guidance of a tissue-engineered soft-robotic ray. Science, 2016, 353( 6295): 158– 162

8 WehnerM, TrubyR L, FitzgeraldD J, MosadeghB, WhitesidesG M, LewisJ A, WoodR J. An integrated design and fabrication strategy for entirely soft, autonomous robots. Nature, 2016, 536( 7617): 451– 455

9 ZhuM J, DoT N, HawkesE, VisellY. Fluidic fabric muscle sheets for wearable and soft robotics. Soft Robotics, 2020, 7( 2): 179– 197

10 CappelloL, GallowayK C, SananS, WagnerD A, GranberryR, EngelhardtS, HaufeF L, PeisnerJ D, WalshC J. Exploiting textile mechanical anisotropy for fabric-based pneumatic actuators. Soft Robotics, 2018, 5( 5): 662– 674

11 TolleyM T, ShepherdR F, MosadeghB, GallowayK C, WehnerM, KarpelsonM, WoodR J, WhitesidesG M. A resilient, untethered soft robot. Soft Robotics, 2014, 1( 3): 213– 223

12 KatzschmannR K, MarcheseA D, RusD. Autonomous object manipulation using a soft planar grasping manipulator. Soft Robotics, 2015, 2( 4): 155– 164

13 PeeleB N, WallinT J, ZhaoH C, ShepherdR F. 3D printing antagonistic systems of artificial muscle using projection stereolithography. Bioinspiration & Biomimetics, 2015, 10( 5): 055003

14 TerrynS, BrancartJ, LefeberD, Van AsscheG, VanderborghtB. Self-healing soft pneumatic robots. Science Robotics, 2017, 2( 9): eaan4628

15 LiS Y, WangK W. Fluidic origami: a plant-inspired adaptive structure with shape morphing and stiffness tuning. Smart Materials and Structures, 2015, 24( 10): 105031

16 RobertsonM A, PaikJ. New soft robots really suck: vacuum-powered systems empower diverse capabilities. Science Robotics, 2017, 2( 9): eaan6357

17 GorissenB, ReynaertsD, KonishiS, YoshidaK, KimJ W, De VolderM. Elastic inflatable actuators for soft robotic applications. Advanced Materials, 2017, 29( 43): 1604977

18 KimW, ByunJ, KimJ K, ChoiW Y, JakobsenK, JakobsenJ, LeeD Y, ChoK J. Bioinspired dual-morphing stretchable origami. Science Robotics, 2019, 4( 36): eaay3493

19 MartinezR V, FishC R, ChenX, WhitesidesG M. Elastomeric origami: programmable paper-elastomer composites as pneumatic actuators. Advanced Functional Materials, 2012, 22( 7): 1376– 1384

20 YangD, VermaM S, SoJ H, MosadeghB, KeplingerC, LeeB, KhashaiF, LossnerE, SuoZ G, WhitesidesG M. Buckling pneumatic linear actuators inspired by muscle. Advanced Materials Technologies, 2016, 1( 3): 1600055

21 JiaoZ D, JiC, ZouJ, YangH Y, PanM. Vacuum-powered soft pneumatic twisting actuators to empower new capabilities for soft robots. Advanced Materials Technologies, 2019, 4( 1): 1800429

22 LiS G, VogtD M, RusD, WoodR J. Fluid-driven origami-inspired artificial muscles. PNAS, 2017, 114( 50): 13132– 13137

23 DeshpandeA R, TseZ T H, RenH L. Origami-inspired bi-directional soft pneumatic actuator with integrated variable stiffness mechanism. In: Proceedings of the 2017 18th International Conference on Advanced Robotics (ICAR). IEEE, 2017, 417– 421

24 LeeJ G, RodrigueH. Origami-based vacuum pneumatic artificial muscles with large contraction ratios. Soft Robotics, 2019, 6( 1): 109– 117

25 RendaF, BoyerF, DiasJ, SeneviratneL. Discrete Cosserat approach for multisection soft manipulator dynamics. IEEE Transactions on Robotics, 2018, 34( 6): 1518– 1533

26 Burgner-KahrsJ, RuckerD C, ChosetH. Continuum robots for medical applications: a survey. IEEE Transactions on Robotics, 2015, 31( 6): 1261– 1280

27 WebsterR J III, JonesB A. Design and kinematic modeling of constant curvature continuum robots: a review. The International Journal of Robotics Research, 2010, 29( 13): 1661– 1683

28 SoflaM S, SadighM J, ZareinejadM. Design and dynamic modeling of a continuum and compliant manipulator with large workspace. Mechanism and Machine Theory, 2021, 164 : 104413

29 TanN, GuX Y, RenH L. Design, characterization and applications of a novel soft actuator driven by flexible shafts. Mechanism and Machine Theory, 2018, 122 : 197– 218

30 Garriga-CasanovasA, Rodriguez y BaenaF. Complete follow-the-leader kinematics using concentric tube robots. The International Journal of Robotics Research, 2018, 37( 1): 197– 222

31 SchillerL, SeibelA, SchlattmannJ. A lightweight simulation model for soft robot’s locomotion and its application to trajectory optimization. IEEE Robotics and Automation Letters, 2020, 5( 2): 1199– 1206

32 YangC H, GengS N, WalkerI, BransonD T, LiuJ G, DaiJ S, KangR J. Geometric constraint-based modeling and analysis of a novel continuum robot with shape memory alloy initiated variable stiffness. International Journal of Robotics Research, 2020, 39( 14): 1620– 1634

33 ZengW H, YanJ Y, YanK, HuangX, WangX F, ChengS S. Modeling a symmetrically-notched continuum neurosurgical robot with non-constant curvature and superelastic property. IEEE Robotics and Automation Letters, 2021, 6( 4): 6489– 6496

34 SinghI, AmaraY, MelinguiA, Mani PathakP, MerzoukiR. Modeling of continuum manipulators using pythagorean hodograph curves. Soft Robotics, 2018, 5( 4): 425– 442

35 GodageI S, Medrano-CerdaG A, BransonD T, GuglielminoE, CaldwellD G. Modal kinematics for multisection continuum arms. Bioinspiration & Biomimetics, 2015, 10( 3): 035002

36 YangJ Z, PengH J, ZhouW Y, ZhangJ, WuZ G. A modular approach for dynamic modeling of multisegment continuum robots. Mechanism and Machine Theory, 2021, 165 : 104429

37 YuanH, ZhouL L, XuW F. A comprehensive static model of cable-driven multi-section continuum robots considering friction effect. Mechanism and Machine Theory, 2019, 135 : 130– 149

38 HuangX J, ZouJ, GuG Y. Kinematic modeling and control of variable curvature soft continuum robots. IEEE/ASME Transactions on Mechatronics, 2021, 26( 6): 3175– 3185

39 BiezeT M, LargilliereF, KruszewskiA, ZhangZ K, MerzoukiR, DuriezC. Finite element method-based kinematics and closed-loop control of soft, continuum manipulators. Soft Robotics, 2018, 5( 3): 348– 364

40 SadatiS M H, ShivaA, RensonL, RuckerC, AlthoeferK, NanayakkaraT, BergelesC, HauserH, WalkerI D. Reduced order vs. discretized lumped system models with absolute and relative states for continuum manipulators. In: Proceedings of Royal Statistics Society International Conference. Belfast, 2019, 1– 10

41 GodageI S, WirzR, WalkerI D, WebsterIIIR J. Accurate and efficient dynamics for variable-length continuum arms: a center of gravity approach. Soft Robotics, 2015, 2( 3): 96– 106