Modeling and Trajectory Tracking Control for a Multi-Section Continuum Manipulator

This paper presents a dynamic modeling and trajectory tracking control design for a hyper-redundant continuum manipulator (HRM), which is able to continuously bend along its length. The HRM has the advantages of high flexibility, large workspace, and low inertia, but also poses challenges of complex kinematics, and dynamic coupling. To address these challenges, the kinematics model was developed based on the piecewise constant curvature (PCC) assumption that captures the relationship between the cable lengths, joint angles, and end-effector pose. Inverse kinematics is presented geometrically by solving trigonometric equations. The dynamic model for the continuum robot (CR) is derived using the Euler–Lagrange representation, which incorporates the effects of gravity and elasticity. Accordingly, three control strategies are developed and applied to a two-section continuum robot, which are the inverse dynamic proportional integral derivative controller (PID), fuzzy logic controller (FLC), and sliding mode controller (SMC). We validate our proposed methods through simulations on various 3D trajectories, utilized MATLAB symbolic math toolbox in conjunction with Simulink simulation; demonstrate and evaluate different close-loop dynamic control responses for different scenarios. This simulation used in the creation of a unique animated simulation for a whole 3D continuum robot while tracking its desired trajectories. Finally, the simulation results demonstrate that our proposed methods can effectively control the HRM with high tracking performance, exhibiting enhanced responses in terms of settling time, rising time, and steady-state error. While the developed SMC outperforms Poth's PID and FLC in terms of robustness and settling time (51% and 10.3%), respectively.


Introduction
Continuum robots are a type of under actuated hyper-redundant manipulators (HRMs), which are typically Inspired by natural biological features, they emulate nature in their movement, such as snakes, elephant trunks, and squid tentacles [1], hyper-redundant is a term which denote robots that have large or infinite degrees of freedom [2].
Continuum Robots (CRs) are fundamentally built based on a flexible backbone on which is mounted a series of rigid discs; their structures continuously bend through elastic deformation along their length, in contrast to conventional robots, which typically consist of a number of rigid links and rotating joints. They can be built with multiple sections; thus, they are theoretically capable of having an endless degree of freedom. Where it shows a significant compliance and great operational capabilities for manipulating objects and interacting with the environment, there for it is a most suitable for minimal invasive surgery [3], medical application [4], and operating in a complex, unstructured environment.
Various techniques can be used to actuate CR, including motor-driven tendons [5], pneumatically operated CRs [6], and an intelligent material, shape memory alloy (SMA) [7], which can alter its shape by passing current throw, where the first two actuation methods use significantly large driving hardware. At the same time, the SMA is best-suitable for miniature CRs.
Continuum robots are complex systems that are challenging to model due to their nonlinearity and infinite degrees of freedom. However, On the basis of various simplifications and assumptions, several research investigations have been attempted to model the kinematics of soft continuum robots.
To the best of our knowledge, the HRM was introduced in 1995 by Chirikjian, G.S. [8]. from then on various methodologies and theories, such as Serret-Frenet frames, Euler-Bernoulli beam equation, and Denavit-Hartenberg parameters (DH), have been utilized to develop a kinematic model [9][10][11]. while the commonly used methods for the dynamic model: finite element method [12] where it used to solve the boundary value problem of a partial differential equation by discretizing the domain of interest into a finite number of subdomains called finite elements and applying a variational principle to minimize the error function and obtain a stable solution. This method is applicable to a variety of complex structures, however obtaining numerical solutions when internal or external motion constraints are present poses a challenge due to the high computational complexity. To obtain a more extensive description of soft robotic arms, the Cosserat theory was proposed The Cosserat technique has been used to control and manipulate flexible robotic arm [13]. Despite the high accuracy of the model created by this technique, the partial differential equations in this model cannot be used in effective control design.
Kazemipour et al. [14] presents a continuum robot dynamic model based on the Euler-Lagrange approach which does not include oversimplified assumptions. In addition they introduces an adaptive sliding mode controller that can adjust its parameters online based on the tracking error. The results shows that the controller can achieve accurate task-space trajectory tracking under different payloads and external forces. Thuruthel et al. [15] proposes a method for modeling and controlling soft robots using first-order dynamical systems based on Cosserat-rod mechanics. in addition they presents two closed-loop control strategies based on the reduced-order state feedback.
The results indicate that controllers developed on this assumption can compensate the errors in modeling accuracy with the increased control frequency. Several studies have explored different techniques for controlling and stabilizing robotic systems using fuzzy logic, neural networks, passivity and flatness [16][17][18][19][20].
We adopted Piece wise constant curvature (PCC) approximation as first step in modeling the kinematic of HRMs, where the bending section of the CR is defined as an arc of a circle and the bending plane is capable of rotating along a fixed axis [21]. While [22,23] covers the inverse kinematics of the PCC model for CRs.
By solving the forward and inverse system kinematic, we developed a system dynamics model using the Euler-Lagrange representation that incorporates the effects of gravity, and elasticity. Where symbolic Math Toolbox used to derive and apply inverse kinematics to robot arms by using symbolic variables for solving system partial derivatives.  Eventually, we developed and applied three various control strategies to control the position and trajectory tracking of the robot. (a) The inverse dynamic PID controller, which utilized the system inverse dynamics to cope with system nonlinearity. (b) The FLC combined with the system's inverse dynamics, which addressed the nonlinearity and coupling effect of the two-link robot arm. (c) The nonlinear SMC with a designed boundary layer δ based on the system stability analysis, which eliminate the chattering phenomena.
By the end, MATLAB Simulink used to simulate and assess different dynamic control loop responses for a given configuration space values, and desired trajectories. This simulation used in the creation of a unique animated simulation for whole 3D continuum robot while tracking its desired trajectories (Online Resource 1). The remainder of this paper is organized as follows. In Section 2, the robot mechanical structure and modelling assumptions are introduced. Section 3 presents the kinematic model for two section CR based on constant curvature (CC) assumption, followed by the dynamic model developed using classic Lagrangian representation, which are shown in Section 4. Three control strategies are developed and investigated in Section 5. Dynamic model validation presented in Section 6. Section 7 presents the simulation study for step and trajectory responses, conducted using MATLAB Simulink. Finally, a conclusion and further research plans are stated in Section 8. Figure 1 shows the continuum robot structure. Each section is comprised of three main parts: a series of disks, a flexible backbone (primary backbone) (PB), and three flexible wires (secondary backbones) (SBs). Each section is made up of a 50 cm long flexible backbone formed Fig. 3 Single section diagram, a The arc is rotated away from the x-z plane by the angle φ, b When φ equals zero, the arc is located on the x-z plane Fig. 4 Two section continuum robot workspaces of composite materials and five disks that are uniformly distributed and rigidly mounted. Each disk features three holes for three driving wires, distributed in a circular pattern 120° apart from one another. In order to control the behavior of the CR, the wires for each section are utilized to apply a moment to the flexible backbone tip. It is crucial to emphasize that the flexible backbone deflection can control the robot's spatial mobility when adequate tension forces are applied to one or two cables at once, which achieves a 2-DOF bending motion of the continuum robot. By using a series of 2-DOF links, a multi-DOF continuum robot arm can be achieved.

The Two-section Continuum Robot Kinematics and Modeling Assumption
Constant curvature has frequently been viewed as a desirable modeling characteristic in hyper-redundant continuum robots due to its simplifications in kinematic modeling.
According to the PCC assumption on modelling CRs presented in [21]. The bending of each section of the robot is considered an arc. The arc parameters are as follows: -, the angle corresponding to the bending of each section. -, the angle of the arc plane. Where i = 1, 2 . in the first section when i = 1 or in the second section when i = 2 . Where, as illustrated in Fig. 1 each subsection has a corresponding bending angle s i .
In Which, ∈ R 4 can be used to express the robot configuration space. Where = 1 1 2 2 T . As shown in Fig. 2, there are three frames with three axes each x i , y i , and z i axis for a two sections CR, with all parameters defined.

Forward Position Kinematics
A continuum robot's geometry provides a method for figuring out how the points along it are posed [21,24,25]. As shown in Also, note that this motion includes a rotation R y ( ) ∈ R 3×3 about the positive y-axis. And rotating the entire arc by about the positive z-axis by moves the robot out of the (x, z) plane, producing a transformation from arc base to tip. (2) and (3): The following matrix, express the end position of each subsection with respect to [Frame i − 1 ], where s i substitute i in the position vector of (4):  By substituting s 2 = l 2 in (6), the end position of the second section with respect to [Frame 0] 0 0 2 can be obtained.

Inverse Position Kinematics
The configuration angles 1 1 1 2 2 2 for two section CR, can geometrically calculated by knowing the end position of each section with respect to [frame 0] shown in Fig. 2. as follows: < 0 Figure 4 shows the CR workspace for each link based on the introduced kinematics model.

Forward Velocity Kinematics
the velocity of each section is obtained by providing the angular velocities of the configuration space ̇ = ̇1̇1̇1̇2 ].
Let the task space of the robot O ∈ R 6 to be: Therefore, by differentiating (10) by the time, the velocity of end point for second section can be obtained: where, the Jacobian matrix is = 0 q ∈ R 6×4 .

Inverse Velocity Kinematics
The configuration space velocities may be acquired in the manner shown below by supplying the task space velocities: Since the system is underactuated, the inverse Jacobian matrix, J −1 ∈ R 4×6 , is constructed using the pseudo-inverse approach.

Forward Acceleration Kinematics
By differentiating (11), the task space acceleration ( ̈ ) can be obtained, as the configuration space angular acceleration is provided: where, ̈ = ̈1̈1̈2̈2 T , and ̇ =̇o q is the differentiation of the Jacobian matrix.

Dynamics of The Two-Section Continuum Robot
In this section, the Eular-Lagrange equation is used to represent whole-body dynamics [10,26,27]. whereby kinetic and potential energy needs to be examined. As depicted in Fig. 1, the robot is divided into sections, each of which is made up of a primary backbone (PB), three secondary backbones (SBs) or actuators, and several discs. Where the MATLAB symbolic toolbox [28,29] is utilized to derive and solve dynamic system

Kinetic Energy Analysis
The three basic components of the kinetic energy analysis are the kinetic energies of the primary backbone (PB), the secondary backbones (SBs), and the kinetic energies of each disc. The kinetic energy for a differential section of the backbone, with mass of dm and inertia moment of dI , The kinetic energy of the PBs of the two sections expressed by KE Pb can be obtained as follows: where, P is the density, A P is the cross-section area, and I P is the second moment of cross-sectional area of the PB.
Calculating the derivative of (5), the velocity of point P can be expressed as follows: The lengths of the SBs are obtained as follows: where, is the distance between the PB and each SB. and each section has its own independent actuators from all other sections. i = 1, 2, l 1,1 , l 1,2 , and l 1,3 are lengths of the first section's actuators (link 1 secondary backbone). And l 2,1 , l 2,2 , and l 2,3 are lengths of the second section's actuators (link 2 secondary backbone). The SBs kinetic energy is divided into two components. The first part KE Sb 1 can be calculated similarly to (14) by substituting each 1 2 with 3 2 (because of the 3 SBs), and each subscript "P" with "S" in which, I S , A S , and S are the second moment of cross-sectional area, the cross-section area, and the density of each SB, respectively.
By the differentiation of the actuators' lengths the second part can be constructed as follows: where, m S is the mass of each SB.
As mentioned in [5], the kinetic energy of the discs located on the first and second sections can be calculated by developing the velocity of each disk, in which all discs are connected to the PB,where the velocity of each disk represents the positions of each disc in the first section with respect to frame 0 when i = 1, and i = 2 for each disc position in the second section with respect to frame 0, where h is the distance between each disk andk i = 1, 2, 3, … , n i . The number of discs in the first and second sections of the CR is indicated by the letters n 1 and n 2 , respectively. n 1 = n 2 = 5 In this study. Where, where, KE D denotes the kinetic energy of the disks placed on the first and the second section. The mass and the mass moment of inertia of each disk are denoted by m D and I D , respectively. The total kinetic energy KE is:

Potential Energy Analysis
The gravitational potential energy and the elastic potential energy form the CR's total potential energy.
The two components of the gravitational energy analysis are the PB and the disc gravitational energies.
In order to compute the gravitational potential energy of the PB GPE Pb , the gravitational acceleration (g) is assumed to be directed in the direction of the positive x 0 -axis. Where, The first and second sections discs' combined gravitational potential energy, denoted by the symbol GPE D , are as follows: The elastic potential energy is given by: where, E P and E S are the modules of elasticity of the PB and the SBs, respectively. (25) represent the total potential energy PE , where it is the summation of (22) and (23), and (24). Where,

Control System Design
Using the dynamic model given in (27) we observe that the continuum manipulator is a complex nonlinear system. In this section, we attempt to apply three nonlinear control algorithms to the manipulator system, including the inverse dynamic control PID, fuzzy logic controller and sliding mode controller.
is the desired configuration space angles that the robot try to reach.
The system coordinates q(t) should follow a reference time profile described by the vector q d (t) , with the error e(t) = q d (t) − q(t) decaying to zero. The desired error equation should be of second order in this case. It can be chosen in the following linear form: From (29), the second order time derivative of the coordinate vector q can be computed as: Substituting expression (30) into the system dynamic model (27) produces the developed control law: where, k p is the proportional gain, k d is the derivative gain, and K i is the integral gain, (t) defined in (28) and ̇ (t) is the differentiation of (28).while k p , k d and K i values are chosen to minimize (t). Figure 5 shows the inverse dynamic PID controller block diagram, where the system dynamic model is solved by MAT-LAB symbolic math toolbox, and the values of vector q are used as runtime feedback to the controller.

Fuzzy Logic Controller
Fuzzy logic control (FLC) is a well-known control approach based on artificial intelligence. It exploits the functionary's previous knowledge of the system to be controlled. The primary Fig. 11 Fuzzy logic controller dynamic responses to step input responsibility of the functionary is to establish decision-based rules by studying the system's behavior and the linguistic variables within the framework of the controlled system. We proposed an FLC controller in conjunction with system inverse dynamics, compensating for system nonlinearity and coupling effects: where, ̈ (f ) = ̈1 Considering the error from (28) (d) ] which means that the value of i is low and decreasing. Figure 6 show the member ship function for designed FLC, where Fig. 7 show the control surface for ̈i (f ) andφ i (f ) . Table 1 present the rule base table for

Sliding Mode Controller
SMC is basically a nonlinear adaptive control algorithm that gives a robust performance by using discontinuous control signals to slide the dynamic system along a hypersurface. The surface of the sliding mode s(t) ∈ R 4 is: where must satisfy Hurwitz condition, > 0 . The tracking error and its derivative value is: where q d (t) is the desired position signal.
From (32), we can see that if s(t) = 0 , then e(t) +̇e(t) = 0 , and we can get e(t) = e(0)exp(− t) . That is, when t → ∞ , Increasing big (IB) position tracking error and speed tracking error will tend to zero exponentially with value. and, To guarantee ṡs < 0 , and avoiding Chattering phenomenon due to high frequency switching of SMC, the proposed control signal designed as: while M (q) μ sat s δ , is a term responsible for the boundary layer. here, ∈ R 4×4 is the gain of the controller which is a diagonal matrix. Where, while the value of provided by solving (38) online with model.
where, δ is the boundary layer in which once the state trajectory reaches this layer, then it remains there.

Validation
Our proposed multi-section continuum robot dynamic model was validated using previously published data from Amouri et al. [30]. Amouri et al. present a 2-DOF cable-driven continuum robot dynamic model using system kinetic and potential energy. A comparison study was conducted using the PID controller developed by Amouri et al. on our proposed dynamic model. Table 2 displays the PID controller constants ( k p ,k i , andk d ) chosen by Amouri et al. Figure 8 illustrates a comparison between the dynamic responses for the bending angle ( i ) developed by our model and those of Amouri et al. Figure 8a, extracted using a data export program. The comparison reveals only slight differences between the two outputs, confirming the validity of our proposed model. Additionally, we compared ( i ) fluctuation response between our proposed work and that of Amouri   Figure 8b. The response comparison demonstrates similar fluctuations between the two systems, further validating our proposed model.

Simulation and Results
In order to simulate the continuum manipulator behavior with two bending sections, the dynamics model has been implemented using MATLAB 20201b. The material and geometric properties of the continuum manipulator under consideration are given in Table 3.

Step Response
The following control strategies try to move the two-section CR into the desired configuration space q(d), which is equivalent to the intended task space 0 (d) = 0 0 For SMC, the values for designed = 3I 4×4 , and = 9I 4×4 . The saturation function limit , are provided by solving (38) online with simulation. Where 4×4 is the identity matrix that belongs to R 4×4 . Figures 9, 10, 11 and 12 show complete dynamic responses for predesigned controllers, including error e (t) and   Table 4 shows the values of overshoot, rise-time, and settling-time.
Where Fig. 9 reveals that the PID controller with designed gain values, despite the exceptional settling-time t s achieved, the PID produces a relatively high overshoot followed by a small and short oscillation, which is caused by the lowest rising time t r of the controller, in its endeavors to achieve zero steady-state error. Due to the overshoot and oscillation, we must emphasize that using the PID in real-time applications is not the ideal choice.
As opposed to that, the FLC shows no overshoot and oscillation. And produces an acceptable rise time t r , while achieving zero steady-state error.
The robust nonlinear sliding mode controller shows an astounding response. where it produces zero overshoot and no oscillation, while achieving the lowest rising time and zero steady-state errors.

Trajectory Tracking Response
This subsection covers the inverse dynamic responses of the designed controller for a given trajectory. To this point, three simulation examples are performed in two and three-dimensional spaces, tracking a circular and sinusoidal trajectory in a 2D plane and a conical helix in a 3D plane.
The first example is a circular trajectory in the x-z plan, where CR tracks a perimeter of a circular path with a diameter of 50 mm in 6 s in the curvature planes 1 = 2 = 0. Figure 14a shows the different controller responses for tracking the circular trajectory. Figure 14b represents a magnified view, where PID shows an overshoot from the desired set point (start point for the circular path), where FLC shows more accurate tracking for trajectory and SMC achieves the efficient tracking with minimal tracking error.   Figure 15 shows the response of the configuration space parameter = 1 1 2 2 T for each controller.
For the second path, a sinusoidal trajectory consisting of two full waves with an inclined angle of 45 • from the x-axis in the x-z plane, where: Figure 16a shows the different controller responses for tracking the sinusoidal trajectory, where Fig. 16b represents a magnified view, where SMC shows great tracking for the trajectory with an error of less than 0.4 mm. Figure 17 Shows the configuration space response for sinusoidal trajectory.
For the three-dimensional space, a conical helix trajectory was simulated for 12 seconds path. Figures 18 and 19 shows the dynamic response for different controller for conical helix trajectory, where it reveals that the PID controller response suffers from a relatively high overshoot of 28.4% which leads to a high peak error; however, it rapidly handled the error, due to its lower settling time. The FLC provides a more precise response, with no overshoot and a reduced settling time of 1.56 s, resulting in an accurate tracking. The SMC shows the most responsive and effective controller, where it produces a zero percent overshoot, a rising time of less than 0.59 s, a settling time of less than 1.4 s the lowest among them all, and zero steady-state error, which leads to the most accurate tracking of trajectories.
Using the data developed during MATLAB Simulink simulation, which constituted from solving the kinematics and dynamics for each and every part and component of CR, a highly accurate dynamic animated simulation program developed using MATLAB environment, where it used to present a three-dimensional space simulation for each previously presented trajectory, plotting every configuration space variable along with CR dynamics and its workspace as shown in Fig. 20.
(Online Resource 1) represent created animated simulation for each and every trajectory, where developed MAT-LAB application Hyper-redundant Continuum robot simulator (HRCR simulator) was created to facilitate the simulation and dynamics plotting graphs.

Conclusion
For the two-link hyper-redundant continuum manipulator consisting of a primary backbone, a secondary backbone, and a series of disks, a kinematic model is derived from the piece-wise constant curvature assumption, and its dynamic model is derived using the Euler Lagrange representation. Three control strategies (inverse dynamic PID, intelligent FLC, and nonlinear SMC) were designed and applied for the dynamic model of the system in order to control the position and orientation of the CR for desired configuration values or along a given trajectory.
The inverse dynamics control rule, in particular, changes the closed-loop system into a new second-order linear system, while the fuzzy logic controller converts human experience and a thorough understanding of system dynamics into an intelligent, Robust controller, whereas the sliding mode control employs a discontinuous controller to compel the system to slide on the hypersurface with stability. All three close-loop dynamic control systems were implemented using MATLAB Simulink environment.
The first control strategy PID, which depends on its gains, is able to track the reference configuration space with high acceleration leading to a low rising time, which causes a relatively high overshoot followed by a small oscillation around desired configuration space values. Despite the fact that inverse dynamics control transforms the system into a decoupled linear second-order system where several linear control strategies can be investigated further to improve performance, the inverse dynamics calculation adds computational complexity and delay, which may affect the robot response in real-time applications. The second control strategy, FLC, on the other hand, eliminates both overshoot and oscillation, resulting in a relatively longer rising time. The third control strategy, SMC, achieves the best performance among them all, where it tracks the reference configuration space trajectory by an acceleration depending on the error Fig. 18 Tracking performance of different controllers for 3D conical helix signal, which leads to a low rise time, no oscillation, and zero % overshoot, resulting in the most accurate trajectory tracking performance. By switching the system between two states, it gives significant robustness to the system. This provides the system with a significant advantage in dealing with uncertainty and disturbance, which can be demonstrated in further research.
Future work focuses on designing and implementing a disturbance observer for handling model uncertainty and external disturbances. Experiments will be carried out using the specified dynamic control techniques and the disturbance observer.