Modeling the Static Force of a Festo Pneumatic Muscle Actuator : A New Approach and a Comparison to Existing Models

In this paper, a new approach for modeling the static force characteristic of Festo pneumatic muscle actuators (PMAs) will be presented. The model is physically motivated and therefore gives a deeper understanding of the Festo PMA. After introducing the new model, it will be validated through a comparison to a measured force map of a Festo DMSP-10-250 and a DMSP-20-300, respectively. It will be shown that the error between the new model and the measured data is below 4.4% for the DMSP-10-250 and below 2.35% for the DMSP-20-300. In addition, the quality of the presented model will be compared to the quality of existing models by comparing the maximum error. It can be seen that the newly introduced model is closer to the measured force characteristic of a Festo PMA than any existing model.


Introduction
In recent years, pneumatic muscle actuators (PMAs, Figure 1) were integrated in many robotic systems [1][2][3][4] and are still particularly favored, especially in applications that have to be lightweight and powerful.As can be seen in [6], most robotic systems driven by PMAs, but also driven by other actuators, require an underlying torque controller.While integrating PMAs in robotic systems, the control of PMAs and especially the precise torque control of PMA-driven joints consequently becomes an essential feature.The modeling of PMAs, strongly coupled to model-based controller designs, is therefore still a much discussed topic.
On the one hand, PMAs have their own dynamics, like hysteresis [7], some thermodynamic effects [8,9] and can even be interpreted as a spring-damper combination [10].On the other hand, it has been shown in [11][12][13] that controlling PMAs, with only a static force map approach, can lead to accurate and high-performance controller designs.The reason is that the PMA dynamics are dominated by the fluid dynamics of the air inside the PMA and the stream of mass characterizing the valve that controls the PMA [14,15].Nevertheless, a model describing the static PMA force precisely is crucial for the accuracy of the torque controller.Finding the most accurate model to describe the static force characteristic is therefore one goal of this paper.
Different approaches for modeling the static force characteristic of PMAs can be found in literature.The most popular approach is based only on air compression [16].Although this approach is valid for McKibben PMAs, it is shown by measurements in [14] that this approach does not hold for Festo PMAs.To improve the accuracy of the model presented in [16], the model has been extended in many different ways [8,9,17,18].It has been shown that these extensions are leading to fewer errors between measured static force maps and the forces predicted by the models.Another approach for the specific use on Festo PMAs has been presented in [12].The basic idea is the approximation of a Festo PMA as a piston with a virtual, pressure-dependent piston area and a spring that counteracts the expansion.A second approach for the specific use on Festo PMAs is presented in [19], where the static force characteristic is supposed to behave like a mechanical spring with variable stiffness.Because of uncertain parameters, all models, except the first presented in [16], have in common that they must be identified by minimization of the error between a measured static force map and the force defined by the model.By comparing the maximum error, it will be shown in this paper how the accuracy of existing models varies.Furthermore, a new approach for modeling the static force characteristic of Festo PMAs will be presented.The model is physically motivated and therefore gives a deeper understanding of the Festo PMA.After introducing the new model in Section 3 and an explicit discussion of existing models in Section 4.1, all models will be compared to a measured force map of a Festo DMSP-10-250 (DMSP-<initial inner diameter in millimeter>-<initial length in millimeter>) and a DMSP-20-300 in Section 4.3.It will be shown that the error between the new model and the measured data is below 4.4% for the DMSP-10-250 and below 2.35% for the DMSP-20-300.In addition, the quality of the presented model will be compared to the quality of existing models by comparing the maximum error.

The Static Force Characteristic
A PMA is characterized by its force map.The PMA force F PMA is a function of the muscle pressure p and its length L, as can be found in [5,17,20,21].The force map (Figure 2) can be read in the following way: Starting with an unfastened PMA at initial conditions, which means initial length and atmospheric pressure, the PMA does not exert any force.Now putting pressure inside, the PMA gets shorter until it is fully contracted (about 25% [5]).The exerted force is still zero.While pulling the pressurized PMA back to initial length, the PMA will react with the maximum force.The exerted force can be varied between zero and a length-dependent upper limit by varying the PMA pressure [13].

A New Model of the Static Force Characteristic
PMAs are a combination of a flexible tube and two stiff aluminium connectors (see Figure 1).The membrane that the tube is made of is a combination of a stiff aramid fiber mesh and a flexible rubber that encloses the air inside the PMA.While putting pressure inside the PMA, the membrane expands and, due to the stiff aramid fibers, the PMA gets shorter.According to [17], the approach for the presented PMA model also is that the fiber length L Fiber stays constant and therefore only the membrane rubber deforms.

PMA Volume
Following the ideas of [14,15,17], the PMA volume can be approximated by a cylinder.The volume of a cylinder is a function of length L and diameter D. Cutting the PMA membrane open and flattening it to a plain (Figure 3), the dependency of the diameter on the length can be calculated by using the Pythagoras theorem [17].Supposing that L Fiber is constant, the diameter equation holds, whereby The index 0 indicates the initial state of the PMA, with atmospheric pressure p 0 , initial length L 0 , inner diameter D 0 , membrane thickness H 0 and fiber angle Θ 0 .In [14], it is proven by measurements that the initial fiber angle is Θ 0 = 28.6 • .This is the only value that can be found in literature.The initial membrane thickness, H 0 = 1.8 mm for both the DMSP-10-250 and the DMSP-20-300, can easily be calculated from the given initial inner diameter D 0 and the measurable outer initial diameter of the PMA rubber tube.Inserting the functional dependency on the diameter back into the approximated volume (1), the PMA volume loses its dependency on the diameter and is only a function of the PMA length: (4)

Energy-Based Modeling
While pulling a contracted PMA, the muscle reacts with a force F PMA against the pulling direction.The virtual work of the PMA W PMA is given through Furthermore, the virtual work of the PMA can be separated into two parts.On the one hand, the virtual work W VAE has to be done to change the included air volume.On the other hand, additional virtual work W Elast is necessary to change the potential energy of the elastic membrane rubber.Therefore, always holds.

Virtual Work of the Changing Air Volume
The virtual work needed to change the included air volume inside is given by

Elastic Energy of the Membrane
While putting pressure inside the PMA, the actuator gets shorter and expands.The deformation of the membrane is a plane state strain and can be described by the strain in direction of the PMA length ε L and the PMA perimeter ε PE .Rotating the coordinate system by the membrane fiber angle Θ, the deformation is given by the strain in fiber direction ε AF and the direction of pure rubber ε RU that is perpendicular to the fiber direction (see Figure 4).Following the approach of a constant fiber length, the strain ε AF is always zero.This means that the PMA membrane only expands perpendicularly to the fibers.
For this paper, a one-dimensional state of stress is assumed and, therefore, according to Hooke's law, the tension inside the rubber σ RU will be approximated by where the modulus of elasticity is supposed to be a function of the PMA length.Following the approach of a one-dimensional membrane deformation, the strain is given by the Pythagoras theorem.
Rotating the tension σ RU back to initial coordinates, it is possible to calculate the tension in length direction and the tension in perimeter direction of the PMA Multiplying the tension with the edge surfaces of the unreeled membrane, the virtual work to deform the membrane in length direction is given through This approach neglects any effect of lateral contraction.The virtual work that is necessary to deform the membrane in perimeter direction can be calculated in an analog way: The negative sign in Equation ( 13) is necessary because an increase in length of the PMA has to result in a decreasing perimeter, and this is defined to be positive for this paper.

Summarizing the Energy
Inserting Equations ( 5), ( 7), ( 12) and ( 13) in Equation ( 6), the force of the PMA is given by It must be noted that Equation ( 14) is positive, in case the PMA exerts a pulling force.
Because the modulus of elasticity of the membrane rubber and the fiber angle cannot be measured directly, both terms are identified by minimization.The modulus of elasticity E RU used in Equations ( 10) and ( 11) is approximated by a polynomial of the third order: This order is chosen by experiment.It can be seen that higher order polynomials do not lead to better results.However, for lower order polynomials, the calculated error, defined later in this paper (Equation (23)), is much higher.Furthermore, the approximately known fiber angle Θ 0 is corrected by a constant d 0 : Θ corr 0 is now used as the new, corrected fiber angle.In this paper, all parameters c j and d j (j ∈ N 0 ) have to be identified through solving the minimization problem (22).The parameters c j and d j are calculated in such a way that the quadratic error between a measured force map and the forces predicted by the models is the smallest.A more detailed discussion of the optimization will be given later in this paper.

Existing Models
Defining a valid model describing the static force characteristic of a PMA is still a much discussed topic.The first model that was introduced to describe a McKibben Muscle [16] was only based on the energy that is needed to change the inner PMA air volume.Effects of elasticity of the membrane material were fully neglected: It is shown in [14] that the remaining error between the measured static force characteristic of a Festo PMA and ( 17) is not negligible.The reason for this is that the Festo PMA stores potential energy in its deformed membrane and the McKibben PMA does not.
A modified version of Equation ( 17) is presented in [8,9], where Equation ( 17) has been corrected by two additional factors c 0 and c 1 : Although the parameters c 0 and c 1 were calculated directly in [8,9], in this paper, both parameters will be identified by solving the minimization (22).A third improved, inspired by Equation ( 17) model variant is presented in [18]: A model for the specific use on a Festo PMA is presented in [12].This model is a combination of the pressure-virtual-piston-area product p • A(L) and a length-dependent counter force F c (L).The idea is that the PMA behaves like a combination of a pneumatic piston with a variable piston area and a mechanical spring that counteracts the expansion of the PMA: In [19], the static force characteristic of a Festo PMA is supposed to be equivalent to a mechanical spring, with a displacement-pressure-dependent spring stiffness k(p, ∆L): •∆L , with ∆L = L − min(L). (21)

Test Rig
The static force maps of a Festo DMSP-10-250 and a DMSP-20-300 have been measured with the test rig depicted in Figure 5.During the measurement, all screws are fastened and the PMA length is fixed.The pressure is varied between the length-dependent lower limit and the maximum pressure.While varying the pressure, the PMA force is measured by a 1D-forcesensor KD9363-0.5tand a GSV-1A measurement amplifier, both from ME-Meßsysteme GmbH (Hennigsdorf, Germany).Combining the pressure-dependent force characteristics for different lengths, it is easy to determine the static force characteristic as shown in Figure 2. To reduce the influence of measurement errors, the measurement process has been repeated ten times.The depicted force map therefore shows the mean values.Both measured force maps, for the Festo DMSP-10-250 and the DMSP-20-300, can be found in the appendix in Tables A1 and A2, respectively.

Results
The measured static force maps are characterizing the force behavior of the Festo DMSP-10-250 and DMSP-20-300, respectively.First, the presented PMA models ( 14), ( 18)-( 21) are identified by minimizing the quadratic error between the measured force map and the force map calculated by the model.The optimization has been solved by using the MATLAB (R2015b, The MathWorks, Inc., Natick, MA, USA) function fminsearch: Because the optimization problem ( 22) is nonlinear, only a local, start point-dependent minimum can be found.The start point for all models was identified by iterative testing.It can be seen that good results-in the sense of a small error-can be achieved if the start points are chosen with respect to the physical meaning and within a proper range.According to the SI units-1 m for lengths and 1 Pa = 1 × 10 −5 bar for pressure-any parameter c j and d j (j ∈ N 0 ) is set to 1 if it is only multiplied to a length-dependent factor and to 1 × 10 −5 if it is multiplied by the pressure.The modulus of elasticity of rubber, necessary for Equation (14), is supposed to be in the area of 1 MPa, and therefore the chosen start point is set to 1 × 10 6 .The initial fiber angle Θ 0 in Equation ( 14) is supposed to be correct.The chosen parameter estimation's start points for all models are given without units in Table 1.
Table 1.Start points of parameter estimation for all models without units.
To demonstrate the validity of the the estimated parameters c 0−3 and d 0 defining the new model ( 14), they are given in Table 2.As assumed, the values for the modulus of elasticity are in the range of MPa = 1 × 10 6 Pa, and the initial fiber angle is only slightly corrected by −5.73 • for the DMSP-10-250 and −4.18 • for the DMSP-20-300, respectively.All identified models can now be used for calculating their own specific force maps.The error between the calculated and the measured force map is defined as the maximum difference between a measured and a calculated force point, normalized to the maximum measured force.Multiplied by 100, the error is given by percentage and is shown in Table 3.It can be seen in Table 3 that the presented model is closer to the measured static force map of the Festo DMSP-10-250 and DMSP-20-300 than any existing model.The error of Equation ( 14) for the Festo DMSP-10-250 is smaller than 4.4% and, for the Festo DMSP-20-300, smaller than 2.35%.Furthermore, it can be seen that Equation (17) leads to a maximum error for both PMAs and does not seem to be accurate enough to describe the static force characteristic of a Festo PMA.

Conclusions
In this paper, a new approach for modeling the static force characteristic of Festo PMAs is presented.After a detailed derivation of the new model ( 14), the validity of the model is demonstrated by a comparison to a measured static force map of a Festo DMSP-10-250 and DMSP-20-300, respectively.It is shown that the maximum error between the"true" force and the one predicted with Equation (14) is smaller than 4.4% for the Festo DMSP-10-250 and smaller than 2.35% for the Festo DMSP-20-300.Furthermore, it is shown that the presented model is closer to the measured data than any existing model that can be found in literature.

Figure 1 .
Figure 1.Different Festo pneumatic muscle actuators (PMAs).The depicted connector form is called DMSP.See the data sheet [5] for further information.

Figure 2 .
Figure 2. Measured static force map of a Festo DMSP-10-250 with an initial length of 250 mm and an initial inner diameter of 10 mm.

Figure 4 .
Figure 4. Schematic of the PMA membrane and the fiber angle.

Figure 5 .
Figure 5. Photo of the test rig that the characteristic force maps have been measured with.

Table 2 .
Estimated parameters for F Martens ; optimization start point is given in Table1.

Table 3 .
Force error of each model by percentage.