Modeling and Compensation of Stiffness-Dependent Hysteresis Coupling Behavior for Parallel Pneumatic Artificial Muscle-Driven Soft Manipulator

Abstract

The parallel driving soft manipulator with multiple extensors and contractile pneumatic artificial muscles (PAMs) is able to operate continuously and has varying stiffness, achieving smooth movements and a fundamental trade-off between flexibility and stiffness. Owing to the hysteresis of PAMs and actuator couplings, the manipulator outputs display coupled hysteresis behaviors with stiffness dependence, causing significant positioning errors. For precise positioning control, this paper takes the lead in proposing a comprehensive model aimed at accurately predicting the coupled hysteresis behavior with the stiffness dependence of the soft manipulator. The model consists of an inherent hysteresis submodule, an actuator coupling submodule, and a stiffness-dependent submodule in series. The asymmetrical hysteresis nonlinearity of the PAM is established by the generalized Prandtl–Ishlinskii model in the inherent hysteresis submodule. The serial actuator coupling submodule is dedicated to modeling the actuator couplings, and the stiffness-dependent submodule is implemented with a fuzzy neural network to characterize the stiffness dependence and other system nonlinearities. In addition, an inverse compensator on the basis of the proposed model is conducted. Experiments demonstrate that this model possesses high accuracy and good generalization, and its compensator is effective in decoupling and mitigating hysteresis coupling of the manipulator. The proposed model and control methods significantly improve the positioning accuracy of the pneumatic soft manipulator.

1. Introduction

Pneumatic artificial muscles (PAMs) are mainly made up of braided woven meshes, rubber tubes, and sealed connectors [1]. As a kind of soft actuator, the PAM features a high power density ratio and a large driving force, breaking through the flexibility limitations of the rigid actuator to widely serve as a bionic actuator in robotic fields for underwater tasks [2], elderly care services [3,4], rescues [5], rehabilitation [6,7,8], and other flexible manipulations [9]. By comparison with rigid manipulators, PAM-driven soft manipulators have a stronger ability to adapt to unstructured environments and produce safer interactions with human beings. Whether contractile PAMs (initial braided angle is less than 54°44′) or extensor PAMs (initial braided angle is more than 54°44′) are used, they can only provide linear motions. PAM-driven soft manipulators usually move in 3D spaces through combined deformations of multiple parallel PAMs. For example, as shown in Figure 1, Walker et al. [10] presented an OCTARM continuum manipulator design (Figure 1a) with three sections, with each section containing different numbers of parallel extensor PAMs, achieving nine total degrees of freedom. Similarly, Falkenhahn et al. proposed a BIONIC HANDLING ASSISTANT (Figure 1b) [11], which has three connections in series, and a total of nine bellow-based pneumatic actuators are arranged and connected in parallel. Reference [12] designed a cylindrical manipulator (Figure 1c) with each section having three radially symmetric pneumatic chambers. In addition, a dual-plane constraint design (Figure 1d) for the backbone structure was proposed in reference [13] to integrate the parallel cables and pneumatic artificial muscles into the continuum robot. A soft manipulator controlled by antagonistically arranged extensor and contractile PAMs was proposed to enhance their anti-interference capability and load-bearing capacity in reference [14] (Figure 1e), and our previous work also designed a soft manipulator driven by three contractile PAMs and one extensor PAM in parallel [15,16], as shown in Figure 1f. Based on opposing configurations of extensor PAMs and contractile PAMs, the soft manipulator can output different stiffnesses varying independently from its position.

Despite the fact that PAM possesses excellent application prospects in soft robots, because of friction loss, rubber tube elasticity, and the unique geometric structure, it presents a severe asymmetric hysteresis phenomenon, preventing PAM-driven soft manipulators from positioning accurately [17,18,19]. Further complicating matters, parallel PAMs contained in the soft manipulator are constrained to interact with each other; thus, the system is negatively affected by the coupled hysteresis among different PAMs, increasing difficulties in precise modeling and control of the PAM-driven soft manipulators to a higher degree. It is noted that the hysteresis coupling effect almost exists in any multiple PAMs-actuated parallel system. In order to compensate and decouple the hysteresis coupling effect, an effective modeling and control strategy should be investigated. Currently, hysteresis modeling and compensation control for smart materials have significantly attracted the attention of researchers [20,21,22]. Now, available hysteresis models can be broadly classified into physics-based models and phenomenology-based models. Physics-based models, such as the Maxwell-slip model, are derived from the working principles of systems and material properties [23], meeting difficulties in establishing hysteresis models with high precision for some simplifications. Therefore, phenomenology-based hysteresis models appear, and plenty of work has been conducted. The models based on differential equations and the models based on operators are dominant phenomenology-based hysteresis models. For example, the Bouc–Wen model, as a representation of the differential-based models, adopts nonlinear differential equations to calculate the PAM hysteresis [24]. Krasnosel’skii–Pokrovskii (KP) models [25,26], modified Prandtl–Ishlinskii (MPI) models [27], Preisach models [28,29], generalized Prandtl–Ishlinskii (GPI) [30,31] models, and some other operator-based models fit hysteresis by weighting and summing multiple hysteresis operators. The above existing hysteresis models are mostly concerned with the inherent hysteresis (compared with coupled hysteresis effects) in PAMs, and different PAM couplings are rarely concerned in hysteresis models. For compensation control of parallel PAMs-based manipulators, mechanical couplings are usually regarded as unmodeled disturbances, and they are solved by introducing feedback controllers parallel with individual inverse hysteresis compensators, complicating the controller design. As an example, reference [32] employed the Preisach model to characterize the hysteresis behavior of PAMs and integrated a feedback controller with the inverse Preisach compensator to address the interference among different PAMs. Although reference [33] proposed a lumped model considering mechanical coupling and simple friction with constant curvature assumption for a pneumatic manipulator, it did not accurately describe the hysteresis behavior of the soft manipulator. Particularly, the model calculation is too complicated to achieve control; hence, it needs to be linearized, and a closed-loop controller is employed simultaneously to enhance control accuracy. To avoid time-consuming calculation and model simplification, reference [34] derived a recurrent neural network-based model for a pneumatic manipulator, focusing on learning the coupled relationship between actuator length variations and the manipulator’s posture. However, this model lacks sufficient calculations for hysteresis effects between pressure and the actuator length responses, and is not suitable for low-level controller design.

Based on the above discussions, the research work on modeling of the coupled hysteresis effect in the multiple PAMs-based manipulator is still open, and successful modeling methods can greatly promote the further development of pneumatic manipulator control. This paper conducts research on modeling and compensation of the coupled hysteresis behavior of the soft manipulator with variable stiffness. The main contributions of this paper are illustrated as follows:

(1)

This paper firstly presents a series of experiments to characterize the coupled hysteresis behavior of the soft manipulator with variable stiffness based on PAMs.

(2)

According to existing literature reviews, now, available hysteresis models mostly concentrate on modeling individual hysteresis effects (compared with coupled hysteresis effects) in the PAM, and few hysteresis models have been established involving hysteresis couplings of the PAM-based manipulator, especially with variable stiffness. Based on the input–output behavior characteristics, a coupled hysteresis model that comprehensively considers stiffness influence is proposed. The newly proposed model is constructed of three serial parts, including an inherent hysteresis part, an actuator coupling part, and a stiffness-dependent part. To our knowledge, the model in this paper innovatively predicts the hysteresis coupling of the parallel PAMs-based robot with variable stiffness, exhibiting outstanding generalization ability and precision. We aim to demonstrate that the proposed model with high precision on pneumatic manipulators takes actuator coupling, inherent hysteresis, variable stiffness mechanism, and other nonlinearities in complicated systems into account as much as possible. Experimental validations show that the proposed model can predict the hysteresis coupling of the manipulator under various operation conditions.

(3)

The model is readily applicable to compensation for parallel PAMs-driven robots. Through directly learning the compensator, decoupling and compensation for the manipulator have been successfully achieved, prominently enhancing its positioning precision.

The rest of the paper is as follows. Section 2 tests and analyzes quantitatively coupled hysteresis responses of a PAM-driven manipulator with variable stiffness. Section 3 establishes a coupled hysteresis model. Section 4 derives a decoupled inverse compensator using the model proposed in Section 3. Additionally, Section 5 of this paper conducts identifications and validations of the proposed model and compensator, accompanied by relevant discussions. Finally, Section 6 of the paper summarizes the conclusions.

2. Coupled Hysteresis Characteristics of the Soft Manipulator

2.1. Experiment Setup

The experimental system for quantitatively investigating and analyzing the input–output responses of the pneumatic manipulator is completed in our laboratory, as shown in the Figure 2, which contains the following elements.

(1)

Soft manipulator: Figure 3a illustrates the structural design of the soft manipulator manufactured in the laboratory. Between the two parallel disks, one extensor PAM and three contractile PAMs are installed. To ensure continuous interaction between contractile and extensor PAMs during movements, the braided woven meshes of the extensor PAM and the contractile PAM are tied up with nylon ties along the axis direction. Figure 3b shows the cross-direction view of the manipulator. It can be seen that contractile PAM I, PAM II, and PAM III are placed evenly around the extensor PAM, indicating an angle of $\angle O_{1}O{O}_2=\angle O_{1}O{O}_3=\angle O_{2}O{O}_3=120°$. To be customized easily and reduce cost, PAMs adopted in the manipulator are manufactured by hand in the laboratory. Geometrical parameters are shown in

Table 1. The structural parameters of the soft manipulator.

	Initial Length	Initial Braid Angel	Initial Diameter	Withstand Voltage Range
Extensor PAM	600 mm	60°	45 mm	2.5 bar
Contractile PAM	600 mm	35°	30 mm	2.5 bar

. Due to parallel and symmetrical structure design, after inflation, the length variation of the extensor PAM can be calculated by the average length variations of different contractile PAMs; thus, the end position of the manipulator is determined by the outputs of contractile PAMs. The extensor PAM has an effect on regulating the stiffness of the manipulator. The manipulator with higher pressure in the extensor PAM has higher stiffness than that of the manipulator with lower pressure under the same pressure in contractile PAMs [15]. In addition, based on the antagonistic driving mechanism of extensor and contractile PAMs, the soft manipulator can achieve stiffness varying independently from its position by simultaneously inflating or deflating two types of PAMs [15,16].

(2)

Driving system: To actuate the manipulator, the driving system consists of an air compressor, four proportional valves (SMC, ITV2050-212L), a pneumatic triplet, a host PC, and a dSPACE. A host PC with Simulink is used to communicate with the dSPACE. The analog signals designed in the Simulink are downloaded to the dSPACE and excite proportional valves; thus, the desired inflation pressure is applied to the soft manipulator.

(3)

Perception system: The experimental perception system mainly includes four pressure sensors and three displacement sensors. During experiments, in order to measure the length variations of contractile PAM I, PAM II, and PAM III, displacement sensor wires are threaded through cable guides located along every contractile PAM. Measurements from displacement sensors and pressure sensors are received by the dSPACE and saved on the PC for further investigation.

2.2. Hysteresis Coupling Effects and Stiffness-Dependence

In order to provide an understanding of the coupled hysteresis behaviors of the parallel PAMs-driven manipulator system, each of three contractile PAMs in the soft manipulator is inflated with $p*=0.65+0.65\sin \left(0.04\pi t\right)$ bar alternately, as well as the other two with zero pressure. Pressure p₀ in the extensor PAM keeps 1 bar. Experimental results on the pressure–length variations of contractile PAMs are shown in the first three lines of Figure 4. And the sampling frequency is set to 1000 Hz. Then, all contractile PAMs are inflated together. Figure 4j–l show the relationships between $p_1$ (the pressure of PAM I) and $y_1$ (the displacement of PAM I), $p_2$ (the pressure of PAM II) and $y_2$ (the displacement of PAM II), and $p_3$ (the pressure of PAM III) and $y_3$ (the displacement of PAM III), respectively. As introduced in the above Section 2.1, when the pressure within contractile PAMs remains unchanged, the manipulator stiffness increases as the pressure of the extensor PAM rises [15]. Therefore, this paper employs the pressure $p_0$ within the extensor PAM to represent the manipulator stiffness. To reveal the influence of manipulator stiffness on its output behavior, relationships between $p_1=0.35\sin (0.04\pi t)+0.35\mathrm{bar}$ and output variations of PAM I, PAM II, and PAM III corresponding to $p_0=0,0.5,1.0,1.5,2.0\mathrm{and}2.5$ bar and $p_2=p_3=0\mathrm{bar}$ are investigated, respectively. Figure 5 shows partial results with $p_0=0,0.5,1.0$ bar. In view of Figure 4 and Figure 5, the hysteresis and mechanical coupling effects in the manipulator system can be clearly seen. (1) Hysteresis effects: There are obvious hysteresis characteristics in PAMs, and the hysteresis is asymmetric and possesses a dead zone. The maximum width ratio of the hysteresis loop calculated in the form of $e_w=(\left|y_{up}-y_{down}\right|/\left|y_{\max }-y_{\min }\right|)\times 100\%$, ($y_{up}^i$ and $y_{down}^i$ are PAM output values corresponding to the middle input value, $y_{\max }^i$ and $y_{\min }^i$ define the maximal and minimal PAM output, respectively) indicates that $e_w$ is within the range of [28.5%, 50.1%]. The hysteresis will cause serious inaccuracy and even oscillation during the manipulator’s motion. (2) Actuator coupling effects: It can be seen that there are coupling behaviors among three contractile PAMs. Taking the first line of Figure 4 as an example, if no coupling exists, the displacements of PAM II ($y_2$) and PAM III ($y_3$) should be equal to zero, with $p_2=p_3=0$, while the outputs of PAM II and PAM III show obvious responses to the input excitation $p_1$. And their hysteresis curves exhibit different appearances although $p_2=p_3=0$. This is because there exist handmade errors in the structure parameters of pneumatic artificial muscles, causing more intricate nonlinearity with coupled hysteresis. (3) Variable stiffness effects: The comprehensive experimental results in Figure 5 suggest that the stiffness variation produces non-negligible effects on coupled hysteresis loops, especially on the outputs of PAMs with low pressure. It is very intuitive to see that the pressure $p_0$ has a nonlinear effect on the outputs of contractile PAMs, which means that the stiffness state acts on the coupled outputs of the soft manipulator system. The output amplitudes of contractile PAM II and PAM III sharply decrease with increasing stiffness, while the output amplitude of the contractile PAM I has a small change compared to PAM II and PAM III. The reason is that pressure in the PAM II and PAM III remains zero, which is much lower than that in the PAM I, and they have weaker anti-interference ability. To our best knowledge, there are few studies comprehensively considering the above-mentioned hysteresis effects and mechanical coupling, especially with variable stiffness effects, when modeling the PAM-based soft manipulator with high precision. A comprehensive coupled hysteresis model should be developed first for better positioning.

3. Modeling of Coupled Hysteresis

The experimental results in Section 2 show that not only do PAMs have inherent hysteresis properties, outputs of different PAMs distributed parallelly in the soft manipulator, but they also depend on each other and show strongly stiffness-dependent characteristics. To describe this complex behavior, a stiffness-dependent and coupled hysteresis model (SDCHM) consisting of three parts is established, as shown in Figure 6. Generalized Prandtl–Ishlinskii (GPI) models GPI₁, GPI₂, and GPI₃ are used to model the inherent hysteresis behavior of each contractile PAM, respectively. The GPI model features a simple structure and high flexibility, suitable for addressing PAM hysteresis with dead zones and asymmetry. Furthermore, the GPI model has an analytical inverse solution, providing convenience for the inverse compensator derivation. A linear coupling structure is cascaded in series in the actuator coupling submodule. $G_{ij}(i=1,2,3;j=1,2,3)$ are gains that describe the input–output coupling in the manipulator system. For circumventing simplifying assumptions and facilitating the completion of compensators, inspired by data-driven modeling approaches, the stiffness dependence of the output behavior of the soft manipulator and other unmodeled errors in the first two subparts are fitted in the stiffness-dependent submodule $g(\cdot )$ using a T-S fuzzy neural network (TSFNN) [35]. The inflation pressure of $p_1$, $p_2$, and $p_3$ in contractile PAMs firstly acts upon the inherent hysteresis submodule, and outputs $v_1$, $v_2$, and $v_3$ from the inherent hysteresis submodule are subsequently input into the actuator coupling submodule. Before the final calculated displacement of ${\widehat{y}}_k(k=1,2,3)$ is generated, the actuator coupling submodule outputs $h_k(k=1,2,3)$ and inflation pressure $p_0$ in the extensor PAM flow into the stiffness-dependent submodule for further calculation.

3.1. Inherent Hysteresis Submodel

Among many kinds of available inherent hysteresis models, the generalized P-I hysteresis (GPI) model [36] is applied for its simplicity, good generalization ability for symmetric and asymmetric hysteresis, and the availability of its analytical inverse model. It should be noted that many other asymmetric hysteresis models for the PAM can also be used, such as the EUPI model [37], MPI model, etc.

The GPI model is defined by using weighted generalized play operators [33] to establish asymmetric and saturated hysteresis. The generalized play operator can be derived by a discontinuous function switching between envelop functions of rising and falling edges. The expression of the generalized play operator is shown in Equation (1):

$$\{\begin{array}{c}w(0)=F_{r,\gamma }[p_k](0)=f_{r,\gamma }(p_k(0),0),\\ w(t)=F_{r,\gamma }[p_k](t)=f_{r,\gamma }(p_k(t),F_{r,\gamma }[p_k](t-1)).\end{array}$$

(1)

with

$$f_{r,\gamma }(p_k(t),w(t-1))=\{\begin{array}{l}\mathsf{\Gamma }\circ {\gamma }_R\circ p_k(t),\mathrm{if}{p}_k(t)>p_k(t-1),\\ \mathsf{\Gamma }\circ {\gamma }_L\circ p_k(t),\mathrm{if}{p}_k(t)<p_k(t-1),\\ w(t-1),\mathrm{if}{p}_k(t)=p_k(t-1).\end{array}$$

(2)

where $\mathsf{\Gamma }$ defines the classical play operator, and ‘$\circ$’ represents the composition of functions. $w(t)$ is the generalized play operator output, and $p_k(t),(k=1,2,3)$ is the input value. r denotes the dead zone threshold. ${\gamma }_L(\cdot )$ and ${\gamma }_R(\cdot )$ are envelop functions satisfying ${\gamma }_R(p_k(t))-r\le {\gamma }_L(p_k(t))+r$. Then, according to the play operator [38], Equation (2) can be expanded as

$$f_{r,\gamma }(p_k(t),w(t-1))=\{\begin{array}{l}\max ({\gamma }_R(p_k(t))-r,w(t-1)),\mathrm{if}{p}_k(t)>p_k(t-1),\\ \min ({\gamma }_L(p_k(t))+r,w(t-1)),\mathrm{if}{p}_k(t)<p_k(t-1),\\ w(t-1),\mathrm{if}{p}_k(t)=p_k(t-1).\end{array}$$

(3)

In order to describe and compensate for the PAM hysteresis with output saturation characteristics, hyperbolic tangent envelope functions are adopted in Equation (4). It is worth noting that envelop functions can be defined as other function forms, and we only chose one form to discuss in this paper.

$$\{\begin{array}{c}{\gamma }_R(p_k(t))=\tanh (a_1^k{p}_k(t)+a_2^k)+a_3^k,\\ {\gamma }_L(p_k(t))=\tanh (b_1^k{p}_k(t)+b_2^k)+b_3^k.\end{array}$$

(4)

Then, the output of the GPI model can be defined as:

$$GPI[p_k](t)={\displaystyle {\int }_0^R{x}^k(r^k)F_{r,\gamma }[p_k](t)}dr$$

(5)

where $x^k(r^k)$ is the weight coefficient of the $GPI[p_k](t)$. In practical application, a finite number of GPI operators are sufficient to characterize the hysteresis; thus, the GPI model in Equation (5) can be discretized as Equation (6):

$$GPI[p_k](t)={\displaystyle \sum _{i=0}^N{x}^k(r_i^k)F_{r_i,\gamma }[p_k](k)}$$

(6)

In Equation (6), N denotes the number of operators. $x^k\left(r_i^k\right)$ denotes the weight coefficient of the ith generalized play operator, and the value of $r_i^k$ is given as

$$r_i^k={\displaystyle \frac{i}{N}}{p}_{k\max }(i=0,1,2,\cdots ,N)$$

(7)

In Equation (7), $p_{k\max }$ represents the maximum input value. Parameters $x_i^k(r_i^k)$, $[a_1^k,a_2^k,a_3^k]$, and $[b_1^k,b_2^k,b_3^k]$ should be identified based on experimental data.

3.2. Actuator Coupling Submodule

Outputs $[h_1, h_2, h_3]$ of the actuator coupling submodule are regarded as the combination of individual outputs from different PAMs represented by $v_1$, $v_2$, and $v_3$. The mathematical formulation of the mechanical coupling submodule is as shown in Equation (8):

$$\{\begin{array}{c}{h}_1(t)=G_{11}{v}_1+G_{21}{v}_2+G_{31}{v}_3,\\ h_2(t)=G_{12}{v}_1+G_{22}{v}_2+G_{32}{v}_3,\\ h_3(t)=G_{13}{v}_1+G_{23}{v}_2+G_{33}{v}_3.\end{array}$$

(8)

where $v_1=GP{I}_1[p_1](t)$, $v_2=GPI[p_2](t)$, and $v_3=GPI[p_3](t)$. $G_{ij}(i=1,2,3;j=1,2,3)$ are coupling coefficients. The actuator coupling submodule is able to extend to a combination of n GPI models, contingent on the number of actuators coupled with each other.

3.3. Stiffness-Dependent Submodule

Neural networks have strong nonlinear mapping abilities to comprehensively and accurately describe highly nonlinear systems. In this paper, the stiffness-dependent submodule $g(\cdot )$ is fitted by the T-S fuzzy neural network (TSFNN) that has good generalization capacity for modeling the nonlinear systems. The structure of $g(\cdot )$ is shown in the following Figure 7.

This structure of $g(\cdot )$ in Figure 7 is referred to the TSFNN structure [39] based on if–then rules. The fuzzification process is achieved in the first layer. Each input element corresponds to three fuzzy variables, and membership relations are represented by the Gaussian membership function calculated in Equation (9):

$$\{\begin{array}{c}{L}_{ij}^1(t)={\mu }_{ij}(t)=\exp \text{{}-{(I_i(t)-m_{ij})}^2/2{\sigma }_{ij}^2\text{}},(i=1,2,\dots ,n;j=1,2,3),\\ I_i(t)\in \text{{}{h}_1(t),h_2(t),h_3(t),p_0(t)\text{}}.\end{array}$$

(9)

where input elements $I_1(t),I_2(t),\cdots ,I_n(t)$ are composed of $\text{{}{h}_1(t),h_2(t),h_3(t),p_0(t)\text{}}$, and there are n = 4 input nodes in total. Figure 7 shows that every input matches with three membership functions ${\mu }_{ij}(i=1,2,\dots ,n;j=1,2,3)$. $L_{ij}^1(t)={\mu }_{ij}(t)$ shows outputs of the first layer. ${\sigma }_{ij}$ and $m_{ij}$ are antecedent parameters representing the width and center of Gaussian function, respectively. There are $N_1=3n$ nodes in the first layer. The complexity of the stiffness-dependent submodule is mainly determined by the number of fuzzy variables. Fewer fuzzy variables will reduce model accuracy, but more fuzzy variables will reduce the real-time performance of the model, and more training samples will be required. The output $L_d^2(t)$ of the second layer is calculated as shown in Equation (10):

$$L_d^2(t)=\prod _{i=1}^n{\mu }_{ij}(t),(j\in \{1,2,3\};d=1,2,\cdots ,N_2)$$

(10)

Then, $L_d^2(t)$ is normalized in the third layer with the output $L_d^3(t)$:

$$L_d^3(t)={\overline{\omega }}_h=L_d^2(t)/{\displaystyle \sum _{d=1}^{3^n}{L}_d^2(t)},(d=1,2,\cdots ,N_3)$$

(11)

The node number in the third layer is equal to that in the second layer, meaning that $N_2=N_3=3^n$.

The fourth layer and the fifth layer in this paper constitute the consequent network part with three parallel subnetwork outputs corresponding to estimated length variations of different contractile PAMs. Each node in the consequent network part carries out a selection of inputs of this layer according to the if–then fuzzy logic rule. Thus, each subnetwork has the same number of nodes as the third layer ($N_4=N_3$). As shown in Equation (12), $L_{dk}^4(t)(d=1,2,\cdots ,N_4)$ defines the output of the kth subnetwork in the fourth layer, where $p_{di}^k(t)$ represents the consequent parameter and $k$ indicates different subnetworks, respectively.

$$L_{dk}^4(t)={\overline{\omega }}_h{\displaystyle \sum _{i=1}^n{p}_{di}^k{A}_i(t)},(d=1,2,\cdots ,N_4;k=1,2,3)$$

(12)

The final outputs $L_k^5(t)$ of the SDCHM are obtained in the fifth layer with the calculation shown in Equation (13):

$$L_k^5(t)={\widehat{y}}_k(t)={\displaystyle \sum _{d=1}^{N_4}{L}_d^4(t)},(k=1,2,3)$$

(13)

In the above stiffness-dependent submodule $g(\cdot )$, $m_{ij}$, ${\sigma }_{ij}$, and $p_{di}^k$ should be identified by experimental data.

4. Decoupling Inverse Compensation

On the basis of the newly developed SDCHM model, an inverse and decoupling compensator is constructed to mitigate coupled hysteresis. In order to compensate for the SDCHM model with three submodules, the stiffness-dependent and decoupled inverse compensator (SDIHC) also includes three parts: the inverse hysteresis part, the decoupling part, and the inverse stiffness-dependent part, as shown in Figure 8.

First, the inverse stiffness-dependent submodule is constructed. $g^{-1}(\cdot )$ in Figure 8 defines the inversion of $g(\cdot )$. Because $g(\cdot )$ is calculated by the TSFNN, $g^{-1}(\cdot )$ is also constructed by the same structure, as shown in Figure 7. The inverse modeling based on neural networks has been utilized for a valid method for identifying the inversion of a known model. Particularly, analytical solutions of the inverse model are challenging to determine.

Second, in the decoupling submodule, the inverse operation of Equation (8) is derived as follows:

$$\{\begin{array}{c}{v}_{1d}(t)={\displaystyle \frac{1}{G_{11}}}(h_{1d}-G_{21}{v}_{2d}(t)-G_{31}{v}_{3d}(t)),\\ v_{2d}(t)={\displaystyle \frac{1}{G_{22}}}(h_{2d}-G_{12}{v}_{1d}(t)-G_{32}{v}_{3d}(t)),\\ v_{3d}(t)={\displaystyle \frac{1}{G_{33}}}(h_{3d}-G_{13}{v}_{1d}(t)-G_{23}{v}_{2d}(t)).\end{array}$$

(14)

Thirdly, the inverse hysteresis component is constructed. The inverse compensation method for the GPI model is applied as the following Equation (15):

$$\{\begin{array}{c}{\widehat{p}}_1(t)=GP{I}^{-1}[v_{1d}](t),\\ {\widehat{p}}_2(t)=GP{I}^{-1}[v_{2d}](t),\\ {\widehat{p}}_3(t)=GP{I}^{-1}[v_{3d}](t).\end{array}$$

(15)

Because the envelop functions ${\gamma }_R$ and ${\gamma }_L$ are invertible, the inverse model of the GPI model is expressed as

$${\widehat{p}}_k(t)=GP{I}^{-1}[{\nu }_{kd}](t)=\{\begin{array}{l}{\gamma }_R^{-1}\circ {\mathsf{\Gamma }}^{-1}\circ {\nu }_{kd}(t),\mathrm{if}{\nu }_{kd}(t)>{\nu }_{kd}(t-1),\\ {\gamma }_L^{-1}\circ {\mathsf{\Gamma }}^{-1}\circ {\nu }_{kd}(t),\mathrm{if}{\nu }_{kd}(t)<{\nu }_{kd}(t-1),\\ {\widehat{p}}_k(t-1),\mathrm{if}{\nu }_{kd}(t)={\nu }_{kd}(t-1).\end{array}(k=1,2,3)$$

(16)

where:

$$\begin{array}{ll}{\mathsf{\Gamma }}^{-1}\circ {\nu }_{kd}(t)& ={\mathsf{\Gamma }}^{-1}[{\nu }_{kd}](t)\\ & ={\widehat{x}}^k(r_0^k){\nu }_{kd}(t)+{\displaystyle \sum _{i=1}^N{\widehat{x}}^k({\widehat{r}}_i^k)}{F}_{{\widehat{r}}_i^k}[{\nu }_{kd}](t)\end{array}$$

(17)

with

$${\widehat{x}}^k(r_0^k)={\displaystyle \frac{1}{x(r_0^k)}}$$

(18)

$${\widehat{r}}_i^k=x(r_0^k)r_i^k+{\displaystyle \sum _{j=1}^{i-1}{x}^k(r_j^k)(r_i^k-r_j^k)}$$

(19)

$${\widehat{x}}^k({\widehat{r}}_i^k)=-{\displaystyle \frac{x^k(r_i^k)}{(x^k(r_0^k)+{\displaystyle \sum _{j=1}^i{x}^k(r_j^k)})(x^k(r_0^k)+{\displaystyle \sum _{j=1}^{i-1}{x}^k(r_j^k)})}}$$

(20)

5. Experimental Results and Discussion

In this section, experiments are performed to verify the proposed SDCHM and its inverse compensator, SDIHC. The SDCHM identification is conducted firstly to identify parameters. Then, the model validation is completed to show the prediction accuracy of the SDCHM. Finally, the SDIHC is implemented for compensating and decoupling the outputs of the soft manipulator.

5.1. SDCHM Model Identification and Validation

Since the SDCHM includes three components, identification procedures for the inherent hysteresis submodule, actuator coupling submodule, and stiffness-dependent submodule are conducted in three steps in order.

In the first step, outputs of the actuator coupling submodule and stiffness-dependent submodule are set to 1. Parameters including $x_i^k(r_i^k)$, $[a_1^k,a_2^k,a_3^k]$, and $[b_1^k,b_2^k,b_3^k]$ in three GPI models are identified. Since the GPI is utilized to describe the inherent hysteresis behavior of the single PAM, pressure-displacement experimental data of each contractile PAM used for identification are collected separately. Taking the identification of ${\mathrm{GPI}}_1$ for instance, the input pressure signal for the manipulator is set as Equations (21) and (22):

$$\{\begin{array}{c}{p}_1(t)={\displaystyle \frac{A_s(t)}{2}}\sin (0.04\pi t)+{\displaystyle \frac{A_s(t)}{2}},\\ p_2(t)=p_3(t)=p_0(t)=0.\end{array}$$

(21)

and

$$A_s(t)=\{\begin{array}{c}1.0\mathrm{bar},\\ 1.5\mathrm{bar},\\ 2.0\mathrm{bar},\\ 2.5\mathrm{bar},\end{array}\begin{array}{c}t\in [0s,50s),\\ t\in [51s,100s),\\ t\in [101s,150s),\\ t\in [151s,200s).\end{array}$$

(22)

Input signal $p_1\left(t\right)$ and corresponding output $y_1\left(t\right)$ of the PAM I, experimentally measured with the experimental setup in Figure 2 with the sampling frequency of 50 Hz, are adopted to identify parameters in ${\mathrm{GPI}}_1$, and they can be found using least square (LS) optimization. The number of generalized play operators is chosen to be N = 10, which decides the model’s accuracy and complexity.

Table 2. Influence of the number of generalized play operators.

N	5	10	15	20
accuracy of ${\mathrm{GPI}}_1$ (Err)	2.65%	1.36%	1.52%	1.63%
number of identified parameters	12	17	22	27

reflects the influence of N on model accuracy and complexity, where $Err=\left({\displaystyle \sum _{t=1}^H\left|y_1(t)-{\widehat{y}}_1\left(t\right)\right|}/{\displaystyle \sum _{t=1}^H\left|y_1(t)-y_{1,mean}(t)\right|}\right)\times 100\%$, and $y_{1,mean}$ computes the mean value of $y_1\left(t\right)$. This indicates that as N increases, the Err value decreases, and it tends to saturate when $N\ge 10$. By repeating the above process, other parameters to be identified in ${\mathrm{GPI}}_2$ and ${\mathrm{GPI}}_3$ can be found. The final identification results of the inherent hysteresis model are presented in Figure 9.

In the second step, parameter identifications of the actuator coupling submodule focus on coupling coefficients $G_{ij}(i=1,2,3;j=1,2,3)$. The stiffness-dependent submodule is still set to 1. $G_{11}=G_{22}=G_{33}=1$ are set to consider the cascade connection between ${\mathrm{GPI}}_i$ and $G_{ii}$, (i = 1, 2, 3). $G_{ij}(i\ne j)$ stands for the effect of the PAM i on the PAM j. Identifications for each parameter $G_{ij}(i\ne j)$ are the same, so we only state the identification procedure of parameter $G_{12}$ in detail. First, we apply the sinusoidal signal in Equation (21) into the identified model ${\mathrm{GPI}}_1$ to obtain the output ${\nu }_1={\mathrm{GPI}}_1[p_1](t)$; then, $G_{12}$ can be determined by the LS method to obtain the minimum result of the objective function $\min \{{[C_{21}{G}_{21}-d_{12}]}^T[C_{21}{G}_{21}-d_{12}]\}$, in which $C_{21}=[{\nu }_1[1],{\nu }_1[2],\cdots ,{\nu }_1[N_c]]$, $d_{12}=[y_2[1],y_2[2],\cdots ,y_2[N_c]]$. $N_c$ indicates the number of experimental data involved in optimization and relies on signal period and sampling frequency. All optimized coefficients $G_{ij}(i=1,2,3;j=1,2,3)$ are listed in

Table 3. Optimized coefficients G_ij.

${\mathit{G}}_{\mathit{i}\mathit{j}}$	i = 1	i = 2	i = 3
j = 1	1	0.1521	0.5368
j = 2	0.1124	1	0.6321
j = 3	0.0356	0.2153	1

.

In the third step, identification of the stiffness-dependent submodule is, in fact, intended to train the TSFNN in Figure 7. The stiffness-dependent submodule is employed to model the output stiffness dependence of the pneumatic soft manipulator and other nonlinearities unable to be completely established in the first two submodels. As mentioned in Figure 6, inputs to the stiffness-dependent submodule consist of virtual displacements generated from the actuator coupling submodule and $p_0(t)$ in the extensor PAM. Therefore, the sample set $S=\{INPUT,OUTPUT\}$ is constructed with the following structure:

$$INPUT=\left[\begin{array}{cccc}{H}_1[t]& 0& 0& H_1[t]\\ 0& H_2[t]& 0& 0\\ 0& 0& H_3[t]& 0\\ p_0[t]& p_0[t]& p_0[t]& p_0[t]\end{array}\right]$$

(23)

with

$$\{\begin{array}{c}{H}_1[t]={\displaystyle \sum _{i=1}^3[G_{i1}\cdot {\mathrm{GPI}}_i[p_i][t]],}\\ H_2[t]={\displaystyle \sum _{i=1}^3[G_{i2}\cdot {\mathrm{GPI}}_i[p_i][t]}],\\ H_3[t]={\displaystyle \sum _{i=1}^3[G_{i3}\cdot {\mathrm{GPI}}_i[p_i][t]}].\end{array}\mathrm{with}t=1,2,\cdots ,N_1$$

(24)

where $H_i,(i=1,2,3)$ is obtained from the output of the actuator coupling submodule. Contractile PAMs are inflated with pressure $p_{\mathrm{i}}[t]={\displaystyle \frac{A_s(t)}{2}}\sin (0.04\pi t)+{\displaystyle \frac{A_s(t)}{2}}, (i=1, 2, 3)$ in Equation (21), and $N_1=f_r\cdot t_{f1}\cdot T_s$, where the sampling frequency is $f_r=20\mathrm{Hz}$, the number of cycles is $t_{f1}=2$, and the period is $T_s=50s$. The output sample set, $OUTPUT$, is shown in the following Equation (25):

$$OUTPUT=\left[\begin{array}{cccc}{y}_1^1[t]& y_1^2[t]& y_1^3[t]& y_1^4[t]\\ y_2^1[t]& y_2^2[t]& y_2^3[t]& y_2^4[t]\\ y_3^1[t]& y_3^2[t]& y_3^3[t]& y_3^4[t]\end{array}\right],\mathrm{with}t=1,2,\cdots ,N_1$$

(25)

In Equation (25), $y_1[t]$, $y_2[t]$, and $y_3[t]$ are the experimental displacements of contractile PAM I, PAM II, and PAM III, respectively. Each column in matrix (25) represents experimentally collected outputs of contractile PAMs in the soft manipulator corresponding to every column input of Equation (24), and superscripts mark the number of columns.

In order to highlight the influence of manipulator stiffness on the manipulator outputs, experiments are repeated under different inflation pressures $p_0$, and $p_0$ increases from 0 to 2.5 bar, with an increment of 0.1 bar. In order to sufficiently train the stiffness-dependent submodel, among the experimental data collected above, the data corresponding to pressure $p_0\in [0 \mathrm{bar},0.2 \mathrm{bar},0.4 \mathrm{bar},0.6 \mathrm{bar},0.8 \mathrm{bar},1.0 \mathrm{bar},1.2 \mathrm{bar},1.4 \mathrm{bar},1.6 \mathrm{bar},1.8 \mathrm{bar},2.0 \mathrm{bar},2.2 \mathrm{bar},2.5 \mathrm{bar}]$ are selected as the training set. Meanwhile, data corresponding to pressure $p_0\in [0.1 \mathrm{bar},0.5 \mathrm{bar},0.7 \mathrm{bar},1.1 \mathrm{bar},1.5 \mathrm{bar},1.7 \mathrm{bar},2.3 \mathrm{bar}]$ are selected as the validation set, and the rest of the samplings are utilized as the test set. Hence, 104,000 training samples in total are employed to train the stiffness-dependent submodule, and 56,000 samples and 48,000 samples act as the validation and test set, respectively. By providing the TSFNN with a diverse dataset that comprehensively covers the normal operating pressure ranges of PAMs, its generalization ability can be enhanced. Furthermore, through training and optimization with large-scale sample data, the computational efficiency and response speed of the TSFNN can be improved.

During the training process of the stiffness-dependent submodule, the particle swarm optimization method is used for antecedent parameter identifications due to its excellent global and local search capacity, and the LS method is utilized to complete consequent parameter identifications to improve the training efficiency of the stiffness-dependent submodule. The loss function is shown in Equation (26):

$$loss={\displaystyle \frac{{\displaystyle \sum _{t=1}^{N_1}[{(y_1[t]-{\widehat{y}}_1[t])}^2+{(y_2[t]-{\widehat{y}}_2[t])}^2+{(y_3[t]-{\widehat{y}}_3[t])}^2]}}{N_1}}$$

(26)

where ${\widehat{y}}_1[t]$, ${\widehat{y}}_2[t]$, and ${\widehat{y}}_3[t]$ are model results regarding the length variations of PAM I, PAM II, and PAM III.

Based on the above analysis, Algorithm 1 gives the identification process of the proposed SDCHM. The neural network training and experiments discussed in this paper are all conducted on a computer equipped with a 3.40 GHz × 64 Intel Core i7 processor, 24 GB of RAM, and a GeForce RTX 3050 GPU.

Algorithm 1: Training process of the SDCHM model

Inputs: The training datasets of pressure and displacements of the soft manipulator

Outputs: The displacement outputs of the manipulator from the SDCHM model Set: number of generalized play operators N = 10; maximum iterations $ite=50$; the population size of particles $L=20$; learning factors $c_1=c_2=0.5$; Initialize parameters: Parameters in the inherent hysteresis submodule such as $x^k(r_i^k)=0,(i=1,2,\cdots ,N;k=1,2,3)$, $[a_1^k,a_2^k,a_3^k]=[0,0,0],(k=1,2,3)$, $[b_1^k,b_2^k,b_3^k]=[0,0,0],(k=1,2,3)$; parameters in the actuator coupling submodule such as $G_{ij}=1,(i=1,2,3;j=1,2,3)$; parameters in the stiffness-dependent submodule such as $m_{ij}~N(0,1)$, ${\sigma }_{ij}~N(0,1)$, and $p_{di}^k~N(0,1)$; 1: Optimize parameters of $x^k(r_i^k)$, $[a_1^k,a_2^k,a_3^k]$, and $[b_1^k,b_2^k,b_3^k]$ in the inherent hysteresis submodule by LS method with input and output data of each contractile PAM in the soft manipulator, and outputs of actuator coupling submodule and stiffness-dependent submodule are set as 1. 2: Obtain the identified inherent hysteresis submodule output and input into the actuator coupling submodule; 3: Calculate the actuator coupling submodule output according to Equation (8); 4: Optimize parameters of $G_{ij},(i=1,2,3;j=1,2,3)$ in the actuator coupling submodule by LS method with objective function $\min \{{[C_{ij}{G}_{ij}-d_{ij}]}^T[C_{ij}{G}_{ij}-d_{ij}]\}$, and output of the stiffness-dependent submodule is set as 1; 5: Obtain the identified actuator coupling submodule output and input into the stiffness-dependent submodule; 6: While $i\le ite$,   $i\leftarrow i+1$   Calculate the stiffness-dependent submodule output according to Equations (9)–(13);   Compute loss function $loss={\displaystyle \frac{{\displaystyle \sum _{k=1}^{N_1}[{(y_1[t]-{\widehat{y}}_1[t])}^2+{(y_2[t]-{\widehat{y}}_2[t])}^2+{(y_3[t]-{\widehat{y}}_3[t])}^2]}}{N_1}}$;   Optimize the consequent parameters by the LS method;   Compute loss function again $loss={\displaystyle \frac{{\displaystyle \sum _{k=1}^{N_1}[{(y_1[t]-{\widehat{y}}_1[t])}^2+{(y_2[t]-{\widehat{y}}_2[t])}^2+{(y_3[t]-{\widehat{y}}_3[t])}^2]}}{N_1}}$;   Optimize the antecedent parameters by the particle swarm optimization method;   Update the parameters. End while Return training results of the SDCHM model

To validate the identified SDCHM model, the stiffness-dependent and coupled hysteresis outputs of the pneumatic manipulator corresponding to pressure in the contractile PAM $p^{*}=0.7\sin \left(0.04\pi t\right)+0.7$ bar and $p_0=0.1,1.5,1.7$ bar are predicted, respectively. To further evaluate the performance of the SDCHM, the results of the TSFNN model whose structure is the same to that of the stiffness-dependent submodule are employed as control groups. For example, Figure 10 displays the comparisons among the proposed model results, TSFNN model results, and the experimental data with $p_0=0.1$ bar. Figures in the first line of Figure 10 are the manipulator outputs corresponding to pressures of $[p_1,p_2,p_3]=[p^{*},0,0]$, and the figures in the middle two rows of Figure 10 are manipulator outputs under inputs of $[p_1,p_2,p_3]=[0,p^{*},0]$ and $[p_1,p_2,p_3]=[0,0,p^{*}]$, respectively. The remaining figures in Figure 10 are manipulator outputs with $[p_1,p_2,p_3]=[p^{*},p^{*},p^{*}]$. The results intuitively indicate that the SDCHM models are able to predict the stiffness-dependent hysteresis coupling behavior of the soft manipulator with high accuracy, surpassing the precision of the TSFNN model. This can be attributed to the limited number of fuzzy sets in the TSFNN used for interpolation and its simple structure when not combined with the inherent hysteresis submodule and stiffness-dependent model, making it difficult to meet the demands of high-precision modeling. Increasing the number of fuzzy sets in the TSFNN would significantly prolong the training time, and it would require collecting more experimental datasets to prevent overfitting in the training process. To quantitatively calculate the model precision, Err, and root mean square error (RMSE), the RMSE is defined as:

$$RMS{E}_i=\sqrt{{\displaystyle \sum _{j=1}^{num}{(y_{ij}-{\widehat{y}}_{ij})}^2}/num},(i=1,2,3)$$

(27)

where num denotes the number of data involved in the calculation and ${\widehat{y}}_{ij}(i=1,2,3)$ and $y_{ij}(i=1,2,3)$ represent SDCHM model outputs and real outputs of contractile PAMs.

Table 4. Comparisons of prediction accuracy between different models.

	SDCHM Model		TSFNN Model
Group	Err (%)	RMSE (mm)	Err (%)	RMSE (mm)
(a)	1.74	1.17	16.45	11.27
(b)	1.03	0.45	10.59	5.49
(c)	1.02	0.40	13.21	7.20
(d)	1.02	0.43	10.35	4.80
(e)	1.70	1.13	17.37	12.98
(f)	1.02	0.40	7.83	2.60
(g)	0.98	0.37	6.32	2.24
(h)	0.96	0.32	6.53	2.69
(i)	1.71	1.15	10.26	4.84
(j)	2.88	2.14	10.10	4.01
(k)	1.53	1.03	16.12	10.99
(l)	1.55	1.04	18.11	13.64

depicts comparisons of prediction accuracy between different models based on the results of Figure 10. Figure 11 shows the RMSE of the model results with $p_0=0.1,1.5,1.7$ bar, intuitively indicating that the SDCHM model has better prediction ability and demonstrates good generalization under various working conditions of the soft manipulator. As mentioned in Section 2, contractile PAMs with the same pressure show inconsistent output behavior because of manual fabrication errors. However, the SDCHM model is capable of capturing complex nonlinear effects.

5.2. Implementation of the Decoupling Inverse Compensator SDIHC

According to Figure 8, the implementation of the SDIHC also requires three parts: the inverse stiffness-dependent submodule identification, decoupling submodule identification, and the inverse hysteresis submodule identification.

Part 1: The inverse stiffness-dependent submodel is identified. The network structure of the $g^{-1}(\cdot )$ depicted in the Figure 8 is identical to the TSFNN shown in Figure 6. To identify the inverse stiffness-dependent submodule, inputs and outputs of $S=\{INPUT,OUTPUT\}$ are taken as the corresponding outputs and inputs of $g^{-1}(\cdot )$. During training of the inverse stiffness-dependent submodule, uniform optimization methods and training parameters like the iteration number and the population size of particles are employed for training the stiffness-dependent submodule.

Part 2: The decoupling submodule is identified. According to the identification of the actuator coupling submodule, parameters in the decoupling submodule have been confirmed. We adopted the Jacobi iterative method [40] to solve the linear simultaneous Equation (15) for programming and computation simplification. Thus, the following iterative formula is used, and m represents the following iterations:

$${\nu }_d^{(m+1)}=B_J{\nu }_d^{(m)}+D^{-1}b$$

(28)

with

$$\{\begin{array}{c}{\nu }_d=[{\nu }_{1d},{\nu }_{2d},{\nu }_{3d}],\\ B_J=D^{-1}(D-A),\\ D=\mathrm{diag}(G_{11},G_{22},G_{33}),\\ A=\left[\begin{array}{ccc}{G}_{11}& G_{21}& G_{31}\\ G_{12}& G_{22}& G_{32}\\ G_{13}& G_{23}& G_{33}\end{array}\right],\\ b=\left[h_{1d},h_{2d},h_{3d}\right].\end{array}$$

(29)

When the error of two adjacent iteration results $\left|{\nu }_{id}^{(m+1)}-{\nu }_{id}^{(m)}\right|\le 0.01,(i=1,2,3)$, the iteration is suspended, and outputs of the decoupling submodel are finally confirmed.

Proof of Convergence for Equation (28).

The Jacobi iterative method declared that, assuming the Jacobi iterative matrix is $\Vert B_J\Vert <1$, the Jacobi iterative method converges to any initial value. According to the identified actuator coupling submodule,

$$\{\begin{array}{c}D=\mathrm{diag}(1,1,1),\\ A=[\begin{array}{ccc}1& 0.1521& 0.5368\end{array};\begin{array}{ccc}0.1124& 1& 0.6321\end{array};\begin{array}{ccc}0.0356& 0.2153& 1\end{array}].\end{array}$$

(30)

Thus, $B_J=D^{-1}(D-A)=[\begin{array}{ccc}0& -0.1521& -0.5368\end{array};\begin{array}{ccc}-0.1124& 0& -0.6321\end{array};\begin{array}{ccc}-0.0356& -0.2153& 0\end{array}]$ and $\Vert B_J\Vert =-0.016<1$. The Jacobi iterative method is convergent for Equation (27). □

Part 3: The inverse hysteresis submodule is identified. The GPI model has an analytical inverse solution, and based on the determined parameters of the SDCHM model in the above section, the inverse hysteresis submodule can be completed following Equations (15)–(20).

Figure 12 illustrates the relationship between desired outputs and compensation results with the SDIHC and its control groups. In the first line of Figure 12, the desired outputs of $y_{1d}=6\sin \left(0.04\pi t\right)+33\mathrm{mm}$, $y_{2d}=y_{3d}=33\mathrm{mm}$ with $p_0=0.7\mathrm{bar}$ are tracked. In the second line of Figure 12, the desired outputs of $y_{2d}=6\sin \left(0.04\pi t\right)+33\mathrm{mm}$, $y_{1d}=y_{3d}=33\mathrm{mm}$ with $p_0=1.5\mathrm{bar}$ are tracked; and in the last line, the desired outputs of $y_{3d}=6\sin \left(0.04\pi t\right)+33\mathrm{mm}$, $y_{1d}=y_{2d}=33\mathrm{mm}$ with $p_0=2.3\mathrm{bar}$ are tracked. Yellow lines in Figure 12 indicate the SDIHC compensator results. As control groups, experiments of individual hysteresis compensator without the decoupling and inverse stiffness-dependent submodules are also conducted, and these are represented by blue lines. Tangerine lines denote the ideal linearity between the desired outputs and reference inputs, and black lines show PID control results.

Table 5. RMSE and $e_w$ comparisons of different compensation methods.

	SDIHC Compensator		Individual Hysteresis Compensator		PID-Based Compensator
Group	${\mathit{e}}_{\mathit{w}}$ (%)	RMSE (mm)	${\mathit{e}}_{\mathit{w}}$ (%)	RMSE (mm)	(%)	RMSE (mm)
(a)	0.58	0.21	6.12	2.31	2.89	3.01
(b)	0.87	1.09	7.13	7.55	3.51	4.33
(c)	1.02	1.21	6.78	8.53	4.03	4.96
(d)	0.83	1.06	7.02	9.12	3.92	4.01
(e)	0.62	0.34	3.58	3.55	1.75	3.36
(f)	0.75	1.02	6.88	11.10	3.42	3.89
(g)	1.86	1.33	10.23	13.86	6.21	5.87
(h)	2.33	1.56	8.89	10.21	3.89	5.12
(i)	0.61	0.29	1.89	2.03	2.05	1.96

shows the RMSE and $e_w$ of different compensators. It can be found that the results without decoupling and inverse stiffness-dependent submodules lead to the worst accuracy because the coupling between parallel PAMs play a significant role in the system behaviors. The SDIHC algorithm successfully completed hysteresis compensation and decoupling, and the RMSE of the compensation results was no more than 1.56 mm. As for the hysteresis compensation effect of the SDIHC compensator, the maximum width ratio $e_w$ of the hysteresis loops with different stiffness states after applying the SDIHC compensator was within the range of [0.58%, 2.33%], clearly lower than the experimental results before compensation in Section 2.

To further validate the SDIHC, a desired signal $y_{1d}=y_{2d}=y_{3d}=abs[40\sin \left(0.04\pi t\right)+5]\mathrm{mm}$ with $p_0=0.6\mathrm{bar}$ is tracked in Figure 13. The steady-state error with the proposed SDIHC was within the range of [−1.33 mm, 1.25 mm]. It can be observed that the proposed compensator provided an effective method for precise end position control of the pneumatic soft manipulator.

6. Conclusions

Soft manipulators with different PAMs parallel to each other process asymmetric hysteresis, actuator couplings, and stiffness-dependent nonlinearities, leading to inaccuracy positioning. To compensate for system nonlinearities and to decouple the actuator coupling effects, a SDCHM model for the parallel PAMs-driven manipulator is proposed. The SDCHM presents a cascaded model architecture with an inherent hysteresis submodule, an actuator coupling submodule, and a stiffness-dependent submodule. Based on the established model, its compensator SDIHC can be achieved easily without any simplification, allowing it to compensate and decouple for actuator hysteresis couplings at different stiffness states. The experimental results suggest that the SDCHM performs well and possesses good accuracy in predicting inputs–outputs of the soft manipulator in different working conditions. In addition, the compensation effect of the inverse compensator SDIHC is significantly better than that of the controller only with the inverse hysteresis compensator, successfully mitigating the hysteresis coupling phenomenon in the PAM-based soft manipulator. Although the proposed modeling and compensating methods are showcased on a manipulator driven by three parallel PAMs, they are believed to be easily extended to parallel driving systems based on more PAMs. However, the proposed model has not yet considered the impacts of high driving frequencies and external interferences, such as changes in external loads, on the hysteresis couplings of the PAM-based manipulator. Additionally, extending the proposed method to soft manipulators with multiple sections is worth further study, and a closed-loop strategy based on the proposed model will be the next key point for more robust control.