<figure><img src="image_1.png" /></figure>

HAL
open science

# Piecewise Linear Strain Cosserat Model for Soft Slender Manipulator

Haihong Li, Lingxiao Xun, Gang Zheng

► To cite this version:

Haihong Li, Lingxiao Xun, Gang Zheng. Piecewise Linear Strain Cosserat Model for Soft Slender Manipulator. IEEE Transactions on Robotics, 2023, 39 (3), pp.2342-2359. 10.1109/TRO.2023.3236942. hal-03913430

HAL Id: hal-03913430

https://inria.hal.science/hal-03913430

Submitted on 5 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Piecewise Linear Strain Cosserat Model for Soft Slender Manipulator

Haihong Li, Lingxiao Xun, Gang Zheng

***Abstract***—Recently soft robotics has rapidly become a novel and promising area of research with many designs and applications due to their flexible and compliant structure. However, it is more difficult to derive the nonlinear dynamic model of such soft robots. The differential kinematics and dynamics of the soft manipulator can be formulated as a set of highly nonlinear partial differential equations (PDEs) via the classic Cosserat rod theory. In this work, we propose a discrete modeling technique named piecewise linear strain (PLS) to solve the PDEs of Cosserat-based models, based on which the associated analytic models are deduced. To validate the accuracy of the proposed Cosserat model, the static model of the conical cantilever rod under gravity as a simple example is simulated by using different discretization methods. Results indicate that PLS Cosserat model is comparable to the mechanical deformation behavior of real-world soft manipulator. Finally, a parameters identification scheme for this model is established, and the simulation as well as experimental validation demonstrate that using this method can identify the model physical parameters with high accuracy.

***Index Terms***—Soft manipulator, dynamics, piecewise linear strain (PLS), Cosserat model, parameters identification

## I. INTRODUCTION

### A. Review of relevant literature

Soft robots are originally inspired by the structure and behavior of animal species such as octopus tentacles, elephant trunks, snakes and caterpillars [1], [2], which have become an emerging type of robots made of continuous slender elastic elements that can be idealized as one-dimensional objects. Slender elastic rods demonstrating large deflections are increasingly prevalent not only in continuum robots, but also for the applications of the soft robots [3], [4]. Despite their potential strengths, the deformability of soft robots produces an infinite degree-of-freedom and highly coupled nonlinear system that is arduous to model. The difficulty in modeling soft robots obviously results from falling short of tools in the field of continuum mechanics.

The beam/rod theory, a subclass of continuum mechanics, provides a theoretical guidance to model numerous problems in engineering due to its generality and simplicity. The classical Euler-Bernoulli beam theory has been a widely used approach for the mechanical modeling of robotic structures in the past years. For example, [5] developed a quasi-static bending model from a geometrically exact Euler-Bernoulli

formulation that generalized a sequence of soft arm designs to predict the result of design changes. Timoshenko beam model allows for the effect of shear deformation, which has been investigated and applied in soft robot modeling [6], [7]. As a geometrically nonlinear extension of the Euler-Bernoulli beam, the Kirchhoff rods conveniently describe bending and torsion deformations for general models of soft continuum robots [8], [9]. In applications with large deflections, classical rod theories in geometric nonlinearities are required. Many models employ the approximate nature of constant curvature (CC) representations which results in computationally simple mathematical models and motivates the use of generalized elastic models such as Cosserat beam theory [10], [11], [12], [13], [14]. As a matter of fact, the Cosserat rod is the geometrically nonlinear generalization of the Timoshenko-Reissner beam, and it can model not only bending and torsion of soft continuum robots, but as well shearing and extension. The dissertation about the problems of statics, dynamics, and stability for continuum robots with slender elastic rod on the basis of Cosserat was exhaustively presented in [15]. Based on the work of [15], the authors put forward a new implicit dynamics framework for the purpose of solving Cosserat partial differential equation (PDE) models and addressing the issue of computational difficulties [16]. In [17], a Cosserat geometrically exact dynamic model of a soft manipulator driven by cables was developed. [18] introduced a method based on the assembly of heterogeneous, active and passive Cosserat rods to model dynamic musculoskeletal architectures that can withstand all modes of deformation. For the work of [19], the authors proposed a dynamic Cosserat modeling approach based on the strain nonlinear parameterization. It is numerically simple with a good accuracy and less degrees of freedom (DoFs). However, this technique leads to computational complexities owing to high-order basis functions required when modeling quite complicated deformation such as global buckling behaviors. The pseudo-rigid body (PRB) 3R model was first established to account for the large deformation of the flexible beam subject to the tip load, and utilized to analyze the compliant mechanisms with high computational efficiency [20]. On the basis of the PRB 3R model, [21] proposed a three-dimensional (3D) static modeling method with high model precision for cable-driven continuum robots, with the aim of the extension from the PRB 3R model to 3D applications.

The rod modeling frameworks to describe soft continuum robots have been laid down in light of the above reviewed beam theories. However, the continuous rod models have infinite-dimension states, and as such, they are rather difficult

Haihong Li, Lingxiao Xun are in Defrost team, Inria, university of Lille, Centrale Lille, CRISTAL - Centre de Recherche en Informatique Signal et Automatique de Lille - UMR 9189, France, and Gang Zheng is with School of Mathematics and Big Data, Foshan University, Foshan 528000 China, and University of Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000 Lille, France (e-mail: haihong.li@inria.fr; lingxiao.xun@inria.fr; gang.zheng@inria.fr).
to model and control than the discrete systems with finite
DoFs. Specifically, modeling and control of soft robots for
some practical applications call for accurate and efficient
models that will enable simulation of mechanical behavior,
improve structural design and development of control systems.
The existing alternative to a finite-dimensional model of the
robot can be useful for control purposes. The mathematical
foundations from the rod theories facilitate an assortment
of discrete computational models. In practice, these discrete
approaches allow to significantly reduce the number of DoFs
of the calculated mechanical model, and simultaneously en-
sure the model precision. There are several typically discrete
modeling techniques for soft continuum robots, which mainly
contains the piecewise constant curvature (PCC) models [22],
finite element models/methods (FEM) [23], [24], [25], and the
piecewise constant strain (PCS) Cosserat models [26].

PCC modeling approach describes soft robots with a finite
number of arcs parameterized by three quantities (i.e., cur-
vature, arc length, and bending plane) [22], resulting in the
reduced-order and relatively simple mathematical models. This
method has been proven to be an excellent technique with
a wide range of applications [27], [28], [29], [30]. Despite
this success, when external loads are considered, the PCC
is inaccurate due to its strong piece-wise constant curvature
assumption for practical applications. Furthermore, the CC
assumption is valid only for the specific design and actuation
manner, and hence it can not be regarded as generic. These
flaws can potentially produce critical behaviors in the real-
world scenario.

FEM-based method which is formulated as ways of ap-
proximate solutions of PDEs has been used to model physical
behavior of soft robots [31], [32]. Nevertheless, to obtain exact
modeling precision, this method demands that the number
of nodes tends to infinity, which increases the dimension of
the system and leads to a higher computation complexity.
Moreover, FEM models for large-deformation 3D nonlinear
elasticity always involve unnecessary computational expense
when modeling long, slender arms like rod-based soft robots
because general deformations of the cross sections of the
model are included.

PCS Cosserat modeling technique [26] extends the PCC
assumption to more general rod kinematics while providing a
dynamic model. It employs a finite set of piecewise constant
strains with discontinuities happening at fixed points along
the rod to model the deformation of the soft robots, and
provides advantages to use with not only a few state variables,
but also a relatively high modeling precision. More recently,
the PCS Cosserat dynamic models have been increasingly
applied to soft robots [33], [34], [35]. With regards to the soft
robots modeled through the PCS approximation, the number of
sections to be divided depends on the designer of the model as
a result of application-specific considerations. This technique
is able to exactly approximate the continuum formulation
when the arm is divided into more sections, yet the dimension
of state variables and computational cost will highly increase.
Consequently, a trade-off between model order reduction and
accuracy should be established. From practical point of view,
the deformation field of any certain section under external

forces is obviously not constant along the soft arm. To this end,
the authors proposed a novel variable-strain parameterization
by the discretization of the continuous Cosserat rod model
onto a finite set of strain basis functions in order to generalize
the PCS method to the case of discrete non-constant strains
[36]. However, this approach has not been extended to the
multisection, dynamic case.

In the past three decades, numerous model parameters identification techniques have been proposed by researchers. Several major contributions for dynamic model identification and its applications in the robotics control have been reported in [37], [38]. These studies assumed that the kinematic models were accurate, however, some typical investigations presented uncertainties of kinematic models [39], [40]. The parameters identification of static model is more obtainable than the dynamic model in practical applications since it just needs information of joint position rather than joint velocity and acceleration. The contributors of [41] established a framework for the special actuator's modeling using Cosserat rods along with experimental data-driven estimates of the actuation and stiffness parameters. In [42], a material characterization approach for identifying the parameters for material models via the target application was proposed to improve their use in the modeling of soft continuum robots. For the work of [43], the authors introduced a procedure of the static model identification towards general serial articulated manipulator and presented an application of the identified static model to validate the proposed methods. In the literature, most of existing parameters identification methods based on Cosserat rod theory were proposed for specific robots, and there is a lack of general approach to identify physical parameters for soft robots based on the discrete Cosserat models.

The Cosserat rod model is the closest to the mechanics
of deformation of the soft robots since it can produce an
exact nonlinearity in the deformations due to bending, torsion,
extension, and shearing. In this paper, we first put forward
a piecewise linear parameterization of the rod shape by its
strain fields, namely piecewise linear strain (PLS) Cosserat
model. This discrete Cosserat technique, striking a balance
between accuracy, robustness and computational complexity,
aims to address the lack of rigorous modeling techniques in
soft robotics and develop a general framework between soft
and rigid robots. Apart from the proposed modeling method,
a systematic parameters identification framework based on the
PLS Cosserat model for soft manipulators with arbitrary actu-
ation manner is constructed, which is necessary and promising
towards the development of model-based controller design.

B. *Contributions*

The soft manipulator is modeled as a continuous assembling
of cross sections moving on the midline in the 3D space
via infinite rigid transformations dependent on the internal
deformations. The geometrical assumption of considering no
cross-sectional deformation allows to describe the soft arm by
using ideas from the Lie group structure of rigid-body motions
and the other concepts from the techniques of differential
geometry. In conclusion, the distinct contributions of this paper
are summarized as follows:
* Propose a discrete modeling approach named PLS for large elastic deformations via the Cosserat rod theory with application to soft manipulators. The analytically integrable models involving geometric, differential kinematic and dynamic models are obtained. Simultaneously, derivation of a novel strain mode selection scheme via PLS Cosserat is provided to model different simplified beam models.

* Present a physical parameter identification strategy for the discrete Cosserat model with arbitrary actuation manner by solving a nonlinear programming (NLP) problem. Experimental setup using a cable-driven soft manipulator made of silicone indicates that the proposed method can effectively get more accurate model physical parameters. The PLS Cosserat model with the identified material parameters is capable of predicting the position of end-effector of the arm with a small error.

* Implement the static simulation of proposed method for soft manipulator actuated by cables, validate the precision of PLS Cosserat model, and compare the systems with different modes.

## C. Outline

The remainder of the paper is structured as follows. In Section II, the continuum models of the soft arm via the Cosserat rod on Lie group are recalled. Section III presents the detailed mathematical derivation of the novel discrete Cosserat model approach and strain mode choice scheme via Cosserat. In Section IV, a model parameters identification method based on the discrete Cosserat model is proposed. In Section V, the comparison of accuracy and computational efficiency of the discrete Cosserat models with those of FEM is performed through a cantilever rod simulation. In Section VI, we design a conical soft silicone manipulator, and use it to implement the parameters identification and model validation. Finally, the conclusion and future work are made.

## II. DERIVATION OF COSSERAT ROD PDES

A Cosserat rod is a one-dimensional slender continuum deformable body, and each cross section of the rod is considered as an infinitesimal material element whose size and shape cannot change under external forces, as illustrated in orange part in Fig. 1. In general, all variables of the rod can be parameterized by the reference arc length $X \in [0, L] \subset \mathbb{R}$ along the undeformed rod and by the time $t \in \mathbb{R}$. In the interest of clarity, we use $(\cdot)'$ and $(\cdot)$ to represent the space and time partial derivatives with respect to $X$ and $t$, respectively. For ease of reference, the nomenclature we selected is summarized in Table I. With the aim to make the discussion of the discretization more concrete, we start by presenting the PDE system describing the Cosserat rod.

As depicted in Fig. 1, the position field of any cross section can be represented by a centerline vector $u(X, t) \in \mathbb{R}^3$ and an orthonormal rotation matrix $R(X, t) \in SO(3)$ with respect to the inertial frame at time $t$. Hence, the homogeneous transformation matrix can be then defined as $g(X, t) = (\begin{matrix} R & u \\ 0^\top & 1 \end{matrix}) \in$

<table><caption>TABLE I
Nomenclature and definitions</caption><thead><tr><th>Symbol</th><th>Unit</th><th>Definition</th></tr></thead><tbody><tr><td>X</td><td>m</td><td>Arc length.</td></tr><tr><td>t</td><td>s</td><td>Time.</td></tr><tr><td>R(X, t)</td><td>—</td><td>Rotation matrix.</td></tr><tr><td>u(X, t)</td><td>m</td><td>Position vector.</td></tr><tr><td>g(X, t)</td><td>—</td><td>The configuration matrix g(X, t) = (R<sub>0</sub><sup>⊤</sup> u<sub>1</sub>)</td></tr><tr><td>K(X, t)</td><td>1/m</td><td>Angular strain in the body frame.</td></tr><tr><td>Q(X, t)</td><td>—</td><td>Linear strain in the body frame.</td></tr><tr><td>Ω(X, t)</td><td>1/s</td><td>Angular velocity in the body frame.</td></tr><tr><td>V(X, t)</td><td>m/s</td><td>Linear velocity in the body frame.</td></tr><tr><td>⟨·⟩</td><td>—</td><td>Mapping from ℝ<sup>3</sup> to so(3), e.g. ã = [0 -a<sub>3</sub> a<sub>2</sub>]<br>Mapping from ℝ<sup>6</sup> to se(3), e.g. ξ̃ = (K̃ Q̃ 0 0), η̃ = (Ω̃ Ṽ 0 0) ∈ se(3)</td></tr><tr><td>dX</td><td>m</td><td>Infinitesimal material element.</td></tr><tr><td>ρ</td><td>kg/m<sup>3</sup></td><td>Density of material.</td></tr><tr><td>R(X)</td><td>m</td><td>Cross-sectional radius.</td></tr><tr><td>A(X)</td><td>m<sup>2</sup></td><td>Cross-sectional area.</td></tr><tr><td>E</td><td>Pa</td><td>Young's modulus.</td></tr><tr><td>ν</td><td>—</td><td>Poisson ratio.</td></tr><tr><td>G</td><td>Pa</td><td>Shear modulus (For the isotropic material, G = E/(2(1 + ν))).</td></tr><tr><td>μ</td><td>Pa·s</td><td>Viscosity modulus.</td></tr><tr><td>I<sub>3</sub></td><td>—</td><td>3 × 3 identity matrix.</td></tr><tr><td>J(X)</td><td>m<sup>4</sup></td><td>Second moment of area tensor in the body frame, J = [J<sub>x</sub> 0 0 0]<br>rod, J<sub>x</sub> = J<sub>y</sub> + J<sub>z</sub>, and J<sub>y</sub> = J<sub>z</sub> = πR<sup>4</sup>/4, J<sub>y</sub>, J<sub>z</sub> are separately the second moments of the area w.r.t. Y- and Z-axis, J<sub>x</sub> is the polar moment of the area around the X-axis). The adjoint map of the Lie algebra, e.g. ad<sub>ξ</sub> = (K̃ Q̃ 0<sub>3×3</sub> K̃) , ad<sub>η</sub> = (Ω̃ Ṽ 0<sub>3×3</sub> Ω̃)</td></tr><tr><td>M</td><td>—</td><td>Cross-sectional mass matrix.</td></tr><tr><td>K<sub>tb</sub></td><td>N·m<sup>2</sup></td><td>Stiffness matrix for torsion and bending, K<sub>tb</sub> = [G 0 0 0]<br>Stiffness matrix for elongation and shearing, K<sub>es</sub> = [E 0 0 0]<br>Damping matrix for torsion and bending, D<sub>tb</sub> = [0 0 0 0]<br>Damping matrix for elongation and shearing, D<sub>es</sub> = [0 0 0 0]<br>Generalized selection matrix.</td></tr><tr><td>Ad<sub>g</sub>(X)</td><td>—</td><td>Mapping from a matrix to a vector. The matrix transforming the velocity or acceleration twist from body frame to inertial frame, i.e., Ad<sub>g</sub>(X) = (R̃ uR̃ 0<sub>3×3</sub> R̃) ∈ ℝ<sup>6×6</sup>. The transformation matrix between the inertial frame and the manipulator base frame. Euclidean norm of a vector or matrix. Total number of sections divided. Number of tests for parameters identification. Lagrangian function. Local distance between the midline of soft rod and the cable. Unit vector tangent to the cable path.</td></tr><tr><td>g<sub>r</sub></td><td>—</td><td></td></tr><tr><td>||·||<sub>2</sub></td><td>—</td><td></td></tr><tr><td>N</td><td>—</td><td></td></tr><tr><td>L</td><td>—</td><td></td></tr><tr><td>d<sub>i</sub>(X)</td><td>m</td><td></td></tr><tr><td>t<sub>ci</sub>(X, t)</td><td>—</td><td></td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Geometric description of a Cosserat rod and its applications.</figcaption></figure>

SE(3). For the cross section at $X \in [0, L]$, we denote its strain by $\xi(X,t) = (K^\top, Q^\top)^\top \in \mathbb{R}^6$ where $K(X,t) \in \mathbb{R}^3$ and $Q(X,t) \in \mathbb{R}^3$ respectively represent the angular strain (bending and torsion) and linear strain (shearing and extension), and denote its velocity by $\eta(X,t) = (\Omega^\top, V^\top)^\top \in \mathbb{R}^6$ where $\Omega(X,t) \in \mathbb{R}^3$ and $V(X,t) \in \mathbb{R}^3$ respectively stand for the angular and linear velocity.

Based on the above notations, the corresponding kinematic model was derived in [35]:

$$g' = g\hat{\xi}, \quad \dot{g} = g\hat{\eta} \qquad (1)$$

$$\eta' = \dot{\xi}(X) - \mathrm{ad}_{\xi(X)}\eta(X) \qquad (2)$$

$$\dot{\eta}' = \ddot{\xi}(X) - \mathrm{ad}_{\dot{\xi}(X)}\eta(X) - \mathrm{ad}_{\xi(X)}\dot{\eta}(X) \qquad (3)$$

as well as the dynamical model defined below [44]:

$$\mathcal{M}\dot{\eta} - \mathrm{ad}_{\bar{\eta}}\mathcal{M}\eta = \mathcal{F}'_i - \mathrm{ad}_{\xi}\mathcal{F}_i + \bar{\mathcal{F}}_e \qquad (4) \\ \mathcal{F}_i = \mathcal{F}_{ie} + \mathcal{F}_{ia}$$

where $\mathcal{M} \in \mathbb{R}^{6 \times 6}$ is the cross-sectional mass matrix, which equals to $\mathcal{M} = \rho[\mathrm{diag}(J_x(X), J_y(X), J_z(X), A(X), A(X), A(X))]$. $\rho$ is the density of material; $J_y(X)$, $J_z(X)$ and $J_x$ are, respectively, the second moment of the area with respect to Y-axis, Z-axis, and the polar moment of the area around the X-axis; $A(X)$ is the cross-sectional area. $\mathrm{ad}_{\cdot}(\cdot)$ is the adjoint map of the Lie algebra; $\mathcal{F}_i \in \mathbb{R}^6$ represents the internal wrench; $\bar{\mathcal{F}}_e \in \mathbb{R}^6$ is distributed external wrench for unit of X; $\mathcal{F}_{ie}$ denotes internally elastic wrench, and $\mathcal{F}_{ia}$ stands for the internal wrench produced by the actuation. The above PDE (4) is defined with the boundary conditions (BCs) of the internal wrenches and configurations of the tip at $X = 0$ and $X = L$

$$\mathcal{F}_i(0) = -\mathcal{F}_{e0}, \quad \mathrm{or} \quad g(0) = g_0 \\ \mathcal{F}_i(L) = \mathcal{F}_{eL}, \quad \mathrm{or} \quad g(L) = g_L \qquad (5)$$

where $\mathcal{F}_{e0}$ and $\mathcal{F}_{eL}$ are tip external wrenches at $X = 0$ and $X = L$, respectively.

As for the internal elastic wrench, the Kelvin-Voigt model [45] can be adopted both for the elastic and viscous members because of the constitutive material behavior of the soft manipulator, i.e.,

$$\mathcal{F}_{ie}(X) = \Sigma(X)(\xi(X) - \xi_0) + \gamma(X)\dot{\xi} \qquad (6)$$

with

$$\Sigma(X) = \begin{bmatrix} K_{tb} & \\ & K_{es} \end{bmatrix}, \quad \gamma(X) = \begin{bmatrix} D_{tb} & \\ & D_{es} \end{bmatrix}$$

where $K_{tb} = \mathrm{diag}(GJ_x(X), EJ_y(X), EJ_z(X)) \in \mathbb{R}^{3 \times 3}$ and $K_{es} = \mathrm{diag}(EA(X), GA(X), GA(X)) \in \mathbb{R}^{3 \times 3}$ are stiffness matrices determined by the material properties and cross-sectional geometry; $D_{tb} = \mathrm{diag}(\mu J_x(X), 3\mu J_y(X), 3\mu J_z(X)) \in \mathbb{R}^{3 \times 3}$ and $D_{es} = \mathrm{diag}(3\mu A(X), \mu A(X), \mu A(X)) \in \mathbb{R}^{3 \times 3}$ are viscosity matrices for Kelvin-Voigt-type viscous damping, in which the scalar 3 represents the Trouton's ratio between shear and extension viscosity for incompressible Newtonian fluids and holds more generally for viscoelastic fluids in the limit of very small strain rates [45], and $\mu$ is the viscosity modulus; $\xi_0$ represents the undeformed configuration of the arm. Note that a rod cross-section does not require to be circular, and it is only needed that the studied manipulator is slender, which implies that its length is much larger than its radius [46].

The actuation wrench $\mathcal{F}_{ia}$ in (4) depends on the type of actuators used, and the most common actuation manners for soft arms are tendon, fluidic and pneumatic actuators. In terms of the external wrench $\bar{\mathcal{F}}_e$ in (4), we can consider distributed load produced by gravity, or point load exerted by external disturbance in accordance with the actual situation.

Generally, the solution of the strong form (4) can be approximated by using many different approaches, such as implicit finite difference method [47], the shooting method [16], assumed mode method [48] and so on. In [26], by introducing the virtual displacement $\delta\phi(X) \in \mathbb{R}^6$, the D'Alembert's principle is used to obtain the weak form of (4), then the length space $[0, L]$ was discretized into $N$ sections, and the strain is assumed to be constant for each section (i.e., PCS: piecewise constant strain) in order to deduce analytic formula. Clearly, such an assumption requires finely spatial discretization which yields relatively high-dimensional system. Such a method can provide enough precision and work well for numerical simulation which might take time. However, PCS will be a big issue when designing model-based controllers.

To reduce the dimension of the deduced system, [19] proposed to globally approximate the strain field as $\xi(X,t) = \Phi(X)q(t)$ where $\Phi(X) = (\Phi_1, \Phi_2, \dots, \Phi_n)$ defines $n$ basis functions to parameterize the strain space (i.e., VS: variable strain). The advantage of VS is that the dimension of the resulted dynamical model is quite low with respect to PCS, but can provide comparable precision if the number of basis function is high enough. However, the choice of the number of basis function is highly dependent on external disturbance. For example, the external disturbance such as concentrated loads will lead to local strain mutation, and a global VS approximation method will be difficult to guarantee the local fitting precision. Moreover, all the kinematics are no more analytically integrable in contrast to PCS, but numerically reconstructed with a quaternion-based integrator. In other words, the deduced model is not anymore analytic, and will not be very friendly for control design.

In this paper, we will combine the advantages of those two methods, i.e., two ideas: small discretization to keep the local approximation precision (PCS), and interpolation to decrease the dimension (VS), so as to propose a piecewise linear strain (PLS) approach which is capable of obtaining the analytic formula of the model and facilitating the control design.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Sketch of comparison among VS, PCS and PLS Cosserat models with real strain.</figcaption></figure>

### III. PIECEWISE LINEAR STRAIN COSSERAT MODEL

In this section, the geometric (or kinematic), differential kinematic and dynamic model of the PLS Cosserat rod will be developed. Afterwards, the rod kinematics is reduced by neglecting any component of the six internal degrees of freedom (DoFs) via the PLS Cosserat model, which can effectively model Euler-Bernoulli (E-B) beam, extensible Kirchhoff (E-K) rod and Timoshenko beam, etc.

#### A. Idea of PLS

Firstly, we would like to give a global picture to illustrate the difference between PCS, VS and PLS Cosserat modeling methods, and highlight the advantage of the PLS method.

Given an arbitrary strain field of soft body subject to external forces, as described in Fig. 2, the PCS, VS and PLS share the common objective: to approximate the real strain field as precise as possible. PCS (proposed in [26]) adopts a local approximation scheme with local constant strain assumption, while VS (proposed in [19]) uses a global approximation manner via the chosen basis functions (polynomials for example). Clearly, PCS will result in high-dimensional system in order to reach a certain precision while VS will suffer from lower precision near the local mutation region even the number of basis function is increased. Nevertheless, PLS combines both advantages of PCS and VS, and proposes to locally use a linear strain to approximate the real one. Such a scheme is able to locally approximate to the strain field of soft body compared to VS, and get much less DoFs than PCS. The sketch described in Fig. 2 intuitively demonstrates that PLS has advantages over those proposed approaches in some engineering applications, and thus provides significance for authors to carry out this work.

In what follows, we will focus on the PLS Cosserat modeling principle. First of all, it is noted that a “section” is defined as a unit block which is able to produce independent mechanical deformation while “segment” is a subset of one section. To put it differently, one section is made up of quite a few segments, as illustrated in Fig. 3. In view of the strain field $\xi(X)$ along the soft manipulator varying with the arc length $X$ at any moment $t$, we divide the whole soft arm into

$N$ sections in the form of $[0, L_1], [L_1, L_2] \cdots [L_{N-1}, L_N]$ (with $L_N = L$). Then, the continuous strain field $\xi(X)$ can be substituted by $N$ continuous strain fields of the form $\{\xi_1(X), \xi_2(X), \cdots, \xi_N(X)\}$. To develop this discrete Cosserat model, for any section $1 \le n \le N$, we make the following two assumptions.

* The strain twists $\bar{\xi}_{n-1}$ and $\bar{\xi}_n$ respectively corresponds to those at the proximal and distal ends of any section $n$, and other strain twists along the section $n$ linearly vary with $X$, as shown in Fig. 3(a). Based on this linear assumption, the principle of the piecewise linear strain (PLS) can be then formulated as

$$ \xi_n(X) = \bar{\xi}_{n-1} \frac{L_n - X}{L_n - L_{n-1}} + \bar{\xi}_n \frac{X - L_{n-1}}{L_n - L_{n-1}} $$

for $X \in [L_{n-1}, L_n]$.

* Since the geometric and differential kinematic models are still linear time-varying systems under the PLS assumption, the section $n$ is subdivided into $k$ ($k \in \mathbb{R}$) infinitesimal segments of the form $[L_{n-1}, L_{n-1} + \Delta X], [L_{n-1} + \Delta X, L_{n-1} + 2\Delta X], \cdots, [L_{n-1} + (k-1)\Delta X, L_{n-1} + k\Delta X]$. In this way, the strain twists $\xi_n(X)$ along the segment $j$ remain constant, i.e., $\xi_n(X) \equiv \xi_n(L_{n-1} + (j-1)\Delta X), X \in [L_{n-1} + (j-1)\Delta X, L_{n-1} + j\Delta X]$. $\Delta X$, called one segment, represents the infinitesimal distance between any adjacent cross sections, as displayed in Fig. 3(b).

**Remark 1:** When developing the proposed discrete Cosserat model, we have introduced the notions of sections $N$ and segments $k$. In general, the number of sections $N$ depends on the deformation variation of the soft manipulator under the external forces. If there exists the large deformation variation in the presence of external loads, $N$ should be larger to yield more precise approximation, which in turn increases the computational cost; small value of $N$ can also be used for the scenarios where the strains are not quickly spatial-varying with respect to $X$. Concerning the number of infinitesimal segments $k$, it is in fact not related to the deformation variation and can take a larger value to improve the integral accuracy which however will increase computational complexity. Thus, a balance between the precision and the complexity needs to be taken into account in practice. In the subsequent derivation of PLS models, it is assumed that these two variables are well pre-chosen according to the above mentioned criteria.

#### B. PLS Cosserat: Geometric Model

Based on the PLS assumptions, for any segment $j$ at time $t$, the system (1), (2) and (3) can be seen as liner time-invariant systems, thus can be analytically solved. In consequence, the initial configuration, velocity and acceleration of any segment $j$ depend on the rightmost counterparts of the previous segment ($j-1$) along the section $n$.

To guarantee the continuity, specifying the rightmost configuration $g(L_{n-1})$ of the section ($n-1$) as the initial value of the system (1) for section $n$, the rightmost configuration of any segment $j$ along the section $n$ at time $t$ can be recursively derived. Taking the rightmost configuration $g(L_{n-1} + j\Delta X)$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Schematic illustration of the PLS Cosserat model. (a) Soft arm divided into several sections. (b) Any section subdivided into quite a few segments.</figcaption></figure>

of the segment j as the initial counterpart of system (1) for segment (j+1), and considering piecewise linear strain along one section, the position and orientation of any cross section at X along the section n at time t analytically yields

$$ 
\begin{aligned}
\boldsymbol{g}(X) &= \boldsymbol{g}(L_{n-1}) \left( \prod_{i=0}^{j-1} e^{\Delta X \Theta_{ni}} \right) e^{(X - L_{n-1} - j \Delta X) \Theta_{nj}} \\
&\triangleq \boldsymbol{g}(L_{n-1}) \boldsymbol{g}_n(X)
\end{aligned}
\quad (7) $$

with $\Theta_{ni} = \alpha_{ni} \hat{\boldsymbol{\xi}}_{n-1} + \beta_{ni} \hat{\boldsymbol{\xi}}_n$, $\Theta_{nj} = \alpha_{nj} \hat{\boldsymbol{\xi}}_{n-1} + \beta_{nj} \hat{\boldsymbol{\xi}}_n$, where $\alpha_{ni} = 1 - \frac{i\Delta X}{L_n - L_{n-1}}$, $\beta_{ni} = \frac{i\Delta X}{L_n - L_{n-1}}$, $\alpha_{nj} = 1 - \frac{j\Delta X}{L_n - L_{n-1}}$, $\beta_{nj} = \frac{j\Delta X}{L_n - L_{n-1}}$, $g_n(X)$ stands for the position and orientation of any cross section at X along the section n w.r.t. the rightmost counterparts of the section (n-1).

Intuitively, compared to the PCS modeling method where $g(X)$ is only a function of one strain field, the model deduced from PLS in (7) depends on both $\hat{\boldsymbol{\xi}}_{n-1}$ and $\hat{\boldsymbol{\xi}}_n$ of a certain section n.

### C. PLS Cosserat: Differential Kinematic Model

The differential kinematics aims to find the mapping between the velocity twist along the manipulator and the time derivative of the deformation twists of the soft arm. Under PLS assumption, by integrating the system (2), we can then analytically obtain the velocity of each cross section at time t (see the detailed calculation of $\eta(X)$ in Appendix A).

According to Proposition 2.25 in [49] that $e^{ad\Theta^\vee}$ is equivalent to $Ad_e\Theta$, we know that the coefficient matrix of the velocity $\eta(L_{n-1})$ of the tip at $X = L_{n-1}$ for the section n is the exponential in the adjoint representation of the Lie group transformation $g_n(X)$, i.e.,

$$ e^{-(X-L_{n-1}-j\Delta X)\mathrm{ad}\Theta_{nj}^{\vee}} \left( \prod_{i=0}^{j-1} e^{-\mathrm{ad}\Delta X \Theta_{ni}^{\vee}} \right) = \mathrm{Ad}^{-1}_{g_n(X)} $$

Defining the coefficient matrices of $\dot{\hat{\boldsymbol{\xi}}}_{n-1}$ and $\dot{\hat{\boldsymbol{\xi}}}_n$ separately as $T_{g_n(X)}$ and $T_{g_n(X)}$, the analytical solution of the velocity in Appendix A is then given by

$$ \eta(X) = \mathrm{Ad}_{g_n(X)}^{-1} \eta(L_{n-1}) + T_{g_n(X)} \dot{\hat{\boldsymbol{\xi}}}_{n-1} + T_{g_n(X)} \dot{\hat{\boldsymbol{\xi}}}_n \quad (8) $$

with

$$ T_{g_{n1}(X)} = e^{-(X-L_{n-1}-j\Delta X)\mathrm{ad}\Theta_{nj}^{\vee}} \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X \Theta_{n\tau}^{\vee}} \right) \right. \\ \left. \qquad \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \alpha_{n(i-1)} \mathrm{d}s \right] \\ + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{nj}^{\vee}} \alpha_{nj} \mathrm{d}s, $$

$$ T_{g_{n2}(X)} = e^{-(X-L_{n-1}-j\Delta X)\mathrm{ad}\Theta_{nj}^{\vee}} \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X \Theta_{n\tau}^{\vee}} \right) \right. \\ \left. \qquad \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \beta_{n(i-1)} \mathrm{d}s \right] \\ + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{nj}^{\vee}} \beta_{nj} \mathrm{d}s, $$

where the symbol $\vee$ is an operator to map a matrix into a vector.

From (8), if we know the strain twists ($\hat{\boldsymbol{\xi}}_{n-1}$ and $\hat{\boldsymbol{\xi}}_n$) and strain twists rates ($\dot{\hat{\boldsymbol{\xi}}}_{n-1}$ and $\dot{\hat{\boldsymbol{\xi}}}_n$) of tips at $X = L_{n-1}$ and $X = L_n$, then the velocity of any cross section at X and time t along the section n can be recursively derived.

Due to the same reason, the acceleration $\dot{\eta}(X)$ of any cross section at X of the section n along the soft manipulator at time t can be analytically computed. Considering linearly variable strain twists along a certain section n and using the property of Lie algebra that $\mathrm{ad}_m n = -\mathrm{ad}_n m$ holds for any element m and n in the Lie algebra, the system (3) can be analytically solved, and the detailed derivation is given in Appendix A. For the sake of simplicity, this analytical solution is re-formulated below in the compact way as

$$ \dot{\eta}(X) = \mathrm{Ad}_{g_n(X)}^{-1} \dot{\eta}(L_{n-1}) + \mathrm{AD}_{g_{n1}(X)} \dot{\hat{\boldsymbol{\xi}}}_{n-1} \\ + \mathrm{AD}_{g_{n2}(X)} \dot{\hat{\boldsymbol{\xi}}}_n + T_{g_{n1}(X)} \ddot{\hat{\boldsymbol{\xi}}}_{n-1} + T_{g_{n2}(X)} \ddot{\hat{\boldsymbol{\xi}}}_n \quad (9) $$

with
$$
\begin{aligned}
\mathrm{AD}_{g_{n1}(X)} ={}& e^{-(X-L_{n-1}-j\Delta X)\mathrm{ad}_{\Theta_{n\tau}^{\vee}}} \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}_{\Delta X}\Theta_{n\tau}^{\vee}} \right) \right. \\
& \left. \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{(s-L_{n-1}-i\Delta X)\mathrm{ad}_{\Theta_{n(i-1)}^{\vee}}\mathrm{ad}_{\eta(s)}\alpha_{n(i-1)}\mathrm{d}s \right] \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}_{\Theta_{n,j}^{\vee}}\mathrm{ad}_{\eta(s)}\alpha_{nj}\mathrm{d}s,
\end{aligned}
$$

and

$$
\begin{aligned}
\mathrm{AD}_{g_{n2}(X)} ={}& e^{-(X-L_{n-1}-j\Delta X)\mathrm{ad}_{\Theta_{n\tau}^{\vee}}} \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}_{\Delta X}\Theta_{n\tau}^{\vee}} \right) \right. \\
& \left. \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{(s-L_{n-1}-i\Delta X)\mathrm{ad}_{\Theta_{n(i-1)}^{\vee}}\mathrm{ad}_{\eta(s)}\beta_{n(i-1)}\mathrm{d}s \right] \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}_{\Theta_{n,j}^{\vee}}\mathrm{ad}_{\eta(s)}\beta_{nj}\mathrm{d}s.
\end{aligned}
$$

Thus, we can use (9) to calculate the acceleration of all cross sections along the section $n$ if knowing strain twists ($\tilde{\xi}_{n-1}$ and $\tilde{\xi}_n$), strain twist rates ($\tilde{\xi}_{n-1}$ and $\tilde{\xi}_n$) and rates of strain twist rate ($\tilde{\xi}_{n-1}$ and $\tilde{\xi}_n$) of tips at $X = L_{n-1}$ and $X = L_n$.

The relation between the velocity twist $\eta(X)$ along the robot and the strain twists ($\tilde{\xi}_{n-1}$ and $\tilde{\xi}_n$), and another relation between the acceleration twist $\dot{\eta}(X)$ and the strain twists ($\tilde{\xi}_{n-1}$ and $\tilde{\xi}_n$) as well as the rate of strain twists ($\dot{\tilde{\xi}}_{n-1}$ and $\dot{\tilde{\xi}}_n$) of the tips for the section $n$ ought to be illustrated for the purpose of deriving the subsequent PLS Cosserat dynamic model. Applying (8) from the base to the tip for all cross sections along the soft manipulator in the chain, we can obtain the mapping as the geometric Jacobian which is an essential tool to describe the differential kinematics and dynamics of the PLS Cosserat model. Defining $S_{(\cdot)} \in \mathbb{R}^{6 \times 6}$ and $\dot{S}_{(\cdot)} \in \mathbb{R}^{6 \times 6}$ as the components of the Jacobian matrix and its partial derivative in time $t$, the structural form of Jacobian and its derivative over the length of the soft arm can be separately expressed as $J(\tilde{\xi}_0, \tilde{\xi}_1, \dots, \tilde{\xi}_{N-1}, \tilde{\xi}_N, X) = [S_0 \quad S_1 \quad S_2 \quad S_3 \dots S_N] \in \mathbb{R}^{6 \times 6(N+1)}$, and $\dot{J}(\tilde{\xi}_0, \tilde{\xi}_1, \dots, \tilde{\xi}_N, \dot{\tilde{\xi}}_0, \dot{\tilde{\xi}}_1, \dots, \dot{\tilde{\xi}}_{N-1}, \dot{\tilde{\xi}}_N, X) = [\dot{S}_0 \quad \dot{S}_1 \quad \dot{S}_2 \quad \dot{S}_3 \dots \dot{S}_N] \in \mathbb{R}^{6 \times 6(N+1)}$.

The geometric Jacobian represents the relation between the velocity twist along the arm and the time derivative of the joint position of the soft manipulator. The generalized joint position vector $q = [\tilde{\xi}_0^\top \quad \tilde{\xi}_1^\top \quad \tilde{\xi}_2^\top \quad \dots \quad \tilde{\xi}_N^\top]^\top \in \mathbb{R}^{6(N+1)}$ composed of strain twists of all linear interpolation nodes is introduced, and the manipulator like a cantilever rod with high ratio of length to the cross-sectional area selected as a candidate is fixed to a mobile base (i.e., the velocity twist of the fixed end $\eta(0) = 0$). As a result, the discrete model of velocity (8) is globally equivalent to

$$ \eta(X) = J(q, X)\dot{q} \quad (10) $$

where the analytical formula of $J(q, X)$ can be found in (34) of Appendix B. Finally, by taking the derivative of (10) with respect to time $t$, the acceleration twist $\dot{\eta}(X)$ arrives at

$$ \dot{\eta}(X) = J(q, X)\ddot{q} + \dot{J}(q, \dot{q}, X)\dot{q} \quad (11) $$

and the analytical formula of $\dot{J}(q, \dot{q}, X)$ is detailed in (35) of Appendix B. We would like to emphasize that, compared to

PCS method, the deduced Jacobian matrix and its derivative in time are significantly different.

## D. PLS Cosserat: Dynamic model

In order to deduce PLS Cosserat rod dynamic model corresponding to the PDE (4), we introduce the relation between virtual displacement and the state vector $\delta\phi(X) = J(q, X)\mathcal{P}\delta q_a$ based on the fact that the strain twist $\tilde{\xi}_N$ of the free end of the soft arm is constrained by the boundary condition. Substituting differential kinematics models (10) and (11) into (4), a nonlinear ODE can be then obtained. Considering that the resulting ODE holds for $\forall \delta q_a^\top \neq 0$, the dynamics for the PLS Cosserat model yields

$$
\begin{aligned}
& \left( \mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} \mathcal{M} J \mathrm{d}X \right) \dot{q} - \left[ \mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} \left( \mathrm{ad}_{J\dot{q}}^{\top} \mathcal{M} J - \mathcal{M} \dot{J} \right) \mathrm{d}X \right] \dot{q} \\
& = \mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} \left( \mathcal{F}'_{ie} - \mathrm{ad}_{\tilde{\xi}^{\top}} \mathcal{F}_{ie} \right) \mathrm{d}X + \mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} \bar{\mathcal{F}}_{e} \mathrm{d}X \\
& \qquad + \mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} \left( \mathcal{F}'_{ia} - \mathrm{ad}_{\tilde{\xi}^{\top}} \mathcal{F}_{ia} \right) \mathrm{d}X
\end{aligned}
\quad (16)
$$

with $\mathcal{P} = [\mathbf{I}_{6N \times 6N}]$, and $\delta q_a = [\delta \tilde{\xi}_0^\top \quad \delta \tilde{\xi}_1^\top \quad \dots \quad \delta \tilde{\xi}_{N-1}^\top]^\top \in \mathbb{R}^{6N}$.

Let us introduce the following coefficient matrices and wrenches from (16):

* $\mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} \mathcal{M} J \mathrm{d}X = M(q)$, the $6N \times 6(N+1)$ mass matrix.

* $-\mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} (\mathrm{ad}_{J\dot{q}}^{\top} \mathcal{M} J - \mathcal{M} \dot{J}) \mathrm{d}X = C(q, \dot{q})$, the $6N \times 6(N+1)$ Coriolis matrix.

* $\mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} (\mathcal{F}'_{ie} - \mathrm{ad}_{\tilde{\xi}^{\top}} \mathcal{F}_{ie}) \mathrm{d}X = F_i(q, \dot{q})$, the $6N \times 1$ internal wrench.

* $\mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} \bar{\mathcal{F}}_{e} \mathrm{d}X = F_e(q)$, the $6N \times 1$ external wrench.

* $\mathcal{P}^{\top} \int_{0}^{L_N} J^{\top} (\mathcal{F}'_{ia} - \mathrm{ad}_{\tilde{\xi}^{\top}} \mathcal{F}_{ia}) \mathrm{d}X = F_a(q)$, the $6N \times 1$ actuation wrench.

Thanks to the definition of the internal wrench and PLS assumption, the stiffness and viscosity matrices can be decoupled from $F_i(q, \dot{q})$ (see simplification details in Appendix C). The concise formulation can be then written as

$$ F_i(q, \dot{q}) = K(q)(q - q_0) + D(q)\dot{q} \quad (17) $$

where $K(q) \in \mathbb{R}^{6N \times 6(N+1)}$ stands for the stiffness matrix, $D(q) \in \mathbb{R}^{6N \times 6(N+1)}$ is the viscosity matrix, and $q_0 = [\xi_{00}^\top \quad \xi_{10}^\top \quad \xi_{20}^\top \quad \dots \quad \xi_{N0}^\top]^\top \in \mathbb{R}^{6(N+1)}$ represents the undeformed configuration of the manipulator.

Aside from (16), the boundary condition in (5) must be considered. Therefore, the generalized PLS Cosserat dynamic system can be given by

$$
\begin{aligned}
M(q)\ddot{q} + C(q, \dot{q})\dot{q} - K(q)(q - q_0) - D(q)\dot{q} - F_e(q) = F_a(q) \\
\Gamma\dot{q} + \sigma(q - q_0) - \mathcal{F}_e(L_N) = -\mathcal{F}_{ia}(L_N)
\end{aligned}
\quad (18)
$$

where the second equation in (18) is derived from the boundary condition (5) and the constitutive law (6), with $\Gamma = [0_{6 \times 6N}, \gamma(L_N)]$, and $\sigma = [0_{6 \times 6N}, \Sigma(L_N)]$.

We would like to emphasize that the actuation wrench depends on the type of actuators used and the distribution of the chosen actuators on the soft arm. Taking the soft arm
driven by cables for example, the non-parallel (spiral) cables distribution manner can create torsion and rotation motion.

Compared to PCS dynamic model where the mass matrix $M(q)$ is square, the deduced PLS dynamic model has a non-square $M(q)$ which is not positive definite. However, by complementing with the boundary condition (5) and the constitutive law (6), a similar Lagrangian model can be obtained. If the viscosity is not considered in (6), i.e., $\gamma = 0$, then it leads to an algebraic equation in (18), which in fact is a differential-algebraic system. This characteristic, distinguished with the PCS and VS methods, is mainly due to the PLS assumption.

### E. Strain Mode Choice Scheme via the PLS Cosserat

According to the Cosserat rod theory, the strain twist of any interpolation node $\bar{\xi}_i$ could take any value in the six dimensional components. However, the internal rod kinematics describing the motions between the cross sections can be constrained by some restrictions in the usual applications of soft robots. To tackle this restriction, by following the similar idea of [19], we decompose the strain twist of any interpolation node $i$, by preserving the same variable order, as

$$ \bar{\xi}_i = \mathbf{B}_a \xi_{ia}^* + \mathbf{B}_e \xi_{ic}^* \quad (19) $$

where $\xi_{ia}^* \in \mathbb{R}^{n_i}$ determines the vector field of the free strains (i.e., the number of DoFs $n_i$) of the interpolation nodes allowed by the rod kinematics, $\xi_{ic}^* \in \mathbb{R}^{6-n_i}$ represents the vector field of constrained strains, $\mathbf{B}_a$ and $\mathbf{B}_c$ stand for the complementary selection matrix of 1 and 0 such that $\mathbf{B}_a^\top \mathbf{B}_a = \mathbf{I}_{(n_i)\times(n_i)}$, $\mathbf{B}_c^\top \mathbf{B}_c = \mathbf{I}_{(6-n_i)\times(6-n_i)}$, and $\mathbf{B}_a^\top \mathbf{B}_c = 0$. Since the variable order in $\bar{\xi}_i$ is preserved during the decomposition, thus the matrices $\mathbf{B}_a$ and $\mathbf{B}_c$ are uniquely determined.

Substituting (19) into (7), the reduced geometric model can be obtained. For the constrained soft manipulator, it is necessary to consider the strain states of all interpolation nodes. Hence, the generalized joint position vector can be expressed as

$$ q = \bar{\mathbf{B}}_a \bar{\mathbf{q}} + \bar{\mathbf{B}}_c \underline{\mathbf{q}} \quad (20) $$

with $\bar{\mathbf{B}}_a = \mathbf{I}_{(N+1)\times(N+1)} \otimes \mathbf{B}_a$, $\bar{\mathbf{B}}_c = \mathbf{I}_{(N+1)\times(N+1)} \otimes \mathbf{B}_c$, where $\otimes$ represents the Kronecker product. In such a way, $\bar{\mathbf{B}}_a \in \mathbb{R}^{6(N+1)\times[n_i(N+1)]}$ is the generalized selection matrix for the allowed states, $\bar{\mathbf{q}} = [\xi_{0a}^* \ \xi_{1a}^* \ \xi_{2a}^* \ \dots \ \xi_{Na}^*] \in \mathbb{R}^{n_i(N+1)}$ includes the allowed DoFs along the arm.

For the PLS Cosserat model reduction, the internal elastic wrench should be divided into two parts: one is constrained wrench in charge of imposing the internal constraints for prohibited strains, and another is the elastic wrench related to the allowed DoFs. Thus, the reduced internal wrench is given by

$$ \mathcal{F}_{ie} = \underbrace{\Sigma(X)\mathbf{B}_a(\xi_a^*(X) - \xi_{i0}^*) + \gamma(X)\mathbf{B}_a\dot{\xi}_a^*(X)}_{\mathcal{F}_{ie}^*} + \mathbf{B}_c\lambda(X) \quad (21) $$

where $\xi_a^*(X)$ for $X \in [L_{n-1}, L_n]$ can be obtained by the linear interpolation of allowed strains of the adjacent nodes, $\xi_{i0}^*$ is the undeformed states of allowed strains, $\lambda(X) \in \mathbb{R}^{(6-n_i)}$ is the constrained wrench.

Inserting (20) into (10) and (11), the reduced PLS Cosserat differential kinematic models yields

$$ \begin{aligned} \eta(X) &= \bar{\mathbf{J}}(\bar{\mathbf{q}}, X) \dot{\bar{\mathbf{q}}} \\ \dot{\eta}(X) &= \bar{\mathbf{J}}(\bar{\mathbf{q}}, X) \ddot{\bar{\mathbf{q}}} + \dot{\bar{\mathbf{J}}}(\bar{\mathbf{q}}, \dot{\bar{\mathbf{q}}}, X) \dot{\dot{\bar{\mathbf{q}}}} \end{aligned} $$

where $\bar{\mathbf{J}}(\bar{\mathbf{q}}, X) = \mathbf{J}(q, X)\bar{\mathbf{B}}_a \in \mathbb{R}^{6 \times [n_i(N+1)]}$ is the reduced body Jacobian matrix. Using the relation $\delta\phi(X) = \bar{\mathbf{J}}(\bar{\mathbf{q}}, X)\bar{\mathbf{P}}\delta\bar{\mathbf{q}}_a$, and substituting the reduced kinematics relations as well as (21) into (4) lead to

$$ \begin{aligned} & \underbrace{\left(\bar{\mathbf{P}}^{\top} \int_{0}^{L_N} \bar{\mathbf{J}}^{\top} \mathcal{M} \bar{\mathbf{J}} \mathrm{d}X\right)}_{\bar{\mathbf{M}}(\bar{\mathbf{q}})} \ddot{\bar{\mathbf{q}}} - \underbrace{\left[\bar{\mathbf{P}}^{\top} \int_{0}^{L_N} \bar{\mathbf{J}}^{\top} \left(\mathrm{ad}_{\bar{\mathbf{J}}^{\top}} \mathcal{M} \bar{\mathbf{J}} - \mathcal{M} \bar{\mathbf{J}}^{\top}\right) \mathrm{d}X\right]}_{\bar{\mathcal{C}}(\bar{\mathbf{q}}, \dot{\bar{\mathbf{q}}})} \dot{\bar{\mathbf{q}}}= \bar{\mathbf{P}}^{\top} \int_{0}^{L_N} \bar{\mathbf{J}}^{\top} \left(\mathcal{F}_{ie}^{*'} - \mathrm{ad}_{\dot{\xi}^*} \mathcal{F}_{ie}^*\right) \mathrm{d}X + \bar{\mathbf{P}}^{\top} \int_{0}^{L_N} \bar{\mathbf{J}}^{\top} \bar{\mathbf{F}}_e \mathrm{d}X \\ & + \bar{\mathbf{P}}^{\top} \int_{0}^{L_N} \bar{\mathbf{J}}^{\top} \left(\mathcal{F}_{ia}^{*'} - \mathrm{ad}_{\dot{\xi}^*} \mathcal{F}_{ia}^*\right) \mathrm{d}X \\ & + \bar{\mathbf{P}}^{\top} \int_{0}^{L_N} \bar{\mathbf{J}}^{\top} \left(\mathbf{B}_c \lambda' - \mathrm{ad}_{\dot{\xi}^*} \mathbf{B}_c \lambda\right) \mathrm{d}X \end{aligned} \quad (22) $$

where $\bar{\mathbf{P}} \in \mathbb{R}^{n_i(N+1) \times n_i N}$ represents a reduced selection matrix, $\bar{\mathbf{q}}_a = [\xi_{0a}^* \ \xi_{1a}^* \ \xi_{2a}^* \ \dots \ \xi_{(N-1)a}^*] \in \mathbb{R}^{n_i N}$ is composed of allowed DoFs of all strain nodes except those of the free end.

In addition, it is worth noting that the reduced boundary condition should be formulated as

$$ \begin{aligned} & \mathbf{B}_a^\top \Sigma(L_N) \mathbf{B}_a (\xi_{Na}^* - \xi_{N0}^*) + \mathbf{B}_a^\top \gamma(L_N) \mathbf{B}_a \dot{\xi}_{Na}^* \\ & = \mathbf{B}_a^\top (-\mathcal{F}_{ia}(L_N) + \mathcal{F}_e(L_N)) \end{aligned} \quad (23a) $$

$$ \lambda(L_N) = \mathbf{B}_c^\top (-\mathcal{F}_{ia}(L_N) + \mathcal{F}_e(L_N)) \quad (23b) $$

It must be also pointed out that all the items except the last one in (22) can be obtained by replacing the Jacobian matrix $\mathbf{J}$ in (16) with $\mathbf{J}\bar{\mathbf{B}}_a$. The following part concerns the calculation of the item with the constrained wrench.

**Theorem 1:** For the PLS Cosserat model with full modes, if the strain field is re-formulated as $\xi(X) = \xi_0 + \Phi(X)q(t)$, where $\Phi(X) \in \mathbb{R}^{6 \times (N+1)}$ is a generalized matrix comprised of coefficient matrix of strain interpolation nodes via the PLS assumption, then there exists a relationship among three quantities (i.e., $\bar{\mathbf{B}}_a$, $\Phi$, and $\mathbf{B}_c$) satisfying the following equality

$$ \bar{\mathbf{B}}_a^\top \mathbf{\Phi}^\top \mathbf{B}_c = \begin{bmatrix} a_1(X)\mathbf{B}_a^\top \\ (b_1(X) + a_2(X))\mathbf{B}_a^\top \\ \vdots \\ (b_{n-1}(X) + a_n(X))\mathbf{B}_a^\top \\ \vdots \\ b_N(X)\mathbf{B}_a^\top \end{bmatrix} \mathbf{B}_c = 0 \quad (24) $$

with $a_n(X) = \frac{L_n - X}{L_n - L_{n-1}}$, and $b_n(X) = \frac{X - L_{n-1}}{L_n - L_{n-1}}$, then the generalized constrained wrench $\bar{\mathbf{F}}_\lambda$ for the reduced PLS Cosserat relates to the constrained wrench of end cross section.

**Proof 1:** Note that we want to prove that $\bar{\mathbf{F}}_\lambda$ is only dependent on the constrained wrench at $X = L_N$. At this aim, let us insert (10) into (2), and it yields

$$ J'\dot{\underline{q}} = -\mathrm{ad}_{\dot{\xi}} J\dot{\underline{q}} + \dot{\xi}(X) = -\mathrm{ad}_{\dot{\xi}} J\dot{\underline{q}} + \Phi(X)\dot{\underline{q}} $$
which holds for $\forall \dot{q} \neq 0$, and thus

$$ \bar{J}' = -\mathrm{ad}_{\xi} \bar{J} + \Phi(X) \bar{B}_{\alpha} \quad (25) $$

Substituting (25) into $\bar{F}_{\lambda}$ with unknown constrained wrench $\lambda(X)$ in (22), then it arrives at

$$ \begin{aligned} \bar{F}_{\lambda} &= \bar{P}^{\top} \int_{0}^{L_N} [ (\bar{J}^{\top} \mathbf{B}_c \lambda)' - \bar{J}'^{\top} \mathbf{B}_c \lambda - \bar{J}^{\top} \mathrm{ad}_{\xi}^{\top} \mathbf{B}_c \lambda ] \mathrm{d}X \\ &= \bar{P}^{\top} \int_{0}^{L_N} [ (\bar{J}^{\top} \mathbf{B}_c \lambda)' + \bar{J}^{\top} \mathrm{ad}_{\xi}^{\top} \mathbf{B}_c \lambda - \bar{B}_a^{\top} \Phi^{\top} \mathbf{B}_c \lambda - \bar{J}^{\top} \mathrm{ad}_{\xi}^{\top} \mathbf{B}_c \lambda ] \mathrm{d}X \\ &= \bar{P}^{\top} \int_{0}^{L_N} (\bar{J}^{\top} \mathbf{B}_c \lambda)' \mathrm{d}X - \bar{P}^{\top} \int_{0}^{L_N} \bar{B}_a^{\top} \Phi^{\top} \mathbf{B}_c \lambda \mathrm{d}X \end{aligned} $$

Clearly, the second item in above equation can be removed in accordance with (24). Consequently, by using (23b), we obtain

$$ \bar{F}_{\lambda} = \bar{P}^{\top} (\bar{J}^{\top} \mathbf{B}_c \lambda)|_{0}^{L_N} = \bar{P}^{\top} \bar{J}^{\top} (L_N) \mathbf{B}_c \bar{B}_c^{\top} ( - \mathcal{F}_{i\alpha}(L_N) + \mathcal{F}_{e}(L_N) ) \quad \blacksquare $$

By combining (22) and (23), one can easily obtain the reduced model based on PLS Cosserat to model several simplified systems ($n_i \le 6$). An exhaustive reference of the reduced systems with the complementary selection matrices $\mathbf{B}_a$ and $\mathbf{B}_c$ for describing different systems is shown in Appendix D.

## IV. MODEL PARAMETERS IDENTIFICATION

In practice, the exact values of the physical parameters are typically unknown or difficult to derive even for the soft robot manufacturers. Besides, even though there is full knowledge of the model and sufficient data available, an accurate description is most often not possible. Therefore, it is necessary and important to identify the physical parameters for a soft robot conveniently and accurately. To reach this goal, we present an efficient algorithm framework to identify the parameters involved in the deduced PLS dynamic model.

Technologically speaking, the joint position vector $q$ is actually difficult to measure in the experiment. However, we can easily determine the position of the end-effector by using the position sensor. Therefore, the proposed identification scheme is based on the measurement of end-effector position. To this aim, it is assumed that $\bar{N}$ sets of different experiments are effectuated and the objective is to seek optimal parameters to minimize the difference between the real measured end-effector position and that obtained from simulation, by satisfying of course the PLS Cosserat static model. Consequently, the parameters identification algorithm can be formulated as the following nonlinear programming (NLP) problem:

$$ \begin{array}{ll} \underset{\delta=(\theta, q_1, q_2, \dots, q_{\bar{N}})}{\arg\min} & f(\delta) = \sum_{i=1}^{\bar{N}} \|u_i - u_{ei}\|_2^2 \\ \text{s.t.} & \begin{cases} \mathbb{K}(\theta, q_1) - \mathbb{F}_a(q_1) = \mathbf{0} \\ \mathbb{K}(\theta, q_2) - \mathbb{F}_a(q_2) = \mathbf{0} \\ \vdots \\ \mathbb{K}(\theta, q_{\bar{N}}) - \mathbb{F}_a(q_{\bar{N}}) = \mathbf{0} \end{cases} \end{array} \quad (26) $$

with

$$ \begin{aligned} \mathbb{K}(\theta, q_i) &= \begin{bmatrix} -\mathbf{K}(q_i) \\ \sigma \end{bmatrix} (q_i - q_{0i}) - \begin{bmatrix} \mathbf{F}_e(q_i) \\ \mathcal{F}_e(L_N) \end{bmatrix}, \\ \mathbb{F}_a(q_i) &= \begin{bmatrix} \mathbf{F}_a(q_i) \\ -\mathcal{F}_{i\alpha}(L_N) \end{bmatrix} \end{aligned} $$

where $u_i = Wg(q_i, L)\Upsilon$ implies the end-effector position provided by the PLS static model, with $W = [I_3 \quad 0]$ and $\Upsilon = [0_3 \quad 1]^{\top}$ in the $i^{th}$ experiment, $g(q_i, L)$ stands for the position and orientation of the end-effector in the $i^{th}$ experiment, with $q_i$ being the strain vector, $\theta$ represents those parameters to be identified, $u_{ei}$ is the $i^{th}$ experimental measurement of the end-effector position by using the sensor.

**Remark 2:** It is worth noting that the number of sections $N$ is pre-chosen in the proposed optimization algorithm since we are interested in identifying the material-related parameter, such as Young's modulus. However, the value of $N$ might be also treated as model parameter to be optimized for the proposed PLS Cosserat model. This will in fact lead to a nonlinear parameter identification problem, which can be handled via a two-stage method. The basic idea of such an approach is to fix firstly the value of $N$ and to solve the optimization problem (26) in the first stage. Once the objective function can not be minimized within the prescribed tolerance, then in the second stage we can set $N$ to $N+1$ and repeat the first stage until the optimal number of the sections $N$ is found to minimize the cost function by well satisfying the tolerance.

In our work, the parameter vector $\theta$ consists of the Young's modulus $E$, shear modulus $G$ and density of material $\rho$ to be identified, noted as $\theta = [E \quad G \quad \rho]$. To solve the proposed optimization problem (26), the Newton-type method is used by attempting to find the optimal solution $\delta^*$ which can generally satisfy the Karush-Kuhn-Tucker (KKT) conditions that there exists multiplier vector $\bar{\lambda}^* \in \mathbb{R}^{6\bar{N}(N+1)}$ such that the following equations hold:

$$ \begin{aligned} \nabla_{\delta} \mathcal{L}(\delta^*, \bar{\lambda}^*) &= \mathbf{0} \\ \bar{\lambda}^* &\neq \mathbf{0} \\ \mathcal{H}(\delta^*) &= \mathbf{0} \end{aligned} $$

with

$$ \mathcal{L} = f(\delta) + \mathcal{H}(\delta)^{\top} \bar{\lambda}, $$

and

$$ \mathcal{H}(\delta) = [\mathcal{H}_1(\theta, q_1)^{\top} \quad \mathcal{H}_2(\theta, q_2)^{\top} \quad \dots \quad \mathcal{H}_{\bar{N}}(\theta, q_{\bar{N}})^{\top}]^{\top} $$

where $\mathcal{H}_i(\theta, q_i) = \mathbb{K}(\theta, q_i) - \mathbb{F}_a(q_i)$ is the static model in the $i^{th}$ experiment, $\mathcal{L}$ represents the Lagrange function, and $\bar{\lambda}$ is the Lagrangian multiplier vector.

**Remark 3:** The Newton-type algorithm is sensitive to the choice of initial guess. Note that if a minimum exists, it is not necessarily unique. In other words, there may be multiple feasible solutions that meet the KKT conditions, known as local minima. However, regardless of the number of local minima, there always exists a unique optimal solution (if it exists). Therefore, to obtain the optimal material parameters of soft manipulator accurately and efficiently, the initial variable $\delta_0$ composed of material parameters and configuration of the manipulator should be carefully chosen. In general, it is recommended to start the iteration with estimates that are close to the true parameter values. Following this criterion, we refer to the initial guesses of material parameters (i.e., Young's modulus $E$, shear modulus $G$ and density $\rho$) of the soft manipulator provided by manufacturers, and choose the undeformed reference shape of the arm as initial configuration.
<figure><img src="image_6.png" /><figcaption>Fig. 4. The side view of the PLS Cosserat model with three continuum sections actuated by four cables.</figcaption></figure>

When the material parameters are not available from robot manufacturers, we can guess the initial values of these parameters within an acceptable interval in accordance with the properties of the material, and then allow the algorithm to run multiple times in order to determine the optimality of the solution for this optimization problem.

## V. SIMULATION COMPARISON OF DISCRETE COSSERAT MODELS

This section is devoted to validating the precision of the proposed PLS Cosserat models, by comparing it with the result obtained via finite element method (FEM). In addition, since PLS Cosserat shares the same local approximation scheme with PCS Cosserat, we will compare as well the precision of PLS model and PCS model.

### A. Simulation Setup

The comparison is effectuated by simulating a cantilever rod under external forces (for example, under gravity). The simulated rod is of conical shape and actuated by cables (see Fig. 4), with total length $L = 0.20$ m, base radius $R_{max} = 1 \times 10^{-2}$ m, tip radius $R_{min} = 5 \times 10^{-3}$ m, Young's modulus $E = 1.1 \times 10^5$ Pa, shear modulus $G = 3.793 \times 10^4$ Pa, and density of material $\rho = 2000$ kg/m$^3$. $\xi_0 = [0, 0, 0, 1, 0, 0]^T$ represents the undeformed straight configuration when the rod is stress-free. Besides, the rod shares the X-axis, Y-axis and Z-axis with the inertial frame, and thus the map $g_r$ between the base frame of the rod and the inertial frame is a $4 \times 4$ identity matrix.

Several cables are attached at the free end of the rod, and parallel to the surface of soft rod to produce the maximum torque and simultaneously to reduce the cables' friction, as shown in Fig. 4. Thus, the friction of the cables can be neglected. Additionally, the local distance between the midline of soft rod and the cable $i \in [1, 4]$ at the cross section X is defined as $d_i(X) \in \mathbb{R}^3$, we can then obtain the cable position vector in the inertial frame $u_{ci} = u + Rd_i$. By taking the derivative of cable position vector w.r.t. arc length and

normalizing, the unit vector $t_{ci}(X, t) \in \mathbb{R}^3$ tangent to cable path yields [17]:

$$ 
\begin{aligned}
t_{ci}(X,t) &= \frac{\mathbf{R}^{-1} \mathbf{u}'_{ci}}{\|\mathbf{R}^{-1} \mathbf{u}'_{ci}\|} = \frac{\mathbf{R}^{-1}(\mathbf{u}' + \mathbf{R}'\mathbf{d}_i + \mathbf{R}\mathbf{d}'_i)}{\|\mathbf{R}^{-1}(\mathbf{u}' + \mathbf{R}'\mathbf{d}_i + \mathbf{R}\mathbf{d}'_i)\|} \\
&= \frac{\mathbf{Q} + \mathbf{K} \times \mathbf{d}_i + \mathbf{d}'_i}{\|\mathbf{Q} + \mathbf{K} \times \mathbf{d}_i + \mathbf{d}'_i\|_2}
\end{aligned}
$$

The distance $d_i(X)$ is fixed once the soft manipulator is designed. Accordingly, the unit tangent vector $t_{ci}$ hinges on the angular strain $K(X)$ and linear strain $Q(X)$ of the soft manipulator. The actuation wrench $\mathcal{F}_{ia}(X)$ for unit of $X$ can be obtained by calculating the torque and force exerted by the cables, and it is given by

$$ 
\begin{aligned}
\mathcal{F}_{ia}(X) &= \begin{bmatrix} d_1 \times t_{c1} & d_2 \times t_{c2} & d_3 \times t_{c3} & d_4 \times t_{c4} \\ t_{c1} & t_{c2} & t_{c3} & t_{c4} \end{bmatrix} T \\
&= \Lambda(X)T
\end{aligned}
\quad (27) $$

where $\Lambda(X) \in \mathbb{R}^{6 \times 4}$ is a matrix function whose columns are composed of vector functions, and $T \in \mathbb{R}^4$ is the system input vector consisted of magnitude of all cables' tension.

With regards to the external wrenches, the case under gravity is considered in our test for the sake of simplicity, which can be then expressed as

$$ \bar{\mathcal{F}}_e(X) = \mathcal{M} \mathrm{Ad}_{g_r}^{-1} \mathcal{G} \quad (28) $$

where the inverse of adjoint representation of the Lie group is used to transform twists from inertial frame to body frame, $\mathcal{G}$ is the gravity acceleration twist w.r.t. the inertial frame, equal to $\mathcal{G} = [0, 0, 0, 0, 0, -9.81]^T$.

By the selection of cable actuation manner, the strain twist at $X = L_N$ satisfies the following differential equation

$$ \Sigma(L_N)(\bar{\xi}_N - \xi_{N0}) + \gamma(L_N)\dot{\bar{\xi}}_N - \mathcal{F}_e(L_N) = -\Lambda(L_N)T \quad (29) $$

Substituting (28) together with (29) into (18), the PLS Cosserat dynamic system driven by cables can be given by

$$ 
\begin{aligned}
& \begin{bmatrix} M(q) \\ \Gamma \end{bmatrix} \ddot{q} + \begin{bmatrix} C(q, \dot{q}) \\ \sigma \end{bmatrix} \dot{q} - \begin{bmatrix} \mathcal{G}(q) \mathrm{Ad}_{g_r}^{-1} \mathcal{G} + F_i(q, \dot{q}) \\ \mathcal{F}_e(L_N) \end{bmatrix} \\
& = \begin{bmatrix} H(q) \\ 0_{6 \times 4} \end{bmatrix} T + \begin{bmatrix} 0_{6N \times 4} \\ -\Lambda(L_N) \end{bmatrix} \dot{T}
\end{aligned}
\quad (30) $$

with

$$ \mathcal{G}(q) = \mathcal{P}^{\top} \int_{0}^{L_N} \mathcal{J}^{\top} \mathcal{M} \mathrm{Ad}_{g_r}^{-1} \mathcal{G} \mathrm{d}X, $$

and

$$ \mathcal{H}(q) = \mathcal{P}^{\top} \int_{0}^{L_N} \mathcal{J}^{\top} \left( \Lambda' - \mathrm{ad}_{\xi(X)}^{\top} \Lambda \right) \mathrm{d}X $$

Redefining each term and naming the coefficient matrices of (30), the generalized PLS Cosserat dynamic model for numerical simulation can be then formulated as the following general Lagrangian structural form

$$ \underline{M}(q)\ddot{q} + \underline{C}(q, \dot{q})\dot{q} - \underline{K}(q, \dot{q}) = \underline{H}(q)T + \overline{H}\dot{T} \quad (31) $$

with the initial conditions at $t = t_0$

$$ q(t_0) = [\bar{\xi}_0^\top \quad \bar{\xi}_1^\top \quad \dots \quad \bar{\xi}_N^\top]^\top $$

and

$$ \dot{q}(t_0) = [\dot{\bar{\xi}}_0^\top \quad \dots \quad \dot{\bar{\xi}}_{N-1}^\top \quad (\mathcal{F}_e - \Lambda T - \Sigma(\bar{\xi}_N - \xi_{N0}))^\top \gamma^{-1}]^\top $$
where $\underline{M}(q)$ is a symmetric positive-definite mass matrix, $\underline{H}(q)$ represents the generalized actuation matrix. It is worth noting that continuous Cosserat rod dynamic formulation of soft robots leads to a PDE in the form of a boundary value problem (BVP), however, the PLS Cosserat dynamic model actuated by cables takes the form of a nonlinear ODE, which paves the way to design the model-based controllers.

## B. Comparison of Accuracy for Static Models

FEM can be regarded as a reliable tool (and reference) to verify the modeling result by selecting finer mesh, since it discretizes the space in a very generic manner without introducing restrictive assumptions. Hence, this paper used the FEM approach as an alternative of the real soft arm. In short, the main idea of FEM is to discretize the geometric shape of the rod by employing finite number of finer elements. The geometric model of the rod is established in the SolidWorks, we then use the FEM in Comsol which is a general FEM software to obtain the equilibrium position of the cantilever rod. Specifically, in terms of spatial discretization of the studied cantilever rod, quadrilateral mesh elements are used, and the rod is spatially discretized into 650 elements along the X-axis regarding the trade-off between the accuracy and computation expense. The mesh average element quality, equals to 0.8033, indicates the high discretization accuracy of the model.

The rod modeled by PCS and PLS is divided into three sections, and the length of each section is separately $9 \times 10^{-2}$ m, $7 \times 10^{-2}$ m and $4 \times 10^{-2}$ m from base to tip. For a more specific comparison and evaluation, the material and geometric parameters for the discrete Cosserat models with the same sections in MATLAB are in accordance with those of the FEM. From the simulation results, the equilibrium position of the end-effector of the cantilever rod via FEM under gravity is $\underline{u}_e = [5.8479, 0, -17.8395]^T$, as shown in Fig. 5(a). It took around 14 seconds to complete one simulation because the FEM for large deformations always requires heavy computation. Likewise, we can derive the positions of the end-effector of the discrete Cosserat static models by using the Newton method. From the perspective of computation time of discrete Cosserat static models, we observe that the systems can converge in less than 2 seconds (i.e., 1.3 s for PCS, and 1.2 s for PLS) mainly due to the use of the basic idea of order-model reduction. To intuitively demonstrate the whole shape accuracy of the soft rod modeled by different modeling methods, the spatial configurations of the models under gravity have been depicted in Fig. 5(c) which highlights the high deformation similarity of the PLS and FEM. Comprehensively considering the comparison results among them illustrated in Table II and Fig. 5, we come to a conclusion that, compared to the PCS, the model via the PLS Cosserat modeling approach fits much better the results obtained by the FEM, with the relative error of the end-effector position less than 5%. In other words, the PLS Cosserat static model is essentially comparable to the FEM in terms of accuracy, which can be further verified by the purple PLS configuration tendency plotted in Fig. 5(b)-(c), almost the same as the FEM in Fig. 5(a). Logically, this can be explained by the fact that the piecewise linear

<table><caption>TABLE II
COMPARISON RESULTS OF PCS AND PLS COSSERAT STATIC MODELS W.R.T. FEM UNDER GRAVITY IN TERMS OF END-EFFECTOR POSITION AND ITS RELATIVE ERROR</caption><thead><tr><th rowspan="2">Modeling method</th><th colspan="3">Position of end-effector (Unit:cm)</th><th colspan="3">Relative error of end-effector position w.r.t. FEM (Unit: %)</th></tr><tr><th>u<sub>x</sub></th><th>u<sub>y</sub></th><th>u<sub>z</sub></th><th>e<sub>x</sub></th><th>e<sub>y</sub></th><th>e<sub>z</sub></th></tr></thead><tbody><tr><td>FEM</td><td>5.8479</td><td>0.0000</td><td>-17.8395</td><td>×</td><td>×</td><td>×</td></tr><tr><td>PCS [26]</td><td>5.3450</td><td>0.0000</td><td>-17.1693</td><td>-8.60</td><td>0.00</td><td>-16.88</td></tr><tr><td>PLS</td><td>5.7787</td><td>0.0000</td><td>-17.8394</td><td>-1.18</td><td>0.00</td><td>-0.01</td></tr></tbody></table>

<table><caption>TABLE III
COMPARISON RESULTS OF DIFFERENT BEAM MODELS VIA PLS COSSERAT W.R.T. FEM UNDER GRAVITY IN TERMS OF END-EFFECTOR POSITION AND ITS RELATIVE ERROR</caption><thead><tr><th rowspan="2">Modeling method</th><th colspan="3">Position of end-effector (Unit:cm)</th><th colspan="3">Relative error of end-effector position w.r.t. FEM (Unit: %)</th></tr><tr><th>u<sub>x</sub></th><th>u<sub>y</sub></th><th>u<sub>z</sub></th><th>e<sub>x</sub></th><th>e<sub>y</sub></th><th>e<sub>z</sub></th></tr></thead><tbody><tr><td>FEM</td><td>5.8479</td><td>0.0000</td><td>-17.8395</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Cosserat</td><td>5.7787</td><td>0.0000</td><td>-17.8394</td><td>-1.18</td><td>0.00</td><td>-0.01</td></tr><tr><td>E-B</td><td>5.4940</td><td>0.0000</td><td>-17.7527</td><td>-6.05</td><td>0.00</td><td>-0.49</td></tr><tr><td>E-K</td><td>5.6925</td><td>0.0000</td><td>-17.8523</td><td>-2.66</td><td>0.00</td><td>0.07</td></tr><tr><td>Timoshenko</td><td>5.3596</td><td>0.0000</td><td>-17.9002</td><td>-8.35</td><td>0.00</td><td>0.34</td></tr></tbody></table>

interpolation technique applied to all cross sections of the proposed PLS Cosserat model allows to precisely approximate the deformation behavior of the soft arm in the real scenario.

## C. Accuracy Comparison of PLS Model with Different Modes

Based on the aforementioned parameters setting of the soft arm, different models via the PLS Cosserat involving Euler-Bernoulli (E-B), extensible Kirchhoff (E-K) and Timoshenko beams are established by strain mode selection. These beams are fixed at $X = 0$, and subject to gravity as well as an external imposed concentrated load with $F_{\text{tip}} = [0, 0, 0, F, 0, 0]^T$ (i.e., tension along X-axis) at $X = L_N$ in the inertial frame, we can then convert the concentrated load to the body frame with $\mathcal{F}_e(L_N) = \text{Ad}_{g(L_N)}^\top \mathcal{F}_{\text{tip}}$. The three-dimensional static simulation is implemented by increasing load with increment of 0.05 N at a time. Fig. 6(a)-(c) displays a contrast of the evolution of the equilibrium configurations among these beams for several sets of the tip load, and Table III shows different beams' end-effector positions versus that of FEM. The results indicate that it is feasible to remove negligible modes in some particular cases with low-precision requirement, which contributes to the real-time simulation and control.

# VI. PLS COSSERAT MODEL VALIDATION

## A. Studied Soft Manipulator

A prototype of soft manipulator, similar to that used in the simulation, was designed to carry out the material parameters identification of the PLS Cosserat model by the real input-output relationship obtained from the experimental setup. The exact geometric parameters of the manipulator are illustrated in Fig. 7. The investigated soft arm is controlled by 4 cables mounted through it from the base to the tip, and the cables
<figure><img src="image_7.png" /><figcaption>Fig. 5. Simulation comparison of three different modeling methods for the cantilever rod before and after deformation. (a) FEM. (b) PCS and PLS Cosserat models in MATLAB. (c) Configuration of the soft manipulator under gravity in three-dimensional space.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Deformation comparison among beam models via PLS Cosserat under external loads. (a) Cosserat versus Euler-Bernoulli. (b) Cosserat versus Timoshenko. (c) Cosserat versus extensible Kirchhoff.</figcaption></figure>

are respectively driven by different weights, as shown in Fig. 8. To obtain the position of any waypoint along the soft arm in real time, the magnetic sensors are placed on the corresponding position depicted in blue point in Fig. 7, and for this purpose, a long and conical hole was made along the whole length of the arm. In addition, several 3D-printed rigid rings are mounted along the soft manipulator to minimize friction between the cables and a single conical piece of silicone. The casting material of the manipulator we used is an isotropic silicone rubber with unknown material parameters (i.e., Young's modulus, and density of material) which are needed to be identified by using the proposed parameters identification scheme.

<figure><img src="image_9.png" /><figcaption>Fig. 7. Soft manipulator prototype and position sensors distribution.</figcaption></figure>

## B. Validation and Discussion

The model validation which determines whether the model is proper enough for its intended use is implemented. Six sets of experiments are performed to acquire the position information of end-effector under the effect of different loads exerted by weights equivalent to the cables' tension. Table IV provides specific input and output values of the experiments, and the position as well as orientation of the manipulator in several cases are displayed in Fig. 8. Subsequently, the experimental outputs obtained by the position sensor are utilized to realize

the parameters identification algorithm, and thus the material parameter vector $\theta$ can be calculated.

It should be pointed out that the proposed parameters identification scheme for the PLS Cosserat model has been established to identify the material-related parameters. Thus, if we swap one soft manipulator to another with the same material, or simply change the configuration of the same robot with the same material (e.g. increasing length of the arm), then it is not necessary to repeat the experiments.

The material parameters obtained by solving the optimization problem (26) are as follows: Young's modulus $E$ =
<figure><img src="image_10.png" /><figcaption>Fig. 8. Several snapshots of the different experiments for implementing parameters identification.</figcaption></figure>

<table><caption>TABLE IV
EXPERIMENTAL SAMPLES FOR IDENTIFICATION ALGORITHM</caption><thead><tr><td>Order of experiment</td><td>Cables' tension (Unit: N)</td><td>End-effector position measured by the sensor (Unit: cm)</td></tr></thead><tbody><tr><td>1</td><td>[0, 0, 0, 0]<sup>T</sup></td><td>[15.77, 0.03, -10.10]<sup>T</sup></td></tr><tr><td>2</td><td>[0, 0, 0, 0.98]<sup>T</sup></td><td>[17.73, 0.02, -6.44]<sup>T</sup></td></tr><tr><td>3</td><td>[0.98, 0, 0, 4.90]<sup>T</sup></td><td>[6.85, -4.77, 7.20]<sup>T</sup></td></tr><tr><td>4</td><td>[0, 0, 0, 1.96]<sup>T</sup></td><td>[18.31, 0.01, -1.90]<sup>T</sup></td></tr><tr><td>5</td><td>[1.96, 0, 0, 1.96]<sup>T</sup></td><td>[13.09, -9.24, -2.46]<sup>T</sup></td></tr><tr><td>6</td><td>[0.98, 0, 0, 0]<sup>T</sup></td><td>[14.72, -4.20, -9.78]<sup>T</sup></td></tr></tbody></table>

2.563 × 10⁵ Pa, shear modulus G = 8.543 × 10⁴ Pa, and density of material ρ = 1.41 × 10³ kg/m³.

After that, model validation is performed to verify the accuracy of the PLS Cosserat model with the identified parameters. It should be emphasized that the loading conditions (e.g. order/method of applied loads) have a minor or negligible impact on the results. To avoid these uncertain disturbances in the experiments, we repeated 10 times for each input, and then recorded the average of the three representative waypoints at X = 6.4 cm, X = 12.8 cm and X = 19.2 cm along the soft manipulator, respectively. Five different groups of control inputs and average of outputs presented in Table V are selected to compare the position and orientation of the soft arm between the experiments and simulations. The comparison results demonstrate that the position vectors at different waypoints along the PLS Cosserat model are almost the same as those of the experiments in three cases illustrated in Fig. 9. As for the remaining two sets of experiments, there are larger absolute errors of the selected waypoints between the model and manipulator than those from the other three cases, which may be due to the occurrence of tiny pleats on the silicone surface caused by relatively larger cables' tension.

The absolute errors of the three waypoints of the model with respect to those of the experiments are all within ±5 mm, as shown in Fig. 10, further showing the effectiveness of the parameters identification method and the accuracy of the

<figure><img src="image_11.png" /><figcaption>Fig. 9. Configuration comparison between simulations and experiments under gravity and different cables' tension. (a) $T = [0, 0, 0, 2.94]^T$. (b) $T = [0, 0, 1.96, 1.96]^T$. (c) $T = [0.98, 0, 0, 2.94]^T$.</figcaption></figure>

proposed PLS Cosserat model.

## VII. CONCLUSION AND FUTURE WORK

A piecewise linear strain Cosserat model for the soft slender manipulator has been developed for the first time, which combines the strengths of the PCS and VS Cosserat models. This method depends on a rigorous mathematical framework via the Lie group theory which facilitates a natural coupling of the position and the orientation variables, and exhibits an important advantage. The PLS Cosserat model has been compared with the discrete models reported in the published literature, showing comparable to the FEM and even better results than the PCS Cosserat model in terms of accuracy. As proved in Section V, the PLS Cosserat model shows great potential to be universally applied to the modeling of slender rod-like soft arms in a real scenario.

The parameter identification based on the PLS Cosserat static model can be considered as a NLP problem with several nonlinear equality constraints. A general framework has been proposed for the material parameters estimation, which is applicable to identify the parameters of the soft manipulators with arbitrary cross-sectional shape and actuation manner. To carry out the model validation, we designed a manipulator
<figure><img src="image_12.png" /><figcaption>Fig. 10. The illustration of absolute errors of the different waypoints (i.e., X = 6.4 cm, X = 12.8 cm, X = 19.2 cm,) for the PLS Cosserat model and soft prototype under five different sets of control inputs. (a) X-axis. (b) Y-axis. (c) Z-axis.</figcaption></figure>

<table><caption>TABLE V
DIFFERENT EXPERIMENTS FOR PLS COSSERAT MODEL VALIDATION</caption><thead><tr><th rowspan="2">Order of the<br>control input</th><th rowspan="2">Cables' tension<br>(Unit: N)</th><th colspan="3">Average of the waypoints<br>along the soft arm measured by<br>the position sensors (Unit: cm)</th></tr><tr><th>X = 6.4</th><th>X = 12.8</th><th>X = 19.2</th></tr></thead><tbody><tr><td rowspan="2">1 (10 times)</td><td rowspan="2">[0, 0, 0, 2.94]<sup>T</sup></td><td>[6.15]</td><td>[12.29]</td><td>[16.81]</td></tr><tr><td>[0.10<br>-0.60]</td><td>[0.06<br>-0.66]</td><td>[0.03<br>2.87]</td></tr><tr><td rowspan="2">2 (10 times)</td><td rowspan="2">[0, 0, 1.96, 1.96]<sup>T</sup></td><td>[5.91]</td><td>[10.97]</td><td>[12.96]</td></tr><tr><td>[0.93<br>-1.03]</td><td>[3.96<br>-2.56]</td><td>[9.42<br>-2.55]</td></tr><tr><td rowspan="2">3 (10 times)</td><td rowspan="2">[0.98, 0, 0, 2.94]<sup>T</sup></td><td>[6.24]</td><td>[11.83]</td><td>[14.97]</td></tr><tr><td>-0.57<br>-0.55</td><td>-2.40<br>-0.58</td><td>-6.03<br>2.53</td></tr><tr><td rowspan="2">4 (10 times)</td><td rowspan="2">[0, 0, 2.94, 6.86]<sup>T</sup></td><td>[4.89]</td><td>[3.96]</td><td>[1.50]</td></tr><tr><td>[1.75<br>1.95]</td><td>[5.16<br>4.74]</td><td>[3.69<br>2.37]</td></tr><tr><td rowspan="2">5 (10 times)</td><td rowspan="2">[0, 0, 0, 5.88]<sup>T</sup></td><td>[6.30]</td><td>[8.68]</td><td>[4.28]</td></tr><tr><td>[0.06<br>0.84]</td><td>[0.15<br>5.84]</td><td>[0.12<br>6.89]</td></tr></tbody></table>

prototype made of silicone and established the experimental platform. Both simulation and experiment results indicate that the proposed scheme is able to identify the material parameters of soft arm with high precision, and provides a foundation for developing the model-based controllers. In our future work, the standard Lagrangian formulation in (31) for PLS Cosserat dynamic model actuated by cables will be used to design static and dynamic model-based controllers due to the relative parameterization of the soft arm.

## ACKNOWLEDGMENTS

This work is partially supported by the National Natural Science Foundation of China (Grant No. 62073081), by the Project of Department of Education of Guangdong Province (Grant No. 2019KZDXM037), by the fund of Guangdong-

Hong Kong-Macao Joint Laboratory for Intelligent Micro-Nano Optoelectronic Technology (No. 2020B1212030010), by the project ROBOCOP [ANR-19-CE19], by the project COSSEROOTS [ANR-20-CE33]. The author Haihong Li would like to acknowledge the support from the China Scholarship Council (Grant no. 202008440356).

## APPENDIX

### A. Analytic solution of differential kinematics model

Based on the aforementioned two assumptions, the right-most velocity of any segment *j* along the section *n* at time *t* can be calculated by means of the integral of DE (2).

$$ \eta(L_{n-1} + j\Delta X) = \left[ \prod_{i=0}^{j-1} e^{-\Delta X_{ad} \xi_n(L_{n-1} + i\Delta X)} \right] \eta(L_{n-1}) + \sum_{i=1}^{j} \left\{ \left[ \prod_{\tau=i}^{j-1} e^{-\Delta X_{ad} \xi_n(L_{n-1} + \tau\Delta X)} \right] \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{(s-L_{n-1}-i\Delta X)\text{ad}\xi_n(L_{n-1}+(i-1)\Delta X)} \dot{\xi}_n(L_{n-1} + (i-1)\Delta X) ds \right\} $$

After that, the velocity of each cross section along the section *n* at time *t* can be analytically obtained

$$ \eta(X) = e^{-(X-L_{n-1}-j\Delta X)\text{ad}\xi_n(L_{n-1}+j\Delta X)} \left\{ \left[ \prod_{i=0}^{j-1} e^{-\Delta X_{\text{ad}}\xi_n(L_{n-1}+i\Delta X)} \right] \eta(L_{n-1}) + \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\Delta X_{\text{ad}}\xi_n(L_{n-1}+\tau\Delta X)} \right) \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{(s-L_{n-1}-i\Delta X)\text{ad}\xi_n(L_{n-1}+(i-1)\Delta X)} \dot{\xi}_n(L_{n-1}+(i-1)\Delta X) ds \right] + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\text{ad}\xi_n(L_{n-1}+j\Delta X)} \dot{\xi}_n(L_{n-1}+j\Delta X) ds \right\} $$

Under the PLS assumption, the velocity of any cross section at *X* and along the section *n* at time *t* can be re-formulated
as follows

$$
\begin{equation}
\begin{aligned}
\eta(X) ={}& e^{-(X-L_{n-1}-j\Delta X)\mathrm{ad}\Theta_{nj}^{\vee}} \left\{ \left( \prod_{i=0}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{ni}^{\vee}} \right) \eta(L_{n-1}) \right. \\
& \left. + \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{n\tau}^{\vee}} \right) \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \alpha_{n(i-1)} \mathrm{d}s \right] \dot{\xi}_{n-1} \right. \\
& \left. + \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{n\tau}^{\vee}} \right) \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \beta_{n(i-1)} \mathrm{d}s \right] \dot{\xi}_{n} \right\} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{nj}^{\vee}} \alpha_{nj} \mathrm{d}s \dot{\xi}_{n-1} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{nj}^{\vee}} \beta_{nj} \mathrm{d}s \dot{\xi}_{n}
\end{aligned}
\tag{32}
\end{equation}
$$

with $\Theta_{ni} = \alpha_{ni}\hat{\xi}_{n-1} + \beta_{ni}\hat{\xi}_n$, $\Theta_{nj} = \alpha_{nj}\hat{\xi}_{n-1} + \beta_{nj}\hat{\xi}_n$, where $\alpha_{ni} = 1 - \frac{i\Delta X}{L_n - L_{n-1}}$, $\beta_{ni} = \frac{i\Delta X}{L_n - L_{n-1}}$, $\alpha_{nj} = 1 - \frac{j\Delta X}{L_n - L_{n-1}}$, $\beta_{nj} = \frac{j\Delta X}{L_n - L_{n-1}}$.

By the same reason, the rightmost acceleration of any segment *j* of the section *n* along the soft manipulator at time *t* can be computed by the use of the integral of (3):

as follows

$$
\begin{equation}
\begin{aligned}
\dot{\eta}(X) ={}& e^{-(X-L_{n-1}-j\Delta X)\mathrm{ad}\Theta_{nj}^{\vee}} \left\{ \left( \prod_{i=0}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{ni}^{\vee}} \right) \dot{\eta}(L_{n-1}) \right. \\
& + \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{n\tau}^{\vee}} \right) \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \alpha_{n(i-1)} \mathrm{d}s \right] \dot{\xi}_{n-1} \\
& + \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{n\tau}^{\vee}} \right) \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \beta_{n(i-1)} \mathrm{d}s \right] \dot{\xi}_{n} \\
& + \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{n\tau}^{\vee}} \right) \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \alpha_{n(i-1)} \mathrm{d}s \right] \dot{\xi}_{n-1} \\
& + \sum_{i=1}^{j} \left[ \left( \prod_{\tau=i}^{j-1} e^{-\mathrm{ad}\Delta X\Theta_{n\tau}^{\vee}} \right) \int_{L_{n-1}+(i-1)\Delta X}^{L_{n-1}+i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \beta_{n(i-1)} \mathrm{d}s \right] \dot{\xi}_{n} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{nj}^{\vee}} \alpha_{nj} \mathrm{d}s \dot{\xi}_{n-1} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{nj}^{\vee}} \beta_{nj} \mathrm{d}s \dot{\xi}_{n} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{n\tau}^{\vee}} \alpha_{n(i-1)} \mathrm{d}s \dot{\xi}_{n-1} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \beta_{n(i-1)} \mathrm{d}s \dot{\xi}_{n} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{n\tau}^{\vee}} \alpha_{n(i-1)} \mathrm{d}s \dot{\xi}_{n-1} \\
& + \int_{L_{n-1}+j\Delta X}^{X} e^{-(X-s)\mathrm{ad}\Theta_{n(i-1)}^{\vee}} \beta_{n(i-1)} \mathrm{d}s \dot{\xi}_{n}
\end{aligned}
\tag{33}
\end{equation}
$$

## B. Jacobian Matrix and Its Derivative in Time

From recursive use of (32) and (33) for the whole manipulator, the thorough calculations of the Jacobian matrix *J*(*q*, *X*) and its partial derivative in time *J̇*(*q*, *dot q*, *X*) are respectively given by

$$
\begin{equation}
\begin{aligned}
\dot{\eta}(L_{n-1} + j\Delta X) ={}& \left[ \prod_{i=0}^{j-1} e^{-\Delta X \mathrm{ad} \boldsymbol{\xi}_n (L_{n-1} + i \Delta X)} \right] \dot{\eta}(L_{n-1}) \\
& + \sum_{i=1}^{j} \left\{ \left[ \prod_{\tau=i}^{j-1} e^{-\Delta X \mathrm{ad} \boldsymbol{\xi}_n (L_{n-1} + \tau \Delta X)} \right] \int_{L_{n-1} + (i-1)\Delta X}^{L_{n-1} + i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\boldsymbol{\xi}_n (L_{n-1} + (i-1)\Delta X)} \mathrm{d}s \ddot{\boldsymbol{\xi}}_n (L_{n-1} + (i-1)\Delta X) \right\} \\
& + \sum_{i=1}^{j} \left[ \left[ \prod_{\tau=i}^{j-1} e^{-\Delta X \mathrm{ad} \boldsymbol{\xi}_n (L_{n-1} + \tau \Delta X)} \right] \int_{L_{n-1} + (i-1)\Delta X}^{L_{n-1} + i\Delta X} e^{-(s-L_{n-1}-i\Delta X)\mathrm{ad}\boldsymbol{\xi}_n (L_{n-1} + (i-1)\Delta X)} \mathrm{d}s \ddot{\boldsymbol{\xi}}_n \right\}
\end{aligned}
\tag{34}
\end{equation}
$$

## C. Simplification of Generalized Internal Wrench

Considering the PLS hypothesis, the acceleration of any cross section at *X* and along the section *n* at time *t* can be re-written

Based on linear constitutive relationship chosen for both the elastic and viscous members, the generalized internal wrench $F_i(q, \dot{q})$ which is a nonlinear function of the joint positions
D. Selection schemes of $\mathbf{B}_a$ and $\mathbf{B}_c$

* All strain modes but two curvatures on the Y-axis and Z-axis are neglected for Euler-Bernoulli beam in 3-D space:

$$ \mathbf{B}_a = \begin{bmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 0 \\ 0 & 0 \end{bmatrix}, \quad \mathbf{B}_c = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}. $$

* For extensible Kirchhoff rod, bending, twist, and extension modes are considered:

$$ \mathbf{B}_a = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{B}_c = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{bmatrix}. $$

* As for Timoshenko beam, all modes except twist and extension about X-axis are included:

$$ \mathbf{B}_a = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{B}_c = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 1 \\ 0 & 0 \end{bmatrix}. $$

REFERENCES

[1] D. Rus and M. T. Tolley, "Design, fabrication and control of soft robots," *Nature*, vol. 521, no. 7553, pp. 467-475, 2015.

[2] S. Rozen-Levy, W. Messner, and B. A. Trimmer, "The design and development of branch bot: a branch-crawling, caterpillar-inspired, soft robot," *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 24-36, 2021.

[3] J. Till and D. C. Rucker, "Elastic stability of cosserat rods and parallel continuum robots," *IEEE Trans. Robot.*, vol. 33, no. 3, pp. 718-733, Jun. 2017.

[4] F. Chen and M. Y. Wang, "Design optimization of soft robots: A review of the state of the art," *IEEE Robot. Autom. Mag.*, vol. 27, no. 4, pp. 27-43, 2020.

[5] G. Olson, R. L. Hatton, J. A. Adams, and Y. Mengüç, "An eulernernoulli beam model for soft robot arms bent through self-stress and external loads," *Int. J. Solids Struct.*, vol. 207, pp. 113-131, 2020.

[6] L. Lindenroth, J. Back, A. Schoisengeier, Y. Noh, H. Würdemann, K. Althoefer, and H. Liu, "Stiffness-based modelling of a hydraulically-actuated soft robotics manipulator," in *Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst.*, Oct. 2016, pp. 2458-2463.

[7] H. Godaba, F. Putzu, T. Abrar, J. Konstantinova, and K. Althoefer, "Payload capabilities and operational limits of eversion robots," in *Annu. Conf. Towards Auton. Robotic Syst.*, 2019, pp. 383-394.

[8] F. Boyer, G. De Nayer, A. Leroyer, and M. Visonneau, "Geometrically Exact Kirchhoff Beam Theory: Application to Cable Dynamics," *J. Comput. Nonlinear Dyn.*, vol. 6, no. 4, pp. 041004-041014, 2011.

[9] A. Novelia, *Discrete Elastic Rods for Simulating Soft Robot Limbs*. University of California, Berkeley, 2018.

[10] F. Boyer and F. Renda, "Poincaré's equations for cosserat media: Application to shells," *J. Nonlinear Sci.*, vol. 27, no. 1, pp. 1-44, 2016.

[11] C. B. Black, J. Till, and D. C. Rucker, "Parallel continuum robots: Modeling, analysis, and actuation-based force sensing," *IEEE Trans. Robot.*, vol. 34, no. 1, pp. 29-47, Feb. 2018.

[12] D. C. Rucker and R. J. Webster III, "Statics and dynamics of continuum robots with general tendon routing and external loading," *IEEE Trans. Robot.*, vol. 27, no. 6, pp. 1033-1044, Dec. 2011.

[13] D.-Q. Cao and R. W. Tucker, "Nonlinear dynamics of elastic rods using the cosserat theory: Modelling and simulation," *Int. J. Solids Struct.*, vol. 45, no. 2, pp. 460-477, 2008.

[14] Y. Haibin, K. Cheng, L. Junfeng, and Y. Guilin, "Modeling of grasping force for a soft robotic gripper with variable stiffness," *Mech. Mach. Theory.*, vol. 128, pp. 254-274, 2018.

[15] J. D. Till, *On the Statics, Dynamics, and Stability of Continuum Robots: Model Formulations and Efficient Computational Schemes*. University of Tennessee, Knoxville, 2019.

[16] J. Till, V. Aloi, and C. Rucker, "Real-time dynamics of soft and continuum robots based on cosserat rod models," *Int. J. Robot. Res.*, vol. 38, no. 6, pp. 723-746, Apr. 2019.

[17] F. Renda, M. Giorelli, M. Calisti, M. Cianchetti, and C. Laschi, "Dynamic model of a multibending soft robot arm driven by cables," *IEEE Trans. Robot.*, vol. 30, no. 5, pp. 1109-1122, Oct. 2014.

[18] X. Zhang, F. K. Chan, T. Parthasarathy, and M. Gazzola, "Modeling and simulation of complex dynamic musculoskeletal architectures," *Nature Commun.*, vol. 10, no. 1, pp. 1-12, 2019.

[19] F. Boyer, V. Lebastard, F. Candelier, and F. Renda, "Dynamics of continuum and soft robots: A strain parameterization based approach," *IEEE Trans. Robot.*, vol. 37, no. 3, pp. 847-863, Jan. 2021.

[20] H.-J. Su, "A Pseudorigid-Body 3R Model for Determining Large Deflection of Cantilever Beams Subject to Tip Loads," *J. Mech. Robot.*, vol. 1, no. 2, Jan. 2009, 021008-1-021008-9.

[21] S. Huang, D. Meng, X. Wang, B. Liang, and W. Lu, "A 3d static modeling method and experimental verification of continuum robots based on pseudo-rigid body theory," in *Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst.*, Nov. 2019, pp. 4672-4677.

[22] R. J. Webster III and B. A. Jones, "Design and kinematic modeling of constant curvature continuum robots: A review," *Int. J. Robot. Res.*, vol. 29, no. 13, pp. 1661-1683, 2010.

[23] V. Sonneville, A. Cardona, and O. Brüls, "Geometrically exact beam finite element formulated on the special euclidean group se (3)," *Comput. Methods Appl. Mech. Engrg.*, vol. 268, pp. 451-474, 2014.

[24] F. Largilliere, V. Verona, E. Coevoet, M. Sanz-Lopez, J. Dequidt, and C. Duriez, "Real-time control of soft-robots using asynchronous finite element modeling," in *Proc. IEEE Int. Conf. Robot. Autom.*, May 2015, pp. 2550-2555.

[25] G. Zheng, O. Goury, M. Thieffry, A. Kruszewski, and C. Duriez, "Controllability pre-verification of silicone soft robots based on finite-element method," in *Proc. IEEE Robot. Automat., Int. Conf.*, May. 2019, pp. 7395-7400.

[26] F. Renda, F. Boyer, J. Dias, and L. Seneviratne, "Discrete cosserat approach for multisection soft manipulator dynamics," *IEEE Trans. Robot.*, vol. 34, no. 6, pp. 1518-1533, Dec. 2018.

[27] I. S. Godage, G. A. Medrano-Cerda, D. T. Branson, E. Guglielmino, and D. G. Caldwell, "Dynamics for variable length multisection continuum arms," *Int. J. Robot. Res.*, vol. 35, no. 6, pp. 695-722, 2016.

[28] B. A. Jones and I. D. Walker, "Kinematics for multisection continuum robots," *IEEE Trans. Robot.*, vol. 22, no. 1, pp. 43-55, 2006.

[29] V. Falkenhahn, A. Hildebrandt, R. Neumann, and O. Sawodny, "Dynamic control of the bionic handling assistant," *IEEE/ASME Trans. Mechatronics*, vol. 22, no. 1, pp. 6-17, 2016.

[30] C. Della Santina, A. Bicchi, and D. Rus, "On an improved state parametrization for soft robots with piecewise constant curvature and its use in model based control," *IEEE Robot. Autom. Lett.*, vol. 5, no. 2, pp. 1001-1008, 2020.

[31] K. M. De Payrebrune and O. M. O'Reilly, "On the development of rod-based models for pneumatically actuated soft robot arms: a five-parameter constitutive relation," *Int. J. Solids Struct.*, vol. 120, pp. 226-235, 2017.

[32] S. Grazioso, G. Di Gironimo, and B. Siciliano, "A geometrically exact model for soft continuum robots: The finite element deformation space formulation," *Soft robot*, vol. 6, no. 6, pp. 790-811, 2019.

[33] T. G. Thuruthel, E. Falotico, F. Renda, and C. Laschi, "Model-based reinforcement learning for closed-loop dynamic control of soft robotic manipulators," *IEEE Trans. Robot.*, vol. 35, no. 1, pp. 124-134, Feb. 2019.

[34] I. Hussain, M. Malvezzi, D. Gan, Z. Iqbal, L. Seneviratne, D. Pratichizzo, and F. Renda, "Compliant gripper design, prototyping, and modeling using screw theory formulation," *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 55-71, 2021.

[35] C. Armanini, I. Hussain, M. Z. Iqbal, D. Gan, D. Pratichizzo, and F. Renda, "Discrete cosserat approach for closed-chain soft robots: Application to the fin-ray finger," *IEEE Trans. Robot.*, vol. 37, no. 6, pp. 2083-2098, Apr. 2021.

[36] F. Renda, C. Armanini, V. Lebastard, F. Candelier, and F. Boyer, "A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation," *IEEE Robot. Autom. Lett.*, vol. 5, no. 3, pp. 4006-4013, 2020.

[37] J. Wu, J. Wang, and Z. You, "An overview of dynamic parameter identification of robots," *Robot. Comput-Integr. Manuf.*, vol. 26, no. 5, pp. 414-419, 2010.

[38] J. Jovic, A. Escande, K. Ayusawa, E. Yoshida, A. Kheddar, and G. Venture, "Humanoid and human inertia parameter identification using hierarchical optimization," *IEEE Trans. Robot.*, vol. 32, no. 3, pp. 726-735, Jun. 2016.

[39] Y. Zhang, S. Chen, S. Li, and Z. Zhang, "Adaptive projection neural network for kinematic control of redundant manipulators with unknown
physical parameters,” *IEEE Trans. Ind. Electron.*, vol. 65, no. 6, pp. 4909–4920, 2017.

[40] D. Chen, Y. Zhang, and S. Li, “Tracking control of robot manipulators with unknown models: A jacobian-matrix-adaption method,” *IEEE Trans. Ind. Informat.*, vol. 14, no. 7, pp. 3044–3053, 2017.

[41] N. K. Uppalapati, G. Singh, and G. Krishnan, “Parameter estimation and modeling of a pneumatic continuum manipulator with asymmetric building blocks,” in *Proc. IEEE Int. Conf. Soft Robot.*, 2018, pp. 528–533.

[42] M. Di Lecce, O. Onaizah, P. Lloyd, J. Chandler, and P. Valdastri, “Evolutionary inverse material identification: Bespoke characterization of soft materials using a metaheuristic algorithm.” *Front. in Robot. AI.*, vol. 8, pp. 790571–790571, 2021.

[43] Y. Han, J. Wu, C. Liu, and Z. Xiong, “Static model analysis and identification for serial articulated manipulators,” *Robot. Comput. Integr. Manuf.*, vol. 57, pp. 155–165, 2019.

[44] F. Boyer, M. Porez, and A. Leroyer, “Poincaré–cosserat equations for the lighthill three-dimensional large amplitude elongated body theory: application to robotics,” *J. Nonlinear Sci.*, vol. 20, no. 1, pp. 47–79, 2010.

[45] J. Linn, H. Lang, and A. Tuganov, “Geometrically exact cosserat rods with kelvin–voigt type viscous damping,” *Mech. Sci.*, vol. 4, no. 1, pp. 79–96, 2013.

[46] S. S. Antman, “The theory of rods,” in *Linear theories of elasticity and thermoelasticity*. Springer, 1973, pp. 641–703.

[47] J. Till, V. Aloi, K. E. Riojas, P. L. Anderson, R. J. Webster III, and C. Rucker, “A dynamic model for concentric tube robots,” *IEEE Trans. Robot.*, vol. 36, no. 6, pp. 1704–1718, Dec. 2020.

[48] C. Della Santina and D. Rus, “Control oriented modeling of soft robots: the polynomial curvature case,” *IEEE Robot. Autom. Lett.*, vol. 5, no. 2, pp. 290–298, 2019.

[49] B. C. Hall et al., *Lie groups, Lie algebras, and representations: an elementary introduction*, 2nd ed. New York: Springer, 2003.