1. Introduction
Soft Robotics is an emerging branch of bio-inspired robotics fostering the development of biomedical and assistive technologies as well as the design of control methodologies for the optimization of human–robot interaction (HRI) in prosthetics, rehabilitation and surgical robotics.
In the field of Soft Robotics, the safety constraints involved in HRI can be fulfilled through adjustable stiffness mechanisms actuated by compliant actuators [
1]; the control of the adaptable compliance can be achieved through some control methodologies inspired by the human musculoskeletal system.
Among the actuation technologies implemented in the control systems of soft robots, Pneumatic Artificial Muscles (PAMs) [
2] are good candidates as biomimetic, light, safe and efficient actuators for rehabilitation robots, wearable exoskeleton robots and energy-efficient walking humanoids (see [
3,
4,
5,
6]).
Analogously with the biological muscle, a PAM only works in the contraction direction; the “muscular activation” generating the contraction force is obtained by inflating the air chamber of the PAM. The contraction force is a nonlinear function of the internal pressure of the PAM. Moreover, the mechanical response of the actuator exhibits an hysteretic behaviour, since the force vs. length characteristics of a PAM depends on the inflation/deflation phase of the actuator.
Nonlinear and hysteretic effects must not be a deterrent to the adoption of PAMs as compliant actuators for assistive and rehabilitation robots, since control engineers have developed nonlinear compensation control strategies that can achieve robust trajectory tracking control of rehabilitation exoskeletons [
7] and prostheses [
8]. The most efficient compensation strategies rely on gray-box models, such as the phenomenological model in [
9] and the lumped-parameter model in [
10]. The visco-elastic response of a PAM can be reproduced through a parallel configuration of a spring, a damper and an active force element (see [
9,
10,
11]). More recently, Al-Ibadi et al. [
12] have proposed a general model of the PAM force depending on the key parameters of initial length, initial inner diameter, the thickness of both the inner rubber tube and the braided sleeve of the actuator. The aforementioned models are not able to describe the hysteretic component of the actuator output.
The friction in PAMs that is generated by the internal structure of the elastic bladder enclosing the pressurized chamber of the actuator is the main source of hysteresis. Moreover, the force/length hysteresis loops, along the cycles of inflation and deflation of the actuator, are asymmetric.
In the majority of the literature references on the trajectory tracking control of PAMs, the hysteresis component of the actuator response is discarded, see, e.g., [
13], whereas other important results on the hysteresis modelling in PAMs are not concerned with the asymmetry of hysteresis loops (see [
14,
15]). The modelling and identification of asymmetric hysteresis loops have been only recently addressed in [
16,
17].
The work [
15] is concerned with the modelling and compensation of the hysteresis in a robotic joint actuated by PAMs. The hysteretic relationship between force and length of a PAM is approximated through a combination of four elementary visco-elastic blocks of a Maxwell-slip model that yields symmetric loops.
Preisach [
18], Prandt–Ishlinskii [
19], Bouc–Wen [
20,
21] and Maxwell-slip [
22] models that can capture the symmetric hysteresis of general friction phenomena may be applied to the identification of the mechanical response of PAMs. However, these models cannot fit the asymmetry of the hysteresis loops in PAMs and, therefore, extended versions are more suitable to PAM control, like the modified Prandt–Ishlinskii model in [
17].
The identification of hysteresis requires, in general, a large number of single slip-blocks or play operators in Maxwell-slip and Prandt–Ishlinskii models, respectively; the number of model parameters to be identified may increase if the hysteresis is not symmetric. For instance, in [
17], a number of 40 parameters is required to capture the asymmetric length/pressure hysteresis of a PAM. The number of parameters can be reduced to 9 if the generalized Bouc–Wen model is used to fit asymmetric loops in the force vs. length curve of a PAM (see [
16]).
This work proposes a novel multistate model of asymmetric hysteresis, whose structure and number of parameters enable the efficient identification and compensation of hysteresis, e.g., in motion control applications. The structure of the model is comparable to that of other multistate friction models [
23], where a parallel connection of multiple elasto-plastic elements describes, in the presliding regime of motion, the symmetric hysteresis and its properties of nonlocal memory and rate independency.
Hysteresis with nonlocal memory is defined as an input–output map, for which the output at any time instant depends on the output at some time instant in the past, the input since then and the past extremum values of the input or the output [
24,
25]. Rate-independency can be also assumed since the amplitude of the hysteresis loops in the presliding regime depends on displacement rather than velocity.
From the multistate structure, our model inherits the aforementioned properties of friction hysteresis. An important extension provided by our model is the identification of asymmetric hysteresis loops. The model is applied to the experimental characterization of a commercial PAM, by identifying the asymmetric hysteresis observed in the nonlinear relationship force vs. length of the actuator, during the cycles of inflation/deflation of the elastic chamber of the PAM.
The paper is structured as follows:
Section 2 points out the main experimental findings of the hysteresis in PAMs. These findings are obtained through the experimental setup developed at the Biomechatronics Laboratory of the University of Catanzaro “Magna Græcia”, Catanzaro, Italy. The experimental observations are taken into account in the subsequent development of the friction model of asymmetric hysteresis presented in
Section 3.
Section 4 presents and discusses the results of the experimental identification of a commercial PAM on the basis of the multistate friction model, together with a comparison with other models. The paper closes with some concluding remarks.
Notation: The saturation function
, the extended signum function
are defined as follows:
where
and
.
2. Hysteresis Extraction from the Mechanical Response of a Fluidic Muscle Actuator
The force output of a PAM is affected by the hysteresis arising from the friction effects between the braided cords of the actuator shell and the rubber bladder enclosing the pressurized chamber of the PAM, as experimentally observed in [
26,
27].
The friction effects in PAMs look like the adhesive forces in dry friction between two contacting surfaces during the presliding phase of motion. The friction forces at the asperity contacts are a function of displacement; the asperities at the interface of the contacting materials deform elasto-plastically, until the displacement increases and the asperity contacts are broken away. Therefore, elasto-plastic models are suitable to the description of the dynamic friction in the presliding regime of motion, as well as to the characterization of the hysteresis in PAMs as a function of the cycles of inflation and deflation of the actuator.
The hysteresis observed in the mechanical response of PAMs is characterized by nonlocal memory and quasi-rate independency and, therefore, the hysteretic behavior of PAMs is similar to that of the presliding regime of motion. Quasi-rate independency is referred to the invariance of the shape of the hysteretic loops within a certain range of velocity of relative contraction/stretching of a PAM.
Other than the friction associated with the purely dynamic part of the actuator response, the force vs. length curve of a PAM has a static component expressed as a nonlinear function of internal pressure and length of the actuator itself.
Both static and dynamic components can be experimentally identified through the procedure originally proposed in [
15], after writing the total force of the PAM as
where
represents the contribution of the static force and
is the hysteresis force.
The term will be efficiently identified through an appropriate hysteresis function with nonlocal memory, as discussed in the next section.
A commercial PAM FESTO fluidic muscle MAS-10 (FESTO AG, Esslingen, Germany) having a nominal length of 200 mm and a diameter of 10 mm, has been used for the identification procedure.
Firstly, it is necessary to extract hysteresis from the overall force
(
1) as follows:
The test rig depicted in
Figure 1 allows for reproducing two sets of experimental conditions in order to extract hysteresis from (
2); the first one is the isobaric setup for measuring
, the second one is the isometric setup used to identify
as a nonlinear function of length and internal pressure of the PAM.
The muscle force is measured through a load cell 3138 S by Phidgets Inc. (Calgary, AB, Canada). The pressure of the muscle is regulated through the proportional pressure regulator MPPE-3-1/8-6-010-B by FESTO AG (Esslingen, Germany). The measurement of the pressure of the actuator is obtained from a pressure transmitter FESTO SPTW-P16. Pressure and force signals of the fluidic muscle are acquired and processed on a PC, by using an Arduino Uno R3 (8-bit AVR RISC-based microcontroller ATmega 328) as data acquisition board.
In isobaric configuration, the internal pressure of the PAM is kept constant. Then, a controlled deformation is applied to the PAM through a ball screw actuated by a NEMA 17 stepper motor (Quimat, Hong Kong, China) and the corresponding force
is measured.
Figure 2 shows the experimental curves of force vs. length obtained at different values of pressure in the range [1–6] [bar]. These curves show asymmetric hysteresis loops.
The isometric tests are most useful to identify the static force component . The isometric curves do not show hysteresis, since there is no motion of the actuator or stretching of the elastic bladder of the PAM. During the static tests, the PAM is constrained at desired length and initial pressure. Then, the internal pressure of the actuator is increased and the contracting force is measured.
The static force vs. pressure curves, obtained at different actuator lengths, are shown in
Figure 3. The experimental curves allow for formulating the constrained model of the force
as a function of internal pressure
p and length
l of the actuator.
The static force can be approximated by the polynomial function
The structure of the function is chosen to be conveniently linear in the control input p; the polynomial order in l gives a satisfactory trade-off between the goodness-of-fit and the degree of nonlinearity of the model.
The parameters in (
3) have been identified through a multi-variable least squares fitting in MATLAB. The fitting function guarantees a root mean square error (RMSE) of
[N] and a square of the coefficient of determination
. The corresponding values of the parameters are listed in
Table 1.
Figure 4 gives the three-dimensional plot of the fitted
as a function of pressure
p and length
l, together with the experimental isometric curves.
Finally, the hysteresis can be extracted from (
2) by subtracting the output of the polynomial fitting model (
3) to the isobaric measurements of the total force
at given pressure.
In the subsequent development, the extracted hysteresis curves can be used in order to derive a multistate model of the friction-induced hysteresis.
3. A Multistate Friction Model for the Identification of the Asymmetric Hysteresis
A methodology for hysteresis modelling, based on a novel multistate friction model, is proposed here as an efficient approach to capture, through a minimal set of model parameters, the asymmetric hysteresis in the force output of PAMs. The potential applications of the multistate model are not only the characterization of the friction-induced hysteresis in PAMs but also, on a different length scale, the identification of the asymmetric hysteresis in general problems of mechanical friction.
The model is obtained as a parallel connection of
n visco-elastic blocks (see
Figure 5). The resulting nonlinear switched system follows a switching rule
depending on the sign of the velocity input
v. The fitting of asymmetric hysteresis is obtained by introducing into the model two sets of parameters. The first one is mapped to the rising phase of motion
and the second one to the descending phase (
). The components of each elementary block are a spring of stiffness
, a viscous element of damping
and a tip sliding on a rough surface in the presence of a Coulomb friction saturation
. The total force is obtained as sum of the visco-elastic forces of the single elements.
Referring to the physical interpretation of
Figure 5, the friction force acting on the sliding tip is determined by a set-valued Coulomb friction law
according to the relative velocity of the tip
. The dynamic behaviour of the deflection of the tip is captured by the state variable
.
Therefore, the multistate friction model can be written as
where
n is the number of blocks; the model parameters follow the switching rule
The property of nonlocal memory of the hysteresis loops, also observed in PAMs, is recovered by the parallel configuration of
n visco-elastic blocks, as shown in [
23]. The same literature results show that the rate-independency of the shape of the hysteresis loops can be also guaranteed by a parallel setup of multiple visco-elastic elements of
Figure 5.
4. Identification Results and Discussion
The fitting of the extracted hysteresis of the MAS-10 fluidic muscle has been obtained on the basis of the multistate friction model of
Section 3. Through the procedure of
Section 2, the hysteresis force
to be identified has been extracted from the largest isobaric loop
vs. length of the fluidic muscle. The isobaric loop at 4 bar corresponds to the full stretching of the actuator giving the largest hysteresis loop within the pressure range of the fluidic muscle.
An interior-point fmincon solver from the MATLAB Optimization Toolbox has been adopted for minimizing the quadratic cost function weighting the error between the predicted PAM response (model output) and the extracted hysteresis force (experimental data).
The final model used for fitting has three elementary blocks and eleven parameters. A total of six parameters can be associated with the first block; this number is reduced to five by adding the constraint . Another two blocks are assigned to the descending phase and, therefore, they contribute to the PAM dynamics through two additional state variables and six parameters activated under the condition .
The vector of the identified parameters reads
Figure 6 shows the simulation output of the fitted hysteresis model compared to the hysteresis loop extracted from the isobaric experiment at a pressure of 4 bar. The comparison provides a good agreement between the model output and the measured data. The goodness-of-fit has been quantified by
and RMSE = 3.5568 N. The value of the identified parameters are shown in
Table 2. The small error at the upper extremities of the curves can be explained by the nonlinearity originated from the local deformations at the interface between the elastic bladder and the end caps of the actuator during the inversion of motion.
By mapping the model parameters to each phase of motion, it is possible to reduce the model order while reproducing the asymmetry of the measured hysteresis loops. Indeed, the peculiar structure of the multistate friction model allows for using more parameters for fitting the descending branch of the hysteresis loop of the tested PAM. Thus, the number of blocks and parameters can be adapted to the nonlinearity in each phase of motion.
Note that the optimal value of the cost function achieved by the final model with eleven parameters is comparable to that obtained by the eighteen-parameters model (three blocks mapped to the two phases of motion). The higher number of viscoelastic elements (six) in the eighteen-parameters model does not improve the goodness-of-fit. Moreover, in the identification of the final model, the number of optimization iterations decreases and the cost function has a faster convergence to the optimal value, as clear from
Figure 7.
The advantages of the multistate friction model in terms of optimality of the identification performance and reduced number of parameters can be appreciated also in comparison with the related works (see [
16,
17,
28]). Furthermore, the proposed model can be implemented in a feedforward control compensation scheme without the implementation issues of other models.
In [
28], a modified Maxwell-slip (MMS) model is developed for identifying the asymmetric hysteresis of a class of piezo-actuators. The MMS model involves the typical switching structure of the Maxwell-slip form that makes the online parameter identification difficult, as pointed out in [
29]. A total of 42 parameters are involved in the general MMS friction model in [
28], whereas 40 parameters characterize the Prandt–Ishlinskii model of a PAM in [
17]. The approach of Bouc–Wen modelling in [
16] allows for reducing the number of parameters, but it requires the online computation of the internal state variable of the model. Therefore, it is not possible to use the Bouc–Wen model for offline compensation strategies of hysteresis. In this respect, the structure of the multistate friction model provides an important advantage, since a closed-form solution can be derived along the hysteresis loops.
The general applicability of the multistate friction model for compensation of hysteresis in tracking control can be confirmed analogously with other literature references [
30], where a multi-element Maxwell-Slip viscoelastic model is used for implementing the feedforward compensation of hysteresis into a cascade PI–PI tracking control architecture.