Dynamic modeling of soft continuum manipulators using lie group variational integration

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.


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][2][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 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 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][9][10][11]. Compared to other actuation mechanisms, magnetic actuation benefits from high dexterity, versatility, and wireless actuation [12][13][14][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.

References Modeling Approach and its Properties
Robot type/ Application [16] Static analysis: Cosserat rod model. 3D elesticity Surgical suture/ strands [17] Beam mechanics based on elastic energy Concentric tubes/ General MIS [18] Static analysis based on screw theory and a virtual-work model Multiple parallel backbones/ General MIS [10,19] Linear elasticity theory Single/Redundant tendons [20] 3D Static analysis with loads: Cosserat rod model General purpose CRs [21] Beam mechanics based on elastic energy (includes both bending and torsion) Concentric tubes/ General MIS [22] Bernoulli-Euler elastica theory: statics, 2D Multibackbone [23] Static analysis based on a virtual-work model Serial Segments/ General surgical end-effectors [24,25] Static analysis: Cosserat rod theory Concentric tubes with and without external loads [26] Static analysis: Cosserat rod theory Magnetic Catheter/ General purposes [27] Loaded static analysis: Cosserat rod theory General MIS [28] Dynamic analysis: Cosserat rod model. 3D elesticity Guidewire/ Interventional Radiology procedures [29] FEM: large deformation and inflation Simulations on general medical robots [30] Lumped-parameter model Multiple parallel shafts/ general Magnetic resonance imaging (MRI)-compatible medical manipulators [31] Pseudo-rigid-body model Multiple parallel shafts/ cardiac robotic catheter [32] 2D static analysis: rigid-link modeling Planar tendon-driven continuum manipulator/ general medical robots 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.

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 R 3 toward deriving the dynamic equations of motion of geometrically exact continuum manipulator by applying Hamilton's principle to the Lagrangian of the system.

Kinematics
The manifold of configuration space of a continuum manipulator considering Boundary Conditions (BCs) is defined as Q ¼ fðO; PÞ 2 C 1 ð�Þ : ½0; L� ! SOð3Þ � R 3 jBCs are satisfiedg in which L is the length of the undeformed continuum manipulator, P : ½0; L� ! R 3 maps the line of continuum manipulator's centroids (i.e. center of mass) to Euclidean space R 3 and the orthogonal transformation O : ½0; L� ! SOð3Þ determines the orientation of moving crosssections at points PðsÞ in the terms of a fixed basis fE 1 ðsÞ; E 2 ðsÞ; E 3 ðsÞg. Therefore, the orientation of each cross-section which is denoted by directors or moving basis fD 1 ðsÞ; D 2 ðsÞ; D 3 ðsÞg can be written as The BCs imply that the clamped cross section is orthogonal to the plane defined by E 1 and E 2 . In addition, a curve qðs; tÞ ¼ ðOðs; tÞ; Pðs; tÞÞ 2 Q characterizes a time-evolved configuration space of the continuum manipulator. The family of tangent vectors to the curve q(t) is defined as _ qðs; tÞ ¼ dqðs; tÞ dt ¼ _ Oðs; tÞ; _ Pðs; tÞ which characterize tangent bundle T q Q to Q at the manifold q(s, t).

Lagrangian and equation of motion
To derive the equations of motion, we first need to introduce the Lagrangian L : T q Q of the system which can be written as � ds |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } Kinetic energy � ds |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } where the matrices C 1 and C 2 are defined as C 1 ≔ diag(GA GA EA) and C 2 ≔ diag(EI 1 EI 2 GJ p ). For brevity, other parameters are defined in Table 2. where each cross section is given by a compact set A ¼ fðx; yÞjx; y 2 Rg, Lie algebra soð3Þ is associated with the Lie group SO (3), and Hat map/ operator^: R 3 ! soð3Þ which is a one-to-one invertible map, i.e., an isomorphism, is defined as The Euler-Lagrange equations are obtained by applying by the Lagrange-d'Alembert principle to the action functional H associated to L, namely By employing the Lagrange-d'Alembert principle, one computes https://doi.org/10.1371/journal.pone.0236121.g001 Table 2. Definition of parameters in Lagrangian (1) as described in [35].
M = ρ 0 × A ρ 0 and A are the body constant mass density and cross section's area.
Taking into account the expressions for δω, δO, and δΓ in Eq (3) and using integration by parts in space and time, we obtain Euler-Lagrange equations with non-conservative force in which we could consider 6 × 1 representations of general non-conservative force vector where F and N are force and moment vectors in R 3 , respectively. Also, T � q Q denotes the cotangent bundle of Q. 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].

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.

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.
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 For an element I i , a i and a i+1 denote its left and right nodes. The configuration of the continuum manipulator at the node a i is given by O a i :¼ Oðs a i Þ and p a i :¼ Pðs a i Þ. Also, o a i denotes the angular velocity of a node a i . Given a node a i , the discrete time evolution of this node is given by the discrete curve ðO j a i ; p j a i Þ 2 SEð3Þ ¼ SOð3Þ � R 3 ; j ¼ 0; � � � ; V and is based on the discrete Euler-Lagrange equations on Lie group SE(3), The discrete variables F j a i associated with the node a i are In time discretization of the continuum manipulator, we have Dp j a i :¼ p jþ1 a i À p j a i . By identifying the configuration space Q of the continuum manipulator with the infinite dimensional Lie group G ¼ C 1 ð½0; L�; SOð3Þ � R 3 , we consider the trivialized Lagrangian L : G � g ! R, where g is a Lie algebra associated with the Lie group G. A spatial discretization of the trivialized Lagrangian for an element I i 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 I i : where V I i is conservative potential energy of an element I i due to the gravity and elasticity and given by " For the whole continuum manipulator: The temporal discretized Lagrangian L j I i approximates the Lagrangian L I i in Eq (5) during the time step Δt is therefore The discrete action sum over the discretized time interval [0, T] = {t 0 , � � �, t j |t j = t j−1 + Δt, t 0 = 0, t V = T}, is computed as follows.
By applying the discrete Lagrange-d'Alembert principle (8), we get the discrete Euler-Lagrange equations for a node a i in a compact form as 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.

Discrete Euler-Lagrange equations for rotations.
• For the left node of the continuum manipulator (a i=0 ) • For the interior nodes of the continuum manipulator • For the right node of the continuum manipulator (a i=N ) where the variable c a i is defined asĉ The Cayley transformation and its inverse are defined in the following form for [35,42]). In addition, Dp a i j t¼t j ¼ p jþ1 a i À p j a i . For discrete Euler-Lagrange equations for rotations, Eqs (10)- (12), one has to solve an implicit expression of the form In order to solve Eq (13) for F 2 SO (3), (the vector U or the right hand sides of Eqs (10)-(12) and the symmetric non-standard inertia matrix J d are known), a Newton iteration method based on the Cayley transformation is applied (as described in [39], Section 3:3:8]).

Discrete Euler-Lagrange equations for translations.
• For the left node of the continuum manipulator (a i=0 ) Dp a 0 • For the interior nodes of the continuum manipulator Dp a iÀ 1

PLOS ONE
Dynamic modeling of soft continuum manipulators • For the right node of the continuum manipulator (a i=N )

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. For simplicity, we assume that each node's position is included in its state vector. Therefore, given node a i , i = {0, � � �, N}, the time-varying dynamic equations based on Eqs (14)-(16) can be written as S jþ1 where

PLOS ONE
Dynamic modeling of soft continuum manipulators is a modified viscous model dissipation force or hysteretic damping force in the form of K j a i � Dt À 1 ðp j a i À p jÀ 1 a i Þ, where � denotes Hadamard product and K j a i 2 R 3 is damping capacity that is independent of frequency of motion and needs to be estimated, H jþ1 a i 2 R 3�6 is the output matrix, and Y jþ1 a i 2 R 3�1 is the output vector. By substituting F j Using commutative property of Hadamard product, we can write G j Then, Hadamard product can be converted to matrix multiplication by the corresponding diagonal matrix of the vector G j a i ðS j a i ÞDt À 1 ðp j a i À p jÀ 1 a i Þ which is denoted by G j a i ¼ G j a i ðS j a i ÞDt À 1 diagðp j a i À p jÀ 1 a i Þ and G j a i 2 R 6�3 . Therefore, Eq (18) may be written as S jþ1 Estimation of the state and output vector is given bŷ whereŜ j a i ¼ ½p jÀ 1 a ip j a i � T 2 R 6 is the estimation of the state vector,K j a i 2 R 3 is the model dissipation error estimates,Ŷ jþ1 a i 2 R 3�1 is the output vector estimates. Finally,Ỹ jþ1 a i denotes the measurement. The block diagram of the filter integrated with the model is shown in Fig 4. To findK j a i for the node a i 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 K jþ1 a i . The cost function for each node a i is given as where e jþ2 a i ¼Ŷ jþ2 a i ÀỸ jþ2 a i , W 2 R 3�3 , and R 2 R 3�3 are positive semi-definite and positive definite matrices, respectively.
In order to derive an optimal estimation law, we need to approximate the output estimation vectorŶ jþ2 a i at the next sampling time j + 2, which is given by its Taylor   Solving Eq (19) forK jþ1 a i by considering Eq (20) yieldŝ Stability and convergence analysis of the filters can be found in [42]. Here we skip the analysis for brevity.

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.

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 kg m 3 , 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 E 1 E 3 -plane with the initial deflection y E 1 E 3 ¼ 3:69. y E 1 E 3 denotes the rotation of the tip around E 1 in the E 1 E 3 -plane. It is worth pointing out that the nondissipative force is the gravity in the E 3 -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.
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 y E 1 E 2 ¼ À 5:53 (i.e., the tip distance is 8 mm from E 2 -axis) and y E 1 E 3 ¼ 6:52 (i.e., the tip distance is 10 mm from E 3 -axis) in the E 3 E 2 and E 1 E 2 planes, respectively. The results of the experiment, simulation, and error are depicted in Fig 7. Maximum and mean absolute error in both E 3 and E 1 axes are 0.15 mm (i.e., 2.12%), and 0.05 mm, respectively. The simulation parameters are summarized in Table 4.

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 r ¼ 1101 kg m 3 , 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 E 2 , Fig 8(a). Also, maximum and mean absolute error in E 2 -axis are 0.56 mm (i.e., 4.87%), and 0.05 mm, respectively. In E 3 -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.
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 D m = 2 mm, height L m = 4 mm, mass M m = 9.6 × 10 −5 kg. The rod moves around E 1 with the initial deflection y E 2 E 3 ¼ À 41:74 in the E 2 E 3 -plane. Also, maximum and mean absolute error in E 3 -axis are 0.55 mm (i.e., 1.16%), and 0.13 mm, respectively. In E 2 -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.
Hereafter, a magnetic field generation setup is employed to actuate the manipulators. The following section introduces magnetic field generation setup and the related background.

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 E 1 -axis. The second pair of smaller coils are placed inside the first pair to produce a field along the E 2 -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 B u = 45 mT. 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 B g = 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 E 2 E 3 -plane using a rotating magnetic field of magnitude B u = 14.5 mT.
The magnetic field produces force and torque F rot and τ rot , respectively, given by where m is the dipole moment of the tip's magnet. The dipole moment can be computed as m ¼ 1 m 0 B r V in which residual magnetism B r 2 [1.32, 1.37] mT, μ 0 is the permeability of vacuum, and the volume of the magnet, V = 4π mm 3 . 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.
As a last experiment, we fabricated a PDMS continuum manipulator with a square crosssection 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 D m = 2 mm, height L m = 3 mm, weight M m = 7.2 × 10 −5 kg. The embedded magnets are induced to pursue a prescribed motion in the E 2 E 3 -plane using a varying magnetic field of initial and final magnitude B u = 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 E 2 -axis direction, respectively. In E 3 -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 E 2 -axis direction, respectively. In E 3 -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. 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.

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.  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   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 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 preexisting 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.

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

PLOS ONE
Dynamic modeling of soft continuum manipulators 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.