Dynamic modeling of soft continuum manipulators using lie group variational integration

Abstract

This paper presents the derivation and experimental validation of algorithms for modeling and estimation of soft continuum manipulators using Lie group variational integration. Existing approaches are generally limited to static and quasi-static analyses, and are not sufficiently validated for dynamic motion. However, in several applications, models need to consider the dynamical behavior of the continuum manipulators. The proposed modeling and estimation formulation is obtained from a discrete variational principle, and therefore grants outstanding conservation properties to the continuum mechanical model. The main contribution of this article is the experimental validation of the dynamic model of soft continuum manipulators, including external torques and forces (e.g., generated by magnetic fields, friction, and the gravity), by carrying out different experiments with metal rods and polymer-based soft rods. To consider dissipative forces in the validation process, distributed estimation filters are proposed. The experimental and numerical tests also illustrate the algorithm’s performance on a magnetically-actuated soft continuum manipulator. The model demonstrates good agreement with dynamic experiments in estimating the tip position of a Polydimethylsiloxane (PDMS) rod. The experimental results show an average absolute error and maximum error in tip position estimation of 0.13 mm and 0.58 mm, respectively, for a manipulator length of 60.55 mm.

About the Authors:

Abbas Tariverdi

Roles Formal analysis, Visualization, Writing – original draft

* E-mail: [email protected]

Affiliation: Department of Physics, University of Oslo, Oslo, Norway

ORCID logo http://orcid.org/0000-0001-5024-9224

Venkatasubramanian Kalpathy Venkiteswaran

Roles Writing – review & editing

Affiliation: Department of Biomechanical Engineering, University of Twente, Enschede, The Netherlands

Ørjan Grøttem Martinsen

Roles Funding acquisition

Affiliations Department of Physics, University of Oslo, Oslo, Norway, Department of Clinical and Biomedical Engineering, Oslo University Hospital, Oslo, Norway

Ole Jacob Elle

Roles Writing – review & editing

Affiliations The Intervention Centre, Oslo University Hospital, Oslo, Norway, Department of Informatics, University of Oslo, Oslo, Norway

Jim Tørresen

Roles Writing – review & editing

Affiliation: Department of Informatics, University of Oslo, Oslo, Norway

Sarthak Misra

Roles Resources

Affiliations Department of Biomechanical Engineering, University of Twente, Enschede, The Netherlands, Department of Biomedical Engineering, University of Groningen and University Medical Centre Groningen, Groningen, The Netherlands

Abstract

This paper presents the derivation and experimental validation of algorithms for modeling and estimation of soft continuum manipulators using Lie group variational integration. Existing approaches are generally limited to static and quasi-static analyses, and are not sufficiently validated for dynamic motion. However, in several applications, models need to consider the dynamical behavior of the continuum manipulators. The proposed modeling and estimation formulation is obtained from a discrete variational principle, and therefore grants outstanding conservation properties to the continuum mechanical model. The main contribution of this article is the experimental validation of the dynamic model of soft continuum manipulators, including external torques and forces (e.g., generated by magnetic fields, friction, and the gravity), by carrying out different experiments with metal rods and polymer-based soft rods. To consider dissipative forces in the validation process, distributed estimation filters are proposed. The experimental and numerical tests also illustrate the algorithm’s performance on a magnetically-actuated soft continuum manipulator. The model demonstrates good agreement with dynamic experiments in estimating the tip position of a Polydimethylsiloxane (PDMS) rod. The experimental results show an average absolute error and maximum error in tip position estimation of 0.13 mm and 0.58 mm, respectively, for a manipulator length of 60.55 mm.

Figures

Table 7

Table 8

Fig 16

Table 1

Fig 1

Table 2

Fig 2

Fig 3

Fig 4

Fig 5

Fig 6

Table 3

Fig 7

Table 4

Fig 8

Fig 9

Table 5

Fig 10

Fig 11

Table 6

Fig 12

Fig 13

Fig 14

Fig 15

Table 7

Table 8

Fig 16

Table 1

Fig 1

Table 2

Citation: Tariverdi A, Venkiteswaran VK, Martinsen ØG, Elle OJ, Tørresen J, Misra S (2020) Dynamic modeling of soft continuum manipulators using lie group variational integration. PLoS ONE 15(7): e0236121. https://doi.org/10.1371/journal.pone.0236121

Editor: Ning Cai, Beijing University of Posts and Telecommunications, CHINA

Received: March 13, 2020; Accepted: June 29, 2020; Published: July 22, 2020

Copyright: © 2020 Tariverdi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the manuscript.

Funding: The author(s) received no specific funding for this work.

Competing interests: The authors have declared that no competing interests exist.

1 Introduction

Reachability, high level of dexterity, and large elastic deformability are the primary driving factors behind the growth of research in the design, modeling, and control of continuum manipulators. Flexible continuum manipulators have recently generated interest in several fields [1–3], especially in minimally invasive surgical robotics and interventional medicine, such as catheter-based endovascular intervention [4, 5] and cardiac surgeries [6, 7]. In contrast to conventional rigid link manipulators, soft manipulators are able to reshape their configurations to allow for redundancies in path planning, and are capable of precise and delicate manipulation of objects in complex and varying environments.

There are numerous candidate actuation mechanisms for continuum manipulators such as tendon-drives and concentric tubes [8–11]. Compared to other actuation mechanisms, magnetic actuation benefits from high dexterity, versatility, and wireless actuation [12–15]. By applying remote magnetic torques on permanent magnets or coils which are embedded inside the body of a manipulator and/or at its tip, one can control the motion and configuration of the manipulator.

This paper aims to develop a computational model for analyzing the dynamics of soft continuum manipulators, which is one of the key challenges in soft robotics. In many tasks, dynamic models of manipulators are essential for control, trajectory planning, and optimal design purposes, especially in Minimally Invasive Surgeries (MIS) for operation in unknown and unstructured environments such as inside the human body. Due to elastic characteristics and geometric nonlinearities (i.e., bending, torsion, shear, elongation, including large deformation) of continuum manipulators, their dynamics have highly nonlinear behavior and are expressed as partial differential equations. Some recent modeling approaches of soft continuum manipulators/robots, which have been employed in the surgical robotics field, are summarized in Table 1.

[Figure omitted. See PDF.]

Table 1. References on dynamics/static analysis of soft continuum manipulators in surgical robotics field.

https://doi.org/10.1371/journal.pone.0236121.t001

Soft continuum manipulators are analogous to specific Cosserat continuums. Therefore, Lie group synchronous variational integrators [35, 36], a novel time and space integration scheme, is employed in this paper to model geometrically exact beams based on the Simo beam model [37] and Hamiltonian formulation. The core idea of this algorithm is to obtain the dynamic behavior of the system while conserving the invariants (energy, momentum maps) of the system, especially for long-time simulations. The distinguishing characteristic of variational integrators is that they define the equations of motion based on the discretized variational principle of the system. Combining the integrators with Lie-group/algebraic techniques enables the algorithm to respect not only the structure of the dynamics and its properties but also preserve the structure of the configuration space. The advantages of employing the Lie group variational integration method compared to other modeling strategies is that the proposed solver is applicable to exact nonlinear dynamic models of continuum manipulators subject to large deformations. The algorithm preserves the symplectic structure, i.e., the invariants of mechanical systems. Also, it allows the usage of different time steps at different points in a given finite element for the geometry of soft manipulators. These properties are investigated in previous work (e.g., [35, 38, 39]), while the main focus of this paper is the experimental validation of the method on magnetically-actuated soft continuum manipulators.

Investigation of previous work in modeling of the continuum manipulators suggests that existing literature focuses primarily on static or quasi-static approaches, or does not provide sufficient experimental validation in realistic application scenarios. By contrast, the main contribution of this article compared with the existing work in literature is the validation of an accurate dynamic model of a soft continuum manipulator, considering spatial motion. Also, it should be noted that the model accounts for the geometric nonlinearities (e.g., large deformation) and respects conservation of dynamical properties of the system (e.g., energy and momentum maps conservation), and structures of configuration space simultaneously. Besides, it should be pointed out that three-dimensional internal and external dissipation forces act on the continuum manipulator and hence affect the dynamics. Therefore, it is necessary to consider these friction/ dissipation forces in the validation process. To this end, distributed prediction filters have been proposed.

In summary, this article’s contributions can be stated as follows.

* Existing studies on the modeling of continuum manipulators primarily consider static or quasi-static approaches. However, in numerous applications, the fully spatial dynamics of manipulators need to be considered for accurate control and design purposes. The primary contribution of this article is the derivation and experimental validation of a dynamic model for forced continuum manipulators with soft materials undergoing spatial deformation. The model accounts for the nonlinearities of the continuum manipulator; loading resulted from magnetic fields, the gravity, and internal and external dissipation forces generated by friction.

* Due to the difficulty in obtaining knowledge about the internal and external dissipation forces, distributed estimation filters have been designed to take these forces into account and capture their behavior.

The rest of the paper is organized as follows. In Section 2, mathematical preliminaries, including the system description and notation, are discussed. Next, Section 3 addresses the algorithm and numerical results. The experimental framework and implementation results are described in Section 4, demonstrating the effectiveness of the theoretical formulation. In addition, Section 5 provides a discussion on the implementation of the modeling algorithm. Finally, Section 6 summarizes the results of this work and draws conclusions and posits directions for future work.

2 Continuum manipulator dynamics

This section is devoted to describing kinematics and full three-dimensional dynamics for continuum manipulators undergoing large deflections (for detailed explanations, refer to the reference [35]). We review the static description of a continuum manipulator in three-dimensional space toward deriving the dynamic equations of motion of geometrically exact continuum manipulator by applying Hamilton’s principle to the Lagrangian of the system.

2.1 Kinematics

The manifold of configuration space of a continuum manipulator considering Boundary Conditions (BCs) is defined asin which L is the length of the undeformed continuum manipulator, maps the line of continuum manipulator’s centroids (i.e. center of mass) to Euclidean space and the orthogonal transformation determines the orientation of moving cross-sections at points in the terms of a fixed basis . Therefore, the orientation of each cross-section which is denoted by directors or moving basis can be written as Fig 1 shows initial and a time-evolved configuration of the continuum manipulator with the free right tip and clamped left end, i.e, BCs: , . The BCs imply that the clamped cross section is orthogonal to the plane defined by and . In addition, a curve characterizes a time-evolved configuration space of the continuum manipulator. The family of tangent vectors to the curve q(t) is defined aswhich characterize tangent bundle to at the manifold q(s, t).

[Figure omitted. See PDF.]

Fig 1. Initial and time-evolved configurations of the continuum manipulator.

The highlighted frames depict cross-sections at discretization points. Fixed bases or material frame are also shown at the fixed end of the manipulators. Also, moving bases are attached to the cross section at the centroid s and the tip of manipulators.

https://doi.org/10.1371/journal.pone.0236121.g001

2.2 Lagrangian and equation of motion

To derive the equations of motion, we first need to introduce the Lagrangian of the system which can be written as(1)where the matrices C1 and C2 are defined as C1 ≔ diag(GA GA EA) and C2 ≔ diag(EI1 EI2 GJp). For brevity, other parameters are defined in Table 2.

[Figure omitted. See PDF.]

Table 2. Definition of parameters in Lagrangian (1) as described in [35].

https://doi.org/10.1371/journal.pone.0236121.t002

where each cross section is given by a compact set , Lie algebra is associated with the Lie group SO(3), and Hat map/ operator which is a one-to-one invertible map, i.e., an isomorphism, is defined as(2)

The Euler-Lagrange equations are obtained by applying by the Lagrange-d’Alembert principle to the action functional associated to , namelyBy employing the Lagrange-d’Alembert principle, one computesThe terms δω, δΩ, and δΓ are defined ([38]) as follows:(3)where .

Taking into account the expressions for δω, δΩ, and δΓ in Eq (3) and using integration by parts in space and time, we obtain Euler-Lagrange equations with non-conservative force as(4)in which we could consider 6 × 1 representations of general non-conservative force vector where and are force and moment vectors in , respectively. Also, denotes the cotangent bundle of . For simplicity, one may think of the cotangent boundle as the space of positions and momenta. For the exact definitions, refer to [40] or [41].

3 Lie group variational integrators for the forced continuum manipulator

In this section, the focus is on analyzing a Lie group variational integrator for continuum manipulators with conservative (e.g., the gravity) and non-conservative forces (e.g., friction and loads inserted by actuators). In the following subsection, the discretized version of the forced Euler-Lagrange Eq (4) is given (for further details and stability analysis, see [35]), and afterward, the estimation process is discussed.

3.1 Modeling

This section is devoted to introducing a Lie group variational integration scheme for continuum manipulators with external loading. First, one needs to consider the spatial discretization of Lagrangian introduced in the previous section. Afterward, discrete Lagrange-d’Alembert equations need to be expressed on Lie group SE(3). These equations are employed to propose a model-based distributed estimation scheme. Fig 2 depicts the modeling procedure in this section.

[Figure omitted. See PDF.]

Fig 2. Steps 1 through 5 toward deriving continuum manipulator discrete dynamics.

https://doi.org/10.1371/journal.pone.0236121.g002

Here notations of the paper are provided. Additionally, concepts and definitions on Lie groups and Lie algebra are presented in Appendix 6.

Notations: The undeformed continuum manipulator’s length [0, L] is spatially discretized into N subsets of length . For an element , ai and ai+1 denote its left and right nodes. The configuration of the continuum manipulator at the node ai is given by and . Also, denotes the angular velocity of a node ai. Given a node ai, the discrete time evolution of this node is given by the discrete curve and is based on the discrete Euler-Lagrange equations on Lie group SE(3), The discrete variables associated with the node ai are defined as . We denote the fixed time step by Δt = tj − tj+1, j = 0, ⋯, V. In time discretization of the continuum manipulator, we have .

By identifying the configuration space of the continuum manipulator with the infinite dimensional Lie group , we consider the trivialized Lagrangian , where is a Lie algebra associated with the Lie group G. A spatial discretization of the trivialized Lagrangian for an element and the total system are computed as follows, respectfully. It should be noted that the evaluation of Lagrangian at midpoints of nodes is employed. Other evaluations of the Lagrangian depending on a different number or combinations of nodes are possible (see [38]).

For an element :(5)where is conservative potential energy of an element due to the gravity and elasticity and given by(6)For the whole continuum manipulator:

The temporal discretized Lagrangian approximates the Lagrangian in Eq (5) during the time step Δt is therefore(7)where and .

The discrete action sum over the discretized time interval [0, T] = {t0, ⋯, tj|tj = tj−1 + Δt, t0 = 0, tV = T}, is computed as follows.The discrete Lagrange–d’Alembert principle is(8)

By applying the discrete Lagrange–d’Alembert principle (8), we get the discrete Euler–Lagrange equations for a node ai in a compact form as(9)

Finally, using the definitions of adjoint and coadjoint actions, and cotangent lift of left translation which are presented in Appendix 6, Eqs (7) and (9) yields Eqs (10)–(12) and (14)–(16) to update rotations and positions of each node.

3.1.1 Discrete Euler-Lagrange equations for rotations.

* For the left node of the continuum manipulator (ai=0)(10)

* For the interior nodes of the continuum manipulator ∀ai, i ∈ {1, ⋯, N − 1}(11)

* For the right node of the continuum manipulator (ai=N)(12)

where the variable is defined as which is approximated by the Cayley transformation as , where The Cayley transformation and its inverse are defined in the following form for convenience with inverse (see, [35, 42]). In addition, .

For discrete Euler-Lagrange equations for rotations, Eqs (10)–(12), one has to solve an implicit expression of the form(13)

In order to solve Eq (13) for F ∈ SO(3), (the vector or the right hand sides of Eqs (10)–(12) and the symmetric non-standard inertia matrix Jd are known), a Newton iteration method based on the Cayley transformation is applied (as described in [39], Section 3:3:8]).

3.1.2 Discrete Euler-Lagrange equations for translations.

* For the left node of the continuum manipulator (ai=0)(14)

* For the interior nodes of the continuum manipulator ∀ai, i ∈ {1, ⋯, N − 1}(15)

* For the right node of the continuum manipulator (ai=N)(16)

Remark 1 For magnetic actuation, we fabricate manipulators with embedded permanent magnets. Consider Magnet i with weight mi in an interval in which ai and ai+1, i ∈ {0, ⋯, N − 1} are considered as left and right nodes of the interval. Therefore, the distributed load per unit length for Nodes ai and ai+1 are approximately considered as . In addition, if Magnet i is embedded at the tip, replaces the distributed load per unit length of Node aN+1 while the distributed load per unit length of Node aN is unchanged.

3.2 Estimation

In this section, online distributed estimation algorithms are developed to predict the model dissipation error. The structure of the estimator mimics the model’s structure, as explained in Section 3.1. To design the estimation protocol, we follow the same line of ideas as in [42] but in distributed multi-systems configuration. We consider each node as an individual system coupled with the other adjacent nodes, i.e., neighbors, in succession. In other words, each node exchanges its local pose (position and orientation) with its neighbors. It should be noted that the estimation filters are designed and implemented for each node. Fig 3 shows the configuration of the distributed filters and nodes.

[Figure omitted. See PDF.]

Fig 3. Configuration of the nodes of the model and the corresponding distributed filters.

Filter ai and Node ai are coupled with the adjacent nodes in succession.

https://doi.org/10.1371/journal.pone.0236121.g003

For simplicity, we assume that each node’s position is included in its state vector. Therefore, given node ai, i = {0, ⋯, N}, the time-varying dynamic equations based on Eqs (14)–(16) can be written as(17)where j = {1, 2, ⋯}, i = {0, 2, ⋯, N}, is the true state vector, is a known input vector, is sufficiently differentiable, is the model dissipation error matrix and in the considered systems is time-varying, is a modified viscous model dissipation force or hysteretic damping force in the form of , where ∘ denotes Hadamard product and is damping capacity that is independent of frequency of motion and needs to be estimated, is the output matrix, and is the output vector.

By substituting into Eq (17), one has(18)Using commutative property of Hadamard product, we can write . Then, Hadamard product can be converted to matrix multiplication by the corresponding diagonal matrix of the vector which is denoted by and . Therefore, Eq (18) may be written asEstimation of the state and output vector is given bywhere is the estimation of the state vector, is the model dissipation error estimates, is the output vector estimates. Finally, denotes the measurement. The block diagram of the filter integrated with the model is shown in Fig 4. To find for the node ai at the time j, we consider a pointwise cost function that penalizes and minimizes the estimation error vector (the error between measurement and the output estimation) at the next sampling time j + 1, and estimated damping capacity . The cost function for each node ai is given as(19)where , , and are positive semi-definite and positive definite matrices, respectively.

[Figure omitted. See PDF.]

Fig 4. Block diagram of the proposed prediction filter ai coupled with the node ai’s model.

The filter employs the model’s output to perform a pointwise optimization problem to predict the damping capacity .

https://doi.org/10.1371/journal.pone.0236121.g004

In order to derive an optimal estimation law, we need to approximate the output estimation vector at the next sampling time j + 2, which is given by its Taylor series expansion as follows(20)where

Similarly, we may expand the ith component of in an first-order Taylor series so thatwhere the hth component of isSolving Eq (19) for by considering Eq (20) yields(21)

Stability and convergence analysis of the filters can be found in [42]. Here we skip the analysis for brevity.

4 Simulation and experimental results

In this section, we investigate and analyze the solver’s performance with different continuum manipulators through experiments. The experiments here are expected to provide validation of the theoretical formulation for a variety of scenarios. As discussed earlier, it is worth remembering that the dynamic equations for translation and rotation are decoupled. Eqs (14)–(16) can be solved explicitly to update nodes translation, while an iterative method—as it is discussed in Section 3.1.1—is necessary to solve Eqs (10)–(12) for updating the rotations. It should be pointed out that the estimation law (21) is implemented for every node to estimate conservative forces. The required parameters for the simulation will be discussed for each experiment.

4.1 Flexible metal rods

As a first case, we consider a cylindrical rod made of aluminum (Al4043/ AlSi5) with diameter of 2 mm, length 200 mm, mass density 2690 , Young’s modulus 75 GPa and Poisson’s ratio 0.33. As a first example, we suppose a planar motion of the rod in the -plane with the initial deflection . denotes the rotation of the tip around in the -plane. It is worth pointing out that the nondissipative force is the gravity in the -axis direction. In addition, we run a simulation with the given specifications with N = 15 discretization nodes. These points are depicted in Fig 5 together with some time-evolved configurations of the rod. For simplicity, only tip positions are used for the comparison with the simulation results. The maximum and mean absolute error are 0.15 mm (i.e., 2.5% of displacement), and 0.05 mm, respectively. The error, simulation and experiment results are shown in Fig 6 and the simulation parameters are summarized in Table 3.

[Figure omitted. See PDF.]

Fig 5. Sample of grabbed images of flexible rod (AlSi05) configurations, in-plane experiment (-plane).

15 discretization nodes, depicted in blue squares, are superimposed on the flexible rod.

https://doi.org/10.1371/journal.pone.0236121.g005

[Figure omitted. See PDF.]

Fig 6.

Simulation and experiment results for flexible rod (AlSi05), in-plane experiment: (a) Tip position in -axis direction. Inset highlights the results in a small time range. (b) Error in -axis direction.

https://doi.org/10.1371/journal.pone.0236121.g006

[Figure omitted. See PDF.]

Table 3. Simulation parameters in Eqs (10)–(16) and (21) for in-plane experiment of flexible cylindrical rod (AlSi05).

https://doi.org/10.1371/journal.pone.0236121.t003

Next, we consider a three-dimensional motion for a rod with the same material as the first case but with diameter d = 1 mm with initial deflections (i.e., the tip distance is 8 mm from -axis) and (i.e., the tip distance is 10 mm from -axis) in the and planes, respectively. The results of the experiment, simulation, and error are depicted in Fig 7. Maximum and mean absolute error in both and axes are 0.15 mm (i.e., 2.12%), and 0.05 mm, respectively. The simulation parameters are summarized in Table 4.

[Figure omitted. See PDF.]

Fig 7.

Simulation and experiment results for flexible rod (AlSi05), out-of-plane experiment: (a) Tip position in -axis direction. (b) Tip position in -axis direction. (c) Error in -axis direction. (d) Error in -axis direction. (e) Tip 3D position: non-planar experiment. (f) Tip 3D position: simulation.

https://doi.org/10.1371/journal.pone.0236121.g007

[Figure omitted. See PDF.]

Table 4. Simulation parameters in Eqs (10)–(16) and (21) for out-of-plane experiment of flexible cylindrical rod (AlSi05).

https://doi.org/10.1371/journal.pone.0236121.t004

4.2 Polymer-based rods

In the second experiment, a cylindrical Polydimethylsiloxane (PDMS) rod is considered. Fig 8 depicts the rod, which has diameter D = 5 mm, length L = 60.5 mm. In addition, for the rod, mass density , Young’s modulus E = 365.12 MPa, and Poisson’s ratio ν = 0.5. In this experiment, the rod is kept straight initially, with the gravity acting along , Fig 8(a). Experiment and simulation results for the tip position and the error are depicted in Fig 9. Also, maximum and mean absolute error in -axis are 0.56 mm (i.e., 4.87%), and 0.05 mm, respectively. In -axis, maximum and mean absolute error are 0.28 mm (i.e., 4.89%), and 0.05 mm, respectively. The simulation parameters are summarized in Table 5.

[Figure omitted. See PDF.]

Fig 8. Sample of grabbed images for Polydimethylsiloxane (PDMS) rod without any embedded magnet in 2D experiment: In-plane motion.

Also, 10 discretization points are superimposed on the soft rod.

https://doi.org/10.1371/journal.pone.0236121.g008

[Figure omitted. See PDF.]

Fig 9.

Simulation and experiment results of Polydimethylsiloxane (PDMS) rod, planar motion: (a) Tip Position in -axis direction. Inset magnifies the results in a small time range. (b) Tip Position in -axis direction. (c) Error in -axis direction. (d) Error in -axis direction.

https://doi.org/10.1371/journal.pone.0236121.g009

[Figure omitted. See PDF.]

Table 5. Simulation parameters in Eqs (10)–(16) and (21) for in-plane experiment of PDMS rod (without magnet).

https://doi.org/10.1371/journal.pone.0236121.t005

For the next experiment, we fabricated a cylindrical PDMS manipulator with a permanent magnet at the tip. The initial and some time-evolved configurations of the rod are depicted in Fig 10. The specifications of the rod are as follows: diameter D = 4 mm, length L = 60.55 mm. In addition, the embedded neodymium magnet is a cylindrical magnet with diameter Dm = 2 mm, height Lm = 4 mm, mass Mm = 9.6 × 10−5 kg. The rod moves around with the initial deflection in the -plane. Also, maximum and mean absolute error in -axis are 0.55 mm (i.e., 1.16%), and 0.13 mm, respectively. In -axis, maximum and mean absolute error are 0.61 mm (i.e., 3.35%), and 0.14 mm, respectively. Tip positions in the experiment, simulation and the error are shown in Fig 11. The simulation parameters are summarized in Table 6.

[Figure omitted. See PDF.]

Fig 10. Sample of grabbed images: Polydimethylsiloxane (PDMS) rod with an embedded magnet at the tip in 2D experiment: Motion in a plane.

Also, 10 discretization points are shown on the soft rod.

https://doi.org/10.1371/journal.pone.0236121.g010

[Figure omitted. See PDF.]

Fig 11.

Simulation and experiment results of Polydimethylsiloxane (PDMS) rod with an embedded magnet at the tip, planar experiment: (a) Tip Position in -axis direction. Inset highlights the results in a small time range. (b) Tip Position in -axis direction. (c) Error in -axis direction. (d) Error in -axis direction.

https://doi.org/10.1371/journal.pone.0236121.g011

[Figure omitted. See PDF.]

Table 6. Simulation parameters in Eqs (10)–(16) and (21) for in-plane and circular-motion experiments of PDMS rod (with tip magnet).

https://doi.org/10.1371/journal.pone.0236121.t006

Hereafter, a magnetic field generation setup is employed to actuate the manipulators. The following section introduces magnetic field generation setup and the related background.

4.3 Magnetic field generation

The setup used here consists of two pairs of Helmholtz coils to generate magnetic fields. Each pair consists of two identical electromagnetic coils, as shown in Fig 12. The first pair of coils generates a uniform magnetic field along the -axis. The second pair of smaller coils are placed inside the first pair to produce a field along the -axis. Two cameras are placed next to the setup to monitor the side view of the workspace. For image acquisition, we use both cameras in a stereo vision setup to reconstruct 3D views of the manipulator’s motion. The setup produces a maximum magnetic field Bu = 45 mT.

[Figure omitted. See PDF.]

Fig 12. Magnetic field generation setup with the stereo vision cameras.

Two nested pairs of Helmholtz coils generate uniform magnetic fields in and -axes direction.

https://doi.org/10.1371/journal.pone.0236121.g012

For the first experiment using magnetic actuation, we use the rod with a neodymium magnet with diameter 2 mm, height 4 mm, and magnetisation N45.

First, a magnetic field Bg = 7.75 mT is applied to compensate for the gravity. Then, the tip of the manipulator is induced to rotate in a circle in the -plane using a rotating magnetic field of magnitude Bu = 14.5 mT.

The magnetic field produces force and torque Frot and τrot, respectively, given bywhere m is the dipole moment of the tip’s magnet. The dipole moment can be computed as in which residual magnetism Br ∈ [1.32, 1.37] mT, μ0 is the permeability of vacuum, and the volume of the magnet, V = 4π mm3. Experiment, simulation results and the error are shown in Fig 13. It should be noted that the error plot shows Euclidean norm of the tip position in the experiment and simulation. Also, the maximum and mean absolute errors are 1.20 mm (i.e., 1.43%) and 0.59 mm, respectively. Since we use the same manipulator as the previous experiment, simulation parameters can be found in Table 6.

[Figure omitted. See PDF.]

Fig 13.

Reconstruction of the scene for circular motion of Polydimethylsiloxane (PDMS) rod with an embedded magnet at the tip: (a,b) 3D experiment and simulation results. (c) Euclidean distance/error of experiment and simulation results in tip position.

https://doi.org/10.1371/journal.pone.0236121.g013

As a last experiment, we fabricated a PDMS continuum manipulator with a square cross-section and two embedded permanent magnets, one at the tip and another in the middle—36.1 mm from the tip—of the manipulator. The embedded neodymium magnets are identical cylindrical magnets with different dipole moment’s directions and diameter Dm = 2 mm, height Lm = 3 mm, weight Mm = 7.2 × 10−5 kg. The embedded magnets are induced to pursue a prescribed motion in the -plane using a varying magnetic field of initial and final magnitude Bu = 20 mT and 19.85 mT. The initial and some time-evolved configurations of the rod are depicted in Fig 14. It should be noted that analysis of magnetic force and torque follows the same procedure as described above. The specifications of the manipulator are as follows: edge length a = 2 mm, length L = 85.5 mm. The maximum and mean absolute error for the tip magnet are 1.00 mm (i.e., 2.24%) and 0.15 mm, in -axis direction, respectively. In -axis direction, the maximum and mean absolute error for the tip magnet are 1.40 mm (i.e., 5.13%) and 0.33 mm, respectively. The maximum and mean absolute error for the middle magnet are 0.47 mm (i.e., 2.05%) and 0.08 mm, in -axis direction, respectively. In -axis direction, the maximum and mean absolute error for the middle magnet are 0.40 mm (i.e., 3.68%) and 0.10 mm, respectively.

[Figure omitted. See PDF.]

Fig 14. Sample of grabbed images: Polydimethylsiloxane (PDMS) rod with two embedded magnets, 2D experiment: Motion in a plane.

Also, 10 discretization points, (blue squares,) are shown on the rod.

https://doi.org/10.1371/journal.pone.0236121.g014

The position of the tip and middle magnets in the experiment and simulation and also the error is shown in Fig 15. For this experiment, the simulation parameters are summarized in Table 7.

[Figure omitted. See PDF.]

Fig 15.

Simulation and experiment results of Polydimethylsiloxane (PDMS) rod with two embedded magnets, 2D experiment: (a) Tip magnet’s position in -axis direction. (b) Error of tip magnet’s position in -axis direction. (c) Tip magnet’s position in -axis direction. (d) Error of tip magnet’s position in -axis direction. (e) Middle magnet’s position in -axis direction. (f) Error of middle magnet’s position in -axis direction. (g) Middle magnet’s position in -axis direction. (h) Error of middle magnet’s position in -axis direction.

https://doi.org/10.1371/journal.pone.0236121.g015

[Figure omitted. See PDF.]

Table 7. Simulation parameters in Eqs (10)–(16) and (21) for in-plane experiment of square cross-section PDMS rod (with 2 magnets).

https://doi.org/10.1371/journal.pone.0236121.t007

5 Discussion

We validate our approach by designing and carrying out different experiments with flexible metal rods and polymer-based soft rods. The results are summarized in Table 8.

[Figure omitted. See PDF.]

Table 8. Maximum and mean absolute error in the experiments using flexible metal rods (AlSi05) and Polydimethylsiloxane (PDMS) rods.

https://doi.org/10.1371/journal.pone.0236121.t008

Table 8 demonstrates the maximum and the mean absolute values of the errors. As we observe from this table, the simulation results closely match the experimental responses, i.e., for Experiments 1 and 2 in which the flexible metal rods (AlSi05) are employed, the worst-case errors are < 0.01% of the manipulators’ length. For dynamic Experiments 3 and 4 in which the PDMS rods are used, maximum of errors respectively are 0.95% and 1% of the manipulator’s length. For the polymer rods, higher errors are due to the uncertainties in fabrication and nonlinear elastic properties. In quasi-static Experiments 5 and 6, the manipulators experience large deformations and external loads; the worst-case errors are less than 2% and 1% of the manipulators’ length. It should be pointed out that compared to the manipulators’ length, the mean absolute deviations are small, which reflect the model’s performance.

During the implementation of the modeling approach, it was observed that the number of nodes affects the frequency of motion. Increasing the number of nodes provides a more accurate solution for the frequency of the system. However, the computation time increases significantly with the number of nodes. Therefore, to be able to run the simulations in a reasonable amount of time and with a small number of nodes, frequency shaping was necessary to be able to match the results.

The following example shows the motivation behind the frequency shaping of the motion. Consider a manipulator with the following specifications: Length L = 0.5 m, mass density , square cross-section with edge length a = 5 cm, Poisson’s ratio ν = 0.35, and Young’s modulus Ehf = 5 × 104 KPa in the high frequency case and Elf = 500 KPa. Fig 16 compares the position between high and low-frequency cases. The base of the manipulator is fixed at the origin. From Fig 16, it is observed that by changing the Young Modulus from Ehf to Elf, tip motion is preserved but in a scaled frequency.

[Figure omitted. See PDF.]

Fig 16.

Simulation results for high and low frequency comparison: (a) High frequency. (b) Low frequency. (c) Error of in high and low frequency cases.

https://doi.org/10.1371/journal.pone.0236121.g016

The frequency of the continuum motion is only dependant on parameters such as length, the moment of inertia of the cross-section, Young modulus, and material density. Then, the natural frequency of the continuum manipulator with a fixed end and free tip can be written as

Consider a manipulator with an equivalent spatial discretization of its central line by elements. denotes the natural frequency of the manipulator in the simulation with elements, and one has

The correction factor which needs to be multiplied by Young Modulus in the simulation is obtained as

By employing this correction factor in the simulation, the effect of the number of discretization elements on the frequency of motion can be eliminated.

High fidelity models are helpful for explaining and predicting the behavior of a system with complex dynamics. However, due to computational constraints, these models may not be employed for closed-loop control purposes in a real-time implementation of robotic applications. Additionally, recent developments in computer simulations demand superior, robust, and efficient numerical frameworks compared to traditional approaches. Discrete geometric mechanics, which are employed in this paper, provides a systematic method to cope with the complexity of continuum manipulators’ dynamics. The necessity of guaranteeing robots’ performance in sensitive applications such as minimally invasive surgeries requires the use of pre-existing knowledge or a model in control architecture to obtain guaranteed and reliable behavior in the presence of disturbances and uncertainties. Although model-free control approaches are easy to implement, they do not provide and ensure any performance level and high control-loop bandwidths.

6 Conclusions and future work

This article studies the estimation and model validation problem of continuum manipulators’ dynamics using Lie group variational integrators. Using magnetic actuation, dynamic and static experiments were conducted on manipulators with rigid and soft materials (e.g., Aluminum and PDMS) to illustrate the validity of the presented algorithm for a wide range of experiments.

Due to the lack of knowledge about friction/damping, distributed predictive filters were designed to provide information about the unknown signals. Therefore, the dynamical model equipped with the estimation algorithm is a self-contained generic model for continuum manipulator integration, which provides us with a systematic approach to employ optimal control theory for realistic trajectory planning in the presence of user/environment-specified constraints. The designing of a controller and the parallel variational integration algorithm are to be investigated as future work.

Supporting information

[Figure omitted. See PDF.]

S1 Appendix. Some concepts and definitions on Lie groups and Lie algebra are presented (References [35] and [43]).

https://doi.org/10.1371/journal.pone.0236121.s001

(PDF)

Citation: Tariverdi A, Venkiteswaran VK, Martinsen ØG, Elle OJ, Tørresen J, Misra S (2020) Dynamic modeling of soft continuum manipulators using lie group variational integration. PLoS ONE 15(7): e0236121. https://doi.org/10.1371/journal.pone.0236121

References

1. Polygerinos P, Correll N, Morin SA, Mosadegh B, Onal CD, Petersen K, et al. Soft robotics: review of fluid-driven intrinsically soft devices; manufacturing, sensing, control, and applications in human-robot interaction: review of fluid-driven intrinsically soft robots. Advanced Engineering Materials. 2017;19(12):1700016.

2. Rus D, Tolley MT. Design, fabrication and control of soft robots. Nature. 2015;521(7553):467–475. pmid:26017446

3. Laschi C, Mazzolai B, Cianchetti M. Soft robotics: Technologies and systems pushing the boundaries of robot abilities. Science Robotics. 2016;1(1):eaah3690.

4. Grady MS, Howard MA, Dacey RG, Blume W, Lawson M, Werp P, et al. Experimental study of the magnetic stereotaxis system for catheter manipulation within the brain. Journal of Neurosurgery. 2000;93(2):282–288. pmid:10930015

5. Burgner J, Swaney PJ, Lathrop RA, Weaver KD, Webster RJ. Debulking from within: a robotic steerable cannula for intracerebral hemorrhage evacuation. IEEE Transactions on Biomedical Engineering. 2013;60(9):2567–2575. pmid:23649131

6. Carpi F, Pappone C. Stereotaxis Niobe® magnetic navigation system for endocardial catheter ablation and gastrointestinal capsule endoscopy. Expert Review of Medical Devices. 2009;6(5):487–498. pmid:19751121

7. Kesner SB, Howe RD. Position control of motion compensation cardiac catheters. IEEE Transactions on Robotics. 2011;27(6):1045–1055.

8. Guo S, Fukuda T, Nakamura T, Arai F, Oguro K, Negoro M. Micro active guide wire catheter system-Characteristic evaluation, electrical model and operability evaluation of micro active catheter. In: Proceedings of IEEE International Conference on Robotics and Automation. vol. 3. Minneapolis, MN, USA: IEEE; 1996. p. 2226–2231.

9. Webster RJ, Jones BA. Design and kinematic modeling of constant curvature continuum robots: a review. The International Journal of Robotics Research. 2010;29(13):1661–1683.

10. Camarillo DB, Milne CF, Carlson CR, Zinn MR, Salisbury JK. Mechanics modeling of tendon-driven continuum manipulators. IEEE Transactions on Robotics. 2008;24(6):1262–1273.

11. Cianchetti M, Ranzani T, Gerboni G, De Falco I, Laschi C, Menciassi A. STIFF-FLOP surgical manipulator: Mechanical design and experimental characterization of the single module. In: 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems; 2013. p. 3576–3581.

12. Heunis C, Sikorski J, Misra S. Flexible instruments for endovascular interventions: improved magnetic steering, actuation, and image-guided surgical instruments. IEEE Robotics & Automation Magazine. 2018;25(3):71–82.

13. Sikorski J, Heunis CM, Franco F, Misra S. The armm system: an optimized mobile electromagnetic coil for non-linear actuation of flexible surgical instruments. IEEE Transactions on Magnetics. 2019;55(9):1–9.

14. Sikorski J, Dawson I, Denasi A, Hekman EEG, Misra S. Introducing BigMag—A novel system for 3D magnetic actuation of flexible surgical manipulators. In: 2017 IEEE International Conference on Robotics and Automation (ICRA). Singapore, Singapore: IEEE; 2017. p. 3594–3599.

15. Thomas T, Kalpathy Venkiteswaran V, Ananthasuresh GK, Misra S. A Monolithic Compliant Continuum Manipulator: a Proof-of-Concept Study. Journal of Mechanisms and Robotics. 2020; p. 1–11.

16. Pai DK. Strands: interactive simulation of thin solids using cosserat models. Computer Graphics Forum. 2002;21(3):347–352.

17. Webster RJ, Romano JM, Cowan NJ. Mechanics of precurved-tube continuum robots. IEEE Transactions on Robotics. 2009;25(1):67–78.

18. Xu K, Simaan N. An investigation of the intrinsic force sensing capabilities of continuum robots. IEEE Transactions on Robotics. 2008;24(3):576–587.

19. Camarillo DB, Carlson CR, Salisbury JK. Configuration tracking for continuum manipulators with coupled tendon drive. IEEE Transactions on Robotics. 2009;25(4):798–808.

20. Jones BA, Gray RL, Turlapati K. Three dimensional statics for continuum robotics. In: 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems; 2009. p. 2659–2664.

21. Rucker DC, Webster RJ, Chirikjian GS, Cowan NJ. Equilibrium conformations of concentric-tube continuum robots. The International Journal of Robotics Research. 2010;29(10):1263–1280. pmid:25125773

22. Xu K, Simaan N. Analytic formulation for kinematics, statics, and shape restoration of multibackbone continuum robots via elliptic integrals. Journal of Mechanisms and Robotics. 2010;2(1):011006.

23. Kai Xu, Simaan N. Intrinsic wrench estimation and its performance index for multisegment continuum robots. IEEE Transactions on Robotics. 2010;26(3):555–561.

24. Rucker DC, Jones BA, Webster RJ III. A geometrically exact model for externally loaded concentric-tube continuum robots. IEEE Transactions on Robotics. 2010;26(5):769–780. pmid:21566688

25. Lock J, Laing G, Mahvash M, Dupont PE. Quasistatic modeling of concentric tube robots with external loads. Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems IEEE/RSJ International Conference on Intelligent Robots and Systems. 2010;2010:2325–2332.

26. Tunay I. Distributed parameter statics of magnetic catheters. Conference proceedings: Annual International Conference of the IEEE Engineering in Medicine and Biology Society IEEE Engineering in Medicine and Biology Society Annual Conference. 2011;2011:8344–8347.

27. Mahvash M, Dupont PE. Stiffness control of surgical continuum manipulators. IEEE Transactions on Robotics. 2011;27(2):334–345.

28. Wen Tang, Tao Ruan Wan, Gould DA, Thien How, John NW. A stable and real-time nonlinear elastic approach to simulating guidewire and catheter insertions based on cosserat rod. IEEE Transactions on Biomedical Engineering. 2012;59(8):2211–2218.

29. Tunay I. Spatial continuum models of rods undergoing large deformation and inflation. IEEE Transactions on Robotics. 2013;29(2):297–307.

30. Jung J, Penning RS, Zinn MR. A modeling approach for robotic catheters: effects of nonlinear internal device friction. Advanced Robotics. 2014;28(8):557–572.

31. Greigarn T, Cavusoglu MC. Pseudo-rigid-body model and kinematic analysis of MRI-actuated catheters. In: 2015 IEEE International Conference on Robotics and Automation (ICRA). Seattle, WA, USA: IEEE; 2015. p. 2236–2243.

32. Roesthuis RJ, Misra S. Steering of multisegment continuum manipulators using rigid-link modeling and fbg-based shape sensing. IEEE Transactions on Robotics. 2016;32(2):372–382.

33. Grazioso S, Di Gironimo G, Siciliano B. A geometrically exact model for soft continuum robots: the finite element deformation space formulation. Soft Robotics. 2019;6(6):790–811. pmid:30481112

34. Venkiteswaran VK, Sikorski J, Misra S. Shape and contact force estimation of continuum manipulators using pseudo rigid body models. Mechanism and Machine Theory. 2019;139:34–45.

35. Demoures F, Gay-Balmaz F, Leyendecker S, Ober-Blöbaum S, Ratiu TS, Weinand Y. Discrete variational Lie group formulation of geometrically exact beam dynamics. Numerische Mathematik. 2015;130(1):73–123.

36. Bobenko AI, Suris YB. Discrete time lagrangian mechanics on lie groups, with an application to the lagrange top. Communications in Mathematical Physics. 1999;204(1):147–188.

37. Simo JC. A finite strain beam formulation. The three-dimensional dynamic problem. Part I. Computer Methods in Applied Mechanics and Engineering. 1985;49(1):55–70.

38. Leitz T, Ober-Blöbaum S, Leyendecker S. Variational Lie Group Formulation of Geometrically Exact Beam Dynamics: Synchronous and Asynchronous Integration. In: Terze Z, editor. Multibody Dynamics: Computational Methods and Applications. Computational Methods in Applied Sciences. Cham: Springer International Publishing; 2014. p. 175–203.

39. Lee T. Computational geometric mechanics and control of rigid bodies.; 2008. Available from: https://deepblue.lib.umich.edu/handle/2027.42/60804.

40. Marsden JE, Ratiu TS. Introduction to mechanics and symmetry. vol. 17 of Texts in Applied Mathematics. New York, NY: Springer New York; 1994.

41. Abraham R, Marsden JE. Foundations of mechanics. 2nd ed. Reading, Mass: Benjamin/Cummings Pub. Co; 1978.

42. Lu P. Optimal predictive control of continuous nonlinear systems. International Journal of Control. 1995;62(3):633–649.

43. Holm DD. Geometric Mechanics: Part II: Rotating, Translating and Rolling. 2nd ed. IMPERIAL COLLEGE PRESS; 2011.