Stiffness Analysis of a Pneumatic Soft Manipulator Based on Bending Shape Prediction

Soft manipulators can perform continuous operations due to their inherent compliance and dexterity, thus enabling safe interactions and smooth movements in confined environments. However, high compliance usually means low load capacity. It is important for a soft manipulator to possess proper flexibility while maintaining an acceptable stiffness to widen its applications. This paper has hence devoted efforts to a kind of variable stiffness mechanism for a soft manipulator actuated by pneumatic artificial muscles (PAMs). Due to the combination of contractile and extensor PAMs, the manipulator is able to vary its stiffness independently from the configuration. The stiffness characteristics of the soft manipulator are quantitatively analyzed by bending shape prediction under different loading and inflation conditions, and the prediction is built upon a nonlinear statics model coupled with PAM nonlinearity and the Cosserat theory. In addition, experimental measurements are conducted to further validate the expected performance of the manipulator design. The experimental and verified theoretical analysis results indicate that the manipulator shape and stiffness are greatly affected by the pressure variation of PAMs, realizing a large bending space with a high output force. The variable stiffness design obviously increases the manipulator’s ability to resist additional interference at the same position.


I. INTRODUCTION
Soft robots (continuum robots made from soft materials) inspired by biological features [1]- [4], such as elephant noses, octopuses, and worms, possess the advantages of adaptability, flexibility, and safety [5]- [8], thus meeting the growing demand for dexterous and human-friendly manipulation. In addition to their increased use in the academic community, soft robots are widely used in industrial operations [9], medicine [8], bionic robots [10], etc. The actuation strategies employed in soft robots mainly consist of pneumatic actuator-driven [11], [12], tendon-driven [13], shape memory alloy (SMA)-driven [14], and electroactive polymer (EAP)-driven mechanisms [15], which provide the merit of high dexterity and compliance [16] along with the extra drawbacks of low carrying capacity and poor anti-interference The associate editor coordinating the review of this manuscript and approving it for publication was Hamid Mohammad-Sedighi . ability. Robots with low stiffness exhibit difficulties in precise positioning. These issues have motivated a recent surge in the development of variable stiffness mechanisms for soft robots to achieve higher stiffness and ameliorate the fundamental trade-off between flexibility and stiffness.
In general, existing mechanisms designed to achieve variable stiffness can be roughly divided into 3 categories dependent on different principles: analytical-based mechanisms, material-based mechanisms, and structural-based mechanisms. In analytical-based mechanisms, a stiffness controller can be realized by sensing disturbing force, and understanding the relationship between the force and robot deformation [17]. However, there are notable challenges in deriving stiffness models for soft robots with high nonlinearity and researching feasible solutions for force perception. Many efforts have recently been devoted to material-based strategies to change the stiffness of continuum robots. Magnetorheological substances [18], electrorheological VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ fluids [19], and other intelligent materials have shown outstanding performance in terms of variable stiffness, but high magnetic fields or voltage fields are needed. Thermally activated materials, such as low melting point alloys [20], shape memory polymers (SMPs) [21], and shape memory alloys (SMAs) [22], can obtain a wide range of changes in stiffness under different temperatures. Nevertheless, most heat-activated materials require a relatively long activation time to transform from heated to cooled operational modes. Although some composites can significantly change their stiffness by modulating their elasticity modulus, simplifying the process and reducing costs remain challenging issues [23].
To avoid the limitations of material-based mechanisms, structural-based mechanisms have been favored by researchers. A popular example of a structural-based approach that allows for variable stiffness is the granular jamming mechanism [24], This method has been researched by modulating the relative movement between particles and ultimately meets the requirement of a large variable stiffness range. Unfortunately, the granular jamming mechanism inevitably increases the mass of the robot. In addition, Yong-Jae Kim presented a novel layer jamming technology that exploits the friction between thin material layers controlled by a confining pressure [25] to vary stiffness. However, the friction-based mechanism depends heavily upon the selection of an elastomeric membrane. It is also possible to use a variable neutral line mechanism to achieve adjustable stiffness. An asymmetric arrangement of tendons and links enabled continuous stiffness modulation in reference [26]. In reference [27], [28], the authors used extending and contracting hydraulic actuators imitating PAMs to attain variable stiffness operations, where one contracting actuator is surrounded by five extending actuators. Such approach is categorized as a technology that use active actuators arranged antagonistically to achieve variable stiffness for soft robots [29]. However, this structure not only increases the mass of the manipulator because of the high bulk modulus fluid but also limits its load capacity because the output force of the extending actuator is smaller than that of the contracting actuator.
In order to avoid excessive manipulator mass and increase the capacity of the loading operations, a soft manipulator actuated by both contractile and extensor PAMs has been proposed in reference [30]. The PAM can be categorized as a soft actuator consisting of a rubber bladder surrounded by braided shells, and it can serve as a contractile actuator when the braided angle is less than 54 • 44'; otherwise, it serves as an extensor actuator [31]. A soft manipulator consisting of PAMs responding to inflation pressure has the advantages of high compliance, low production cost, and a high powerto-weight ratio. In addition, unlike the OctArm [32], which consists only of extensor PAMs, the novel manipulator can achieve an output stiffness that can be varied independently of the position. The same strategy in reference [30] for variable stiffness was also applied in a soft gripper [33]. To reduce the size of the fingers, the gripper design used different types of PAMs, but it relocated the contractor muscles and transmitted their force to the fingers through tendons. However, references [30], [33] both investigated variable stiffness characteristics by some experimental measurements, and their kinematics were analyzed based on the hypothesis of constant curvature under no loading conditions. Inspired by the work in reference [34], this paper analyzes the manipulator stiffness using the Cosserat theory. This theory can derive the relationship among the manipulator motion shapes, structural properties, and loading capacities based on the basic principle of continuum medium mechanics. For modeling the continuum manipulator, the Cosserat theory is more accurate than techniques based on constant curvature assumptions [35]- [39] or pseudo-rigid body models [13] and is more efficient than the finite element approach [40]. Different from the work of Haibin et al. [34], who modeled the grasping force of an SMA-driven soft manipulator with a hypothesis of neglecting elongation strain and tangential strain, we break the above limitations and predict the manipulator motion profile under different pressure and external loading conditions by comprehensively considering the geometric and material nonlinearities of the PAM, the unique structural properties and the gravity of the manipulator. Then, the manipulator stiffness is derived by calculating the ratio coefficients between end-effector movements and external point loads.
The contributions of this paper hence lie on the structural strategy analysis for the variable stiffness of the pneumatic manipulator. Stiffness characteristic analyses based on manipulator shape prediction are conducted. By means of mathematical analysis and experimental validation, the soft manipulator is proven to possess the advantages of high loading capacity and variable stiffness independent of the position.
The rest of this paper is organized as follows: In Section II, we summarize the structure strategy of the soft manipulator. The procedure of theoretical analysis based on the modified output force model of the PAM and the Cosserat theory is conducted in Section III. Section IV presents the manipulator performance under different loading conditions, and further its stiffness characteristics are analyzed. Finally, the conclusions are outlined in Section V.

II. DESIGN DESCRIPTIONS OF THE MANIPULATOR A. THE STRUCTURE DESIGN
The pneumatic manipulator with variable stiffness is exhibited in Fig. 1(a). The connected extensor and contractile PAMs constitute the main structure of the manipulator and offer high flexibility and compliance. For implementation simplicity and cost savings, the PAM used in the manipulator is made in house and is mainly consists of PET braided shells and an elastomeric rubber bladder inside the braided shells. In the structural design of the manipulator, one central extensor PAM is evenly surrounded by three contractile PAMs. All PAMs are mounted to two mounting plates with a diameter  of 150 mm at both ends of the manipulator, and each contractile PAM is 35 mm away from the center of the manipulator. In particular, to ensure that the contractile PAMs are always in contact with the extensor PAM, nylon ties are utilized to pass through two adjacent crossing points in the braided shells and are located approximately every 25 mm along the length of the manipulator. The circles in Fig. 2 represent the location of the nylon ties. In addition, similar ties are located along the length of the external side of each contractile PAM and serve as cable guides for the displacement sensor wires. The displacement sensors are mounted to measure the contraction amount of the contractile PAMs. To eliminate the extra preload applied to the manipulator by the displacement sensor wire, the wire is wound around variable diameter pulleys, as displayed in Fig.1

B. WORKING PRINCIPLE
Due to the symmetrical distribution of the contractile PAMs relative to the center extensor PAM, the manipulator length along the neutral axis (always coincident with the extensor muscle center) equals an average of the length of two kinds of muscles and varies with each muscle motion. In addition, the manipulator movement can be decomposed into bending and extending motions since each PAM can be actuated independently.
It is worth noting that the output force of a contractile PAM is much higher than that produced by an extensor PAM with the same geometric size. Hence, the contractile PAMs located outside the manipulator ensure that the manipulator possesses a higher bearing capacity. The desired maximization of the output force and required payload is the decisive factor for the number and placement of contractile PAMs. At the same time, the center extensor PAM producing a larger deformation in length allows the manipulator to reach higher curvatures and a larger workspace, compared to a manipulator made purely of contractile PAMs. In addition, the combination of two kinds of PAMs makes it possible for the manipulator to achieve stiffness variation decoupled with its end-effector position variation. When the manipulator reaches a certain position, the stiffness can be tuned by simultaneously adjusting the inflation pressure of the extensor PAM and contractile PAMs. For example, in Fig. 3(a), the manipulator is inflated, while the extensor and contractile PAMs are all inflated in Fig. 3(b). Because the actuating forces generated by contractile and extensor PAMs are in opposite directions, the pressure in the extensor PAM in Fig. 3(b) must be higher than that in Fig. 3(a) to keep the initial lengths the same before applying loads. When the manipulators are deflected by the same payload, the manipulator in the higher pressure mode moves a shorter distance L 2 than that L 1 in the lower pressure operating mode, meaning that the manipulator stiffness becomes higher as the total pressure in the structure is increased.

III. MODEL ANALYSIS
Unlike traditional rigid-body robots, soft manipulators are subjected to a wide range of continuous deformations such as bending, torqueing and stretching. First, the output force VOLUME 8, 2020 of the PAM is modeled, comprehensively considering the structural properties and nonlinear elasticity of the material. Then the Cosserat theory is utilized to establish a static model of the manipulator, which provides a theoretical basis for predicting the performance of the soft manipulator in the following section.

A. OUTPUT FORCE MODEL OF THE PAM
An accurate output force model of the PAM can affect the performance prediction of the pneumatic manipulator. According to previous research reported in the literatures [41]- [45], the material compliance and unique structural characteristics of the PAM significantly complicate the modeling process. However, there are few studies on the modeling of the output force model of the PAM that consider both the nonlinear material elasticity [46] and irregular cylindrical shape [43] of the actuator. Thus, we are motivated to comprehensively incorporate the nonlinear PAM elasticity and structural characteristics into an ideal model based on the virtual work theory [42]. Fig. 4 shows the geometrical structure of the PAM. The basic structural parameters of the PAM include the current braid angle θ (from the axial plane of the PAM to the braid fiber), the current length L, and the current outer diameter D. In addition, the variable l in Fig. 4(b) represents the fiber length, and n is the number of turns that the fiber makes around the rubber tube. Since the braided shells are much stiffer than the rubber tube, it can be assumed that the length of the braided shell remains constant during the movement of the PAM. Thus, according to Fig. 4(b), the braided shell parameters and contraction ratio can be derived via geometric constraints: where ξ defines the contraction ratio of the PAM, and θ 0 , D 0 , and L 0 represent the initial angle, initial braid diameter, and initial length, respectively. The output force f ideal of a contractile PAM calculated in reference [42] is displayed as follows: where p represents the pressure in the PAM. When calculating the extension force of the extensor PAM, the force direction is opposite to that of the contractile PAM, and ξ represents the extension ratio. Fig. 5 shows the experimental and model results of Eq. (4) of force/pressure curves for a contractile PAM (the initial length: 600 mm, the initial diameter: 30 mm, and the initial braid angle: 35 • ). It can be seen that the PAM behaves with a slight hysteresis at certain contraction lengths in the isometric experiments, and there is an obvious discrepancy between the ideal model prediction measurements. In this paper, the ideal model is modified from two aspects. A loss of elastic energy is first considered using the Mooney-Rivlin theory, which can describe the nonlinear elasticity behavior of almost all rubber materials. As defined in reference [47], the Mooney-Rivlin strain energy function is given by Eq. (5), where C 10 and C 20 are Rivlin coefficients [48], and I 1 , I 2 and I 3 are Cauchy-Green strain tensors, which are expressed of the three principal stretch ratios λ 1 , λ 2 , and λ 3 , respectively, where, λ 1 = L L 0 = λ describes the axial stretch ratio along the longitudinal axis of the PAM (λ < 1 means the contraction state, and λ > 1 means the extension state), λ 2 represents the circumferential deformation, and λ 3 is the radial deformation. By applying the bladder incompressibility assumption, I 3 = 1 in Eq. (6). Thus, the elastic force produced by the bladder is expressed as follows: L 0 calculates the bladder volume by using the geometry size displayed in Fig. 4. Limited by existing experimental conditions, we set C 10 = 610000 and C 20 = −22000 according to reference [49] in the process of calculation. Using conservation of energy, the modified output force F of the PAM is shown in Eq. (8) by considering strain energy in the bladder, Second, the influence of the noncylindrical shape of the PAM on the modeling accuracy is taken into account. When the PAM works, a conical shape appears near the end fittings, and the proportion of the conical part is relevant to the inflation pressure and contraction ratio. Consequently, a nonlinear polynomial depending on both the air pressure and the actuator length is adopted to supply compensation for the ideal model, as shown in Eq. (9) where, [a 1 , a 2 , a 3 , b 1 , b 2 , b 3 ] are parameters needed to be identified by experimental datasets in Fig. 5. Additional attention should be paid to a special phenomenon in Fig. 5 before identifying Eq. (9). Due to the rubber tube elasticity and the space left between the braided sleeve and the inner tube in the initial state, there is an active actuating pressure p a for the actuator. When the inflation pressure is less than p a , the output force of the PAM equals zero. Therefore, the final improved output force model of the PAM is modeled as Eq. (10) pressure p a is a function of the contraction ratio, where k 1 and k 2 are identified using experiment data of different contraction ratios and active actuating pressure in  Fig. 6 shows the output force of the PAM for the theoretical and experimental data. The results of the ideal model serve as control groups. It is obvious that the results of the modified model agree well with the experimental data. To further validate the modified model in this paper, the output force of the PAM with a contraction ratio ξ = 0.2 is calculated in Fig. 7(a). The results show satisfactory precision, and its maximum absolute error ratio is 2.5%. Fig. 7(b) suggests the output force plot of another PAM  with an initial length of 580 mm when its contraction ratio is ξ = 0.1 and the maximum absolute error ratio is 2.7%. Thus, the improved model is used in the following section to consider the structural and material nonlinearities of the PAM in modeling the manipulator performance.

B. STATICS FORMULATION
To model the relationship between the load exerted at the end effector and the position, the Cosserat theory is employed to predict the continuous deformation of the manipulator in this paper under different loading conditions. The Cosserat theory has recently shown influential prospects in analyzing continuum robots. However, most researchers provided kinematic and static models for tendon-driven [50] or SMA-driven [34] manipulators and usually neglected the axial strain of the manipulator. In reference [51], although a pneumatic continuum robot was studied based on the Cosserat theory, the geometric characteristics of the actuators were not considered in the model, resulting in larger errors. Moreover, the structure in reference [51] consisted of extensor PAMs without the ability to change stiffness when decoupled of position, and there was a lack of systematic stiffness analysis. In our research, the nonlinear geometric property of the PAM, the manipulator gravity, and external loads are taken into account to calculate extensions, large curvatures, and shear deformations in the manipulator. Fig. 8(a) shows the deformation of the manipulator backbone under the actuating force, external load and gravity. In Fig. 8(a), o-xyz is defined as a global frame, which is stationary relative to the manipulator base, and O − d 1 d 2 d 3 is a local frame located at any point of the manipulator. Fig. 8(b) displays a cross-sectional diagram of the manipulator. The red circles represent contractile PAM I, PAM II, and PAM III, and the blue circle represents the extensor PAM. The output force produced by contractile PAMs are represented by F ci (i = 1, 2, 3), and F e is used to characterize the force produced by the extensor PAM. The force values can be calculated by Eq. (10). This section uses the conventions as followed in [50]: The spatial derivative of position vector r(s) can be described as Thus, the kinematic function of the manipulator in the global frame are defined specifically in the following way: where ' denotes a derivative with respect to s, and^denotes that the u(s) is converted into a skew-symmetric matrix. As shown in Fig. 8(a), the distributed force is defined as f (s), which is typically equal to the gravity. The load at the end of the manipulator is defined as W . The force F p and moment M are produced by output forces of the PAMs. We denote n(s) = n 1 (s)d 1 + n 2 (s)d 2 + n 3 (s)d 3 and m(s) = m 1 (s)d 1 + m 2 (s)d 2 + m 3 (s)d 3 as internal force and moment vectors. Therefore, the static equilibrium equations in the space can be obtained as follows: In this paper, extension, large curvature, and shear deformations are all considered. In order to simplify the model and improve the calculation efficiency, we assume that the internal forces n(s) and internal moments m(s) and the kinematic vectors υ(s) and u(s) follow the constitutive equations [52]: where K 1 = K 2 = GA, K 3 = EA, S 1 = EI * 1 , S 2 = EI * 2 , and S 3 = GJ . E is the Young's modulus, G is the shear modulus, and A is the area of the cross section, which can be calculated according to the geometric parameters in TABLE 1. The second moments of the I * 1 and I * 2 , and the polar moment of inertia of the cross section J are along d 1 , d 2 , and d 3 , respectively: Eq. (15) is decomposed with respect to d 1 , d 2 , and d 3 to obtain the statics governing the manipulator deformation, as shown in Eq. (18), where, The force F p is produced by the PAMs as shown in Fig. 8(b), which can be calculated according to Eq. (10). Eq. (14) and Eq. (18) form a set of ordinary differential equations, which need boundary conditions for solution. At the base of the manipulator, the kinematic boundary conditions are as follows: At the end of the manipulator, its rotation and displacement are related to the moment and total forces. Due to the forces produced by PAMs, a moment M applied to the end of the manipulator can be calculated as follows: where R c1 = R c2 = R c3 represent radiuses of the contractile PAMs, and R E represents radius of the extensor PAM, which can be found in TABLE 1.
The concentrated forces at the end of the manipulator are consisted by F p and W , thus, In this paper, the nonlinear problem is solved via the bvp4c provided in MATLAB. Algorithm 1 shows the specific algorithm for solving the model in Eq. (14) and Eq. (18)

C. EXPERIMENTAL VALIDATION
The above model is able to predict the bending shapes and thus the position of the manipulator end when the inflation pressure and the external forces are given. To verify the modeling method, a series of experiments are conducted with different inflation methods and external loads.
Figs. 9(a) and 9(b) present a schematic diagram of the shape configuration tests and an experimental setup. The experimental setup includes the manipulator system described in Fig. 1(a), a PC, an air compressor, an air source triplet, four air pressure sensors, a 9-axis sensor, an Arduino Mega 2560, a sheet of graph paper, a laser pointer,  Under the influence of external loads and actuating pressure, the soft manipulator bends into different shapes. The shape configurations of the manipulator are first tested with zero external load. During the inflation process, the extensor PAM is first inflated to a desired value, and then three contractile PAMs are inflated to 15 kPa. Finally the pressure of PAM I is subsequently increased the desired value. Thus the manipulator bends with different shapes, mainly by depending on the pressure in the extensor PAM and PAM I. To describe the influence of payloads, other experiments are performed with the external weights attached to the end plate of the manipulator. With the test device shown in Fig. 9, the proposed model above can be validated with different combinations of the actuating pressure and external loads. Due to the same pressure in two of the contractile PAMs, the manipulator bending with the variables in Eq. (18) υ 2 = 0, u 1 = 0, and u 3 = 0. TABLE 3 shows part of the experimental configurations for model validation. PC I, PC II and PC III represent the pressure values in contractile PAM I, PAM II and PAM III, respectively. PE represents the pressure value in the extensor PAM, and Load represents the external load attached to the end of the manipulator. Fig. 10 illustrates the experimental and simulation results. In order to emphasize the effects of structural and material properties of the PAM along with the distributed weights and external loads on the performance estimation of the manipulator, previous research [53] on the constant curvature model is used as control results. In Fig. 10, the red points represent the end positions of the manipulator measured by experiments, and the green solid lines are used to illustrate the bending direction in experiments; thus, their angles away from the vertical direction represent the bending angles of the In order to quantitatively evaluate the difference between the experimental results and the simulation results, the absolute error (e p ) along the distance direction between the experimental and simulation results is defined by where (x p , z p ) represents the end position of the manipulator calculated by the proposed model, and (x e , z e ) is the corresponding experimental results. e px = x p − x e and e pz = z p − z e represent the absolute errors of the manipulator's position in the x and z directions, respectively. In addition, we define Eq. (24) to calculate the absolute error of bending angle, where θ p and θ e represent the theoretical and experimental bending angles of the manipulator end, respectively. TABLE 3 shows the above errors between the simulation and experimental results. The average e px , e pz and e p are 19.2 mm, 19.1mm and 27.5 mm, respectively. Compared with the total length of the manipulator (600 mm), this error accounts for a small ratio, i.e., 3.2%, 3.2%, and 4.6% of the total manipulator length, respectively. The average e θ is 2.7 • , which is 3% of the maximum bending angle (90 • ).
These results all indicate that the simulation results of the model in this paper are consistent with the experimental values. By contrast, the results calculated by the constant curvature model are also compared with the experimental results, where e p between the experiment and model results is 123.8 mm, which accounts for 20.6% of the total manipulator length. It is worth noting that the manipulator undergoes bending and extension as a result of the inflating pressure and external loads due to its compliant nature. Therefore, the influences of structural and material nonlinearities, gravity, and external loads are important for predicting the bending shape of the soft manipulator more accurately.

IV. STIFFNESS CHARACTERISTICS
The stiffness characteristic reflects the relationship between the deformation and the external force applied to the soft manipulator. Consequently, to characterize the manipulator stiffness, we first analyze the shape of the manipulator under loading and unloading conditions. In this section, the influences of the inflation methods and loads on the bending profiles of the manipulator are first predicted based on the verified model in Section III. Then, the stiffness of the manipulator is characterized quantitatively based on the shape prediction under different inflation cases, thus providing a reference for the manipulator application.

A. BENDING SHAPE PREDICTION
According to the experiments in Fig. 9, it is obvious that pressure in the actuators and the weights of the payload affect the manipulator shape due to its compliant nature. By using the model in Section III, the bending shapes of the manipulator are predicted under the following conditions. In Fig. 11, the pressure in the extensor PAM is 10 kPa, 20 kPa, and 40 kPa from left to right, and the load applied to the manipulator end is 0 kg, 0.5 kg, and 1.0 kg from top to bottom. For each combination of pressure values in the extensor PAM and the load, the pressure in PAM I varies from 20 to 70 kPa with 10 kPa increment, while the other two contractile PAMs maintain a constant pressure of 15 kPa. The bending angle of the manipulator end increases with increasing pressure in PAM I and reaches a maximum with a pressure of 70 kPa in PAM I. We find that depending upon the inflation pressure and the attached load, the bending angle of the manipulator end changes, and its end moves to a range of different locations. The bending angle trends against the air pressure in PAM I under different loading conditions are illustrated in Fig. 12.
To fully illustrate the effects of inflation pressure and payload variations on the bending behavior of the soft manipulator, TABLE 4 summarizes the analysis according to Fig. 12. In the first column, the effects of the pressure change in the extensor PAM are observed. These results are collected by comparing the maximum bending angle corresponding to 10, 20, and 40 kPa of pressure in the extensor PAM, which confirms that the higher pressure in the extensor PAM achieves a higher bending angle of the manipulator end. The result shows a 22.4% increase on average among the variations in the first column. In the remaining columns, the effects of the pressure change in the contractile PAM are the focus. In the second column, the pressure in the extensor PAM equals 10 kPa, and the bending angle variations are collected with the pressure in the contractile PAM I varying from 20 to 70 kPa. By analogy, the data in the third and fourth columns are obtained with the pressure in the extensor PAM set to 20 kPa and 40 kPa, respectively. The average percentage increases of bending angles are 145%, 155%, and 134%. They are all more than 100%, meaning that the novel   manipulator has the ability to attain a large increase in curvature. To further present the load capacity of the manipulator, Fig. 13 shows the comparative bar charts of the percentage decrease of the bending angle when the load varies from 0.5 kg to 1 kg. The figure clearly shows that the increase in pressure PE results in less reduction of the bending angle. By collecting the bending angle reductions corresponding to 20 kPa, 30 kPa, 40 kPa, 50 kPa, 60 kPa, and 70 kPa of pressure in PAM I, the average percentage decrease of the bending angle decreases from 15.9% to 10% when the pressure in the extensor PAM increases from 10 kPa to 40 kPa. These data indicate that when the load is doubled, the bending angle of the manipulator reduces by no more than 16%, which proves that the novel structure attains large output force to lift the medium load without a considerable impact on its range of motion.

B. STIFFNESS ACCESSMENT
It is known that the novel manipulator stiffness is inversely proportional to its displacement under a certain load condition. In other words, if the bending shape of the manipulator without load is defined as A, on the basis of shape A, VOLUME 8, 2020 the bending shape of the manipulator applied a load to the end is defined as B. Therefore, the proportion between the load and movement between the ends of shape A and shape B can be used to calculate the stiffness of the manipulator. Based on the shape prediction of the novel manipulator shown in Fig. 11, the positions of the manipulator end with zero load and a load of 1 kg applied to the end are summarized in Fig. 14. Fig. 14(a) shows the end positions in the x direction under different inflation configurations, and Fig. 14(b) displays the end positions in the z direction. The arrows indicate the displacements of the manipulator end under loading conditions. The stiffness results against the pressure in contractile PAM I are shown in Fig. 15 such that the stiffness characteristics of the manipulator are obtained. In order to verify the validity of the calculation, the experimental data are obtained by the same experimental setup in Fig. 9. After each PAM is inflated to the desired value, the positions of the manipulator end are marked before and after the weight (1 kg) is vertiaclly arranged at the center of the end disk. Based on performing load-unload tests five times in each case, the average stiffness can be gained under different inflation conditions. In Fig. 15, diamond symbols represent the experimental data, and solid lines represent the theoretical data. Moreover, the maximum model errors in the x direction and y direction are estimated as 3.26% and 4.3%, respectively.
As shown in Fig. 15, the stiffness in the z direction is higher than that in the x direction in the same case. The stiffness of the novel manipulator increases with increasing pressure in the middle extensor PAM, and a higher stiffness increase arises at a lower pressure of the contractile PAM. TABLE 5 reports the effects of pressure changes on the stiffness of the manipulator quantitatively. The first column  shows the percentage increase of stiffness when the pressure of the extensor PAM increases from 10 kPa to 20 kPa. The stiffness increases in the x direction range between 28% and 54% for a range of stiffness related to data collected for 10, 20, 30, 40, 50, 60 and 70 kPa of pressure in the contractile PAM; analogously, the stiffness increases in the z direction range between 5% and 40%. In the second column, the stiffness increases between 10% and 46% in the x direction and between 2% and 11% in the z direction with the change of extensor pressure from 20 to 40 kPa.
The above analyses illustrate the variation trends of the manipulator stiffness response to the inflation pressure. The stiffness variation is coupled with the end position of the manipulator. From Fig. 15, we can also see that the stiffness decreases with the increase of the bending angle of the manipulator, which causes a poor anti-interference ability. This problem can be solved by the variable stiffness mechanism decoupled from its position. Fig. 16 shows the movement of the manipulator applying the same load with two kinds of stiffness states. When there is no load applied to the manipulator, the blue solid line calculates the manipulator shape with , and x = 371.1mm and z = 427.8mm, respectively. After a load of 1 kg is applied to the manipulator, the blue solid line moves to the blue dotted line, while the red solid line moves to the red dotted line. Comparing the manipulator with a totally lower pressure (represented by the solid blue line), it can be calculated that the displacement of the manipulator represented by the red line decreases by 9% when moving in the x direction and 14.8% when moving in the z direction, which means that a higher stiffness is obtained when the total pressure is increased. In fact, there are infinite combinations that can actuate the manipulator to the same target location. Fig.17 exhibits the variable stiffness of the manipulator at the same unloaded position (374.2 mm, 427.5 mm) shown in Fig. 16 theoretically and experimentally when the inflation pressure of the extensor PAM increases from 0 kPa to 40 kPa. It can be clearly seen that the stiffness of the manipulator can be higher with a wholly higher pressure, meaning that the stiffness performance of the manipulator can be varied independent from its position. The results indicate the efficacy of the variable stiffness mechanism in this paper.

V. CONCLUSION AND FUTURE WORK
Due to the combination of contractile and extensor PAMs, this paper has introduced a soft manipulator possessing high compliance and a variable stiffness mechanism. Compliance for the manipulator is preferable for tasks that require adaptability and flexibility to the environment, and the variable stiffness mechanism provides higher stiffness without position changes when strong anti-interference ability is required. Then, the Cosserat theory, coupled with both the material nonlinear elasticity and shape irregularity of the PAM, is utilized to comprehensively characterize the performance of the novel manipulator for the first time. Based on the shape prediction under different loading conditions, the influence of inflation pressure on the stiffness of the manipulator is analyzed. The theoretical and experimental results reveal that 1) the bending space of the manipulator increases with the inflation pressure up in the extensor PAM, which allows the manipulator to attain a lager bending angle than a manipulator constructed only by contractile PAMs; 2) the contractile PAMs evenly distributed around the extensor PAM allow the manipulator to achieve a high load capacity; 3) the manipulator stiffness increases with the increase of inflation pressure in the extensor PAM; and 4) the variable stiffness of the manipulator is realized by adjusting the pressure in both the contractile and extensor PAMs simultaneously while leaving the shape invariant.
The stiffness analysis in this paper provides a guideline for soft manipulator design, motion and control. Improving stiffness performance based on optimizing the structure will be the focus of future work. Moreover, ongoing work will aim to further replicate one link to allow multiple modules. Each module will be able to provide the similar performances in terms of motion capability and stiffening capability.