Use of McKibben Muscle in a Haptic Interface †

Dynamic Modeling and Experimental Validation of a Haptic Finger Based


Introduction
Teleoperation consists of the control of a remote machine or device (slave), using an appropriate interface, called a master.The master is haptic if it is able to transmit the sense of touch, for example generating a force reflection to the operator.In general, haptics refers to the study of touch feedback in different applications, including virtual reality (such as for example in medical simulators), gaming, and tele-robots (such as surgical robots or remote manipulators for space exploration).
Among the different haptic devices that are able to generate haptic feedback on different sites of human body, the haptic glove is the most used, because hands play a dominant role in perception and manipulation tasks when grasping virtual objects [1].
When designing a haptic glove, the architecture of the mechanism that is able to transmit the force to the single fingers, as well as the type of sensors and the actuation, must be chosen.Generally, the problem has to be solved for a single finger, and replicated to all the fingers of a hand.
Starting from the performance specifications in terms of degrees of freedom, weight, size, dexterous capabilities, and safety level, multiple mechanisms for the transmission of the reflection force have been developed [2].In particular, in order to assure the coincidence between the center of the relative rotation of the links of the mechanism and the mean physiological rotational axes of the joints of the fingers, different solutions have been designed.Heo et al. [2] classified five mechanical architectures: serial linkages with rotational joints that are coincident with the physiological rotational axes, linkages with a remote center of rotation, redundant linkages [3], tendon-driven mechanisms, and serial linkages attached to distal segment.
As regards the sensing of the user's intended motion, the reflection force on the finger and the position of the joints should generally be measured.Heo et al. [2] described some different solutions that have been adopted in the literature.For the measurement of the contact force, force sensing resistors, pneumatic pressure sensors, and strain gauge sensors are frequently used.For the measurement of motion sensing, a bending sensor or a rotary encoder can be used.
To drive the mechanism that is devoted to transmitting the reflection force to the finger, different types of actuators have been used.Conventional electric motors have been efficaciously employed for the availability, reliability, and simplicity of control [2,4].Due to their intrinsic self-adaptability to the compliance of the biological tissues, non-conventional deformable actuators, such as shape memory alloy actuators [5], may be employed for this purpose.Pneumatic actuators provide some advantages over electric ones in terms of a high power-to-weight ratio, simplicity, safety, low cost, and easy maintenance.However, non-linearities due to air compressibility and actuator friction make it difficult to control pneumatic actuators [1].Between the pneumatic actuators, the soft actuators present different advantages, due to very high force/weight ratio and low friction [6].For all these reasons, among the pneumatic actuators, the McKibben pneumatic artificial muscle has been employed in several haptic devices, manipulators, and gloves [7][8][9][10].
The authors developed a prototype of a haptic finger, with a mechanism with one degree of freedom for force feedback actuated by McKibben's muscle.The device has been studied both as a rehabilitation system [11] and as a haptic system.In particular, a mathematical model of the device was developed, which was validated in a specific condition [12].
In the paper, the prototype of a haptic finger is presented, and the non-linear model developed is described.Therefore, further different validation tests of the model are presented.Finally, several simulations are presented, which are aimed at highlighting the dynamic behavior of the entire device.

The Prototype
Figure 1 shows the general scheme of the developed device.As regards the mechanism for the transmission of the reflection force on the fingertip, a one degree-of-freedom planar four-bar linkage has been chosen.The operator rests his hand on the fixed frame of the device, while his fingertip is held by a specially designed support integral with the beam force sensor mounted on the coupler of the four-bar mechanism (Figure 2).The operator, manipulating a virtual object or actuating a remote device, moves the finger between the position of maximum extension (Figure 2a) and that of maximum flexion (Figure 2b).For each position of the fingertip imposed by the operator, a corresponding angle α 1 of the rocker (Figure 1) is measured by a rotary encoder.The control unit, solving the forward kinematics equations of the mechanism, calculates the position of the fingertip (coordinates u, w in the sagittal plane, Figure 3a).Depending on the geometry and the stiffness of the remote or virtual object, the control unit evaluates the reflection force reference F ref that must be applied to the fingertip.The control unit also calculates the error e between the force reference F ref and the force F measured by the sensor, and provides the command signal for the pressure control proportional valve (FESTO ® MPPE-3-1/4-10-010B) V ref , which is used to regulate the upstream pressure p 1 of a fluidic resistance R. Regarding the actuation of the feedback force mechanism, a McKibben pneumatic muscle (Shadow Robot Company Ltd., London, UK) has been chosen.The pressure control proportional valve through the fluidic resistance R supplies the muscle and consequently generates a force F mu that, in parallel with a spring force F m , pulls a tendon that is connected to the coupler of the four-bar mechanism (Figures 1 and 2).The dimensional synthesis of the four-bar linkage was carried out by ensuring that the path generated by the coupler point P (fingertip holder) should be included into the natural workspace of a finger, which was calculated by taking into account the dependency among the rotations of each link of the finger (Figure 3a).The selected lengths of the proximal, middle, and distal phalanxes are respectively 4 cm, 2.5 cm, and 2 cm, while the rotation of the distal interphalangeal joint (DIP) joint has been considered to be equal to 2/3 of the proximal interphalangeal joint (PIP) joint rotation [11].
The dimensions of the links of the four-bar mechanism, whose scheme is shown in Figure 3b, are reported in Table 1.The operator's fingertip is positioned at point P, while the metacarpophalangeal joint (MCP) is coincident with the fixed hinge O. Tendon 5 is connected to coupler 2 in X1; it passes through the low-friction support point X and is pulled by the McKibben muscle.The dimensional synthesis of the four-bar linkage was carried out by ensuring that the path generated by the coupler point P (fingertip holder) should be included into the natural workspace of a finger, which was calculated by taking into account the dependency among the rotations of each link of the finger (Figure 3a).The selected lengths of the proximal, middle, and distal phalanxes are respectively 4 cm, 2.5 cm, and 2 cm, while the rotation of the distal interphalangeal joint (DIP) joint has been considered to be equal to 2/3 of the proximal interphalangeal joint (PIP) joint rotation [11].
The dimensions of the links of the four-bar mechanism, whose scheme is shown in Figure 3b, are reported in Table 1.The operator's fingertip is positioned at point P, while the metacarpophalangeal joint (MCP) is coincident with the fixed hinge O. Tendon 5 is connected to coupler 2 in X1; it passes through the low-friction support point X and is pulled by the McKibben muscle.The dimensional synthesis of the four-bar linkage was carried out by ensuring that the path generated by the coupler point P (fingertip holder) should be included into the natural workspace of a finger, which was calculated by taking into account the dependency among the rotations of each link of the finger (Figure 3a).The selected lengths of the proximal, middle, and distal phalanxes are respectively 4 cm, 2.5 cm, and 2 cm, while the rotation of the distal interphalangeal joint (DIP) joint has been considered to be equal to 2/3 of the proximal interphalangeal joint (PIP) joint rotation [11].
The dimensions of the links of the four-bar mechanism, whose scheme is shown in Figure 3b, are reported in Table 1.The operator's fingertip is positioned at point P, while the metacarpophalangeal joint (MCP) is coincident with the fixed hinge O. Tendon 5 is connected to coupler 2 in X 1 ; it passes through the low-friction support point X and is pulled by the McKibben muscle.

The Dynamic Model
The dynamic model of the haptic finger has been developed according to the logic represented in the block diagram of Figure 5.

The Dynamic Model
The dynamic model of the haptic finger has been developed according to the logic represented in the block diagram of Figure 5.The positioning in the sagittal plane (u, w) of each link of the mechanism can be calculated solving the non-linear system of Equations ( 1)-( 3), thus obtaining ϑ1, ϑ2, and β: sin sin cos tan cos cos sin Then, the trajectory of the fingertip in the (u, w) sagittal plane (Figure 3a) is calculated.
The stroke CAtt of the actuator depends on the length of the tendon l5 and the preload stroke of the spring, CAtt1, and is calculated by solving Equation (4): ( ) where l5max corresponds to the maximum length XX1 of the tendon (Figure 3b), which occurs at maximum finger flexion.The force F applied to the fingertip is calculated as a function of α1, imposing the equilibrium of the coupler 2 (Figure 6).The direction of F is considered to be perpendicular to the coupler 2 in point P, corresponding to the fingertip support, due to low friction and consequent sliding between the fingertip and the same support.The equivalent mass of the system mc has been considered as centered at point G belonging to the link 2. The static analysis yields the following Equations ( 5)-( 7):  The positioning in the sagittal plane (u, w) of each link of the mechanism can be calculated solving the non-linear system of Equations ( 1)-( 3), thus obtaining ϑ 1 , ϑ 2 , and β: Then, the trajectory of the fingertip in the (u, w) sagittal plane (Figure 3a) is calculated.
The stroke C Att of the actuator depends on the length of the tendon l 5 and the preload stroke of the spring, C Att1 , and is calculated by solving Equation (4): where l 5max corresponds to the maximum length XX 1 of the tendon (Figure 3b), which occurs at maximum finger flexion.The force F applied to the fingertip is calculated as a function of α 1 , imposing the equilibrium of the coupler 2 (Figure 6).The direction of F is considered to be perpendicular to the coupler 2 in point P, corresponding to the fingertip support, due to low friction and consequent sliding between the fingertip and the same support.The equivalent mass of the system m c has been considered as centered at point G belonging to the link 2. The static analysis yields the following Equations ( 5)-( 7): Once the reflection force F has been calculated, and considering that: (i) the force reference Fref is a function of the position (u, w) of the fingertip, depending on the geometry and the mechanical characteristic of the virtual/remote object, (ii) the fingertip position (u, w) is a function of the rocker rotation α1, it is possible to evaluate the force error e: ( ) For the control, a proportional-integral-derivative (PID) algorithm has been adopted: where KP, KI, and KD are, respectively, the proportional, integral, and derivative gains.After an iterative tuning approach, a simple proportional control came out to be sufficiently accurate and stable for the application, with a proportional coefficient of KP = 0.28.
A second-order transfer function has been used to model the dynamic behavior of the pressure control proportional valve (Pressure Control Valve block): where KV is the static gain, σn is the natural frequency, and ζ is the damping factor of the valve.The pneumatic resistance (Fluidic Resistance block) has been modeled according to the ISO 6358 [13].The mass air flow passing through the resistance can be calculated in sonic or subsonic condition, depending on the ratio between the downstream and upstream pressures: where P1 is the upstream absolute pressure, Pmu is the downstream absolute pressure, C is the sonic conductance, b is the critical ratio, and ρ0 = 1.18 kg/m 3 is the air density in normal conditions.The C and b parameters are reported in Table 1.
Since R = 287.2J/(kgK) is the air constant, and T is the absolute air temperature, assuming the variation of the internal volume V of the pneumatic muscle as negligible, the time derivative of the absolute internal pressure of the muscle for an isothermal transformation can be expressed as: Once the reflection force F has been calculated, and considering that: (i) the force reference F ref is a function of the position (u, w) of the fingertip, depending on the geometry and the mechanical characteristic of the virtual/remote object, (ii) the fingertip position (u, w) is a function of the rocker rotation α 1 , it is possible to evaluate the force error e: For the control, a proportional-integral-derivative (PID) algorithm has been adopted: where K P , K I , and K D are, respectively, the proportional, integral, and derivative gains.After an iterative tuning approach, a simple proportional control came out to be sufficiently accurate and stable for the application, with a proportional coefficient of K P = 0.28.A second-order transfer function has been used to model the dynamic behavior of the pressure control proportional valve (Pressure Control Valve block): where K V is the static gain, σ n is the natural frequency, and ζ is the damping factor of the valve.The pneumatic resistance (Fluidic Resistance block) has been modeled according to the ISO 6358 [13].The mass air flow passing through the resistance can be calculated in sonic or subsonic condition, depending on the ratio between the downstream and upstream pressures: where P 1 is the upstream absolute pressure, P mu is the downstream absolute pressure, C is the sonic conductance, b is the critical ratio, and ρ 0 = 1.18 kg/m 3 is the air density in normal conditions.The C and b parameters are reported in Table 1.
Since R = 287.2J/(kgK) is the air constant, and T is the absolute air temperature, assuming the variation of the internal volume V of the pneumatic muscle as negligible, the time derivative of the absolute internal pressure of the muscle for an isothermal transformation can be expressed as: with the actual length l and the radius r of the actuator given by: where, for the employed McKibben muscle (with a 20 mm diameter), C att is the stroke of the pneumatic muscle, l 0 is the initial preloaded length, and n t is the number of turns of the fibers, whose length h is 178 mm.The internal relative pressure of the muscle p mu can be calculated integrating Equation (12).In order to estimate the force exerted by the pneumatic actuator as a function of the internal pressure p mu , a model of the McKibben muscle must be implemented.Various approaches have been proposed in the literature [14].Ideal simple models consider only the kinematic relationship between the braided sheath and the inner tube, but unfortunately, they suffer when comparing theoretical and experimental results [15].To overcome this issue, different models that account for the elastic energy that is stored in the bladder were developed based on the principle of virtual work [16][17][18].Another approach to derive an analytical expression of the actuation force consists of writing balance equations of a free body diagram of the muscle [15].Among these types of models, Ferraresi et al. [19] developed and validated a McKibben model, which is able to take into account the effects of the thickness and elasticity of the inner tube.Due to a good compromise between the simplicity and accurate numerical results, the latter model was used to calculate the force exerted by the muscle as a function of the internal pressure of the actuator [19]: where r i and l i are the initial radius and length of the muscle at rest, respectively, s i is the initial thickness of the inner chamber, and E is the Young modulus of the chamber.Finally, the force F x is directly connected to the force exerted by the muscle F mu and the spring F m by the following relations: where k m is the spring constant, and F m0 is its preload.

Experimental Validation of the Model
In order to verify the capability of the device to generate correct force feedback on the fingertip and validate the model, several experimental tests have been performed on the prototype presented in Section 2.
The angular rotation of the rocker 4 (Figure 3b) has been measured by a rotational encoder (BDK series, 1024 pulses, Baumer electric, Frauenfeld, Switzerland), which was processed and acquired by the incremental encoder interface board dSPACE ® DS3002.The measurement of the force applied to the fingertip was made through a planar beam resistive sensor (Futek FR1020, capacity 89 N, combined error 0.25% rated output, Irvine, CA, USA) conditioned by a full bridge strain gauge module Meco, model MecoStrain.The acquisition of the force sensor signal was provided by a dSPACE ® Multi-Channel A/D Board DS2002.The force sensor is mounted on the coupler by special grasping (Figure 7a).The fingertip support has been made by low-friction plastic material, as shown in Figure 7b, to ensure that the exchanged force can be considered perpendicular to the sensor itself.
The control algorithm was implemented in MATLAB-Simulink ® environment, dSpace Control Desk.The command signal of the pressure control proportional valve V ref was generated by a dSPACE ® D/A board DS2101.In order to evaluate the dynamic performance of the haptic interface, several experimental tests were conducted.An operator, after positioning his hand on the device with a fully extended index finger, as shown in Figure 4b, is invited to freely perform flexions and extensions of the index finger (Figure 8b), as if he were virtually manipulating an object.The virtual object that was chosen for the experimental tests has a constant stiffness Kobj, so that the feedback force is proportional to the rotation of the rocker α1 according to the equation Fref = Kobj Δα1.During the manipulation, both the angular rotation of the rocker α1, which is freely imposed by the operator, and the reflection force generated on the fingertip by the haptic device were acquired.In all of the tests, a purely proportional control (KP = 0.28 V/N) was implemented, while the stiffness of the virtual object Kobj has been changed in the different tests, and fixed respectively equal to 8 N/rad, 10 N/rad, 20 N/rad, 30 N/rad, 40 N/rad, and 50 N/rad.Then, each test is characterized both by different temporal rotations of the rocker α1=α1(t) imposed by the operator, and by different stiffness Kobj assigned to the virtual object.
In order to validate the model, simulations were performed by imposing, as an input, the same rotation law of rocker α1 = α1(t) recorded in each experimental test, and calculating the reflection force generated by the system with the same control (purely proportional control with KP = 0.28 V/N), and with same stiffness of the virtual object Kobj.
Figure 9 shows the experimental results obtained in the different tests described above, compared with the results of the simulation conducted under the same conditions.The figures show the trends of the experimental and simulated feedback force as a function of time, for different stiffnesses of the object, between 8 N/rad and 50 N/rad.Figure 8a shows an image of the entire experimental setup during the static experimental tests.In order to evaluate the dynamic performance of the haptic interface, several experimental tests were conducted.An operator, after positioning his hand on the device with a fully extended index finger, as shown in Figure 4b, is invited to freely perform flexions and extensions of the index finger (Figure 8b), as if he were virtually manipulating an object.The virtual object that was chosen for the experimental tests has a constant stiffness K obj , so that the feedback force is proportional to the rotation of the rocker α 1 according to the equation F ref = K obj ∆α 1 .During the manipulation, both the angular rotation of the rocker α 1 , which is freely imposed by the operator, and the reflection force generated on the fingertip by the haptic device were acquired.In all of the tests, a purely proportional control (K P = 0.28 V/N) was implemented, while the stiffness of the virtual object K obj has been changed in the different tests, and fixed respectively equal to 8 N/rad, 10 N/rad, 20 N/rad, 30 N/rad, 40 N/rad, and 50 N/rad.Then, each test is characterized both by different temporal rotations of the rocker α 1 = α 1 (t) imposed by the operator, and by different stiffness K obj assigned to the virtual object.In order to evaluate the dynamic performance of the haptic interface, several experimental tests were conducted.An operator, after positioning his hand on the device with a fully extended index finger, as shown in Figure 4b, is invited to freely perform flexions and extensions of the index finger (Figure 8b), as if he were virtually manipulating an object.The virtual object that was chosen for the experimental tests has a constant stiffness Kobj, so that the feedback force is proportional to the rotation of the rocker α1 according to the equation Fref = Kobj Δα1.During the manipulation, both the angular rotation of the rocker α1, which is freely imposed by the operator, and the reflection force generated on the fingertip by the haptic device were acquired.In all of the tests, a purely proportional control (KP = 0.28 V/N) was implemented, while the stiffness of the virtual object Kobj has been changed in the different tests, and fixed respectively equal to 8 N/rad, 10 N/rad, 20 N/rad, 30 N/rad, 40 N/rad, and 50 N/rad.Then, each test is characterized both by different temporal rotations of the rocker α1=α1(t) imposed by the operator, and by different stiffness Kobj assigned to the virtual object.In order to validate the model, simulations were performed by imposing, as an input, the same rotation law of rocker α 1 = α 1 (t) recorded in each experimental test, and calculating the reflection force generated by the system with the same control (purely proportional control with K P = 0.28 V/N), and with same stiffness of the virtual object K obj .
Figure 9 shows the experimental results obtained in the different tests described above, compared with the results of the simulation conducted under the same conditions.The figures show the trends of the experimental and simulated feedback force as a function of time, for different stiffnesses of the object, between 8 N/rad and 50 N/rad.
The most relevant differences between the simulated and the experimental results occurred for more rigid virtual objects, in cases where the operator imposed low-rocker rotation values (low reaction forces) and quickly reversed the direction of rotation.
Nevertheless, the model highlights the ability to predict with good accuracy the trend of the force generated by the device in dynamic conditions.
Based on these results, the model was considered validated and sufficiently reliable to predict the dynamic performance of the device, even under conditions different from those tested.This analysis is presented and discussed in the next section.The most relevant differences between the simulated and the experimental results occurred for more rigid virtual objects, in cases where the operator imposed low-rocker rotation values (low reaction forces) and quickly reversed the direction of rotation.
Nevertheless, the model highlights the ability to predict with good accuracy the trend of the force generated by the device in dynamic conditions.
Based on these results, the model was considered validated and sufficiently reliable to predict the dynamic performance of the device, even under conditions different from those tested.This analysis is presented and discussed in the next section.

Dynamic Assessment of the Device
The validated model was used to evaluate the dynamic performance of the device.A first series of simulations concerned the study of the response to step input, for different proportional gain K P values of the pure proportional control (Figure 10).Due to pure proportional control, a steady state error is observed.As expected, when increasing the proportional gain K P , the steady-state error decreases, but for proportional gain values close to 2 V/N, the system tends to instability, and in fact becomes unstable for K P values equal to 3 V/N.The validated model was used to evaluate the dynamic performance of the device.A first series of simulations concerned the study of the response to step input, for different proportional gain KP values of the pure proportional control (Figure 10).Due to pure proportional control, a steady state error is observed.As expected, when increasing the proportional gain KP, the steady-state error decreases, but for proportional gain values close to 2 V/N, the system tends to instability, and in fact becomes unstable for KP values equal to 3 V/N.

Figure 1 .
Figure 1.Scheme of the device.

Figure 1 .
Figure 1.Scheme of the device.

Figure 2 .
Figure 2. Mechanism for the transmission of the reflection force to the fingertip: (a) Finger at maximum extension; (b) finger at maximum flexion.

Figure 2 .Figure 3 .
Figure 2. Mechanism for the transmission of the reflection force to the fingertip: (a) Finger at maximum extension; (b) finger at maximum flexion.Robotics 2019, 8, x FOR PEER REVIEW 4 of 12

Figure 4
Figure4shows the prototype of the haptic finger that was realized and used for experimental tests.The labels refer to the Figures2 and 3.The frame and the links of the mechanism were made of steel.Plain bearings were used to realize the revolute joints O, O 1 , A, and B.

Figure 4
Figure 4 shows the prototype of the haptic finger that was realized and used for experimental tests.The labels refer to the figures 2 and 3.The frame and the links of the mechanism were made of steel.Plain bearings were used to realize the revolute joints O, O1, A, and B.

Figure 4 .
Figure 4.The prototype of the device.(a) Detail of the four-bar linkage and of the sensors.(b) Positioning of the operator's finger.

Figure 4 .
Figure 4.The prototype of the device.(a) Detail of the four-bar linkage and of the sensors.(b) Positioning of the operator's finger.

Figure 5 .
Figure 5. Block diagram of the dynamic model of the haptic finger.Being known the rotation angle of the rocker α1, depending on the position of the operator finger, and the force Fx generated by the McKibben muscle and by a parallel spring, the Device Dynamics block calculates the feedback force F that is applied to the fingertip and the stroke of the McKibben muscle CAtt.The positioning in the sagittal plane (u, w) of each link of the mechanism can be calculated

Figure 5 .
Figure 5. Block diagram of the dynamic model of the haptic finger.Being known the rotation angle of the rocker α 1 , depending on the position of the operator finger, and the force F x generated by the McKibben muscle and by a parallel spring, the Device Dynamics block calculates the feedback force F that is applied to the fingertip and the stroke of the McKibben muscle C Att .The positioning in the sagittal plane (u, w) of each link of the mechanism can be calculated solving the non-linear system of Equations (1)-(3), thus obtaining ϑ 1 , ϑ 2 , and β:

Figure 6 .
Figure 6.The free body diagram of the coupler (represented in the thick blue line).

Figure 6 .
Figure 6.The free body diagram of the coupler (represented in the thick blue line).

Figure 7 .Figure 8 .
Figure 7. Force sensor.(a) The sensor and the mounting gripper.(b) The sensor in static calibration operation.

Figure 7 .
Figure 7. Force sensor.(a) The sensor and the mounting gripper.(b) The sensor in static calibration operation.

Figure 7 .Figure 8 .
Figure 7. Force sensor.(a) The sensor and the mounting gripper.(b) The sensor in static calibration operation.

Table 1 .
Parameters of the model.

Table 1 .
Parameters of the model.

Table 1 .
Parameters of the model.