Variable Curvature Modelling Method of Continuum Robots with Contraints

: The inherent compliance of continuum robots holds great promise in the fields of soft manipulation and safe human-robot interaction. This compliance reduces the risk of damage to the manipulated object and the surroundings. However, continuum robots have theoretically infinite degrees of freedom, and this high flexibility usually leads to complex deformations with external forces and positional constraints. How to describe this complex deformation is the main challenge for modelling continuum robots. In this study, we investigated a novel variable curvature modeling method for continuum robots, considering external forces and positional constraints. The robot configuration curve is described by the developed mechanics model, and then the robot is fitted to the curve. To validate the model, a 10-section continuum robot prototype with a length of 1 m was developed. The ability of the robot to reach the target points and track complex trajectories with load verified the feasibility and accuracy of the model. The ratio of the average position error of the robot endpoint to the robot length was less than 2.38%. This work may serve a new perspective for design analysis and motion control of continuum robots.


Introduction
Owing to intrinsic compliance, environmental adaptability, and operational safety, continuum robots have become an area of great interest in robotics [1]. A large number of researchers have conducted related research [2], [3], and many new continuum robots have been developed. Renda et al. [4] and Xu et al. [5] designed cable-driven continuum robots, Marchese et al. [6] and Tutcu et al. [7] developed pneumatic continuum robots. Gu et al. [8], [9] developed dielectric elastomer-driven continuum robots.
These continuum robots show great potential in a wide range of applications, such as medical equipment, unstructured environment exploration and soft manipulation [10]. However, continuum robots have a theoretically infinite number of degrees of freedom (DOF), and this high flexibility allows complex deformations of the robot in response to external forces and positional constraints. Therefore, the accurate and efficient modelling of continuum robots with external force remains challenging [11].
Different from the traditional rigid rod robot, the continuum robots achieved movement by deformating themselve. Therefore, the kinematics of continuum robots can be replaced by mechanical analysis. The continuum robots are usually discretized into a series of points when solving the model. A spatial curve can be determined by these points position mathematically. The continuum robots fitted to these spatial curves as closely as possible [12], [13]. The calculation methods of spatial curves mainly include constant curvature approach and variable curvature approach.
The constant curvature approach is a simplified approach to modelling continuum robots, which assumes that the curvature of curve is same between discrete points. The advantages of closed kinematics and easily solution make this method widely used of continuum robots modelling. In order to improve the accuracy and dexterity, Jones and Walker [14] developed a geometrical approach for modeling a constant curvature continuum robot, which was used to obtain a closed-form model. Freixedes et al. [15] established an optimization framework for continuum robots based on the assumption of constant curvature and deduced optimized structural parameters. They proposed the kinematic model to describe the deflection characteristics of the contact-assisted continuum robot, and carried out the experimental verification. Webster et al. [16] unified the kinematics and differential kinematics results of singlesegment and multi-segment continuous robots with constant curvature into the common coordinate system and symbol setting. Della Santina et al. [17] analyzed the flaws of constant models, considering an alternative state representation to solve these issues. Simulation cases are used to support the theoretical analysis. In order to simulate contact with the environment, Schiller et al. [18] [20] proposed a quantitative modeling method to obtain three-dimensional reconstructions of the configurations of a continuum robot with a variable curvature.
Gonthina et al [21] proposed a cross-sectional modeling method for variable curvature continuum robots based on Eulerian spiral curves. They compared the simulation results with the constant curvature method, proving that the proposed method is a significantly better match for various configurations of robot hardware. The modeling approach described above provides a good description of the forward and inverse kinematics of a soft continuous robot with variable curvature, but focuses mainly on the geometric description of the robot, without considering the effect of external forces. Therefore, additional visual or displacement sensors are needed to measure the robot configuration in real time. Some scholars have considered the influence of external forces and proposed some meaningful modeling methods.
Godage et al. [22] simulates the transient and steady-state dynamics of the continuum robot prototype based on the lumped parameter model. Renda et al. [23], [24]   The rest of this paper is organized as follows. A variable curvature model of continuum robots is developed and solved by the optimization algorithm in Section 2. In section 3, two sets of experiments (moving to the target points with load, and complex trajectory following with load) are performed to validate the model. Finally, the conclusion is presented in Section 4.

Main Model
The model proposed in this study is applicable to many drive types continuum robots. Because we used a pneumatic continuum robot for the experiments, so we modeled the pneumatic continuum robot as an example. In general, the design of pneumatic continuum robots has a central backbone, and several sections [26], as shown in Fig. 1(a). Each section is composed of three PAMs and a constraint disk, allowing for a 2-DOF bending motion. The motion of the robot can be achieved by deformation of the backbone, which can be bent in three-dimensional space by varying the length of the drive PAMs.
Therefore, the configurations of the backbone can be taken as the configurations of the robot, and the robot model established in this paper mainly takes the backbone as the object to analysis.
The backbone is discrete into n elements, and i N represents the ith nodes. There are m constraint disks, where the segment between analysis was performed on the micro element i Seg of the robot in the inertial coordinate system O   , as shown in Fig. 1 where  denotes the integral variable.
Suppose the rate of section angular displacement  with respect to arc coordinate s is  , Using the infinitesimal rotation theory of a rigid body, the relationship between the  and quaternions ( 1 q , 2 q , 3 q , 4 q ) can be derived: , and  is denoted as where L denotes the length of the robot, and ( , , ) x y z k k k represents the bending stiffness around the coordinate axis, and where d is the diameter of backbone, and E and G are the Young's modulus and shear modulus, respectively.
The variation of Eq. (7) is Writing Eq. (9) as matrix form, The external force potential energy p E can be expressed by the internal force F as In the O   , the force balance of i Seg is (14) Dividing the sides of Eq. (14) by s  , and considering the condition In the spindle coordinate system p xyz  , Eq. (15) needs to be rewritten as

Discretization of the equations
The quaternions at node Performing linear interpolation in The quaternions at nodes 1 i N  and i N are combined into an 8-order array, which is denoted as Then, Eq. (17) can be expressed as 4 1, The average and derivative of  in i Seg are denoted as i  and i   , respectively.
The average and derivative of  in i Seg are denoted as i  and i   , respectively.
Substituting (19) into (23), we can obtain 4 1, Eq. (12) can be discretized as Eq. (25) can be simplified as where i A is arranged diagonally, and each matrix is moved 4 rows and 4 columns to the upper left, then overlapping elements are added, and matrix A is formed. In the same way, i B is arranged horizontally, and each matrix is moved four columns to the left, then overlapping elements are added, and the matrix B is formed. Eq. (24) can be expressed as (28) where The projection of T in O   can be expressed by a quaternions as Then, i T can be expressed by i S as follows: The average and derivative of F in i Seg are denoted as i F and i  F , respectively.
Eq. (13) can be discretized as U are arranged horizontally, where each matrix is moved by four columns to the left to form matrix U. Then, Eq. (33) can be abbreviated as follows: Substituting Eq. (28) and (39) into Eq. (6), the variation of t E is as follows: Eqs. (41) and (42) are the discrete balance equations.

Boundary conditions
According to Eq. (1), the coordinate of any point on the backbone is When continuum robot operating, the endpoint usually needs to be controlled to move to a specified position, so the position of the endpoint needs to be limited. Let the include point of the robot be at the origin O, we obtain: is the desired coordinate of the endpoint.
In addition, the pose of the robot at the initial and end point also needs to be constrained:

Solving algorithm
The continuum robot model composed of Eqs.
By using the least square method, Eq. (46) can be transformed into an optimization problem: where u is the dimension of h, The particle swarm optimization (PSO) algorithm and Levenberg-Marquardt (LM) algorithm are used to solve (36). The PSO algorithm has strong global but weak local optimization ability, and LM algorithm is the opposite. Therefore, we first use PSO algorithm to find the appropriate iteration initial value, and then bring it into LM algorithm to optimize the solution. When finding the iteration initial value, there are two steps. In the first step the PSO algorithm is used, and in the second step the gradient descent method is used to modify the particle velocity term in PSO algorithm. In this way, a better initial iteration value can be obtained. The process is shown in Algorithm 1.

Experimental preparation
A 10 section prototype of a continuum robot driven by PAMs was developed, as shown in Fig. 3(a).
The structure of the prototype is shown in Fig. 3(b). The backbone is an elastic rod, and each section is driven by three PAMs. By changing the length of the PAMs, which can be achieved by inflating and deflating [27], [28], the robot can be controlled to bend in a three-dimensional space. The control system  The analysis of the robot is performed in four main spaces: actuator space, joint space, configuration space and task space. For the pneumatic continuum robot, actuator space refers to the air pressure of the PAMs, joint space refers to the length of the PAMs, configuration space is determined by the discrete point quaternions, and task space refers to the coordinates of the robot endpoint. The mapping from joint space to task space is forward kinematics ( for f ) and vice versa is inverse kinematics ( inv f ), as shown in Fig. 3(c). The length of the PAMs is proportional to the the internal air pressure, but the air pressure is difficult to control precisely. When the inflation velocity is certain, the length of the PAMs is proportional to the inflation/deflation time. Therefore, we control the length of the PAMs by controlling the inflation/deflation time, which can be achieved by controlling the switching time of the solenoid valves, and the relationship between them is shown in Fig. 4. Fig. 4 The relationship between the length of PAMs and the time of inflation/deflation. Different from traditional link robots, the model of continuum robot can usually obtain numerical solutions rathr than analytical solutions. When solving continuum robot model, the difference format and optimization algorithm will inevitably introduce calculation errors, and the number of discrete elements usually has a greater impact on the calculation errors. We compared the calculation errors of different discrete elements and found that when the number of elements increases from 5 to 20, the calculation errors gradually decreases, and when it increases from 20 to 50, the calculation error is basically stable.
At the same time, the calculation time is directly proportional to the number of elements. When the number of elements is 20, the calculation time is 5ms, as shown in Fig. 5. Therefore, 20 is choosed as the appropriate number of discrete elements.

Arriving at the target points
Continuum robots are often required to operate at certain specific locations. Based on the proposed model, experiments to reach the target points with a 50g load on the endpoint are performed first. Five target points in space are defined, and based on the calculations of the model, the robot is controlled to reach the target point, and the ratio of the Euclidean distance between the ideal and the actual position of the robot endpoint to the robot length is defined as the error.

Complex trajectory following
In practical application scenarios, soft continuum robots often need to move along certain trajectories.
Experiments to control the movement of the robot along a complex trajectory in space with a 50g load are carried out to verify the adaptability of the proposed model. A trajectory can be decomposed into a collection of discrete points, and only need to control the robot to move to these discrete points in sequence. A circular trajectory is defined, and it is discretized into 8 discrete nodes. The functional of the trajectory is 2 2 2 2 ( 75) 30 x y z     . Based on the calculation results of the model, the robot is controlled to move to these discrete nodes in turn.

Conclusion
Based on the principle of virtual work and vector mechanics, a continuum robot modeling method is proposed, and the exteral payload and position constraints are considered. The finite difference method are used to discretize the equations, and the PSO and LM optimization algorithm is used to numerically solve the model.