Stiffness modelling and analysis of soft fluidic-driven robots using Lie theory

Abstract

Soft robots have been investigated for various applications due to their inherently superior deformability and flexibility compared to rigid-link robots. However, these robots struggle to perform tasks that require on-demand stiffness, that is, exerting sufficient forces within allowable deflection. In addition, the soft and compliant materials also introduce large deformation and non-negligible nonlinearity, which makes the stiffness analysis and modelling of soft robots fundamentally challenging. This paper proposes a modelling framework to investigate the underlying stiffness and the equivalent compliance properties of soft robots under different configurations. Firstly, a modelling and analysis methodology is described based on Lie theory. Here, we derive two sets (the piecewise constant curvature and Cosserat rod model) of compliance models. Furthermore, the methodology can accommodate the nonlinear responses (e.g., bending angles) resulting from elongation of robots. Using this proposed methodology, the general Cartesian stiffness or compliance matrix can be derived and used for configuration-dependent stiffness analysis. The proposed framework is then instantiated and implemented on fluidic-driven soft continuum robots. The efficacy and modelling accuracy of the proposed methodology are validated using both simulations and experiments.

1. Introduction

The field of robotics has experienced a paradigm shift from traditional robots, made of rigid links and joints, to soft robotic structures (Laschi et al., 2016). Soft robots often distinguish themselves from traditional robots by their intrinsic dexterity, biological morphology, and adaptability to unstructured environments (Shepherd et al., 2011; Arezzo et al., 2017; Rus and Tolley 2015). These versatile capabilities can be (in)directly credited to the compliant properties of the soft materials used to construct the robots. Soft robots have been investigated for safety-related applications in the industrial sector where robots work closely together with humans (Stilli et al., 2017), in minimally invasive interventions (Cianchetti et al., 2014; Abidi et al., 2018), and rehabilitation (Camp et al., 2021).

The stiffness (or its inverse, i.e., the compliance) of a soft robot can be mathematically described by the Young’s modulus defined as the relationship between an applied load and its corresponding deflection. The term soft robot fundamentally implies that the Young’s modulus of such a robot is a few orders less than 10⁹ − 10¹⁰ Pa, a Young’s modulus range for conventional robotic materials such as metals and hard plastics (Majidi 2014). Soft robots, however, are made of materials with a Young’s modulus of approximately 10⁴ − 10⁶ Pa, and they behave in a compliant way when in physical interaction with the environment. As a result, it might be challenging for soft robots to exert sufficient forces on demand. Hence, the investigation of stiffness property and stiffness adjustment is profoundly important to the design, modelling and control of soft robots (Mahvash and Dupont 2011; Naselli and Mazzolai 2021; Komatsu et al., 2021).

1.1. Related work

A number of stiffening mechanisms have been investigated to vary compliance of soft robots. These can be categorised into semi-active techniques, such as actuators that can change their Young’s modulus, and active actuators (Manti et al., 2016). Techniques for achieving stiffness variation in semi-active and active ways include the integration of granular material (Konstantinova et al., 2022), use of an alloy with a low melting point (Peters et al., 2019), and the combination of tendon-driven and air pressurisation through antagonistic actuation (Shiva et al., 2016). Researchers have highlighted challenges in embedding these stiffness mechanisms in soft manipulators for space-restricted applications, such as minimally invasive surgery (Runciman et al., 2019). Apart from the aforementioned mechanism design, the stiffness characterisation is advantageous for understanding the payload capability of the soft robots. An experimental-based stiffness assessment is presented by Ranzani et al. (2016), characterising the planar stiffness of a single pneumatic-driven modular manipulator. A similar approach is followed by Mustaza et al. (2018), where the spatial compliance ellipsoid was experimentally identified.

Once the stiffness properties of soft robots are comprehended, the subsequent stage is to formulate their statics or dynamics models. Various beam theories, such as the Euler–Bernoulli beam theory (Stilli et al., 2018) and the Timoshenko beam (Lindenroth et al., 2016), have been investigated to establish the forward kinematics of soft robotic manipulators. The utilisation of the Euler–Bernoulli beam theory has also demonstrated that the flexural rigidity of tendon-driven robots is subject to variation based on tendon stiffness, thus influencing robot kinematics (Oliver-Butler et al. 2019). The Piecewise Constant Curvature (PCC) model is a well-known method for statics and dynamics modelling of soft robots (Webster and Jones 2010; Caasenbrood et al., 2023). This model assumes that the curvature along the robot’s length is constant. When calculating the resulting strains from applied moment and torque, the PCC model requires consideration of bending and elongation stiffness (Sadati et al., 2017). Common variable curvature models include the Cosserat/Kirchhoff rod theory (Renda et al., 2018; Grazioso et al., 2019). A geometrically exact Cosserat model was proposed in Trivedi et al. (2008) to describe the forward kinematics of multi-segment continuum robots, considering material nonlinearities. Based on the Cosserat rod theory, the forward kinematics are modelled for a three-segment tendon-driven continuum robot with extensible sections (Amanov et al., 2021). Furthermore, a real-time dynamics modelling framework, based on the dynamic Cosserat rod equations, is applied to a tendon-driven continuum robot, parallel robot, and fluidic-driven soft robot (Till et al. 2019). In addition, a Cosserat model-based method is proposed for a pneumatic-driven continuum robot with braided chambers to then further sense environmental stiffness in real-time (Sadati et al., 2020). A new discretised model and a Matlab package (TMTDyn) was then presented for hybrid rigid-continuum systems based on Euler–Bernoulli beam and Cosserat rod model (Sadati et al., 2021). Based on the material properties, numerical analysis can be conducted, for example, using Finite Element Modelling (Elsayed et al., 2014; Koehler et al., 2019; Zupan and Saje 2003). In summary, the aforementioned modelling techniques are built on known stiffness properties of robots. The material properties (e.g., Young’s modulus), together with bending, elongation and shear stiffness, are usually considered to calculate curvatures and strains from actuation forces. The resulting curvatures and strains complete the statics or dynamics models. However, a thorough stiffness description and analysis for soft robots, for example, the Cartesian stiffness, is not investigated in these models.

The importance of understanding configuration-dependent stiffness has been proven for traditional rigid link and continuum robots (Ajoudani et al. 2015, 2017; Gravagne and Walker 2002; Stella et al., 2023). The difference between stiffness controllability in traditional rigid robots compared with soft robots primarily lies in the generation of compliance. In rigid robots, the compliance is provided by finite variable stiffness joints (Wolf et al., 2015). The stiffness matrix in Cartesian space is determined based on the Jacobian projection (Salisbury 1980; Alici and Shirinzadeh 2005). Instead, many soft robots are considered continuum manipulators without physical joints; the compliance is distributed along the robots (Naselli and Mazzolai 2021). A pioneer work for compliance analysis of continuum robots is to use ellipsoid analysis (Gravagne and Walker 2002), analogous to the techniques used in the rigid-link robots. The stiffness properties of soft robots can be characterised by discretising their structures into finite elements with a diagonal stiffness matrix, using the Cosserat rod model for instance (Renda et al., 2018). With this approach, researchers have proposed a Cartesian impedance control method for planar soft robots composed of multiple segments, assuming a PCC model (Della Santina et al., 2020). Moreover, there are efforts towards a detailed analysis of passive stiffness in soft robots, such as the development of an augmented Jacobian-based method using Ordinary Differential Equations (ODEs) to obtain the compliance matrix (Rucker et al., 2011), and the establishment of statics models based on the Cosserat rod theory. Such a concept was then adopted and developed by Black et al. (2017) to analyse the stiffness response of a parallel continuum robot under different configurations. Inspired by the work of Rucker et al. (2011), a new set of compliance differential equations was developed, directly coupled to the solution of the forward kinematics. These equations were employed to derive the compliance matrix of continuum robots (Smoljkic et al., 2014). However, the robotic manipulators were approximated to the slender Kirchhoff rod, whilst large elongation observed in soft robots was not considered (Gong et al., 2021; Caasenbrood et al., 2023).

As for the specific mathematical tool for the compliance modelling in robotic field, Lie theory (Sonneville et al., 2014) and, in particular, screw theory have been utilised in robotics for it concisely interprets rigid-body motions, such as kinematics and statics modelling (Renda et al., 2017), parameter calibration (Sun et al., 2020; Cibicik and Egeland 2021) and identification (Fu et al. 2020), and compliance analysis (Qi et al., 2016; Selig and Ding 2001; Ding and Dai 2010). In addition, an analytical model was presented in Hussain et al. (2021) to model and analyse the stiffness ellipsoids for a tendon-driven soft-rigid gripper.

1.2. Problem statement and motivation

Soft robots are compliant, generate continuous deformation, and are essentially infinite dimensional systems. As presented in Section 1.1, the stiffness properties of soft robots are usually the foundation to establish statics and dynamics models (see Table 1). Then, understanding and modelling the configuration-dependent stiffness are essential to achieve on-demand stiffness control (Stella et al., 2023). Nevertheless, current stiffness modelling techniques primarily apply finite or partial differentiation (e.g., (Mahvash and Dupont 2011; Black et al. 2017; Rucker et al., 2011)) and map virtual joint stiffness to the task space using the Jacobian (e.g., (Renda et al., 2018; Della Santina et al., 2020)). These calculations might be extensive. For instance, the dimensions of the Jacobian matrix are usually large and depend on the number of virtual joints. Moreover, the description of the compliance distribution along the robot is equivalently important, but less investigated in the current literature. Table 1 further supports the gap in exploring a general stiffness analysis.

Table 1.

Summary of soft robotic modelling and stiffness/compliance analysis.

			Material	Cartesian compliance	Configuration-dependent	Compliance
Key authors	Robot type	Model (model type)	Behaviour	modelling	compliance analysis	distribution
Trivedi et al. (2008)	Pneumatic-drivene	Cosserat rod (statics)	Nonlinear	✘	✘	✘
Rucker et al. (2011)	(Concentric) rodi	Cosserat rod (statics)	Linear	Finite differentiation	✘	✘
Smoljkic et al. (2014)
Black et al. (2017)	Parallel roboti	Cosserat rod (statics)	Linear	Finite differentiation	Ellipsoid	✘
Oliver-Butler et al. (2019)	Tendon-driveni	Cosserat rod (statics)	Linear	✘	✘	✔
Della Santina et al. (2020)	Pneumatic-driveni	PCC (dynamics)	Linear	Jacobian projection	✘	✘
Sadati et al. (2021)	Pneumatic-drivene	Reduced order (dynamics)	Linear	✘	✘	✘
Hussain et al. (2021)	Tendon-driveni	Screw theory (dynamics)	Linear	Jacobian projection	Ellipsoid	✘
Amanov et al. (2021)	Tendon-drivene	Cosserat rod (statics)	Linear	✘	✘	✘
Caasenbrood et al. (2023)	Pneumatic-drivene	PCC (dynamics)	Nonlinear	✘	✘	✘
This work	Pneumatic-drivene	Lie theory (statics)	(Non)linear	Compliance integration	Matrix structure analysis	✔

✘ and ✔ denotes the corresponding modelling or analysis is not considered or considered, respectively; The superscript ^e denotes the backbone of the robot can elongate. The superscript ⁱ denotes the elongation of the backbone is negligible or inextensible.

As such, a systematical modelling approach and analysis of configuration-dependent stiffness and compliance distribution remain challenging. Moreover, for elastomer-based soft robots, the elongation of the robots’ backbone often reduces the cross-sectional area. This behaviour introduces nonlinearities and also influences stiffness properties of these robots. Motivated by the aforementioned considerations, a stiffness modelling and analysis framework should make an important contribution to the soft robotics community.

1.3. Contributions and outline

In our previous work (Shi et al., 2021), we proposed a Cartesian stiffness modelling method based on the PCC model for a single fluidic-driven soft robot. The paper included preliminary validation data on the tip stiffness. Going beyond our previous work, this paper proposes a Lie theory–based stiffness modelling and analysis framework for soft fluidic-driven robotic manipulators to investigate their underlying stiffness characteristics and deliver a configuration-dependent stiffness/compliance analysis. Our framework, illustrated in Figure 1, incorporates forward kinematics and stiffness modelling, considering the nonlinearity introduced by the elongation. Hence, the linear and nonlinear material behaviours are explored, embedded in the PCC model and the Cosserat rod model. To validate our framework, we apply it to a single soft fluidic-driven robotic manipulator and a manipulator composed of two segments. This work provides an approach to analyse the general stiffness based on the commonly used static kinematics models. In particular, the contributions made by this work are Table 2:

1. Our proposed Lie theory–based modelling framework (see Section 2) is capable of configuration-dependent stiffness modelling and analysis at different static robot configurations from forward kinematics models (e.g., the PCC and Cosserat rod models). This framework considers nonlinear responses resulting from large longitudinal deformations (e.g., when elongating or bending) as presented in Sections 2.3 and validated in Sections 5.2.

2. Our methodology demonstrates that the robot compliance can be derived by various integration schemes based on the obtained static configurations of soft robots, without using Jacobian projections or finite differentiation (see Section 2.1). We showcase that our approach allows deriving the compliance matrix by direction integration with forward kinematics models, for example, the proposed Cosserat rod-based compliance model (see Section 2.2.2).

3. Our modelling framework can thoroughly reveal the compliance distribution, and for the first time, we achieve a detailed configuration-dependent compliance/stiffness modelling and analysis for extensible fluidic-driven continuum robots (see Section 3 and Section 4). In Section 5, we validate computational results compared to experiments, using soft manipulators made of one and two robotic segments.

Figure 1.

Diagram of the proposed Lie theory–based framework incorporating forward kinematics, stiffness modelling and analysis. This stiffness modelling and analysis framework accommodates the nonlinear behaviours of kinematics, exemplified by the piecewise constant curvature (PCC) model and the Cosserat rod model, and achieves a configuration-dependent stiffness analysis.

Table 2.

Parameters and notations.

Symbol	Definition
$\widehat{(\cdot )}$ or (⋅)ˆ	Mapping from ${\mathbb{R}}^3$ to $\mathfrak{s}\mathfrak{o}\left(3\right)$ , or from ${\mathbb{R}}^6$ to $\mathfrak{s}\mathfrak{e}\left(3\right)$
(⋅)∨	Inverse operation of $\widehat{(\cdot )}$ or (⋅)ˆ, that is, ${\widehat{(\cdot )}}^{\vee }=(\cdot )$
(⋅) b	A variable expressed in the body (material) frame
(⋅) o	A variable expressed in the global (spatial) frame
{⋅}	Representation of a frame, that is, the global and body frame are {O o , x o , y o , z o }, {O b , x b , y b , z b }
(⋅) (s)	A variable has a function of the arc length variable s
(⋅) i	A variable at the position i along the curve s
(⋅),x|y|z	The x, y or z component of a variable
(⋅) s	∂(⋅)/∂s, derivation of a variable with the arc length s
$\stackrel{\sim }{(\cdot )}$	A variable after elongation
δ(⋅)	An infinitesimal variation of a variable
Δ(⋅)	A variation of a variable
(⋅) j	The jth iteration step in the nonlinear compensation
‖ ⋅‖ F	Frobenius norm of a variable
‖ ⋅‖ 2	Euclidean norm of a variable
T	$\in {\mathbb{R}}^6$ , a twist expressed in the Plücker axis coordinate
W	$\in {\mathbb{R}}^6$ , a wrench expressed in Plücker ray coordinate
S a	$\in {\mathbb{R}}^6$ , a screw vector expressed in the frame {a}
g ab	$\in {\mathbb{R}}^{4\times 4}$ , the homogeneous transformation of the frame {b} expressed in the frame {a}
R ab	$\in {\mathbb{R}}^{3\times 3}$ , the translation vector of the frame {b} expressed in the frame {a}
p ab	$\in {\mathbb{R}}^3$ , the rotation matrix of the frame {b} expressed in the frame {a}
Ad ab	$\in {\mathbb{R}}^{6\times 6}$ , adjoint matrix of the frame {b} expressed in the frame {a}
Δ	$\in {\mathbb{R}}^{6\times 6}$ , elliptical polar operator
c(s) i	$\in {\mathbb{R}}^{6\times 6}$ , compliance density matrix at the position i, c(s) i = diag [c se (s) i , c bt (s) i ]
c se (s) i	$\in {\mathbb{R}}^{3\times 3}$ , compliance density matrix containing shear and elongation, its inverse is k se (s) i
c bt (s) i	$\in {\mathbb{R}}^{3\times 3}$ , compliance density matrix containing bending and torsion, its inverse is k bt (s) i
C(s) i	$\in {\mathbb{R}}^{6\times 6}$ , the Cartesian compliance matrix, with the tip compliance matrix denoted by C
C t (s) i	$\in {\mathbb{R}}^{3\times 3}$ , the translational compliance matrix
K(s) i	$\in {\mathbb{R}}^{6\times 6}$ , the Cartesian stiffness matrix, with the tip stiffness matrix denoted by K
κ(s) i	Curvature of an arc along the central backbone
s(s) i	Length of an arc along the central backbone
ϕ(s) i	Angle between the bending plane of the arc and x-z plane
v(s) i	$\in {\mathbb{R}}^3$ , derivative of the position p oi (s) with respect to s
u(s) i	$\in {\mathbb{R}}^3$ , derivative of angular change with respect to s
n i (s)	$\in {\mathbb{R}}^3$ , internal force along the rod
m i (s)	$\in {\mathbb{R}}^3$ , internal moment along the rod
f e (s) i	$\in {\mathbb{R}}^3$ , distributed external force per unit s
f g (s) i	$\in {\mathbb{R}}^3$ , distributed gravitational force per unit s
l e (s) i	$\in {\mathbb{R}}^3$ , external moment per unit s
f P (s) i	$\in {\mathbb{R}}^3$ , distributed force from pressurisation per unit s
l P (s) i	$\in {\mathbb{R}}^3$ , distributed moment from pressurisation per unit s
Ω	A general actuation variable
P k	Actuation pressure in the ith chamber
λ 1|2|3	The principle material stretch ratios
k d	$\in {\mathbb{R}}^3$ , a position vector of the kth chamber
k φ	The angle of the kth chamber related to the +x-axis
α 1	Angle between two adjacent chamber pairs
α 2	Angle between two adjacent chambers
A(s)	Area of the cross-section
I ,x|y|z (s)	Second moment of area about the x-, y- or z-axes
r m1|m2	Radius of the soft robot and the inner lumen
r o	Radius of the central position of the chamber
r P	Radius of the pneumatic chamber
l 0|t|b	length of the chamber, top and bottom sealing part
p w	$\in {\mathbb{R}}^4$ , augmented variable containing position and stiffness
Δf	$\in {\mathbb{R}}^3$ , variation of the applied force
Δx|y|z Δf	$\in {\mathbb{R}}^3$ , variation of the displacement resulting from the applied force Δ f
ΔE x|y|z	The position differences from different simulations

The remainder of this paper is organised as follows: Section 2 proposes the stiffness modelling framework and elaborates its derivations. Section 3 presents the instantiation of the proposed framework based on pneumatic-driven soft continuum robots. The simulation and model-based stiffness analysis are discussed in Section 4. The experimental validation is then presented in Section 5 to demonstrate the efficacy and explore the accuracy of the proposed framework. Finally, the corresponding discussions and conclusions are presented in Section 6.

2. Theoretical framework of the stiffness modelling and analysis

In line with Figure 1, this section presents the stiffness modelling and analysis framework in detail. Compliance matrices are first formulated using Lie theory based on screws. The configuration-dependent Cartesian compliance is built on the PCC and Cosserat models. This is the fundamental difference compared to current approaches that focus on forward kinematics only (see the dotted red rectangle in Figure 1). Furthermore, our proposed framework considers nonlinear responses resulting from large longitudinal deformations of the incompressible material during elongation and bending. The output of the proposed framework includes the forward kinematics and the general stiffness/compliance matrix, resulting in a configuration-dependent stiffness modelling and analysis.

2.1. Cartesian compliance matrix derivation

The basics of screws and Lie theory are summarised in Appendix B. In general, the relationship between a twist T and a wrench W is described by

$\mathbf{T}=C\mathbf{W},\mathbf{W}=K\mathbf{T}$

(1)

C is the 6 × 6 compliance matrix, and K is the 6 × 6 stiffness matrix with the relationship C = K⁻¹. C and K are interchangeable and both used in this paper.

The soft continuum robot can be discretised into finite elements. So in this paper, we use i to denote the ith discretised element and s to denote the continuous curve along the backbone (see Figure 2). The material compliance density $c{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}$ of the robot at position i, written in the body frame (Figure 2(c)) is

$c{\left(s\right)}_i^b=\left[\begin{array}{ll}{c}_{\mathit{se}}{\left(s\right)}_i^b& 0\\ 0& c_{\mathit{bt}}{\left(s\right)}_i^b\end{array}\right]=\left[\begin{array}{ll}{k}_{\text{se}}^{-1}{\left(s\right)}_i^b& 0\\ 0& k_{\mathit{bt}}^{-\mathit{1}}{\left(s\right)}_i^b\end{array}\right]$

(2)

where s is the arc length along the backbone, and the compliance density matrix $c_{\mathit{se}}{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}=\text{diag}\left[\mathit{1}/\mathit{GA}{\left(\mathit{s}\right)}_{\mathit{i}},\mathit{1}/\mathit{GA}{\left(\mathit{s}\right)}_{\mathit{i}},\mathit{1}/\mathit{EA}{\left(\mathit{s}\right)}_{\mathit{i}}\right]$ contains the shear and elongation items, and $c_{\mathit{bt}}{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}=\text{diag}\left[\mathit{1}/\mathit{E}{\mathit{I}}_{\mathit{i},\mathit{x}}\left(\mathit{s}\right),\mathit{1}/\mathit{E}{\mathit{I}}_{\mathit{i},\mathit{y}}\left(\mathit{s}\right),\mathit{1}/\mathit{G}{\mathit{I}}_{\mathit{i},\mathit{z}}\left(\mathit{s}\right)\right]$ contains the bending and torsion items. A(s) _i is the area of the cross-section, which decreases from A(s) _i to $\tilde{\mathit{A}}{\left(\mathit{s}\right)}_{\mathit{i}}$ when the robot elongates, as shown in Figure 2(a). G is the shear modulus, and E is the Young’s modulus with E = 2G (1 + μ), where μ is the Poisson’s ratio. I _i,x (s), I _i,y (s) and I _i,z (s) are the moment of inertia about the local principle axes, varying with the elongation as well.

Figure 2.

Illustration of the deformed robot with a bending angle θ (when the cross-section is circular). (a) The variation of the cross-section is non-negligible under a large longitudinal deformation for the incompressibility of the material. (b) The soft manipulator comprises serially connected compliant elements. (c) An element of the manipulator.

The material compliance density written in the global frame $c{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{o}}$ is described in

$c{\left(s\right)}_i^o={\mathrm{Ad}}_{oi}^{-T}\left(s\right)c{\left(s\right)}_{\mathrm{i}}^{\mathrm{b}}{\mathrm{Ad}}_{oi}^{-1}\left(s\right)$

(3)

where Ad _oi (s) is the adjoint matrix (see Appendix B), defined by rotation matrix R _oi (s) and displacement vector p _oi (s). As the length of the discrete element at position i is ds _i , the material compliance matrix of the element is $c{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}\mathrm{d}{s}_{\mathit{i}}$ , written in the body frame (Selig and Ding 2009).

From Figure 2(b), the whole manipulator is regarded as the combination of flexible elements (Qiu and Dai 2021; Qi et al., 2016). When a wrench W^o is applied at the tip of a manipulator, the wrench is transmitted to each compliant element. Based on (47), the transmitted wrench $\mathbf{W}{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}$ at the ith element satisfies ${\mathbf{W}}^{\mathbf{o}}={\text{Ad}}_{\mathit{oi}}\left(\mathit{s}\right)\mathbf{W}{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}$ . Based on (51), the transformation of a twist from the body frame to the global frame is $\mathbf{T}{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{o}}={\text{Ad}}_{\mathit{oi}}^{-\mathit{T}}\left(\mathit{s}\right)\mathbf{T}{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}$ . Therefore, the total deformation twist T ^o equals the aggregated twists of each flexible element expressed in the global frame resulting in (4) when combined with (3)

$\begin{array}{l}{\mathbf{T}}^{\mathrm{o}}={\int }_0^l\mathrm{d}\mathbf{T}{\left(s\right)}_i^o={\int }_0^l{\text{Ad}}_{\mathit{oi}}^{-T}\left(s\right)\mathrm{d}\mathbf{T}{\left(s\right)}_i^b\\ ={\int }_0^l{\text{Ad}}_{\mathit{oi}}^{-T}\left(s\right)c{\left(s\right)}_i^b\mathbf{W}{\left(s\right)}_i^b\text{d}\mathit{s}\\ ={\int }_0^l{\text{Ad}}_{\mathit{oi}}^{-T}\left(s\right)c{\left(s\right)}_i^b{\text{Ad}}_{\mathit{oi}}^{-1}\left(s\right){\mathbf{W}}^{\mathbf{o}}\text{d}\mathit{s}\\ =\underset{C^o}{\underbrace{\underset{}{{\int }_0^l{\text{Ad}}_{\mathit{oi}}^{-T}\left(s\right)c{\left(s\right)}_i^b{\text{Ad}}_{\mathit{oi}}^{-\mathit{1}}\left(\mathit{s}\right)\text{d}\mathit{s}}}}{\mathbf{W}}^{\mathrm{o}}=C^o{\mathbf{W}}^{\mathbf{o}},\end{array}$

(4)

where l is the integrated length with the arc length s. Equation (4) depends on the material properties and configuration (p _oi (s), R _oi (s)). Equation (4) also denotes that the compliance matrix C ^o can be obtained by integrating the compliance density expressed in a same global coordinate. Therefore, the material compliance matrix $C{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{o}}$ along the manipulator can be derived by varying the integration length, where the 36 entries of this matrix are integrated individually (Selig and Ding 2009; Qi et al., 2016; Qiu and Dai 2021). The calculation of Jacobian, for example, in Hussain et al. (2021), is not required. The stiffness matrix can then be obtained directly by $K^{\mathit{o}}={\left({\mathit{C}}^{\mathit{o}}\right)}^{-\mathit{1}}$ .

Equation (4) represents the general integration form of the Cartesian compliance based on static robot configurations, illustrating that the compliance is configuration-dependent. Various numerical integration schemes can achieve the compliance integration of (4), once the static configurations of robots are obtained from the forward kinematics. For instance, the midpoint or trapezoidal integration schemes could be applied. In this case, the total Cartesian compliance $C{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{o}}$ up to the position i is the discrete representation of (4), which yields (5)

$C{\left(s\right)}_i^o=\sum _{h=1}^i{\mathrm{Ad}}_{oh}^{-T}\left(s\right)c{\left(s\right)}_h^b{\mathrm{Ad}}_{oh}^{-1}\left(s\right)\mathrm{d}{s}_h$

(5)

The number of i is the discrete element numbers. In Section 2.2.1, we introduce the PCC-based compliance modelling approach, utilising either the midpoint or trapezoidal integration scheme with robot configurations derived from the PCC model. The discrete arc numbers in the PCC model are denoted by i, and high-order integration schemes can also be considered. Additionally, in Section 2.2.2, we present the Cosserat rod-based compliance modelling approach, employing a high-order Runge–Kutta integration scheme. We demonstrate the benefits of performing compliance integration simultaneously with the forward kinematics model in this case.

It is important to note that the selection of the compliance integration scheme is not limited by the statics models. Our approach allows for the exploration of various integration schemes to derive compliance based on the static configurations of robots. Therefore, any integration scheme can be considered and evaluated within our framework.

2.2. Configuration-dependent compliance model

As detailed in Section 2.1, (4) needs to be combined with the forward kinematics models. We complete our modelling framework by proposing the PCC-based and the Cosserat rod-based compliance modelling.

2.2.1. PCC-based compliance model

In the PCC model, the spatial backbone is depicted via a set of serially connected discrete arc whose curvature parameters are [κ(s) _i , s(s) _i , ϕ(s) _i ] of the ith element. κ(s) _i is the curvature, s(s) _i is the curve length and ϕ(s) _i is the angle of the arc that rotates from the x-z plane. Using the exponential mapping from Lie theory, the homogeneous transformation g_i(i+1) from ith arc to (i + 1)th arc using screws (Webster and Jones 2010; Shi et al., 2021). The general form of g_i(i+1) is shown in (43) (see Appendix C).

Substituting all the arc parameters [κ(s) _i , s(s) _i , ϕ(s) _i ] into g_i(i+1), the general forward kinematics g _oi (s) is

$g_{\mathit{oi}}\left(\mathit{s}\right)=\prod _{\mathit{h}=\mathit{1}}^{\mathit{i}}{g}_{\left(\mathit{h}-\mathit{1}\right)\mathit{h}}=\left[\begin{array}{ll}{\mathit{R}}_{\mathit{o}\mathit{i}}\left(\mathit{s}\right)& {\mathit{p}}_{\mathit{o}\mathit{i}}\left(\mathit{s}\right)\\ \mathit{0}& \mathit{1}\end{array}\right].$

(6)

Equation (6) can be used to get the adjoint matrix. The PCC-based compliance modelling is derived by combining (5)–(6). The number of i is the discrete arc numbers in the PCC model and influences the accuracy of the calculated compliance $C{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{o}}$ . It is worth mentioning that the piecewise constant curvature assumption is more applicable when the gravitational force is negligible and no external force is applied.

2.2.2. Cosserat rod-based compliance model

The constant curvature assumption is not needed in the Cosserat rod model. However, the traditional Cosserat model focuses on the statics or dynamics, the compliance/stiffness modelling is not included. A Cosserat rod-based compliance model is proposed in this paper, where the stiffness modelling is directly incorporated, without calculating the static robot configurations in advance.

The mechanics of the rod can then be derived by a set of partial differential equations within the Lie group (Till et al. 2019). By convention, the subscript s means the derivative with arc length s, namely, ∂/∂s. The variable $v{\left(\mathit{s}\right)}_{\mathit{i}}^{\mathit{b}}=R_{\mathit{o}\mathit{i}}^{\mathit{T}}\left(\mathit{s}\right)p_{\mathit{o}\mathit{i}}{\left(\mathit{s}\right)}_{\mathit{s}}$ denotes the first derivative of the position in the body frame, and ${u{\left(s\right)}_i^b=\left(R_{oi}^T\left(s\right)R_{oi}{\left(s\right)}_s\right)}^{\vee }$ is the curvature in the body frame. The differentiation of the internal force is described by n _i (s) _s , and the differentiation of the internal moment is m _i (s) _s . Therefore, the derivative of the position vector p(s) and the rotation matrix R(s) can be written as

$p_{oi}{\left(s\right)}_s=R_{oi}\left(s\right)v{\left(s\right)}_i^b,R_{oi}{\left(s\right)}_s=R_{oi}\left(s\right)\widehat{u}{\left(s\right)}_i^b.$

(7)

The distributed forces and moment along the rod are denoted by f _e (s) _i and l _e (s) _i , respectively. The classic equilibrium of the Cosseart rod satisfies

$\{\begin{array}{l}{\mathit{n}}_{\mathit{i}}{\left(\mathit{s}\right)}_{\mathit{s}}=-{\mathit{f}}_{\mathit{e}}{\left(\mathit{s}\right)}_{\mathit{i}},\\ {\mathit{m}}_{\mathit{i}}{\left(\mathit{s}\right)}_{\mathit{s}}=-{\mathit{p}}_{\mathit{oi}}{\left(\mathit{s}\right)}_{\mathit{s}}\times {\mathit{n}}_{\mathit{i}}\left(\mathit{s}\right)-{\mathit{l}}_{\mathit{e}}{\left(\mathit{s}\right)}_{\mathit{i}}.\end{array}$

(8)

The linear constitutive equation relates u(s) _i , v(s) _i to the loading n _i (s), m _i (s), and the linear relationship yields

$\{\begin{array}{ll}{n}_i\left(s\right)& =R_{oi}\left(s\right)k_{se}{\left(s\right)}_i^b\left[v{\left(s\right)}_i^b-v^{\ast }{\left(s\right)}_i^b\right],\\ m_i\left(s\right)& =R_{oi}\left(s\right)k_{bt}{\left(s\right)}_i^b\left[u{\left(s\right)}_i^b-u^{\ast }{\left(s\right)}_i^b\right].\end{array}$

(9)

The variables $v^{\ast }{\left(s\right)}_i^b$ and $u^{\ast }{\left(s\right)}_i^b$ denote the initial strains and curvatures of a rod. In particular, the value of a straight rod has $v^{\ast }{\left(s\right)}_i^b=e_3={\left[\mathrm{0,0,1}\right]}^T$ and $u^{\ast }{\left(s\right)}_i^b={\left[\mathrm{0,0,0}\right]}^T$ . $k_{se}{\left(s\right)}_i^b$ and $k_{bt}{\left(s\right)}_i^b$ are defined in (2).

Recalling (4), the derivative of the compliance matrix, namely the compliance density, can be obtained by

${\left[C{\left(s\right)}_i^o\right]}_s={\mathrm{Ad}}_{oi}^{-T}\left(s\right)c{\left(s\right)}_i^b{\mathrm{Ad}}_{oi}^{-1}\left(s\right).$

(10)

Finally, combining (7) to (10) and choosing the state variables y as ${\left[p,R,m,n,C^o\right]}_{1\times 54}$ , the full set of ODEs of the Cosserat rod-based compliance model can be summarised as

$\{\begin{array}{l}{p}_{oi}{\left(s\right)}_s=R_{oi}\left(s\right)v{\left(s\right)}_i^b,\begin{array}{ll}v{\left(s\right)}_i^b& =c_{se}{\left(s\right)}_i^b{R}_{oi}^T\left(s\right)n_i\left(s\right)+v^{*}{\left(s\right)}_i^b\end{array}\\ R_{oi}{\left(s\right)}_s=R_{oi}\left(s\right)\hat{u}{\left(s\right)}_i^b,\begin{array}{ll}u{\left(s\right)}_i^b& =c_{bt}{\left(s\right)}_i^b{R}_{oi}^T\left(s\right)m_i\left(s\right)+u^{*}{\left(s\right)}_i^b\end{array}\\ n_i{\left(s\right)}_s=-f_e{\left(s\right)}_i\\ m_i{\left(s\right)}_s=-p_{oi}{\left(s\right)}_s\times n_i\left(s\right)-l_e{\left(s\right)}_i\\ {\left[C{\left(s\right)}_i^o\right]}_s={\mathrm{Ad}}_{oi}^{-T}\left(s\right)c{\left(s\right)}_i^b{\mathrm{Ad}}_{oi}^{-1}\left(s\right),c{\left(s\right)}_i^b=\mathrm{diag}{\left[c_{se}\left(s\right),c_{bt}\left(s\right)\right]}_i^b\end{array}$

(11)

In general, modelling the compliance and static robot configurations can be treated as separate entities once the configurations are derived. However, utilising (11), we illustrate that the compliance distribution $C{\left(s\right)}_i^o$ along the robot can be directly obtained by solving the ODEs of the static Cosserat rod within the same integration loop. This integrated approach allows us to efficiently compute the compliance while considering the robot’s static configurations simultaneously without increasing the complexity of solving the ODEs of the Cosserat rod.

Numerical integration methods, such as a high-order Runge–Kutta approach, can be used to solve (11). It is noteworthy that if the integration truncation error of R _oi (s) results in R _oi (s) ∉ SO(3), R _oi (s) might be represented by quaternion when doing the integration (Rucker 2018). If the initial conditions are known, (11) could be integrated directly as an initial value problem. When the initial states are unknown, (11) can be solved as a boundary value problem using the shooting method (Florian and Damien 2020). Equation (11) gives the analytical representation of the distributed compliance based on the conventional static Cosserat rod equations, where the compliance distribution is not included (Renda et al., 2018; Till et al. 2019). Moreover, the compliance matrix in the proposed Cosserat rod-based compliance model is derived directly without the prerequisite of the Jacobian or finite difference approximations (Rucker et al., 2011; Smoljkic et al., 2014).

2.3. Nonlinear responses compensation

As introduced in Section 2.1, the fluidic-driven robot could have large elongation under high pressurisation (Figure 2(a)). This elongation introduces nonlinear kinematic responses (Gazzola et al., 2018; Caasenbrood et al., 2023). To account for this nonlinearity, we introduce a compensation method in this section, inspired by the nonlinear matrix structural analysis (NMSA) method presented in Naselli and Mazzolai (2021). The NMSA method addresses large deformations in soft bodies by applying load in predefined steps. This method enables the use of a linear strain–stress relation to consider large deformations. Therefore, our proposed framework can generally be applied when an overall strain–stress assumption is less effective, without changing the fundamental models in Section 2.1–2.2, such as the method for constructing the compliance density matrix (see (2)). The load/actuation is applied incrementally, and even with a linear material model, the computation is nonlinear because the system parameters, such as the stiffness matrix density, are updated based on the previous step’s results. The three procedures are:

2.3.1. Procedure 1

The actuation is applied in an incremental way ΔΩ, which is either by predefined n steps or an absolute threshold value b_Ω. Using predefined steps, the incremental actuation variable ΔΩ can be described by

$\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Omega }/n.$

(12)

Ω is the maximum actuation variable, which could be the fluidic pressure, fluid volume, or external loads in fluidic-driven robots. System parameters update in each step and ΔΩ is used in every step instead of Ω. In the threshold method, ΔΩ equals

$\mathrm{\Delta }\mathrm{\Omega }=\left\{\mathrm{\Delta }\mathrm{\Omega }\left(t\right):{\Vert \mathrm{\Delta }\mathrm{\Omega }\left(t\right)\Vert }_2\le b_{\mathrm{\Omega }}\right\}.$

(13)

Equation (13) denotes that the system parameters update when the variation of the actuation variable in the 2-norm is larger than b_Ω. ΔΩ(t) is defined as the variation of Ω after the current update starts, which is reset to zero when the next update step begins. Equation (12) can be used if the maximum value of actuation is known, for example, a step input. By contrast, (13) is more convenient if the actuation changes with time, for instance, a ramp or a sine wave input.

2.3.2. Procedure 2

From (1), the increment of the twist at position i, in the jth step in body frame can be described by

${\left(\mathrm{\Delta }{\mathbf{T}}_i^b\left(s\right)\right)}^j={\left[c{\left(s\right)}_i^b\mathrm{\Delta }{\mathbf{W}}_i^b\left(s\right)\right]}^j.$

(14)

where, ${\left(\mathrm{\Delta }{\mathbf{W}}_i^b\right)}^j$ is the wrench at the position i in the jth iteration step. The new deformation twist ${\left(\mathrm{\Delta }{\mathbf{T}}_i^b\left(s\right)\right)}^j$ (e.g., v(s) and u(s)) can be used to construct the variation of the system variables Δy ^j , to update the new state variables y ^j as follows

${\mathbf{y}}^j={\mathbf{y}}^{j-1}+\mathrm{\Delta }{\mathbf{y}}^j,$

(15)

where y ^j is used to model the forward kinematics of robots.

2.3.3. Procedure 3

At the end of the jth step, the compliance matrix density (2) updates based on the results from the (j − 1)th step. When the material stretches, the principle stretches are defined as [λ₁, λ₂, λ₃]. A practical way to update the compliance matrix is treating the material as incompressible (λ₁λ₂λ₃ = 1) and undergoing uniaxial extension. This yields

${\lambda }_2={\lambda }_3=\sqrt{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\lambda }_1.$}\right.}$

(16)

Specifically, when the numerical models including the PCC-based and Cosserat rod-based compliance models build on a linear soft material behaviour, we refer to as the linear Cosserat model (LCM) and linear PCC model (LPCC). Models considering the nonlinear compensation of the soft material are referred to as the nonlinear Cosserat model (NCM) and nonlinear PCC model (NPCC) in this paper.

2.4. Matrix-based stiffness analysis

This proposed method could be used to reveal how stiffness distributes along the robot. The 6 × 6 material compliance or stiffness matrix has the structure of

$C=\left[\begin{array}{ll}{C}_{fv}& C_{mv}\\ C_{f\omega }& C_{m\omega }\end{array}\right]={\left[\begin{array}{ll}{K}_{fv}& K_{mv}\\ K_{f\omega }& K_{m\omega }\end{array}\right]}^{-1}=K^{-1}.$

(17)

The compliance matrix C is symmetrical and consists of four sub-matrices, where C _fv is the force-translational matrix, C _mω is the moment-rotational matrix, C _mv , C _fω are the coupling matrices. The corresponding stiffness matrices are K _fv , K _mω , K _mv , and K _fω .

The compliance matrix $C{\left(s\right)}_i^o$ derived from (5) or (11) is written in the global frame {O _o , x _o , y _o , z _o }, which can be expressed in the body frame {O _i , x _i , y _i , z _i } by

$C{\left(s\right)}_i^b={\text{Ad}}_{oi}{\left(s\right)}^{T}C{\left(s\right)}_i^o{\text{Ad}}_{oi}\left(s\right).$

(18)

Based on a derived Cartesian compliance matrix, that is, C ^o (or C ^b ) at the robot tip, the deflected tip configuration $g_d^o$ of a robot undergoing an applied wrench ΔW ^o (or ΔW ^b ) can be determined. Assuming the current tip configuration is $g_c^o$ (defined by (46)), the new configuration can be calculated from the resulting twists via exponential mapping in (19)

$g_d^o=e^{\left(\mathrm{\Delta }{\mathbf{T}}^o\right)^}{g}_c^o=g_c^o{e}^{\left(\mathrm{\Delta }{\mathbf{T}}^b\right)^}.$

(19)

where ΔT^o = C ^o ΔW ^o and ΔT^b = C ^b ΔW ^b . For instance, the defected position is the difference between the two 3 × 1 translation vectors from $g_d^o$ and $g_c^o$ .

2.4.1. Stiffness/Compliance ellipsoid

Stiffness ellipsoid (SE) or compliance ellipsoid (CE) is the projection between the twist and wrench, where the input twist/wrench is restricted to a hyper-sphere, and the resultant wrench/twist will be a hyper-ellipsoid (Gravagne and Walker 2002). The SE and CE is described by

$S\mathrm{E}=\left\{\delta \mathbf{W}:\delta {\mathbf{T}}^T\delta \mathbf{T}=1\right\},\mathrm{C}\mathrm{E}=\left\{\delta \mathbf{T}:\delta {\mathbf{W}}^T\delta \mathbf{W}=1\right\}.$

(20)

δW and δT are the infinitesimal change of a wrench or twist. Equation (20) also means $\delta {\mathbf{T}}^T{\left(C{C}^T\right)}^{-1}\delta \mathbf{T}=1,\delta {\mathbf{W}}^T{\left(K{K}^T\right)}^{-1}\delta \mathbf{W}=1$ . For CE, there are six principle axis vectors β _i , (i = 1…6) with magnitudes $1/\sqrt{{\gamma }_i},\left(i=1\dots 6\right)$ . β _i and γ _i are the eigenvectors and the eigenvalues of CC ^T . The stiffness matrix has the same principle axis vectors β _i with the compliance matrix but with inverse magnitudes $\sqrt{{\gamma }_i}$ . The dimension of a full stiffness or compliance ellipsoid is $\in {\mathbb{R}}^{6\times 6}$ . A typical way of visualising the compliance ellipsoids in 3D space is by subtracting the 3 × 3 translation and rotation sub-matrices from the 6 × 6 compliance matrix (Black et al. 2017; Hussain et al., 2021). Once the relation between the applied force (or moment) and resulted displacements in the Cartesian space is obtained, the ellipsoids provide a visualisation of the compliance properties.

2.4.2. Eigenscrew decomposition

The ellipsoid is essentially a graphic representation of direct eigenvalue decomposition of the related stiffness or compliance matrix, but this decomposition is dependent on the choices of coordinates. To achieve a coordinate-free analysis, eigenscrew decomposition of (17) can be adopted. The eigenvalue of such decomposition is called eigenstiffness and the corresponding eigenvector is called eigenscrew (Huang and Schimmels 2000). Rewriting (18) in stiffness matrix form and multiplying by the elliptical polar operator Δ yields

$\begin{array}{ll}{K}^0\mathbf{\Delta }& ={\text{Ad}}_{ob}{K}^b{\text{Ad}}_{ob}^T\mathbf{\Delta }={\text{Ad}}_{ob}{K}^b\mathbf{\Delta }{\text{Ad}}_{ob}^{-1}={\mathrm{Ad}}_{ob}\left(K^b\mathbf{\Delta }\right){\text{Ad}}_{ob}^{-1},\end{array}$

(21)

where ${\text{Ad}}_{ab}^T\mathbf{\Delta }=\mathbf{\Delta }{\text{Ad}}_{ab}^{-1}$ . This means that K ^o Δ and K ^b Δ have the same eigenvalues. Thus, KΔ preserves eigenvalues under coordinate transformations. As such, the eigenscrew problem is formulated as the eigenvalue decomposition of (KΔ)ζ _i = η _i ζ _i , (i = 1…6), where η _i is the eigenstiffness and ζ _i is the eigenscrew. η ₁ ∼ η ₆ have three positive values and three negative values.

To describe the configuration-dependent stiffness, the augmented position quantity p _w is proposed as

$p_w=\left[x,y,z,\tau \right],$

(22)

where p _w contains the position and stiffness information. x, y, z denote the position, and τ is the criteria to describe the stiffness, for example, τ can be the derived eigenstiffness η from the eigenscrew decomposition.

In this way, the frame-invariant eigenscrew decomposition can be utilised to analyse the stiffness property of the soft robot. Similarly, the counterpart, eigencompliance decomposition, can also be obtained (Dai and Ding, 2005). For more information, details of the eigenscrew decomposition can be found in Huang and Schimmels (2000); Chen et al. (2015). In summary, the general scheme of the proposed configuration-dependent stiffness/compliance modelling and analysis framework is summarised in Figure 1.

3. Case study: Pneumatic-driven soft robot

3.1. Continuum robot prototype

The soft robot used in this paper follows the large-scale STIFF-FLOP design, which is devised for minimally invasive surgery and made of a highly deformable elastomer with a fluidic-driven principle (Fraś et al., 2015). The robot has six individually reinforced circular chambers. Braided in-extensible threads constrain the radial inflation while allowing longitudinal expansion under pressurisation. In this way, the pressurisation does not change the shape and perimeter of all actuation chambers and the strain in the circumferential direction of the reinforcement layer is negligible (Polygerinos et al., 2015).

The body of an individual cylindrical segment is composed of silicone material [Ecoflex 00-50 Supersoft, SmoothOn] with an overall length of 5.0 cm, and an outer radius of 1.25 cm. The top and bottom of the robot are sealed using rubber material [Dragon Skin 30]. A lumen with a 0.44 cm radius is reserved for feeding through a medical device, such as a camera. The length of the chamber is 4 cm, and its internal radius is 0.225 cm. The adjacent two chambers are paired together. The three chamber pairs are evenly located in the periphery of the silicone cylinder with 120° interval. Actuating one or two chamber pairs results in omni-directional bending, while actuating all three chamber pairs simultaneously generates pure elongation. The detailed dimensions can be found in Figure 3.

Figure 3.

Parameters of the soft robot. l₀, l _t , and l _b are the length of the chamber, the top and the bottom cap of the manipulator.

Following the methodology proposed in Section 2, the specific implementation on the pneumatic-driven manipulator is detailed in the following sections. Section 3.2 derives the compliance density and the updates, and Section 3.3 describes the model implementation.

3.2. Compliance matrix derivation and update

The vector P contains the pressure in six chambers, and P₁ = P₂, P₃ = P₄, P₅ = P₆ because the adjacent two chambers are actuated as one chamber pair. ${}^{k}d\in {\mathbb{R}}^3$ is the position vector of the kth chamber written in the body frame, k = 1⋯6, described by

$d_k=R_z\left({\phi }_{.}^k\right){\left[r_o,\mathrm{0,0}\right]}^T,R_z\left({\phi }_k\right)=\left[\begin{array}{lll}{\cos }^k\phi & -{\sin }^k\phi & 0\\ {\sin }^k\phi & {\cos }^k\phi & 0\\ 0& 0& 1\end{array}\right],$

(23)

where ^k φ is the angle of the kth chamber related to the positive direction of the x-axis, which can be calculated by α₁ and α₂. α₁ is the angle between two chambers in two adjacent chamber pairs, and α₂ is the angle between two adjacent chambers of the same pair, shown in Figure 3.

In line with (16), the element compliance distribution $c{\left(s\right)}_i^b$ can be built on (2) considering the stretch ratio λ₁ via

$\{\begin{array}{l}A{\left(s\right)}_i=\pi \left(\frac{r_{m1}^2-r_{m2}^2}{{\lambda }_1}-6r_p^2\right)\\ I_{,x|y}{\left(s\right)}_i=\frac{\pi r_{m1}^4-\pi r_{m2}^4}{4{\lambda }_1^2}-{\displaystyle \sum }_{k=1}^6\left(\frac{\pi r_p^4}{4}+\frac{{\left({}^k{d}_{,x|y}\right)}^2\pi r_p^2}{{\lambda }_1}\right)\\ I_{,z}{\left(s\right)}_i=\frac{2\left(\pi r_{m1}^4-\pi r_{m2}^4\right)}{4{\lambda }_1^2}-6\left(\pi \frac{r_p^4}{4}+\frac{r_0^2\pi r_p^2}{{\lambda }_1}\right).\end{array}$

(24)

r_m1, r_m2, r _P are the radii of outer, inner manipulator and pressure chamber, which are shown in Figure 3, the subscript x|y denotes around the x- or y-axis. The braided constraint of the chambers can be described as ${\tilde{r}}_P=r_P\sqrt{1-{\lambda }_1^2\cos {\left(\gamma \right)}^2}/\sin \gamma$ (Sadati et al., 2020). When the thread is densely braided, γ ≈ π/2, r _P could be regarded as constant. Other radii would follow the updates as shown in (24).

According to Section 2.3, the system parameters from (24) in the (j + 1)th step can be updated via ${\lambda }_1^j$ derived from the jth step. This can be achieved by substituting $r_{m1}^j,r_{m2}^j,r_1^j,r_0^j,{\left({}_{.}{}^{k}d\right)}^j$ into (24) to obtain $A{\left(s\right)}_i^{j+1},I_{,x|y}{\left(s\right)}_i^{j+1},I_{,z}{\left(s\right)}_i^{j+1}$ .

The actuation variable Ω here is the pressure P in three chamber pairs. According to (12) and (13), the predefined-step method yields

$\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Delta }P=\left[P_1/n,P_3/n,P_5/n\right].$

(25)

The threshold-based method results in

$\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Delta }P=\left[\mathrm{\Delta }{P}_1\left(t\right),\mathrm{\Delta }{P}_3\left(t\right),\mathrm{\Delta }{P}_5\left(t\right)\right].$

(26)

ΔΩ is replaced by ΔP in following sections.

3.3. Instantiation of the configuration-dependent compliance models

The derivation of the compliance matrix (2) and the corresponding update ((24)) are the same in the PCC and Cosserat rod model. In order to derive the forward kinematics in each step, the actuation wrench needs to be calculated first. In each iteration, the force vector F _P written in the body frame equals ${\sum }_{k=1}^6{F}_{P_k}={\sum }_{k=1}^6\left(P_k{A}_P{e}_3\right)$ , where $A_P=\pi r_P^2$ , P _k is the pressure in each chamber, $F_{P_k}$ is the force in each chamber. The totally moment T _P vector written in the body frame equals ${\sum }_{k=1}^6{T}_{P_k}={\sum }_{k=1}^6{(}^{k}d\times F_{P_k})$ .

3.3.1. PCC-based compliance model

As discussed in Section 2.2.1, the set of curvature parameters [κ(s) _i , s(s) _i , ϕ(s) _i ] is used to derive the forward kinematics from the ith to the (i + 1)th element. The length of arc s(s) _i in each iteration is

$s{\left(s\right)}_i={\left(s_0\left(s\right)\right)}_i+\frac{F_{P,z}{\left(s_0\left(s\right)\right)}_i}{EA{\left(s\right)}_i}={\left({\lambda }_1\right)}_i{\left(s_0\left(s\right)\right)}_i.$

(27)

where s ₀ (s) _i is the initial arc length. Curvature κ(s) _i yields

$\begin{array}{l}\kappa {\left(s\right)}_i=\sqrt{{\kappa }_{,x}{\left(s\right)}_i^2+{\kappa }_{,y}{\left(s\right)}_i^2},\\ {\kappa }_{,x}{\left(s\right)}_i=\frac{T_{P,x}}{E{I}_{,x}{\left(s\right)}_i},{\kappa }_{,y}{\left(s\right)}_i=\frac{T_{P,y}}{E{I}_{,y}{\left(s\right)}_i}.\end{array}$

(28)

T _P,x and T _P,y are the components of moments in the x- and y-axes. The angle ϕ(s) _i is calculated by

$\varphi {\left(s\right)}_i=-\mathrm{atan}2\hspace{0.17em}\left({\kappa }_{,x}{\left(s\right)}_i,{\kappa }_{,y}{\left(s\right)}_i\right).$

(29)

Substituting (27)–(29) to (6)–(5), the PCC-based compliance model in one iteration (i.e., the LPCC model) can be realised.

3.3.1.1. Update principle

Equations (27)–(29) are within one iteration, where the iteration number j is dropped. To update the parameters, the arc parameters in the (j + 1)th step can be treated as the superposition on the jth step according to (15). The initial value of the arc length is updated as $s{\left(s\right)}_i^{j+1}={\left({\lambda }_1\right)}_i^{j+1}s{\left(s\right)}_i^j$ . Specially, updates of the arc parameters $\left[\kappa {\left(s\right)}_i^{j+1},s{\left(s\right)}_i^{j+1},\varphi {\left(s\right)}_i^{j+1}\right]$ from the jth step to the (j + 1)th step are detailed as

$$\begin{array}{l}s{\left(s\right)}_i^{j+1}=s{\left(s\right)}_i^j+\frac{F_{\mathrm{\Delta }\mathit{P}}^{j+1}s{\left(s\right)}_i^j}{EA{\left(s\right)}_i^{j+1}},\\ k_{,x}{\left(s\right)}_i^{j+1}=k_{,x}{\left(s\right)}_i^j+\frac{T_{\mathrm{\Delta }\mathit{P},x}^{j+1}}{{EI}_{,x}{\left(s\right)}_i^{j+1}},\\ k_{,y}{\left(s\right)}_i^{j+1}=k_{,y}{\left(s\right)}_i^j+\frac{T_{\mathrm{\Delta }\mathit{P},\mathit{y}}^{j+1}}{{EI}_{,y}{\left(s\right)}_i^{j+1}},\\ \varphi {\left(s\right)}_i^{j+1}=-\mathrm{atan}2(k_{,x}{\left(s\right)}_i^{j+1},k_{,y}{\left(s\right)}_i^{j+1}).\end{array}$$

(30)

ΔP is the discretised actuation value from (25) or (26). The stretch ratio λ₁ in the jth step can be derived and substituted into (24) to update the system parameters. This allows the compliance matrix to be updated in each step.

3.3.2. Cosserat rod-based compliance model

Considering the gravitational force f _g (s) _i = ρA(s) _i g, the force equilibrium in a discretised element σ at the position i is

$\begin{array}{l}{\int }_{i\sigma }^{\left(i+1\right)\sigma }{f}_e{\left(i\sigma \right)}_{i}d\sigma +n_{i+1}\left(i\sigma +\sigma \right)-n_i\left(i\sigma \right)-\\ {\displaystyle \sum }_{k=1}^6\{P_k{A}_P[R_{o\left(i+1\right)}\left(i\sigma +\sigma \right)-R_{oi}\left(i\sigma \right)]e_3\}=0,\end{array}$

(31)

where σ is the discretised length of the arc. The first order Taylor expansion yields n_i+1 (iσ + σ) = n _i (iσ) + n _i (iσ) _s dσ and R_o(i+1) (iσ + σ) = R _oi (iσ) + R _oi (iσ) _s dσ. Equation (31) can then be differentiated with σ and replacing iσ with the arc parameter s, which yields

$\begin{array}{l}{n}_i{\left(s\right)}_s=-f_e{\left(s\right)}_i+{\displaystyle \sum }_{k=1}^6\left[P_k{A}_P{R}_{oi}{\left(s\right)}_s{e}_3\right]\\ =-f_g{\left(s\right)}_i+f_P{\left(s\right)}_i,\end{array}$

(32)

f _e (s) _i is the distributed external force and equals to f _g (s) _i here, and f _P (s) _i is the distributed force resulting from pressurisation. Similarly, the moment balance can be described in (33)

$$\begin{array}{l}{m}_{i+1}\left(i\sigma +\sigma \right)-m_i\left(i\sigma \right)+p_{o\left(i+1\right)}\left(i\sigma +\sigma \right)\times n_{i+1}\left(i\sigma +\sigma \right)-p_{oi}\left(i\sigma \right)\times n_i\left(i\sigma \right)+{{\displaystyle \int }}_{i\sigma }^{i\sigma +\sigma }\left[l_e{\left(i\sigma \right)}_i+p_{oi}\left(i\sigma \right)\times f_e{\left(i\sigma \right)}_i\right]d\sigma \\ -A_P{\displaystyle \sum }_{k=1}^6{P}_k\left[p_{o\left(i+1\right)}\left(i\sigma +\sigma \right)+R_{o\left(i+1\right)}{\left(i\sigma +\sigma \right)}^{k}d\right]\times \left[R_{o\left(i+1\right)}\left(i\sigma +\sigma \right)e_3\right]-\left[p_{oi}\left(i\sigma \right)+R_{oi}{\left(i\sigma \right)}^{k}d\right]\times \left[R_{oi}\left(i\sigma \right)e_3\right]=0.\end{array}$$

(33)

Similar to (32), (33) can be differentiated with σ resulting in

$$\begin{array}{l}{m}_i{\left(s\right)}_s=-l_e{\left(s\right)}_i-{\widehat{p}}_{oi}{\left(s\right)}_{s}n{\left(s\right)}_i+A_P{\displaystyle \sum }_{k=1}^6{P}_k.\\ R_{oi}\left(s\right)\left[\left(v{\left(s\right)}_i^b+u{\left(s\right)}_i^b\times {}^{k}d\right)\times e_3+{}^k\widehat{d}\left(\widehat{u}{\left(s\right)}_i^b{e}_3\right)\right]\\ =-l_e{\left(s\right)}_i-{\widehat{p}}_{oi}{\left(s\right)}_s{n}_i\left(s\right)+l_p{\left(s\right)}_i.\end{array}$$

(34)

The external distributed moment l_e(s) _i is equal to zero and l _P (s) _i is the distributed moment resulting from the pressurisation.

The wrench at the chamber end is balanced by the internal integrated loading [n(l), m(l)], and the exerted external wrench [n _e (l), m _e (l)] including the exerted load and the weight of the tip appendage, which includes the sensor and plastic sealing plate, with a total weight of 8 g in this prototype. The boundary conditions can be described by

$\begin{array}{l}n\left(l\right)=n_e\left(l\right)+R_{ot}\left(l\right)\sum _{k=1}^6{A}_P{P}_k{e}_3,\\ m\left(l\right)=m_e\left(l\right)+R_{ot}\left(l\right){\sum _{k=1}^6}^{}{}^{k}d{A}_P{P}_k{e}_3,\end{array}$

(35)

where R _ot (l) is the rotation matrix at the tip position. In summary, the whole set of ODEs is

$\{\begin{array}{l}{p}_{oi}{\left(s\right)}_s=R_{oi}\left(s\right)v{\left(s\right)}_i^b\\ R_{oi}{\left(s\right)}_s=R_{oi}\left(s\right)\hat{u}{\left(s\right)}_i^b\\ n_i{\left(s\right)}_s=-f_g{\left(s\right)}_i+f_P{\left(s\right)}_i\\ m_i{\left(s\right)}_s=-{\hat{p}}_{oi}\left(s\right)n_i\left(s\right)+l_P{\left(s\right)}_i\\ {\left[C{\left(s\right)}_i^o\right]}_s={\mathrm{Ad}}_{oi}^{-T}\left(s\right)c{\left(s\right)}_i^b{\mathrm{Ad}}_{oi}^{-1}\left(s\right).\end{array}$

(36)

The problem can be solved by (36) with boundary condition (35) (Till et al. 2019). In this paper, (36) is solved using the Fourth-order Runge–Kutta Matlab solver. In this way, the orthonormality of the rotation matrix is maintained during the numerical integration (details are reported in Appendix D). The shooting method is implemented via the embedded optimisation Matlab solver, fsolve (), with the Levenberg–Marquardt method.

3.3.2.1. Update principle

The curvature ${\left[u{\left(s\right)}_i^b\right]}^j$ and strain ${\left[v{\left(s\right)}_i^b\right]}^j$ superposes upon the previous value ${\left[u{\left(s\right)}_i^b\right]}^{j-1}$ and ${\left[v{\left(s\right)}_i^b\right]}^{j-1}$ , so (15) is implemented and updated by

$\begin{array}{l}{\left[v{\left(s\right)}_i^b\right]}^j={\left[c_{se}{\left(s\right)}_i^b{R}_{oi}^T\left(\mathrm{s}\right)n_i\left(s\right)\right]}^j+e_3,\\ {\left[u{\left(s\right)}_i^b\right]}^j={\left[c_{bt}{\left(s\right)}_i^b{R}_{oi}^T\left(s\right)m_i\left(s\right)\right]}^j+{\left[u{\left(s\right)}_i^b\right]}^{j-1},\end{array}$

(37)

with the integration length ${\lambda }_1^j$ updated in each step. Finally, the gravitational force f _g is treated as an individual external load on the effect of pressurisation. In this case, f _g is considered in the last step. The overall implementation procedure of the nonlinear Cosserat model is summarised in Figure 4.

Figure 4.

Implementation scheme of the nonlinear Cosserat method under large elongation. The grey block means it only executes once. The inner loop is to solve the boundary value problem of the set of ODEs in (36). And the outer loop is to deal with the hyper-elasticity resulting from the large elongation and update parameters. The implementation logic is the same for the nonlinear PCC model, by replacing the inner loop (blue rectangle) by the PCC model derived in Section 3.3.1.

3.3.2.2. Modelling a robot made of two segments

For a soft manipulator made of two robotic segments (see Figure 5), our compliance modelling approach will need to consider intermediate conditions. For the second segment, (31)—(37) can be applied by replacing P₁ ∼ P₆ with P₇ ∼ P₁₂. The connection link is considered as rigid. Firstly, the kinematic continuity between the connection link and two segments (i.e., when s = l₁ and s = l₁ + l _k ) need to be satisfied which yields (38)

$\begin{array}{l}{p}_{o{l}_1}\left(l_1^{-}\right)=p_{o{l}_1}\left(l_1^{+}\right),R_{o{l}_1}\left(l_1^{-}\right)=R_{o{l}_1}\left(l_1^{+}\right),\\ p_{o\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{-})=p_{o\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{+}),\\ R_{o\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{-})=R_{o\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{+}).\end{array}$

(38)

where l₁ and l₂ are the backbone lengths for the first and second robotic segment, respectively. The superscripts ⁻ and ⁺ denote the left and right limits (see Figure 5). In addition, the force and moment conditions between the link the two segments are presented in (39)

$\begin{array}{l}{n}_{l_1}\left(l_1^{-}\right)=R_{o{l}_1}\left(l_1^{-}\right)F_P+n_{l_1}\left(l_1^{+}\right),P_k\in \left[P_1,P_6\right],\\ m_{l_1}\left(l_1^{-}\right)=R_{o{l}_1}\left(l_1^{-}\right)T_P+m_{l_1}\left(l_1^{+}\right),P_k\in \left[P_1,P_6\right],\\ n_{\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{-})=n_{\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{+})\\ -R_{o\left(l_1+l_k\right)}\left(l_1+l_k\right)F_P,P_k\in \left[P_7,P_{12}\right],\\ m_{\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{-})=m_{\left(l_1+l_k\right)}{(\left(l_1+l_k\right)}^{+})\\ -R_{o\left(l_1+l_k\right)}\left(l_1+l_k\right)T_P,P_k\in \left[P_7,P_{12}\right],\end{array}$

(39)

where F _P and T _P are the summed force and moment generated from the actuation pressure and written in the body frame (see Section 3.3). The boundary conditions at the tip of the second robotic segment are the same as in (35), with P_k ∈ [P₇, P₁₂] and s = l₁ + l _k + l₂. The resulting ODEs can be solved using the shooting method.

Figure 5.

Illustration of a robot with two connected segments in series. The dimensions of each segment are listed in Figure 3. The length/thickness of the connection part l _k is 6 mm. The chamber pressures of the second segment are denoted from P₇ to P₁₂.

4. Simulation results of the forward kinematics and compliance/stiffness modelling

In this section, computational results are presented based on the established models in Section 3. Simulation 1 presents a comparison of the four forward kinematics models: the linear Cosserat model (LCM), linear PCC model (LPCC), nonlinear Cosserat model (NCM) and nonlinear PCC model (NPCC), and investigates the derived results of the 6 × 6 compliance matrix. The discretised element number in the PCC and Cosserat rod models is chosen as 50 for the interest of compliance integration. Simulation 2 analyses the compliance using the compliance ellipsoid. Simulation 3 presents the configuration-dependent stiffness by the eigenscrew decomposition. In the simulations, the gravity direction is along the −z-axis. In the following sections, we use the terminologies in-plane and out-of-plane bending to denote when the backbone of manipulators bends within a single plane and through different planes, respectively.

4.1. Simulation 1—model implementation and comparison

The first set of simulations compares the LPCC, NPCC, LCM and NCM models. Pressure inputs were in the range of 0–1.5 bar, resulting in the manipulator to bend towards the −x (P₁ = P₂) or +x (P₃ = P₄ = P₅ = P₆) direction. The simulated results when P₁ = P₂ are illustrated in Figure 6(a), showing that the four models have different kinematic responses. In general, the differences are less for low pressure values (e.g., 0.6 bar). Figure 6(b) reports on the results when two chamber pairs are actuated simultaneously, showing that all four models return similar output when the pressure values are less than 0.3 bar. The NCM and NPCC models have larger bending motions compared to the LCM and LPCC, in particular, when the pressure reaches 1.5 bar. This implies that the consideration of the cross-sectional deformation has a significant influence on the forward kinematics calculations. When comparing the PCC and Cosserat models, it is observed that the Cosserat models result in larger bending angles.

Figure 6.

Results of Simulation 1 - The bending comparison of the four models when (a) P₁ = P₂ and (b) P₃ = P₄ = P₅ = P₆, with the maximum pressure of 1.5 bar. The soft robotic manipulator bends in the x-z plane.

The second set of computational experiments investigates the compliance model when P₁ = P₂ = 1.5 bar (ϕ = 0°) and P₁ = P₂ = P₃ = P₄ = 1.2 bar (ϕ = 240°). The applied external forces Δf _x|y|z along three main axes ranged between 0.1 ∼ 0.5 N. By integrating the conventional Cosserat rod model (Till et al. 2019) based on the boundary conditions (35), the robot configuration under external forces can be returned. By solving (36) assuming the external force are zero, the configuration-dependent compliance matrix C ^o can be derived first, the configuration subject to external forces can then be obtained from (19) using C^o, and the tip wrench $\mathrm{\Delta }\mathbf{W}={\left[\mathrm{\Delta }{f}_{x\left|y\right|z}^T,{\left(p\left(l\right)\times \mathrm{\Delta }{f}_{x\left|y\right|z}\right)}^T\right]}^T$ . The differences of the Cartesian tip position ΔE from two methods are compared using (40)

$\mathrm{\Delta }E=\sqrt{{\left(x_{cm}-x_{ecm}\right)}^2+{\left(y_{cm}-y_{ecm}\right)}^2+{\left(z_{cm}-z_{ecm}\right)}^2.}$

(40)

x|y|z _cm are tip positions from solving the Cosserat model, and x|y|z _ecm are calculated tip positions yielding from the compliance matrix. Figures 7(a) and (b) show the comparisons of the tip poses under applied forces from direct solving the Cosserat rod model (coloured in black) and from compliance matrix (the red stars). The extracted tip differences ΔE _x|y|z are shown in Figures 7(c) and (d). It is observed that the differences of the tip position are below 0.15 cm for both in-plane bending cases where ϕ = 0° and ϕ = 240°.

Figure 7.

Results of Simulation 1, showing the deflections from the 6 × 6 Cartesian compliance matrix, where the star marks is the tip position derived from the compliance matrix, based on (19) and (36). The black lines derive from directly solving the Cosserat rod model, assuming external forces Δf_x|y|z are applied at the tip. The results when ϕ = 0° and ϕ = 240° are shown in (a) and (b), respectively. The differences of the tip position ΔE _x , ΔE _y , ΔE _z along the corresponding axes are shown in (c) and (d).

Remark 1

The differences between the nonlinear and linear models become non-negligible especially under high pressure level (see Figure 6(b)). In nonlinear models, high pressure inputs are applied in predefined steps following (25) or (26). Similarly, high external forces can be applied incrementally using the compensation procedures reported in Section 2.3. For Simulation 1, the tip deflection displacements are calculated for small external forces (e.g., below 0.5 N) using two approaches: (1) solving the conventional Cosserat rod model satisfying boundary conditions and (2) calculating the deformed twists from the derived compliance matrix. For in-plane bending with different values of ϕ, both approaches yield similar results. This preliminarily shows the integrated Cartesian compliance matrices are applicable for predicting the deflection displacement under external forces.

4.2. Simulation 2 – Compliance ellipsoid

In line with Section 2.4, the compliance ellipsoid is then adopted to achieve the model-based compliance/stiffness analysis. The nonlinear Cosserat rod model is used in this section. Two sets of simulations were conducted to investigate the compliance distribution along the manipulators under different configurations. The simulation setting is as follows:

• Bending and elongation for a robotic manipulator: For the bending simulation, values for P₁ and P₂ ranged between 0 ∼ 1.5 bar. For the elongation simulation, all chambers were pressurised simultaneously from 0 ∼ 0.9 bar. The difference of pressure values in the two scenarios results in elongation ratios of a similar level.

• Bending for a manipulator made of two robotic segments: In the first simulation, P₁, P₂, P₇, and P₈ were actuated together. In the second simulation, P₁, P₂, P₉, P₁₀, P₁₁, and P₁₂ were actuated. In both simulations, the pressure was between 0 ∼ 1.2 bar.

Figure 8(a) presents the visualisation of the compliance ellipsoids at 0.5l, 0.6l, 0.7l, 0.8l, l along the manipulator when one chamber pair is actuated. The figure indicates that the compliance monotonically increases from the middle to the tip position, from a lower to a higher pressure. This finding can be quantitatively verified from the extracted minimum and maximum eigenvalues of the compliance matrices. For instance, when the pressure increases from 0 to 1.5 bar, the minimum eigenvalues at the middle and tip position increments from 0.06 to 0.11 and 0.13 to 0.36, respectively. Similarly, the corresponding maximum eigenvalues increases from 0.40 to 0.72 and 1.84 to 3.60. Similar findings can be observed from Figure 8(b), that is, the compliance increases under a higher pressure and at a closer distance to the robot tip. Compared to Figures 8(a) and (b) illustrates the robot becomes more compliant for large bending motions. In particular, the overall maximum and minimum eigenvalues in Figure 8(b) are about 2–3 times higher than the values in Figure 8(a). Figure 8(c) shows the results from the elongation scenario. As occurred in the bending test, it is also observed that the overall compliance increases under high pressurisation. When all the chambers are pressurised to 0.9 bar, the maximum eigenvalue at the tip position reaches 9.34 and increases about five times, compared with the value as 1.84 when the pressure is zero.

Figure 8.

Results of Simulation 2—the visualisation of the compliance ellipsoids (scaled by a factor of 0.2) with regards to different configurations for a manipulator with one single segment. The corresponding minimum and maximum eigenvalues of the compliance matrices reported at the bottom. Bending when (a) one chamber pair (P₁ = P₂) is actuated and (b) two chamber pairs (P₃ = P₄ = P₅ = P₆) are actuated: the compliance ellipsoid along the manipulator when the pressure increases from 0 to 1.5 bar in 0.3 bar steps. (c) Elongation: the compliance ellipsoids along the manipulator when P₁ ∼ P₆ increase from 0 to 0.9 bar in 0.3 bar steps.

For a manipulator consisting of two robotic segments, the first and second segments are coloured in blue and red, with a 6 mm thick connection part. Figure 9(a) shows the results when P₁ = P₂ = P₇ = P₈. It indicates that the maximum and minimum eigenvalues at the tips of both robots increase (e.g., the maximum values increase between 90% ∼ 100%) with increasing pressure, which indicates an overall increase in compliance. The maximum eigenvalues from the tip of the second robot are about 6–8 times larger than the values from the tip of the first robot. The eigenvalues also demonstrate that the compliance of the second robot has a greater variation compared to the compliance of the first robot. Figure 9(b) reports the results when P₁ = P₂ = P₉ = P₁₀ = P₁₁ = P₁₂ with similar conclusions, that is, the tip compliance of the second robot is about 6–8 times larger than the tip compliance from the first robot. Compared to Figure 9(a), the minimum eigenvalues from Figure 9(b) show a smaller variation.

Figure 9.

Results of Simulation 2 – Compliance ellipsoids (scaled by a factor of 0.1) at the tip of a robot made of two robotic segments, when (a) P₁ = P₂ = P₇ = P₈ and (b) P₁ = P₂ = P₉ = P₁₀ = P₁₁ = P₁₂. The pressure increases from 0 to 1.2 bar in 0.3 bar steps. The maximum and minimum eigenvalues of the compliance matrices are reported with the magnitude of the actuation pressure.

Remark 2

Simulation 2 shows that the distribution of the compliance along the manipulator can be derived based on the forward kinematics models using our proposed method. This is fundamentally different from, for example, the Jacobian projection (Hussain et al., 2021) or using derivatives (Rucker et al., 2011). In particular, (36) describes how the compliance matrices can be integrated with the static Cosserat rod models simultaneously.

4.3. Simulation 3 – Eigenscrew decomposition

To analyse the coordinate-independent stiffness, the eigenscrew decomposition can be derived based on Section 2.4.2. In the simulation, the pressure was in the range of 0–1 bar and implemented using the nonlinear Cosserat Rod model. Figure 10(a) shows the distribution of p _w , where τ is the maximum positive eigenstiffness. It can be observed the workspace is dome-shaped and enlarged by the elongation. The overall eigenstiffness decreases from 4.37 to 1.49.

Figure 10.

Results of Simulation 3 – Analysis of stiffness at the tip position based on the eigenscrew decomposition. The maximum and minimum positive eigenstiffness values are extracted and then interpolated using natural method from Matlab to assess the configuration-dependent stiffness, which are displayed from the left to the right side in (b), (c), (d). The maximum pressure in each chamber was 1.0 bar. (a) Demonstration of p _w . The colorbar illustrates the eigenstiffness attached to the points of workspace. (b) The stiffness projection onto the x-y plane. (c) The stiffness projection onto the y-z plane. (d) The stiffness projection onto the z-x plane.

Figures 10(b) to (d) show that both maximum and minimum eigenstiffness values are projected onto the x-y, y-z, and z-x planes. Figure 10(b) shows the stiffness decreases from the centre to the periphery of the workspace. For example, the eigenstiffness value is 4.37 at the centre, but that value is 1.93 at the vertex. Figures 10(c) and (d) convey the same message that the stiffness projected onto the y-z plane and z-x plane also decreases when the workspace becomes larger. Please note that only the enclosed area by the red dotted points shown in Figures 10(b) to (d) is the reachable workspace.

Remark 3

Simulation 3 shows that the stiffness of the manipulator varies in the workspace, and the stiffness capability, described by eigenstiffness (independent from the coordinates), decreases under a high pressure level. Large axial elongation makes the robotic manipulator have more compliant responses. This will be further investigated and validated in Section 5.3.

5 Experimental validation of the forward kinematics and stiffness model

In this section, thorough validation procedures of the proposed method and its results are presented comparing our numerical models (LCM, NCM, LPCC and NPCC) from Section 3.3 with experimental data in Experiments 1–3. Experiments 1–2 use a manipulator made of one segment for validation. Experiment 1 reveals the nonlinear kinematic responses are observed and accommodated. Experiment 2 validates the configuration-dependent stiffness/compliance model. Experiment 3 validates the configuration-dependent compliance modelling using manipulators made of one and two robotic segments.

5.1. Experimental setup

The experimental setup is presented in Figure 11. As shown in Figure 11(a), three proportional pressure regulators (Camozzi K8P) are used to modulate pressurised air in each chamber pair. These are controlled by analogue signals (0 ∼ 10 V) with a pressure feedback channel (0 ∼ 3 bar). A compressor (BAMBI MD Range Model 150/500) is chosen as the pneumatic reservoir supplying constant pressurised air to the pressure regulators. NI-DAQ USB-6341 devices with multiple digital and analogue I/Os are used to monitor and control the pressure regulators. The NI-DAQ USB-6341 provides PWM signals, which are converted to analogue signals (0 ∼ 10 V) by Pulse-Width Modulation (PWM) converters, to control each pressure regulator. In Experiment 3, another three pressure regulators are added and controlled in the same way.

Figure 11.

Experimental setup. It comprises two subsystems: (a) The hardware for power and actuation mainly contains the air compressor and electronics, such as pressure regulators, PWM converter (denoted using DAC here) and NI USB-6341 boards, to actuate and control the soft robotic manipulator. A laptop running Matlab software collects data and sends control commands. (b) The forward kinematics and stiffness validation setup contains a 7 DoF Franka Emika Panda robot with the soft robotic manipulator mounted on the end-effector, NDI Aurora magnetic trackers, a 6 DoF force/torque sensor, and a linear rail.

Figure 11(b) shows a 6 Degree-of-Freedom (DoF) Aurora sensor mounted on the tip of the soft robotic manipulator to measure its pose using an NDI Aurora tracking system. The position and orientation tracking errors of the Aurora system (with 6 DoF sensors and a cube volume) are less than 0.48 mm and 0.30°, respectively. Stiffness experiments were conducted by pulling the manipulator’s tip via a 0.19 mm diameter nickel wire connected to an IIT-FT17 6 DoF Force/Torque (F/T) sensor. The F/T sensor was mounted on a motorised linear rail (Zaber X-LSM100 A), which controlled the pulling speed and displacement. A 1.5 m length wire was chosen to guarantee a high confidence in the angle of the applied pulling force. For instance, when Δf _x is applied along the x-axis, a 1 cm displacement error in the z-axis would result in an angular error of about arctan 0.01/1.5 = 0.38°, which is negligible. For the stiffness validation experiments, the soft manipulator was fixed to a 7 DoF robotic arm (Franka Emika Panda) to arrange the desired pose and measure the configuration-dependent stiffness. All hardware was connected to a laptop running Matlab software to capture data, control the soft robot and the 7 DoF robotic arm, and post-process the acquired data.

5.2. Experiment 1 – Validation of the forward kinematics and nonlinear compensation

To validate the forward kinematics and the nonlinearity, three sets of experiments were conducted by pressurising one, two and three chamber pair(s) simultaneously. When one chamber pair was pressurised, the pressure was from 0 to 1.5 bar and followed the sequence: (P₃, P₄ → P₅, P₆ → P₁, P₂). When two chamber pairs were pressurised, the pressure was from 0 to 1.2 bar and followed (P₁, P₂, P₃, P₄ → P₃, P₄, P₅, P₆ → P₁, P₂, P₅, P₆). When three chamber pairs were actuated, the pressure in all six chambers was from 0 ∼ 1.2 bar and increased simultaneously. During these experiments, the tip position of the soft manipulator and chamber pressures were recorded at 40 Hz. The maximum pressure values in both bending experiments were set differently, so the maximum bending angle reached the same value in both experiments. The pressure was incrementally increased by 0.15 bar. In the simulation, the threshold-based method was chosen to discretise the pressure with a threshold value of 0.15 bar, when one chamber pair was pressurised according to (26). The finite-step method was chosen with n = 6, when two chamber pairs were pressurised based on (25).

Figure 12 illustrates extracted static elongation and bending angles. The experimental data are plotted in blue circles and linearly interpolated by curve fitting (see blue curve). The LCM, NCM, LPCC and NPCC are plotted in green, magenta, dotted black, and dotted red, respectively. In addition, the elongation ratio is displayed at the bottom of Figures 12(a) to (c), derived by computing the change of the length of the backbone of the manipulator over the default length. Proportionality constants are 0.4 cm, 52.5°, and 30.0° per 0.1 of the elongation ratio in Figures (a), (b) and (c), respectively. Figures 12(b) and (c) report the elongation ratios for the purpose of illustrating the level of elongation when the robot has a bending motion. From Figure 12(a), it can be seen that the manipulator achieves a maximum elongation of 2.7 cm. The results of all numerical models and experiments are similar for pressures up to 0.45 bar and an elongation ratio of 0.15. Above these values, the curves of the LCM and LPCC behave in a linear way compared to the two nonlinear models and experiments. The overall Root Mean Square Error (RMSE) of two nonlinear models is 0.07 cm, while the RMSE of two linear models is 0.49 cm (see Table 3).

Figure 12.

Results of Experiment 1 – nonlinear kinematic responses: Experiments are compared to results from (non)linear models. Applied pressures are plotted against (a) elongation displacement and elongation ratios when all chamber pairs are actuated, and angles and elongation ratios when (b) one chamber pair and (c) two chamber pairs are pressurised. The orange stars demonstrate the level of elongation when any nonlinearity becomes significant. (a) Reports the elongation ratios from experiments; (b) and (c) report the ratios from the NCM. Comparisons are reported on in Extension 1.

Table 3.

Summary of results for Experiment 1: Comparison of errors for the linear and nonlinear simulation models with experiments.

Pressurisation	Model	Max.	Min.	RMSE	Pressurisation	Max.	Min.	RMSE	Pressurisation	Max.	Min.	RMSE
		[cm]	[cm]	[cm]		[°]	[°]	[°]		[°]	[°]	[°]
Elongation Figure 12(a)	LPCC	1.06	0.05	0.49	Bending Figure 12(b)	28.38	0.74	12.45	Bending Figure 12(c)	38.20	0.19	16.36
	LCM	1.06	0.05	0.49		24.10	1.24	10.38		33.45	1.68	14.10
	NPCC	0.13	0.03	0.07		12.70	0.42	5.91		11.90	0.85	5.61
	NCM	0.13	0.03	0.07		5.43	1.10	3.79		6.62	1.58	4.28

In Figures 12(b) and (c), the maximum bending angles are similar and around 110°. The maximum elongation ratio is about 0.2 and 0.4 when one and two chamber pairs are pressurised. The nonlinearity is notable when the elongation ratio is larger than 0.1. The nonlinear models outperform linear model sets. For example, the NCM shows the smallest RMSE in two bending scenarios, which are 3.79° and 4.28°. By comparison, the linear models are less accurate. Specifically, the LCM, with the RMSE of 10.38° and 14.10°. Moreover, the Cosserat rod models show a better performance that the PCC models. A detailed summary of the errors is shown in Table 3.

Figure 13 shows the results for the forward kinematics validation, with detailed summary of errors in Table 4. The top row plots of Figure 13 show the triangular-shaped pressurisation pattern of the chamber pairs followed by corresponding position measurements. The bottom two rows plot the positional error for the (non)linear Cosserat rod and PCC models compared to the experimental results in the x-, y-, and z-directions.

Figure 13.

Results of Experiment 1 – Forward kinematics validation: (a) One chamber pair actuation with pressure amplitudes of 0.9, 1.2, and 1.5 bar. (b) Two chamber pairs actuation with pressure amplitudes of 0.6, 0.9, and 1.2 bar. x _exp , y _exp , and z _exp are the experiment positions. x _sim , y _sim , and z _sim are the simulation results for the nonlinear Cosserat model. $LC{M}_{e_x}, LC{M}_{e_y}, \text{and}\hspace{0.17em}LC{M}_{e_z}$ are errors from the linear Cosserat model. $NC{M}_{e_x}, NC{M}_{e_y}, \text{and}\hspace{0.17em}\hspace{0.17em}NC{M}_{e_z}$ are errors from the nonlinear Cosserat model. $LPC{C}_{e_x}, LPC{C}_{e_y}, \text{and}\hspace{0.17em}LPC{C}_{e_z}$ are errors from the linear PCC model. $NPC{C}_{e_x}, NPC{C}_{e_y}, \text{and}\hspace{0.17em}NPC{C}_{e_z}$ are errors from the nonlinear PCC model. The comparisons are also reported in the Extension 1.

Table 4.

Summary of results for Experiment 2: Tip position errors.

	Model	$\left|e_x\right|$	$\left|e_y\right|$	$\left|e_z\right|$	$r_{e_x}$	$r_{e_y}$	$r_{e_z}$
		[cm]	[cm]	[cm]	[cm]	[cm]	[cm]
Figure 13(a): One chamber Pair	LPCC	1.27	0.35	0.92	0.35	0.14	0.27
	LCM	1.10	0.41	0.77	0.31	0.15	0.23
	NPCC	1.03	0.38	0.77	0.30	0.15	0.22
	NCM	0.64	0.39	0.62	0.24	0.20	0.21
Figure 13(b): Two chamber	LPCC	0.74	1.27	1.99	0.28	0.39	0.58
	LCM	0.64	1.07	1.74	0.23	0.32	0.51
	NPCC	0.62	0.89	1.57	0.20	0.29	0.46
	NCM	0.48	0.43	0.70	0.18	0.17	0.24

$r_{ex}$ is the RMSE error in the x-axis, $r_{ey}$ in the y-axis, $r_{ez}$ in the z-axis. $\left|e_x\right|, \left|e_y\right|, \left|e_z\right|$ are the maximum absolute errors in the corresponding directions.

Figure 13(a) shows the forward kinematics results when one chamber pair is pressurised with different amplitudes. During 2.5 ∼ 12 s, the errors of the linear and nonlinear models are similar for this value of pressurisation. During 15.5 ∼ 25 s, the amplitude is set to 1.2 bar, resulting in an increase of the error values for two linear models. The maximum error reaches 0.76 cm for the LCM and 0.91 cm for the LPCC. This error further increases during 28.5 ∼ 38 s, when the maximum pressure is 1.5 bar. In this instance, the maximum error in the x-axis is 1.03 cm for the NPCC and 0.64 cm for the NCM. This error is larger for the linear models (1.26 cm for the LPCC and 1.09 cm for the LCM).

Figure 13(b) shows the results when two chamber pairs are pressurised. Similar to Figure 13(a), the error values of the four linear and nonlinear models are almost identical and below 0.5 cm, when the pressure inside each of the two chamber pairs is below 0.6 bar during 2.5 ∼ 12 s. The discrepancy between the linear and nonlinear models again becomes larger with the increase of the pressure. For instance, when the pressure is increased to 1.2 bar, the maximum errors in the x-, y-, and z-directions for the LCM are 0.64, 1.07, 1.74 cm. Conversely, the errors of the NCM are 0.47, 0.43, 0.70 cm in three directions. Overall, the NCM returns the highest accuracy with an RMSE value below 0.24 cm compared to RMSE values of 0.46, 0.51, 0.58 cm for the LCM, NPCC, and LPCC (see Table 4).

5.3. Experiment 2 – Validation of the stiffness and compliance for a one robotic segment

Two types of experiments are used to validate the compliance/stiffness modelling. The translational compliance matrix $C_t^o$ is validated. The experimental compliance ${\left(C_t^o\right)}_e$ is obtained by

${\left(C_t^o\right)}_e=\left[\begin{array}{lll}\frac{\mathrm{\Delta }{x}_{\mathrm{\Delta }{f}_x}}{\mathrm{\Delta }{f}_x}& \frac{\mathrm{\Delta }{x}_{\mathrm{\Delta }{f}_y}}{\mathrm{\Delta }{f}_y}& \frac{\mathrm{\Delta }{x}_{\mathrm{\Delta }{f}_z}}{\mathrm{\Delta }{f}_z}\\ \frac{\mathrm{\Delta }{y}_{\mathrm{\Delta }{f}_x}}{\mathrm{\Delta }{f}_x}& \frac{\mathrm{\Delta }{y}_{\mathrm{\Delta }{f}_y}}{\mathrm{\Delta }{f}_y}& \frac{\mathrm{\Delta }{y}_{\mathrm{\Delta }{f}_z}}{\mathrm{\Delta }{f}_z}\\ \frac{\mathrm{\Delta }{z}_{\mathrm{\Delta }{f}_x}}{\mathrm{\Delta }{f}_x}& \frac{\mathrm{\Delta }{z}_{\mathrm{\Delta }{f}_y}}{\mathrm{\Delta }{f}_y}& \frac{\mathrm{\Delta }{z}_{\mathrm{\Delta }{f}_z}}{\mathrm{\Delta }{f}_z}\end{array}\right]=\left[\begin{array}{lll}{c}_{11}& c_{12}& c_{13}\\ c_{21}& c_{22}& c_{23}\\ c_{31}& c_{32}& c_{33}\end{array}\right],$

(41)

where Δf = (Δf _x , Δf _y , Δf _z ) is the pulling force along the x-, y-, and z-axes and measured by the sensor. The corresponding displacements were recorded by the Aurora system. For instance, $\mathrm{\Delta }{x}_{\mathrm{\Delta }{f}_x}$ is the displacement along the x-axis resulting from Δf _x . The experimental compliance matrix can then be derived column by column, as shown in (41). During all experiments, the pulling displacement values were set to 0.25 cm, 0.3 cm, 0.35 cm, 0.4 cm, respectively, with a pulling speed of 0.02 cm/s.

For planar stiffness validation (Shiva et al., 2016; Ranzani et al., 2016), the manipulator had planar motions, within the x-z plane. The stiffness along the x-, y- and z-axes was investigated. When the manipulator bent to the -x-axis, the pressure of P₁ and P₂ was in the range of 0 ∼ 1.5 bar. When the manipulator bent to the +x-axis, the pressure of P₃, P₄, P₅, and P₆ was in the range of 0 ∼ 1.2 bar. For Cartesian compliance validation, the 3 × 3 compliance matrix $C_t^o$ was validated under six different configurations, where the experiments were conducted pressurising one, two and three chamber pairs in the range of 0 ∼ 1.2 bar. The compliance was investigated when the robot bends with different values of ϕ. Extension 2 reports on the details of this experiment.

Figure 14(a) shows the stiffness response in three directions (i.e., x-, y-, and z-axes) with one chamber pair pressurised. As the bending angle increases due to pressurisation up to about 0.9 bar, the experimentally acquired stiffness values along the x-axis decreases first and then reaches a bending angle of about 90° with a maximum stiffness of 1.21 ± 0.17 N/cm. The four models return different error values when compared with the experimental results. The errors between all models and the experimental data are within 0.10 N/cm when the pressure is less than 0.6 bar. When the pressure further increases to 1.5 bar, the error is 0.28 ± 0.17 N/cm for the NCM, 0.05 ± 0.17 N/cm for the LCM, 0.47 ± 0.17 N/cm for the NPCC, and 0.42 ± 0.17 N/cm for the LPCC. By contrast, the overall stiffness monotonically decreases in the y-direction, from 0.57 ± 0.04 N/cm to 0.27 ± 0.02 N/cm. The error values of the two nonlinear models are within 0.05 N/cm. The error is larger than 0.10 N/cm for the linear models. The overall stiffness also decreases in the z-axis, from 7.33 ± 0.33 N/cm to 0.42 ± 0.03 N/cm.

Figure 14.

Results of Experiment 2 – Planar stiffness along the x-, y-, and z-axes. (a) Stiffness response when P₁ and P₂ are actuated and between 0 ∼ 1.5 bar. (b) Stiffness response when P₃, P₄, P₅, and P₆ are actuated and between 0 ∼ 1.2 bar. The experimental data points are plotted using blue circles. Stiffness is calculated from the measured displacement and force data. The shaded areas in orange colour span from the minimum to maximum stiffness values. The computational results from our models are plotted in a solid pink, green as well as dotted red and black coloured curve.

Figure 14(b) shows the results when two chamber pairs are pressurised. The stiffness response in all directions follows a similar pattern. However, all stiffness values are lower. For example, the stiffness value in the y-direction is 0.27 ± 0.02 N/cm in Figure 14(a); this value is about two times lower (0.12 ± 0.01 N/cm) in Figure 14(b). As mentioned earlier, the pressurisation of one chamber pair with 1.5 bar and two chamber pairs with 1.2 bar each achieve a bending angle of approximately 90°. In line with previous results, the modelling errors of the linear models increase with the pressure. For instance, at a bending angle of about 90°, the maximum stiffness error along the y-axis is 0.08 ± 0.02 N/cm when one chamber pair was pressurised and 0.16 ± 0.01 N/cm when two chamber pairs are actuated.

To quantitatively assess the accuracy, the Frobenius norm is utilised, and the error e _C in (42) is defined as the difference of the Frobenius norm between the simulated and experimental compliance

$e_C=\left|\frac{\Vert {\left(C_t^o\right)}_e{\Vert }_F-\Vert C_t^o{\Vert }_F}{{\Vert {\left(C_t^o\right)}_e\Vert }_F}\right|\times 100\%.$

(42)

Table 5 gives a summary of the results for the spatial compliance validation. Here, the first column presents the actuation pressure. For each of the actuation patterns, average values of the nine compliance elements derived from (41) are displayed from the third to 11th column. The Frobenius norm and the error defined by (42) are shown in the last two columns. Table 5 shows that the overall compliance is configuration-dependent and decreases with the pressurisation, for example, the Frobenius norm of the experimental compliance matrix is 3.38 cm/N when the pressure P is [0.6,0.6,0,0,0,0] bar (see row 1), but that value increases to 12.57 cm/N when the pressure is [0.9,0.9,0.9,0.9,0.9,0.9] bar (see row 6). Furthermore, the accuracy of the linear models decreases with the pressure. For instance, when P = [0.6, 0.6, 0, 0, 0, 0] bar, the four models return a similar performance, with the error as 4.07%, 8.68%, 4.07%, 7.76% for NCM, LCM, NPCC, and LPCC models, respectively. When the pressure increases to P = [0.9, 0.9, 0.9, 0.9, 0.9, 0.9] bar, the corresponding errors are 5.16%, 58.80%, 5.90%, 58.40%. Although the maximum errors of the nonlinear models in Table 5 are larger than 58%, the maximum errors of the nonlinear PCC and the nonlinear Cosserat rod model are below 10% and 15%, respectively.

Table 5.

Summary of the results for the validation of spatial compliance in Experiment 2 using a one-segment robot.

Setup	Model	c 11	c 21	c 31	c 12	c 22	c 32	c 13	c 23	c 33	$\Vert .\Vert$ F	Error (%)
		[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]
P = 0.6 [1, 1, 0, 0, 0, 0] bar	Experiment	2.28	0.15	0.63	0.12	2.23	0.07	0.70	0.03	0.50	3.38	—
	NCM	2.2	0	0.72	0	2.52	0	0.73	0	0.42	3.52	4.07
	LCM	1.92	0	0.63	0	2.21	0	0.64	0	0.37	3.09	8.64
	NPCC	2.2	0	0.67	0	2.51	0	0.68	0	0.38	3.52	4.07
	LPCC	2	0	0.58	0	2.23	0	0.57	0	0.32	3.12	7.76
P = 1.2 [1, 1, 0, 0, 0, 0] bar	Experiment	1.45	0.13	1.26	0.05	3.55	0.14	1.34	0.13	1.65	4.58	—
	NCM	1.67	0	1.47	0	3.62	0	1.48	0	1.81	4.85	6.17
	LCM	1.29	0	1.06	0	2.62	0	1.07	0	1.26	3.52	22.94
	NPCC	1.91	0	1.50	0	3.57	0	1.48	0	1.57	4.83	5.74
	LPCC	1.63	0	1.07	0	2.63	0	1.07	0	0.99	3.58	21.63
P = 0.9 [1, 1, 1, 1, 0, 0] bar	Experiment	3.92	−1.15	1.05	−0.91	3.12	1.72	1.22	1.25	2.14	6.24	—
	NCM	4.58	−1.01	1.07	−1.02	3.40	1.87	1.10	1.90	2.20	7.00	12.17
	LCM	2.72	−0.53	0.62	−0.53	2.11	1.08	0.63	1.09	1.19	4.11	34.14
	NPCC	4.40	−0.78	1.02	−0.77	3.53	1.75	1.00	1.72	1.69	6.63	6.24
	LPCC	2.6	−0.32	0.66	−0.32	2.53	0.82	0.66	0.81	0.70	4.04	35.26
P = 0.9 [1, 1, 0, 0, 1, 1] bar	Experiment	3.87	1.22	0.97	0.95	3.02	−1.46	1.20	−1.56	2.17	6.18	—
	NCM	4.59	1.01	1.07	1.04	3.38	−1.90	1.10	−1.90	2.20	7.00	13.34
	LCM	2.72	0.53	0.62	0.54	2.10	−1.09	0.63	−1.09	1.19	4.11	33.45
	NPCC	4.41	0.78	1.02	0.75	3.55	−1.71	1.00	−1.72	1.69	6.63	7.35
	LPCC	2.83	0.35	0.58	0.34	2.45	−0.98	0.57	−0.99	0.82	4.18	34.59
P = 1.2 [0, 0, 1, 1, 1, 1] bar	Experiment	1.25	0.11	−1.55	0.08	7.12	0.15	−1.92	0.20	4.72	8.88	—
	NCM	1.44	0	−2.12	0	8.00	0	−2.16	0	5.59	10.32	14.88
	LCM	1.07	0	−1.22	0	3.51	0	−1.23	0	2.18	4.60	48.80
	NPCC	2.35	0	−2.71	0	7.39	0	−2.69	0	4.44	9.72	8.20
	LPCC	1.93	0	−1.45	0	3.51	0	−1.46	0	1.50	4.75	47.13
P = 0.9 [1, 1, 1, 1, 1, 1] bar	Experiment	8.85	0.11	0.12	0.10	8.92	0.14	0.07	0.05	0.31	12.57	—
	NCM	9.34	0	0	0	9.34	0	0	0	0.28	13.22	5.16
	LCM	3.66	0	0	0	3.66	0	0	0	0.18	5.18	58.80
	NPCC	8.36	0	0	0	8.36	0	0	0	0.26	11.83	5.90
	LPCC	3.67	0	0	0	3.67	0	0	0	0.16	5.23	58.40

Furthermore, the compliance ellipsoids when P = [0.9, 0.9, 0.9, 0.9, 0, 0] bar are projected in Figure 15. Each row, from left to right, shows the ellipsoids’ projection onto the x-y, y-z, and x-z planes for Δf = [0.5 cos θ, 0.5 sin θ, 0] N, Δf = [0, 0.5 cos θ, 0.5 sin θ] N, and Δf = [0.5 cos θ, 0, 0.5 sin θ] N in Figures 15(a) to (c), respectively. The results demonstrate that both linear models and both nonlinear models return a similar performance. The size of the major and minor axes of the compliance ellipsoids from the two linear models return large errors.

Figure 15.

Results of Experiment 2 – Ellipsoid visualisation of the compliance matrices. NCM, LCM, NPCC, LPCC represent the simulated compliance from the nonlinear Cosserat, linear Cosserat, nonlinear PCC and linear PCC model, separately. The compliance ellipsoid is shown for P = [0.9, 0.9, 0.9, 0.9, 0, 0] bar and projected onto the x-y, y-z, and x-z plane from the left to the right column for (a) Δf = [0.5 cos θ, 0.5 sin θ, 0] N, (b) Δf = [0, 0.5 cos θ, 0.5 sin θ] N and (c) Δf = [0.5 cos θ, 0, 0.5 sin θ] N, with θ = 0 ∼ 2π.

5.4. Experiment 3 – Validation of the compliance for a robot made of two robotic segments

The results of Experiment 3 validate our approach for a manipulator made of two robotic segments, where the robot can generate out-of-plane bending. As reported in Experiments 1–2, the NCM method provides the most accurate forward kinematics and compliance modelling. Therefore, we selected the NCM method for validation in this section. During the experiments, the two robotic segments were actuated achieving six configurations with actuation pressures ranging from 0.6 to 1.2 bar. Similar to Experiment 2, the overall tip compliance was experimentally determined using (41). The pulling forces Δf were applied using calibrated 1 g and 2 g weights, and the corresponding displacements were recorded by an Aurora tracker attached to the overall tip. Compliance errors between simulations and experiments were described using Frobenius norms and defined by (42). Extension 3 reports on the compliance measurement in this experiment.

Figure 16 displays the simulated compliance ellipsoids for a manipulator made of two segments and the captured robot images from experiments under six configurations. The compliance ellipsoids are scaled by a factor of 0.1 for visualisation purposes. The two robots are coloured in blue and red, and the connection part is black. The figure demonstrates that our model is applicable to robots made of two robotic segments and is highly aligned with the experimental results, returning maximum tip position errors of less than 6.0 mm. Table 6 details the validation results for the compliance matrices. The results show that the compliance modelling errors, evaluated by Frobenius norms, are between 10% ∼ 20%. Compared to the Frobenius norms of the compliance matrices from Table 5, the values in Table 6 at the tip of the robot made of two robotic segments are generally 5 ∼ 10 times larger. For instance, forces on the order of 0.1 N typically result in displacements on the order of 1 cm when the robot made of two robotic segments has a bending motion. In contrast, tip forces on the order of 1 N produce displacements on the order of 1 cm for one robotic segment. The tip compliance of the robot made of two robotic segments exhibits significant compliant behaviour in the x- and y-directions, while the compliance along the z-axis is much smaller, as shown in Figure 16(f) and the last row of Table 6.

Figure 16.

Results of Experiment 3 – Simulated compliance ellipsoids at the tip of a manipulator of two robotic segments under six configurations when (a) P₁ = P₂ = P₇ = P₈ = 0.9 bar, (b) P₁ = P₂ = P₉ = P₁₀ = P₁₁ = P₁₂ = 0.9 bar, (c) P₁ = P₂ = P₉ = P₁₀ = 1.2 bar, (d) P₅ = P₆ = P₉ = P₁₀ = 1.2 bar, (e) P₁ = P₂ = P₅ = P₆ = P₇ = P₈ = P₉ = P₁₀ = 0.9 bar, and (f) P₁ ∼ P₁₂ = 0.6 bar. An image of the robot’s configuration from experiments is included with the simulations. The experimental and simulated overall tip positions are presented with the unit of cm.

Table 6.

Results for the validation of spatial compliance in Experiment 3 using a manipulator made of two robotic segments.

Setup	Model	c 11	c 21	c 31	c 12	c 22	c 32	c 13	c 23	c 33	‖ ⋅‖F	Error (%)
		[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]	[cm/N]
P = 0.9 [1,1,0,0,0,0,1,1,0,0,0,0] bar	Experiment	9.74	4.50	10.20	5.50	34.45	3.35	5.84	4.12	19.75	43.45	—
	NCM	7.84	0	11.81	0	36.99	0	11.86	0	23.88	47.75	9.87
P = 0.9 [1,1,0,0,0,0,0,0,1,1,1,1] bar	Experiment	22.02	3.54	9.86	6.35	32.10	2.74	9.87	3.21	4.35	42.51	—
	NCM	28.02	0	12.67	0	35.46	0	12.73	0	6.65	49.08	15.46
P = 1.2 [1,1,0,0,0,0,0,0,1,1,0,0] bar	Experiment	19.87	−7.89	15.21	−7.98	30.05	9.95	14.21	7.32	19.75	48.99	—
	NCM	23.04	−9.99	16.13	−9.96	34.64	8.83	16.21	8.93	23.75	56.34	15.02
P = 1.2 [0,0,0,0,1,1,0,0,1,1,0,0] bar	Experiment	19.38	−8.75	−13.32	−9.70	30.12	−11.10	−13.21	−7.64	19.86	48.80	—
	NCM	22.97	−10.06	−15.83	−9.85	33.96	−9.54	−15.84	−9.57	23.75	56.02	14.80
P = 0.9 [1,1,0,0,1,1,1,1,1,1,0,0] bar	Experiment	39.65	8.97	23.14	15.70	49.65	−21.15	18.72	−12.53	20.50	79.21	—
	NCM	45.75	11.08	26.28	11.70	59.20	−15.79	26.52	−15.63	27.12	92.09	16.26
P = 0.6 [1,1,1,1,1,1,1,1,1,1,1,1] bar	Experiment	56.21	3.21	2.42	2.55	59.32	3.45	0.08	0.05	0.49*	81.93	—
	NCM	67.72	0	0	0	67.72	0	0	0	0.54	95.77	16.88

Note: the value with * was identified using a 0.2 N force, that is, the value of c₃₃ when P = [0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6] bar.

6. Discussions and conclusion

To model the stiffness/compliance of soft robots and achieve the configuration-dependent stiffness analysis, a Lie theory–based stiffness modelling and analysis framework was proposed and experimentally validated using soft manipulators made of one and two robotic segments. In addition, a thorough stiffness/compliance analysis is presented both from simulation and experiment. The main findings of soft manipulators with an elongation capability can be summarised as

1. Large elongations in soft robotic manipulators often introduce nonlinear responses. We proposed two model sets, linear model sets (LPCC and LCM) for small elongations and nonlinear model sets (NPCC and NCM) for large elongations. This improves the modelling accuracy in both forward kinematics and stiffness/compliance.

2. The overall compliance of soft manipulators is influenced by the elongation, for example, the used pneumatic-driven manipulators show compliant behaviours under a high pressure level.

6.1. Nonlinear compensation

Experiment 1 shows the nonlinearity appears for an elongation ratio over 0.15 in cases of elongation (see Figure 12(a)) and 0.1 in cases of bending (see Figures 12(b) to (c)). Below these elongation ratios, the material deforms in a linear manner. The results show that these nonlinearity can be captured by the proposed compensation method, for the nonlinear Cosserat rod and PCC model both show a satisfactory efficacy reproducing the nonlinear responses compared to the linear models. The nonlinearity introduced by elongation further influences the forward kinematics of robots, with the increase of the pressure (see Figure 13). In addition, the results show that both Cosserat models return lower error values compared to the PCC models, as the constant curvature might not be strictly satisfied. Further, the predicted elongation displacement and bending angle of the four models in Figure 12 are slightly higher than the measured results when the pressure is low, which may come from the linearisation of the material model.

Our framework’s nonlinear compensation procedures are based on the theoretical support of the nonlinear matrix structure analysis approach outlined in Gazzola et al. (2018). This allows our framework to consider large deformations using a linear strain–stress relation without the need to change the fundamental models. Similar updates have been reported in (Gazzola et al., 2018; Sadati et al., 2021) to account for the cross-sectional deformations resulting from elongation. It is worth noting that by introducing the stiffness update (see (24)), our compensation approach provides an approximation that is equivalent to the hyper-elastic Neo-Hookean and Mooney–Rivlin models, particularly when the longitudinal strain is below 30% (Gazzola et al., 2018). Figure 12 illustrates that our compensation approach remains effective even when the longitudinal strain is up to approximately 60%.

6.2. Compliance/stiffness modelling and analysis

Experiment 2 corroborates that the compliance is configuration-dependent, in addition, it increases with the pressure. The compliance results of the two nonlinear models show less errors than the two linear models, thanks to the proposed nonlinear compensation method. As concluded in Table 5, there are differences in the compliance matrices when P changes from [0.9, 0.9, 0.9, 0.9, 0, 0] bar to [0.9, 0.9, 0, 0, 0.9, 0.9] bar, but their Frobenius norms are similar. This can be explained by that the configuration determines the shape of the compliance ellipsoids, whereas the pressurisation level determines its size. Experiment 3 provides additional validation results for a robot consisting of two robotic segments. According to Table 6, the compliance modelling errors for the two-segment robot are between 10% ∼ 20%. By contrast, Table 5 indicates that compliance modelling errors for a single robotic segment using the NCM are smaller and between 5% ∼ 15% As Figures 13 and 16 show, errors in forward kinematics modelling can result in compliance errors, since compliance is dependent on the configuration. The increased modelling errors observed in Experiment 3 may be attributed to error propagation, as errors in the first manipulator accumulate in the second manipulator. Additionally, it should be noted that discrepancies between the two manipulators resulting from the fabrication process (e.g., the use of reinforced fibres) could introduce unmodelled errors.

The experimental error sources of the compliance validation mainly come from: firstly, the measured stiffness values are influenced by the pulling displacement under different configurations. For instance, in directional stiffness validation, the measured stiffness in the x-axis is 1.40 N/cm when the tip is pulled by 0.4 cm. This is around 0.35 N/cm larger than the stiffness of 1.05 N/cm under 0.25 cm pulling displacement. Secondly, the pressurised chambers may stiffen the robot to some extent. For instance, the predicted compliance values from the nonlinear models are larger than experimental data when the pressure increases (see Figure 15 and Table 5). In addition, the Aurora tracking system could introduce measurement errors ( $\le 0.5$ mm) and the pressure regulators could add pressure control errors ( $\le 0.03$ bar).

The coordinate-independent stiffness analysis shows that soft robots with elongation capability have an enlarged workspace, however, with a lower stiffness (see Section 4.3). This is in agreement with the analysis in Section 4.2: When the bending angle increases, the overall stiffness decreases. More specifically, the deformation of the cross-section resulting from the elongation reduces the overall stiffness (e.g., see (24)). There is a trade-off between the flexibility introduced by elongation and the stiffness capability.

6.3. Remarks of proposed methodology

The proposed stiffness modelling and analysis framework can account for nonlinear responses resulting from elongation, even when a linear strain–stress relationship is used, thus improving the accuracy of the model. The framework can also directly incorporate well-investigated forward kinematics models in the field of soft robotics, such as the PCC and Cosserat rod model, to deliver a configuration-dependent compliance modelling and analysis, as reported in Section 2, and detailed stiffness/compliance analysis tools are presented. The framework also reveals the compliance at the tip position and the overall distribution of compliance along the manipulator. Furthermore, the proposed framework has been extensively validated through simulations and experiments, with detailed analysis using soft manipulators consisting of one and two robotic segments (see Sections 4 and 5).

Moreover, the Lie theory–based stiffness modelling preserves the inherent symmetrical property of the compliance matrix as the total compliance is built on the exact material compliance. By contrast, the methods in Rucker et al. (2011); Smoljkic et al. (2014); Black et al. (2017) cannot strictly guarantee such symmetry. In addition, the Cartesian compliance tells how the robot deflect under external loads (see ((19))), so it can also be utilised for, for example, stiffness or force control, and perturbation rejection (Mahvash and Dupont 2011; Mustaza et al., 2018; Lai et al., 2021).

6.4. Applicability of the modelling framework

This work proposes a comprehensive stiffness modelling and analysis framework for soft continuum robots, particularly for fluidic-driven soft robots made of silicone material. We chose the STIFF-FLOP manipulator as a test bed for validating our approach due to its nonlinear strain–stress relationship caused by large longitudinal elongation, making it a complex system to model. However, it is important to note that our proposed methodology is not restricted to fluidic-driven robots, but can be applied to other types of manipulators such as tendon-driven and parallel continuum robots that employ forward kinematics models built on stiffness density matrices, as demonstrated in (2) and previous studies (Till et al. (2019); Amanov et al. (2021); Rucker et al. (2011)). These models leverage stiffness properties to calculate curvatures and strains from the actuation space, such as pressure or tendon forces, which are then used to complete the statics and dynamics models. To further demonstrate the applicability of our framework to other robots, we applied it to a parallel continuum robot of a Stewart–Gough configuration (Black et al. (2017)). Although simulation results are only presented in Appendix E as supporting material, they confirm that our framework has potential to be used in various other robotic systems.

To summarise, the applicability of our framework is not restricted by the robot types. Rather, it hinges on the availability of compliance or stiffness properties (either identified or modelled) of the robotic systems. Forward kinematics models based solely on curve fitting (Zhao et al. (2022)) or the length of each bellow (Garbin et al. (2018)) cannot be utilised with our stiffness modelling framework since no quantitative stiffness properties are explicitly incorporated into these statics models.

6.5. Future work

As discussed by Manti et al. (2016), antagonistic actuation principles can influence the stiffness of soft robots while maintaining a position, for instance, tendon-driven robots (Oliver-Butler et al. 2019). As for the fluidic-driven robots in this study, this stiffening effect is compromised because the elongation and configuration play a more important role with regards to changing the overall stiffness. However, this may still influence the accuracy of the modelling. As shown in the Experiment 2, the experimental stiffness can be about 10% ∼ 15% higher than the predicted stiffness from nonlinear models when the pressure increases. This may be tackled by taking the stiffness of the pressurised chambers as a superposition on the stiffness of the main body of robots in future study (Chikhaoui et al., 2019).

Second, the computations in this study are offline. Our purpose is to model and analyse the configuration-dependent compliance, so the real-time performance is beyond the scope of this study. All simulations and experiments were executed on a computer equipped with an Intel i7 processor and 8 GB RAM, running Matlab R2019. The computational time for the simulations conducted in Section 4 and Section 5 is summarised in Table 7. Table 7 demonstrates that the computational time tends to be higher under the maximum actuation pressure level, in particular for the Cosserat rod model sets. In addition, the average computational time is about 0.01–0.02 s, 0.07–0.12 s, 0.8–1.5 s, and 3–6.3 s for the LPCC, NPCC, LCM and NCM models, respectively. As such, both the (non)linear PPC model sets need less computational time than the (non)linear Cosserat rod model sets. Meanwhile, the nonlinear model sets (i.e., the NPCC and NCM) have a higher computation complexity compared to the linear model sets (i.e., the LPCC and LCM), due to the iteration procedures from the nonlinear compensation. Given the real-time computation speed of the Cosserat model was raised around 1 kHz (Till et al., 2015), future studies will also explore how to achieve real-time stiffness modelling and control computation based on the methodology proposed in this study.

Table 7.

Summary of the computational time.

	Model	Average time [s]	Time @ P max [s]
Simulation 1: Figure 6(a)	LPCC	0.0099	0.0087
	LCM	0.89	1.07
	NPCC	0.07	0.12
	NCM	2.98	5.03
Simulation 1: Figure 6(b)	LPCC	0.018	0.019
	LCM	1.06	1.26
	NPCC	0.12	0.20
	NCM	3.21	5.41
Simulation 2: Figure 8(a)	NCM	4.37	8.10
Simulation 2: Figure 8(b)	NCM	4.87	8.82
Simulation 2: Figure 8(c)	NCM	2.83	5.01
Simulation 2: Figure 9(a)	NCM	6.14	10.96
Simulation 2: Figure 9(b)	NCM	6.19	11.02
Simulation 3: Figure 10	NCM	4.87	5.03
Experiment 1: Figure 12(a)	LPCC	0.0095	0.0078
	LCM	1.37	1.72
	NPCC	0.12	0.20
	NCM	3.49	6.66
Experiment 1: Figure 12(b)	LPCC	0.011	0.0076
	LCM	0.92	0.98
	NPCC	0.10	0.18
	NCM	3.05	4.92
Experiment 1: Figure 12(c)	LPCC	0.015	0.012
	LCM	0.97	1.09
	NPCC	0.12	0.22
	NCM	3.08	5.11
Experiment 1: Figure 13(a)	LPCC	0.0092	0.0071
	LCM	0.84	1.04
	NPCC	0.10	0.18
	NCM	2.98	4.77
Experiment 1: Figure 13(b)	LPCC	0.011	0.0095
	LCM	0.89	1.11
	NPCC	0.12	0.20
	NCM	3.31	5.45
Experiment 2: Table 5	NCM	3.94	6.78
Experiment 3: Table 6	NCM	6.32	9.78

Time @ P_max denotes the computation time under the maximum pressure in a simulation or an experiment.

The nonlinear compensation in Section 2.3 is comparable to the nonlinear matrix structural analysis in Naselli and Mazzolai (2021), on the other hand, the iterative updates increase the computation complexity (see Table 7). As such, it is also worthy to incorporate the exact hyper-elastic models into the proposed framework (Marechal et al., 2021) to ease the computation. In addition, it is also encouraging to explore the research in the solid mechanics to study the influences of the cross-sectional deformation (Kumar and Mukherjee 2011).

In summary, the proposed stiffness modelling and analysis method can be used to quantitatively assess the stiffness/compliance response of soft robots under different configurations, which may provide new insights to explore on-demand stiffness analysis and optimisation of soft robots with regards to application requirements. Moreover, in serially connected soft robots, the proposed method could be used to optimise the configuration with high stiffness capability along a desired direction.

Supplemental Material

Supplemental Material

Supplemental Material