A Hybrid Controller for a Soft Pneumatic Manipulator Based on Model Predictive Control and Iterative Learning Control

Due to the outstanding characteristics of the large structural flexibility and strong dexterity of soft robots, they have attracted great attention. However, the dynamic modeling and precise control of soft robots face huge challenges. Traditional model-based and model-free control methods find it difficult to obtain a balance between complexity and accuracy. In this paper, a dynamic model of a three-chamber continuous pneumatic manipulator is established based on the modal method. Moreover, a hybrid controller integrating model predictive control (MPC) and iterative learning control (ILC) is proposed, which can simultaneously perform model parameter learning and trajectory tracking control. Experimental results show that the proposed control method can optimize the parameters of the dynamic model in real time with less iterations than the traditional model-free method and have good control performance in trajectory tracking experiments. In the future, the proposed dynamic model and the hybrid controller should be verified on a multi-section manipulator.


Introduction
Differently from rigid robots, soft continuum manipulators benefit from great compliance, lightweight, safer interaction, and a high power-to-weight ratio [1][2][3]. There has been increasing research on continuum robotics in recent years. They have a promising future and provide new capabilities in some fields, such as surgical applications [4,5] where continuum robots perform better when passing through narrow passages between organs for their flexible mechanisms and dexterous mobility [6]. There have been many attempts at soft robots in bionics design, such as octopus tentacles [7,8], manta [9], caterpillars [10], and elephant trunks [11,12]. Soft continuum manipulators are mainly divided into two categories according to their actuation methods: cable-driven manipulators [13,14] and pneumatically actuated manipulators [15,16]. Compared with general cable-driven continuum robots, pneumatic robots are more lightweight and compliant. This is due to the physical properties of the air used for both actuation and structure. However, the advantages of soft pneumatic continuum robots come with the trade-off of complexity in their modeling and control [17]. Due to the continuity, non-linear flexible material, and the unlimited number of degrees of freedom of the manipulator [18], it is challenging to establish its model.
The common assumptions of traditional rigid multi-segment robots do not apply to soft robots. At this stage, for soft manipulators, kinematics modeling is mainly studied [19,20]. OctArm by Walker et al. [21] is a successful implementation of a multi-section continuum robot using a pneumatic muscle actuator, according to the model proposed by Jones et al. [22]. However, the system dynamics is still a challenge and the using of different materials will make it even more difficult. Soft materials such as silicone, eco flex, rubber, or gel have high hyper-elasticity [23], which needs to be considered in the 1.
Taylor expansion was used to solve kinematics singular solution problems.

2.
Universal dynamic modeling for the three-chamber continuous pneumatic manipulator is presented based on the modal method.

3.
A hybrid controller integrating the model predictive control method and iterative learning control method was proposed.
The scientific flowchart of this paper is shown in Figure 1. The rest of this paper is organized as follows. Section 2 describes the designing and manufacturing method of the multi-chamber continuous soft manipulator. Section 3 describes the kinematics and dynamic model. Section 4 describes the derivation process of the model-based parameter adaptive learning control algorithm. In Section 5, simulations based on different control algorithms are conducted and compared. In Section 6, experiments are carried out and the results are analyzed. Finally, conclusions are made in Section 7.

System Description
The prototype is shown in Figure 2a. It contains chambers, inner layer, fiber, and the outer layer. The three chambers are separated by 120° intervals and their radii is equal. This structure enables the soft pneumatic manipulator to extend and contract, as well as bend in the space. Inspired by the design in [34], each pneumatic chamber is reinforced with a thin nylon fiber, forming a tight spiral on the outside of the cavity in the form of a diagonal. The radial expansion of the soft manipulator is limited and axial stretching is allowed. The purpose of covering the fiber layer with a silica layer is to prevent the fiber from slipping, which will affect the bending angle of the manipulator.
The manufacturing process of the soft pneumatic manipulator is shown in Figure 3. Firstly, silica gel (Ecoflex-0300) is used as the raw material, and the A-liquid and the Bliquid are mixed with the same quality. After stirring, they are put into a vacuum box to remove the dissolved air so as to avoid bubbles forming in the soft arm. Subsequently, the inner molds, as shown in Figure 3a, were coated with a release agent (Vaseline) to facilitate demolding. The silica gel is poured into the mold and placed in the heating box for solidification. After the manipulator is formed, the fiber layer is wound around it. The thickness of the fiber layer is 2 mm. Then it is put into the outer layer mold again, and the evenly mixed silica gel is poured.

System Description
The prototype is shown in Figure 2a. It contains chambers, inner layer, fiber, and the outer layer. The three chambers are separated by 120 • intervals and their radii is equal. This structure enables the soft pneumatic manipulator to extend and contract, as well as bend in the space. Inspired by the design in [34], each pneumatic chamber is reinforced with a thin nylon fiber, forming a tight spiral on the outside of the cavity in the form of a diagonal. The radial expansion of the soft manipulator is limited and axial stretching is allowed. The purpose of covering the fiber layer with a silica layer is to prevent the fiber from slipping, which will affect the bending angle of the manipulator.

length change trajectory Discussion
Problem and Target Methods Verifications and results

System Description
The prototype is shown in Figure 2a. It contains chambers, inner layer, fiber, outer layer. The three chambers are separated by 120° intervals and their radii i This structure enables the soft pneumatic manipulator to extend and contract, as bend in the space. Inspired by the design in [34], each pneumatic chamber is rei with a thin nylon fiber, forming a tight spiral on the outside of the cavity in the fo diagonal. The radial expansion of the soft manipulator is limited and axial stret allowed. The purpose of covering the fiber layer with a silica layer is to prevent t from slipping, which will affect the bending angle of the manipulator.
The manufacturing process of the soft pneumatic manipulator is shown in F Firstly, silica gel (Ecoflex-0300) is used as the raw material, and the A-liquid and liquid are mixed with the same quality. After stirring, they are put into a vacuum remove the dissolved air so as to avoid bubbles forming in the soft arm. Subsequen inner molds, as shown in Figure 3a, were coated with a release agent (Vaseline) to f demolding. The silica gel is poured into the mold and placed in the heating solidification. After the manipulator is formed, the fiber layer is wound around thickness of the fiber layer is 2 mm. Then it is put into the outer layer mold again, evenly mixed silica gel is poured.  The manufacturing process of the soft pneumatic manipulator is shown in Figure 3. Firstly, silica gel (Ecoflex-0300) is used as the raw material, and the A-liquid and the Bliquid are mixed with the same quality. After stirring, they are put into a vacuum box to remove the dissolved air so as to avoid bubbles forming in the soft arm. Subsequently, the inner molds, as shown in Figure 3a, were coated with a release agent (Vaseline) to facilitate demolding. The silica gel is poured into the mold and placed in the heating box for solidification. After the manipulator is formed, the fiber layer is wound around it. The thickness of the fiber layer is 2 mm. Then it is put into the outer layer mold again, and the evenly mixed silica gel is poured.
represents the normalized position parameter. When located at the b 0   . When located at the top, 1 Figure 3. The manufacturing process of the multi-chamber soft pneumatic manipulator. (a) Design of the inner layer mold. It is composed of 3 chambers with a diameter of 6mm and an inner wall with a diameter of 25 mm. The inner wall is spirally convex, which is convenient for subsequent Fiber winding. The top fixing is used for sealing and positioning and the mold is made by 3D printing. (b) The manufacturing process of the soft manipulator, including the inner and outer mold design.

Kinematic Modeling
Based on the assumption of constant curvature, a continuous curved section can be described by three independent parameters {θ, ϕ, λ}. As shown in Figure 4, θ indicates the rotating angle, which is the angle between the bending plane and the x-axis. ϕ indicates the bending angle. λ indicates the bending radius. Note the length variation of the three chambers {l 1 , l 2 , l 3 }, the original length L 0 and the radius r of soft manipulator. The parameter relationship of the bending section can be written as: As for soft pneumatic manipulators, robot mapping describes the coordinate transformation between task space and configuration space. While the traditional Denavit-Hartenberg [35] is not applicable due to the flexibility of the proposed manipulator, the homogeneous transformation matrix can be derived by applying the standard rotational and translational matrices. As shown in Figure 4, according to the bending parameters {θ, ϕ, λ}, the homogeneous transformation matrix shows the relationship between the base coordinate system and the tip end coordinate system of the manipulator. It can be written as: where ξ ∈ [0, 1] represents the normalized position parameter. When located at the base, ξ = 0. When located at the top, ξ = 1. The generalized coordinates q = [l 1 l 2 l 3 ] T . R z and  If the pressure in the three chambers is the same, the manipulator will stretch along the axial direction. At this time, both the bending angle and the bending radius tend to infinity, which leads to singularities in the kinematics solution. On the other hand, from the above Equations (1)-(3), it can be seen that bringing      , respectively, and then the elements in the homogeneous transformation matrix, will produce many triangular and denominator terms. As a result, the calculation will be very complicated. Thus, in this paper, the high-order Taylor expansion of each element will be expressed as a numerically stable multivariate polynomial form. The modal matrix can be obtained as follows: where . It can be found that this model is completely composed of numerically stable multivariate polynomials, which eliminates the problem of singularities and reduces the complexity of calculations, making it easy to find derivation. If the pressure in the three chambers is the same, the manipulator will stretch along the axial direction. At this time, both the bending angle and the bending radius tend to infinity, which leads to singularities in the kinematics solution. On the other hand, from the above Equations (1)-(3), it can be seen that bringing {l 1 , l 2 , l 3 } into {θ, ϕ, λ}, respectively, and then the elements in the homogeneous transformation matrix, will produce many triangular and denominator terms. As a result, the calculation will be very complicated. Thus, in this paper, the high-order Taylor expansion of each element will be expressed as a numerically stable multivariate polynomial form. The modal matrix can be obtained as follows: where u = l 2 1 + l 2 2 + l 2 3 − l 1 l 2 − l 1 l 3 − l 2 l 3 . It can be found that this model is completely composed of numerically stable multivariate polynomials, which eliminates the problem of singularities and reduces the complexity of calculations, making it easy to find derivation.

Dynamic Modeling
When establishing the dynamic model, the curved section is divided into countless small slices. The thickness of each slice is s · ∆ξ, where s = λϕ represents the arc length of the central axis. Note the modal rotation matrix as R(q, ξ) ∈ R 3×3 and the modal position vector as P(q, ξ) ∈ R 3 .
The linear velocity at the center of a slice is: and the angular velocity with respective to the base coordinate axes is: Thus, the velocity vector of a single slice is: where J b ξ (q) is the modal Jacobian matrix. Note the mass of continuum section is m, the inertial matrix of any slice is: According to the velocity vector, the kinetic energy of any slice can be derived by: and with the assumption of equal density, the total kinetic energy is then calculated by integrating all slice kinetic energies as: where the generalized mass matrix is defined as in Equation (13) Due to the particularity of the silicone material, the potential energy of the soft manipulator consists of three parts: gravitational potential energy, elastic potential energy, bending potential energy. The gravitational potential energy can be calculated by modal position vector as: The elastic potential energy is caused by the elongation or contraction of the silicone material chamber. It can be described as: where K e = diag{k e1 , k e2 , k e3 } is the elastic stiffness matrix of the joint space and k ei is the elastic stiffness coefficient of each chamber.
To simplify the calculation of the bending potential energy, the assumption is cited that the bending potential energy is proportional to the square of the bending angle. Thus, the bending potential energy of soft continuum manipulator can be derived by where K b is the bending stiffness coefficient of the arm which is assumed to be constant. In summary, the potential energy of the soft continuum manipulator is The Lagrangian function of the system is defined by L = K − P, and the equation of motion can be obtained by applying the Euler-Lagrange equation formulated in Equation (18).
The newly derived dynamics of the soft pneumatic manipulator is obtained in matrix form as: where C q, . q is the centrifugal and Coriolis forces matrix that can be derived by: and G(q) is the conservative forces vector that can be derived by G(q) = ∂p/∂q i .

Control Algorithm
In the process of establishing the mathematical model, in order to simplify the modeling difficulty, assumptions such as the constant curvature [36] and equal linear density were made. However, due to the non-linearity of the soft material itself, the relationship between the pressure and the length change of the chambers based on the dynamic model is not completely accurate. In order to further improve the accuracy of the model, a parameter adaptive learning control method based on the dynamic model is proposed. The control algorithm diagram is shown in Figure 5 as follows. The Lagrangian function of the system is defined by L = K − P, and the equation o motion can be obtained by applying the Euler-Lagrange equation formulated in Equation (18).
The newly derived dynamics of the soft pneumatic manipulator is obtained in matrix form as: , (19 where   , C q q  is the centrifugal and Coriolis forces matrix that can be derived by:

Control Algorithm
In the process of establishing the mathematical model, in order to simplify the mod eling difficulty, assumptions such as the constant curvature [36] and equal linear density were made. However, due to the non-linearity of the soft material itself, the relationship between the pressure and the length change of the chambers based on the dynamic mode is not completely accurate. In order to further improve the accuracy of the model, a pa rameter adaptive learning control method based on the dynamic model is proposed. The control algorithm diagram is shown in Figure 5 as follows.

Parameters of the Control Model
According to the dynamic model of the continuous section, the motion equation of a single chamber after decoupling can be equivalently expressed as: (21 where J is the equivalent moment of inertia, D is the equivalent damping coefficient, K is the equivalent stiffness, a is the cross-sectional area. Converted to the form of the state space equation: (22 Through the Laplace transform, the system transfer function is obtained as:

Parameters of the Control Model
According to the dynamic model of the continuous section, the motion equation of a single chamber after decoupling can be equivalently expressed as: where J is the equivalent moment of inertia, D is the equivalent damping coefficient, K is the equivalent stiffness, a is the cross-sectional area. Converted to the form of the state space equation: Through the Laplace transform, the system transfer function is obtained as: However, in a real control situation, the system runs discretely. In this paper, a zeroorder holder is used to discretize the continuous system, and the sampling time is set as 0.1 s. The system pulse transfer function is derived through z transformation: Thus, the control model parameters can be expressed as p = p 1 p 2 p 3 p 4 T .

Hybrid Controller
The soft pneumatic manipulator is assumed to be a linear time-invariant system. A linear MPC controller will be used. The cost function of the model predictive control method is: Subject to: where y re f is the reference length of the chamber, Q and R are the output and input weight factors, respectively, H is the prediction horizon. The input pressure of the soft pneumatic manipulator corresponding to the system parameters at the sampling time is optimized through model predictive control. The optimized input air pressure is used as the feedback pressure.
Iterative learning control has no requirements on the model, but the dynamic model parameters of the soft pneumatic manipulator derived previously can be used as the starting point of the iteration, which can greatly accelerate the iteration speed and reduce the number of iterations.
where E r i = y re f i − y i representing error between reference length and actual chamber length; k is sampling instant; i is the number of iterations. The feedforward input coefficient λ is a constant.

Model Parameters Iterative Learning Law
The relationship between model input and output is expressed as: where y (kt) is the model estimated chamber length. Through the inverse z transformation, we can obtain: where k = 2, 3, . . . , N, and N = T/T s , T s is sampling time. Then we have: Input the same pressure, the error between the actual change of chamber length and the change predicted by the model is: where y (k) = ϕ(k) T p i is the estimated length change, y(k) = ϕ(k) T p t is the true length change and p t is the true model parameter.
The estimated model is constantly being updated by minimizing the error E t i , as shown by the following quadratic function. Then the Gauss-Newton method [37] is applied to update the model parameters. min

Parameter-Adaptive-Learning Control Algorithm
A complete algorithm description about parameter-adaptive-learning control based on the soft continuum section dynamic model is presented in Algorithm 1. Note ε r as the error tolerance between the reference length and the actual length. ε t is the error tolerance between the actual length and the model estimated length.

Algorithm 1. Model parameters learning adaptive control
Step1: Initialize the model parameter M 0 , Plan length trajectory y re f ; Step2: Obtain the initial pressure u 0 based on MPC; Step3: Compare with length tolerance ε t and ε r .
Step4: Obtain u m i trough model parameters iterative learning law and MPC Step5: Obtain u l i based trough ILC Step6: Obtain total input pressure u i , real chamber length y t i and model estimated length y p i trough real system and estimated model, respectively. Then go to step 3.
To summarize, in this section, a model-based parameter adaptive learning control algorithm is proposed for the soft pneumatic manipulator. The previously established dynamics model of the manipulator is used. A hybrid controller integrating MPC and ILC is proposed and integrated into the parameter learning control law. Real-time optimization and update of model parameters is realized in the form of online learning. It can improve the accuracy of the dynamic model and achieve trajectory tracking with a small number of iterations. This control algorithm is verified by simulation, and results show that the control proposed algorithm has faster convergence and higher control precision.

Simulation
In order to verify the effect of the control algorithm on the parameter learning of the dynamic model and the desired chamber length change trajectory tracking, a singlechamber soft pneumatic manipulator is used to carried out the verification simulation.

Control Algorithm Simulation
The dynamic parameters of the soft continuous section are shown in Table 1. Then we can obtain the initial continuous time model: The discretization of the control system is realized by using the zero-order holder with the sampling time of 0.1 s. The discretization result is: and the initial parameter of the model is p init = [0.016 0.015 −0.0468 0.837] T . In the simulation, we assume that the real continuous-time system is: Similarly, the discretization result is: The reference length which represents the change of the chamber length is set to be the tracking target. The minimum jerk theory [38] is used to plan the motion trajectory. When jerk movement is minimized, higher-order polynomials will be reduced to fifth-order polynomials. The reference trajectory can be expressed as: The relevant parameters of the control algorithm during the simulation are shown in Table 2. The simulation results are shown in Figure 6. We can see that using the initial model parameters obtained by the dynamic model as the starting point, the model parameters gradually converge to the real model parameters during the adaptive control learning process. This result effectively proves the convergence of the control algorithm. It also shows that the real system parameters can be obtained through the rapid iterative learning of the control algorithm and if we take the dynamic model parameters as the initial values, the trajectory tracking simulation can be achieved with better performance.
In addition, as the number of iterations increases, the input pressure is brought into the real model and the estimated model, respectively. It can be seen from Figure 7a that the chamber length of the estimated model gradually tends to the real chamber length after model parameter iterative learning. The root mean square error (RMSE) between the estimated chamber length and the real length is shown in Figure 7b. After 15 iterations, the RMSE approaches 0. This result indicates that the proposed control algorithm can effectively improve the accuracy of the model parameters. In Figure 7c we can see that the real chamber length gradually tends to the reference length. The RMSE between them is presented in Figure 7d which shows that, after 20 iterations, the RMSE is 0.57 mm.
Compared to the 150 mm long manipulator, this RMSE value is very small. This result indicates the good convergence performance of the control algorithm in length control.
Sensors 2023, 23, 1272 11 of 20 presented in Figure 7d which shows that, after 20 iterations, the RMSE is 0.57 mm. Compared to the 150 mm long manipulator, this RMSE value is very small. This result indicates the good convergence performance of the control algorithm in length control.

Comparison with Traditional Model-Free Algorithms
A diagram of a traditional model-free closed-loop feedback control algorithm is shown in Figure 8 (classic PID controller is selected here). Assume that the rest of the simulation parameters are the same as Table 2 and reference trajectory is also set as Equation (37). The simulation results are shown in Figure 9. It can be seen from Figure 9a that in the first 5 control cycles, the error between the real chamber length and the reference chamber length is large. The real chamber length fluctuates greatly. RMSE is about 5 mm, which is always called the "overshoot" phenomenon. Figure 9b shows the RMSE between the reference chamber length and the real chamber length, which are obtained by the classic PID control algorithm and the proposed control algorithm, respectively. The RMSE value of the proposed parameter adaptive learning control algorithm is significantly smaller than the PID controller in the beginning. Moreover, the proposed control algorithm needs less iterations to achieve good trajectory tracking performance. After the 5th iteration of learning, the RMSE is 0.4 mm, which is lower than the RMSE value of PID method after 25 presented in Figure 7d which shows that, after 20 iterations, the RMSE is 0.57 mm. Compared to the 150 mm long manipulator, this RMSE value is very small. This result indicates the good convergence performance of the control algorithm in length control.

Comparison with Traditional Model-Free Algorithms
A diagram of a traditional model-free closed-loop feedback control algorithm is shown in Figure 8 (classic PID controller is selected here). Assume that the rest of the simulation parameters are the same as Table 2 and reference trajectory is also set as Equation (37). The simulation results are shown in Figure 9. It can be seen from Figure 9a that in the first 5 control cycles, the error between the real chamber length and the reference chamber length is large. The real chamber length fluctuates greatly. RMSE is about 5 mm, which is always called the "overshoot" phenomenon. Figure 9b shows the RMSE between the reference chamber length and the real chamber length, which are obtained by the classic PID control algorithm and the proposed control algorithm, respectively. The RMSE value of the proposed parameter adaptive learning control algorithm is significantly smaller than the PID controller in the beginning. Moreover, the proposed control algorithm needs less iterations to achieve good trajectory tracking performance. After the 5th iteration of learning, the RMSE is 0.4 mm, which is lower than the RMSE value of PID method after 25

Comparison with Traditional Model-Free Algorithms
A diagram of a traditional model-free closed-loop feedback control algorithm is shown in Figure 8 (classic PID controller is selected here). Assume that the rest of the simulation parameters are the same as Table 2 and reference trajectory is also set as Equation (37).
cates the good convergence performance of the control algorithm in length control.

Comparison with Traditional Model-Free Algorithms
A diagram of a traditional model-free closed-loop feedback control algorith shown in Figure 8 (classic PID controller is selected here). Assume that the rest o simulation parameters are the same as Table 2 and reference trajectory is also set as E tion (37). The simulation results are shown in Figure 9. It can be seen from Figure 9a that i first 5 control cycles, the error between the real chamber length and the reference cham length is large. The real chamber length fluctuates greatly. RMSE is about 5 mm, whi always called the "overshoot" phenomenon. Figure 9b shows the RMSE between the erence chamber length and the real chamber length, which are obtained by the classic control algorithm and the proposed control algorithm, respectively. The RMSE valu the proposed parameter adaptive learning control algorithm is significantly smaller the PID controller in the beginning. Moreover, the proposed control algorithm needs iterations to achieve good trajectory tracking performance. After the 5th iteration of le ing, the RMSE is 0.4 mm, which is lower than the RMSE value of PID method aft The simulation results are shown in Figure 9. It can be seen from Figure 9a that in the first 5 control cycles, the error between the real chamber length and the reference chamber length is large. The real chamber length fluctuates greatly. RMSE is about 5 mm, which is always called the "overshoot" phenomenon. Figure 9b shows the RMSE between the reference chamber length and the real chamber length, which are obtained by the classic PID control algorithm and the proposed control algorithm, respectively. The RMSE value of the proposed parameter adaptive learning control algorithm is significantly smaller than the PID controller in the beginning. Moreover, the proposed control algorithm needs less iterations to achieve good trajectory tracking performance. After the 5th iteration of learning, the RMSE is 0.4 mm, which is lower than the RMSE value of PID method after 25 control cycles. Most importantly, the accuracy of the model parameters is gradually adjusted and perfected by the proposed control algorithm, which cannot be achieved by the traditional model-free closed-loop feedback control algorithm.

Effect of Initial Model Parameter Value
To analyze the influence of the initial value of the model parameters on the performance of the control algorithm, simulations withs different initial model parameters are It can be seen from Figure 10 that for different initial values of model parameters, the algorithm can guarantee the convergence to the real model parameters. When init p is used as the initial model value, the convergence speed is the fastest. After the 12th iteration, it has almost converged to the real model parameters and satisfies the tolerance. However, under the other initial values, at least 20 iterations are required to meet the tolerance. This indicates that the establishment of an accurate initial values of model parameters can speed up the convergence of the control algorithm and reduce the amount of iterative learning.

Effect of Initial Model Parameter Value
To analyze the influence of the initial value of the model parameters on the performance of the control algorithm, simulations withs different initial model parameters are carried out. Three different parameters are set, such as p 1 = [1; 1; 1; 1] T , p 2 = [−1; −1; −1; −1] T , and p 3 = p init , where p init is the dynamic model parameter. Taking p 3 as a detailed example, all control parameters are the same as Table 2 except the initial model value.
It can be seen from Figure 10 that for different initial values of model parameters, the algorithm can guarantee the convergence to the real model parameters. When p init is used as the initial model value, the convergence speed is the fastest. After the 12th iteration, it has almost converged to the real model parameters and satisfies the tolerance. However, under the other initial values, at least 20 iterations are required to meet the tolerance. This indicates that the establishment of an accurate initial values of model parameters can speed up the convergence of the control algorithm and reduce the amount of iterative learning.
(a) (b) Figure 9. PID control simulation results. (a) shows the true chamber length, which can also tend to the reference length after more iterations. (b) the RMSE between the PID control algorithm and the proposed control algorithm.

Effect of Initial Model Parameter Value
To analyze the influence of the initial value of the model parameters on the performance of the control algorithm, simulations withs different initial model parameters are It can be seen from Figure 10 that for different initial values of model parameters, the algorithm can guarantee the convergence to the real model parameters. When init p is used as the initial model value, the convergence speed is the fastest. After the 12th iteration, it has almost converged to the real model parameters and satisfies the tolerance. However, under the other initial values, at least 20 iterations are required to meet the tolerance. This indicates that the establishment of an accurate initial values of model parameters can speed up the convergence of the control algorithm and reduce the amount of iterative learning.

Effect of Reference Length Change Trajectory
To analyze the influence of the reference length change trajectory, we set different types of reference trajectories to obtain the length change of the manipulator's chamber. All control parameters are the same as Table 2, except the reference trajectories.
In the simulation results in Figure 11, a trajectory tracking effect under the step length change, the rectangular wave length change, and the sine wave length change are, respectively, shown. As shown in Figure 11b,d,f, the RMSE under the sine wave reference trajectory becomes almost 0 after 15 iterations, and the RMSE under the step signal reference trajectory also changes from 5.3 mm to 1 mm after 15 iterations. Even if the step chamber length change is used as the reference trajectory, the tracking error can still gradually converge towards the decreasing direction with iterations, and the error is reduced by 72%. It is further verified that the proposed control algorithm can maintain convergence under different reference trajectories. model parameters will converge to the real model parameters at different speeds. Taking the dynamic model parameters as the initial values can effectively reduce the number of iterations.

Effect of Reference Length Change Trajectory
To analyze the influence of the reference length change trajectory, we set different types of reference trajectories to obtain the length change of the manipulator's chamber. All control parameters are the same as Table 2, except the reference trajectories.
In the simulation results in Figure 11, a trajectory tracking effect under the step length change, the rectangular wave length change, and the sine wave length change are, respectively, shown. As shown in Figure 11b,d,f, the RMSE under the sine wave reference trajectory becomes almost 0 after 15 iterations, and the RMSE under the step signal reference trajectory also changes from 5.3 mm to 1 mm after 15 iterations. Even if the step chamber length change is used as the reference trajectory, the tracking error can still gradually converge towards the decreasing direction with iterations, and the error is reduced by 72%. It is further verified that the proposed control algorithm can maintain convergence under different reference trajectories.

Experiments
In this section, we conduct some preliminary experiments for the single-chamber and the three-chamber soft pneumatic manipulators to verify our control method. Figure 12 shows the experimental setup of the soft manipulators. The upper computer is used to run control algorithms, dynamic model solving, and other scenarios that require higher computing power. The lower computer uses the advantages of the MCU to control and drive the proportional valve.

Experiments
In this section, we conduct some preliminary experiments for the single-chamber and the three-chamber soft pneumatic manipulators to verify our control method. Figure 12 shows the experimental setup of the soft manipulators. The upper computer is used to run control algorithms, dynamic model solving, and other scenarios that require higher computing power. The lower computer uses the advantages of the MCU to control and drive the proportional valve.

Verification Experiment of the Control Algorithm for Single-Chamber Manipulator
The single-chamber soft manipulator has only one degree of freedom. The change of the manipulator's length is planned. The corresponding expressions for setting the desired length change trajectory are the same as Equation (34). Set the sampling period T = 8 s, and the sampling time T s = 0.2 s. The setting of the control algorithm parameters is consistent with the simulation. Take the single-chamber soft manipulator dynamic model parameter p init = [0.016 0.015 − 0.0468 0.837] T as the initial value of the model for iterative learning control.
Experiment results are shown in Figure 13. As shown in Figure 13a, the manipulator begins to elongate after being pressurized. It can be seen from Figure 13c that in the first sampling period, the maximum chamber length change is 21.3 mm, which is 8.7 mm away from the expected length change peak 30 mm. The relative error is. As the number of iterations of the control algorithm increases, the real chamber length change keeps approaching the desired length change trajectory. After 10 iterations, the maximum length of the chamber change is 28.3 mm, which is 1.7 mm away from the expected peak trajectory, and the error is reduced to 5%. Figure 13d shows the RMSE of the real chamber length and the expected trajectory at each sampling time. It can be seen that RMSE has been decreasing from the initial 7.3 mm, and when the sampling ends after 16 iterations, it becomes 0.9 mm, which verifies the expected trajectory tracking performance of the proposed control algorithm. Similarly, Figure 13e,f show the performance of the control algorithm model parameter learning. The model estimated chamber length continues to tend to the actual chamber length, and the corresponding RMSE continues to decrease.  Figure 12. Experimental setup. The control algorithm is implemented on host computer (64-bi erating system, i5 CPU, and 16 GB RAM) based on ROS operating system. The pressure of the ch ber is controlled by the SMC electric proportional valve (ITV1050). The NDI Aurora electromag positioning system is used to capture the end pose of the actuator and communicate with the through a serial port. The PCB broad is mainly used for hardware control drive.

Verification Experiment of the Control Algorithm for Single-Chamber Manipulator
The single-chamber soft manipulator has only one degree of freedom. The chang the manipulator's length is planned. The corresponding expressions for setting the sired length change trajectory are the same as Equation (34) Experiment results are shown in Figure 13. As shown in Figure 13a, the manipu begins to elongate after being pressurized. It can be seen from Figure 13c that in the sampling period, the maximum chamber length change is 21.3 mm, which is 8.7 mm a from the expected length change peak 30 mm. The relative error is. As the numbe iterations of the control algorithm increases, the real chamber length change keeps proaching the desired length change trajectory. After 10 iterations, the maximum le of the chamber change is 28.3 mm, which is 1.7 mm away from the expected peak tr tory, and the error is reduced to 5%. Figure 13d shows the RMSE of the real cham length and the expected trajectory at each sampling time. It can be seen that RMSE been decreasing from the initial 7.3 mm, and when the sampling ends after 16 iterati it becomes 0.9 mm, which verifies the expected trajectory tracking performance of the posed control algorithm. Similarly, Figure 13e,f show the performance of the contro gorithm model parameter learning. The model estimated chamber length continue tend to the actual chamber length, and the corresponding RMSE continues to decreas

Verification Experiment of the Control Algorithm for Multi-Chamber Manipulator
Compared with the single-chamber actuator, the three-channel soft manipulator can perform linear contraction and bending motions in the space. We plan the trajectory of the three chamber lengths in joint space. As shown in Figure 14a, each chamber uses a differ-

Verification Experiment of the Control Algorithm for Multi-Chamber Manipulator
Compared with the single-chamber actuator, the three-channel soft manipulator can perform linear contraction and bending motions in the space. We plan the trajectory of the three chamber lengths in joint space. As shown in Figure 14a, each chamber uses a different form of step signal as the desired length change trajectory, and the peak change length of each chamber is 30 mm, 20 mm, and 10 mm, respectively. For the convenience of display, the expected length change of each chamber is brought into Equation (2) to obtain the expected bending angle of the soft arm under the assumption of equal curvature segments.
The experimental results are shown in Figure 14. It can be seen from Figure 14c that in the first iteration when the dynamic model is used as the initial value of the soft arm model, the peak bending angle is 102 • . Compared with the expected trajectory peak bending angle of 126 • , the error between the real bend angle and the desired bend angle (obtained through desired chamber lengths) is 20.1%. As shown in Figure 14d, the RMSE of the first iteration in the multi-chamber manipulator is 18.97 • . With the iterative calculation of the control program, the error gradually decreases. This result shows that the three-chamber soft manipulator is complicated, which makes it difficult to model the dynamics model and the accuracy is low. However, with continuous iterative operation, under the effect of the parameter adaptive learning control algorithm, the actual bending angle is basically the same as the expected trajectory trend, and the RMSE in each sampling period is gradually reduced. The RMSE value of the dynamics model changes from 18.97° to 1.92° after 10 iterations. It is reduced by 87% from the original value. This result verified the effectiveness of the model-based parameter adaptive learning control algorithm for the three-chamber soft manipulator.

Trajectory Tracking Experiment for the Multi-Chamber Manipulator
In this experiment, trajectory planning will be carried out in the task space, and the expected trajectory length of each chamber in the joint space will be calculated through the inverse kinematics solution. The experiment is carried out in the x-y plane to plan the equilateral triangle trajectory and the side length is 125 mm.
The experiment process is shown in Figure 15a. The end position is recorded by the NDI electromagnetic positioning system. It can be seen in Figure 15b that the experimental trajectory of the three-chamber continuous soft manipulators basically matches the desired equilateral triangle trajectory, the maximum absolute error is within 15 mm, and the average error is 4.7 mm, compared to the 150-mm-long manipulator. This result indicates This result shows that the three-chamber soft manipulator is complicated, which makes it difficult to model the dynamics model and the accuracy is low. However, with continuous iterative operation, under the effect of the parameter adaptive learning control algorithm, the actual bending angle is basically the same as the expected trajectory trend, and the RMSE in each sampling period is gradually reduced. The RMSE value of the dynamics model changes from 18.97 • to 1.92 • after 10 iterations. It is reduced by 87% from the original value. This result verified the effectiveness of the model-based parameter adaptive learning control algorithm for the three-chamber soft manipulator.

Trajectory Tracking Experiment for the Multi-Chamber Manipulator
In this experiment, trajectory planning will be carried out in the task space, and the expected trajectory length of each chamber in the joint space will be calculated through the inverse kinematics solution. The experiment is carried out in the x-y plane to plan the equilateral triangle trajectory and the side length is 125 mm.
The experiment process is shown in Figure 15a. The end position is recorded by the NDI electromagnetic positioning system. It can be seen in Figure 15b that the experimental trajectory of the three-chamber continuous soft manipulators basically matches the desired equilateral triangle trajectory, the maximum absolute error is within 15 mm, and the average error is 4.7 mm, compared to the 150-mm-long manipulator. This result indicates that the proposed control method has good control accuracy of the soft pneumatic manipulator.

Discussion
From the above experiments we can see that the proposed hybrid controller has good performance in controlling the length change, bend angle, and tip end position of sof pneumatic manipulators. However, the proposed controller needs to combine the MPC and ILC together, which is more complicated than most used model-free controller o model-based controller. In this paper, the proposed method is only verified on a single section pneumatic manipulator, its effect on a multi-section manipulator is unknown.

Conclusions
In this paper, the structural design and manufacturing process of a soft pneumati manipulator is presented and a corresponding dynamic model based on modal method i derived which solves the problem of instability caused by the singular points in the CP kinematics. Moreover, to control the manipulator in real time and with good accuracy, hybrid controller integrating MPC and ILC is proposed and verified through simulation and experiments. The simulation results show that the more accurate initial values o model parameters obtained by the dynamic model, the faster the convergence of the con trol algorithm. The experimental results show that the proposed hybrid control algorithm has good performance in controlling the manipulators' length. During the chamber length control experiment, the RMSE value of the dynamics model changes from 18.97° to 1.92 after 10 iterations, which is an 87% reduction from the original value. In the equilatera triangle trajectory experiment, the maximum tracking error is 15 mm, and the averag tracking error is 4.7 mm. This result indicates that the proposed control method has good control accuracy of the soft pneumatic manipulator. To sum up, the proposed dynami model can be used when there exist kinematics singular solution problems, and the hybrid controller can be used when the convergence of the control algorithm is low.
Although the proposed dynamic model and the hybrid controller have been verified through experiments, the experiments are relatively simple and the prototype is a single section manipulator. In the future, the application of the proposed method on the contro of the soft multi-section manipulator will be verified and more complicated experiment will be conducted.

Discussion
From the above experiments we can see that the proposed hybrid controller has good performance in controlling the length change, bend angle, and tip end position of soft pneumatic manipulators. However, the proposed controller needs to combine the MPC and ILC together, which is more complicated than most used model-free controller or modelbased controller. In this paper, the proposed method is only verified on a single-section pneumatic manipulator, its effect on a multi-section manipulator is unknown.

Conclusions
In this paper, the structural design and manufacturing process of a soft pneumatic manipulator is presented and a corresponding dynamic model based on modal method is derived which solves the problem of instability caused by the singular points in the CP kinematics. Moreover, to control the manipulator in real time and with good accuracy, a hybrid controller integrating MPC and ILC is proposed and verified through simulations and experiments. The simulation results show that the more accurate initial values of model parameters obtained by the dynamic model, the faster the convergence of the control algorithm. The experimental results show that the proposed hybrid control algorithm has good performance in controlling the manipulators' length. During the chamber length control experiment, the RMSE value of the dynamics model changes from 18.97 • to 1.92 • after 10 iterations, which is an 87% reduction from the original value. In the equilateral triangle trajectory experiment, the maximum tracking error is 15 mm, and the average tracking error is 4.7 mm. This result indicates that the proposed control method has good control accuracy of the soft pneumatic manipulator. To sum up, the proposed dynamic model can be used when there exist kinematics singular solution problems, and the hybrid controller can be used when the convergence of the control algorithm is low.
Although the proposed dynamic model and the hybrid controller have been verified through experiments, the experiments are relatively simple and the prototype is a singlesection manipulator. In the future, the application of the proposed method on the control of the soft multi-section manipulator will be verified and more complicated experiments will be conducted.