A New Model Based on Adaptation of the External Loop to Compensate the Hysteresis of Tactile Sensors

This paper presents a novel method to compensate for hysteresis nonlinearities observed in the response of a tactile sensor. The External Loop Adaptation Method (ELAM) performs a piecewise linear mapping of the experimentally measured external curves of the hysteresis loop to obtain all possible internal cycles. The optimal division of the input interval where the curve is approximated is provided by the error minimization algorithm. This process is carried out off line and provides parameters to compute the split point in real time. A different linear transformation is then performed at the left and right of this point and a more precise fitting is achieved. The models obtained with the ELAM method are compared with those obtained from three other approaches. The results show that the ELAM method achieves a more accurate fitting. Moreover, the involved mathematical operations are simpler and therefore easier to implement in devices such as Field Programmable Gate Array (FPGAs) for real time applications. Furthermore, the method needs to identify fewer parameters and requires no previous selection process of operators or functions. Finally, the method can be applied to other sensors or actuators with complex hysteresis loop shapes.


Introduction
Tactile sensors are arrays of force sensing units or taxels. They are used in robotics to detect contact with objects [1][2][3], for instance in handling applications to improve the dexterity of artificial robotic hands [4,5]. Many different approaches have been proposed to make these sensors. Most of them are based on piezoresistive [6][7][8] or capacitive principles [9][10][11], although other transduction principles such as optical [12] or piezoelectric [13] have also been exploited. However, all these sensors show common errors in its operation such as hysteresis, nonlinearity and mismatching [14,15], so it is necessary to compensate these errors to obtain a more precise response. In addition, this compensation has to be carried out in highly demanding real time tasks, so smart tactile sensors with local electronics powerful enough to process the large amount of data from the sensor array in real time are required. The authors have proposed smart tactile sensors with local electronics based on FPGAs [16,17]. These devices consist of logic blocks and dedicated modules that allow the implementation of complex arithmetical operations, and their main advantage is the parallel execution of processes [18]. Therefore, this approach seems suitable to implement the compensation algorithms in real time control tasks. However, the complexity of the involved mathematical operations affect the speed and consumption of power and system resources, so another significant goal is to reduce both as much as possible.
Virtually all sensors and actuators based on smart materials present undesired complex hysteretic nonlinearities when driven with sufficiently high amplitudes. To compensate for these nonlinearities, as in piezoelectric actuators, many efforts have been made including the feed-forward control as the most common approach [19,20]. The main idea of this compensation is to develop a mathematical model of the hysteresis that can be inverted; the inverted model can then be connected in cascade before the actuator input in order to obtain a linearized response. Likewise, in the case of using a tactile sensor, the inverted model can be placed after the sensor output to compensate the hysteretic behavior and obtain a linearized output. In the case of working with a sensor array, it is also possible to reduce the mismatching between different taxels, since all their outputs are equilibrated at the same time as the hysteresis is compensated. Therefore, the goal is to obtain an accurate model of hysteresis of the sensor and compensate for all these errors.
The hysteresis modeling methods commonly used are phenomenological. Other physics-based models require complex differential equations which need high computing power and a long resolution time, so they are more difficult to implement in the local electronics of smart sensors [21]. However, the phenomenological methods are built from experimental data without considering the physical properties of the actuators or sensors [22]. The most common method to compensate for hysteresis in actuators and sensors is the Preisach model [23,24]. This model is difficult to implement due to the large amount of data required to achieve a good approximation. The Prandtl-Ishlinskii model (PI) [25,26] is a subclass of the Preisach model that has become popular because, unlike the Preisach, its inverse can be calculated analytically and its implementation is much simpler. The classical model of Prandtl-Ishlinskii (CPI) approaches the sensor response with a weighted sum of hysteresis play operators [27] but it can only be used to model symmetrical hysteresis curves. In order to model asymmetrical hysteresis loops, different alternatives have been proposed such as the generalized model of Prandtl-Ishlinskii (GPI) [28], which is based on the use of two different envelope functions to combine with play operators. As a special case of this GPI model, a Prandtl-Ishlinskii model that combines two different asymmetrical play operators is proposed in [29] to independently characterize the ascending and descending branches of a piezoelectric actuator hysteresis loop. Another alternative for the modeling and compensation of the asymmetrical hysteresis nonlinearities is a modified Prandtl-Ishlinskii model (MPI) described in [30], which replaces the linear input function of the classical play operator by a generalized input function based on a third-order polynomial. We can also find hysteresis compensation methods based on the use of two dominant continuous functions, one ascending and another descending, which converge to a turning point without memory saturation [31,32]. The dominant functions here are built from high order polynomials, and the whole model is built from these functions through a nonlinear transformation of the coordinate axis (herein we call it POLY model).
This paper presents a new method to model and compensate hysteresis nonlinearities based on the adaptation of the external loop (herein we call it ELAM). It is a phenomenological model, which builds two continuous and monotonic curves, one increasing and the other decreasing, from linear interpolation of the experimentally measured hysteresis external loop. From these two curves, all internal hysteresis loops are approximated using a different procedure depending on whether it is in an increasing or decreasing branch. Moreover, the approximation is made in two intervals defined by a so called split point, with different adaptations of the external loop in each interval. In order to evaluate the effectiveness of the method, the accuracy of the proposed model is compared with that achieved by other methods referred to above. The methods used for this purpose are the generalized model of Prandtl-Ishlinskii (GPI), the modified model of Prandtl-Ishlinskii (MPI) and that based on polynomial fitting [28,30,31].
A discussion about the number of parameters to be estimated in each model, the computational complexity, and the achieved accuracy is undertaken. The ELAM model is shown as the most accurate method, so its inverse is calculated and the hysteresis nonlinearities of the complete tactile sensor array are compensated. The results confirm that it is a very efficient method to be implemented in real-time control systems using smart tactile sensors. The ELAM method is also applicable to all types of hysteresis loops obtained from other sensors or actuators, and it can provide more accurate models than other methods when the hysteresis loops show complex shapes.
The remainder of this paper is organized as follows. Section 2 shows the tactile sensor and the set-up used to obtain the experimental data. Section 3 explains the different modeling hysteresis methods used to compare with that proposed. Section 4 introduces the so called ELAM method. Section 5 deals with the parameter identification of the models. Section 6 shows the results and related discussions. Finally, some concluding remarks are provided in Section 7.

PCB Based Sensor
The tactile sensor used in this article consists of a set of electrodes and addressing tracks fabricated on a flexible printed circuit board (PCB). Atop of these electrodes, a thin film of conductive polymer, such as piezoresistive material is placed. Specifically, a conductive water-based ink of this polymer is deposited by spin-coating technique on a flexible plastic sheet, obtaining a smooth, homogeneous and conductive thin film [33,34]. The most interesting thing about this process is that it is cheap and allows the manufacturing of flexible and low cost tactile sensors. The sensor consists of 16 × 16 taxels and a spatial resolution of 2.54 mm. Figure 1 shows a section view and a top view of a taxel, besides a picture of the complete array of the tactile sensor. The resistance between two electrodes associated to each taxel changes when the exerted force against the taxel changes. The readings of the whole tactile array are registered by means of well-known interface electronics [33] designed to achieve a good static performance, so electro-mechanic relays are used to implement the switches to select the rows as the array is scanned. The data acquisition is achieved with the NI-USB 6259 device from National Instruments (National Instruments Corporation, Austin, TX, USA). Sixteen analog inputs are multiplexed to scan up to 16 × 16 taxels in our testing platform. In order to measure the hysteresis of the tactile sensor based on a PCB, a sequence of uniform pressures on the entire matrix of the tactile sensor is exerted. Six consecutive loading-unloading cycles with different points of return for the same ascending curve (see Figure 2a), and with different starting points to rise from the same descending curve (see Figure 2b), are applied. Thus, the descending and ascending behavior of the sensor are respectively characterized. The cycles are performed with an increase of 2 psi between pressures, and in Table 1 the pressure sequences are shown. These hysteresis curves represent the average of the output produced by all of the tactile sensor taxels after repeating each cycle five times. The time interval between the new pressure level being exerted and voltage output being registered by the acquisition board is 2 s. The platform of pneumatic calibration Tekscan PB100E (Teckscan Inc., South Boston, MA, USA) ( Figure 3a) is used to apply the uniform pressures over the tactile sensor matrix. In order to quantify the hysteresis exhibited by the sensor, the hysteresis error as the difference in sensor output voltage to the same applied pressure is identified when these pressures are exerted on the ascending and descending branches of the cycles. The maximum and average hysteresis errors are referenced to the highest output value to obtain a percentage of the error relative to full scale. The maximum error due to hysteresis is 25.3% of full scale output, while the average error is 10.0% of full scale out. In Figure 2c, the frame obtained at a pressure of 40 PSI is shown as an example to illustrate the mismatching of the output obtained with the taxels of the tactile sensor.

Measurement Setup
In order to analyze the behavior of the tactile sensor output two different measurement platforms are used. The first is based on a pneumatic commercial equilibration/calibration device (Tekscan PB100E Teckscan Inc., South Boston, MA, USA [35]) (see Figure 3a) to obtain readings of the whole tactile matrix under the same uniform pressure. The sensor is placed in a slot of a chamber where one side is rigid and the other is a flexible wall. When the chamber is pressurized the wall exerts an even pressure on the sensor. An electro-valve Pneumax 171E2N.T.D.0009 (Pneumax S.p.a., Milano, Italy) [36], which allows the flow from a compressor, is added to set the desired pressure selected in a computer software.
The second platform is used to register the sensor response to pressure exerted by an object with a known shape on the sensor surface. A piece of fabric between the object and the sensor was added to improve the contact. With the purpose of applying a force on the object, a motorized platform is used, which is composed of a T-NA08A50 linear actuator ( Zaber, Vancouver, BC, Canada) and two T-LA60A actuators of Zaber Technologies (Zaber, Vancouver, BC, Canada) [37] (see Figure 3b). The T-NA08A50 actuator (Zaber, Vancouver, BC, Canada) allows a force to be exerted on the z-axis, while the T-LA60A actuators allow movement of the platform along the x and y axes. One sensor Nano17 from ATI Industrial Automation (ATI Industrial Automation, Apex, NC, USA) [38] was added at the tip of the motor in the vertical direction to register the actual force exerted on the objects and then on the sensor.  Figure 4 illustrates the procedure to compensate the hysteresis followed by the four methods that are compared in this paper. The pressure exerted against a taxel is identified by p(t), and v(t) is its output voltage. H(p(t)) represents the output with the corresponding measured hysteresis at the sensor. The goal is to find a model Hm(p(t)) which fits the experimental data H(p(t)) with the greatest possible accuracy, low computational cost, and which can be inverted to obtain Hm −1 (v(t)). Once this model is inverted, it is possible to obtain pm(t) which ideally is equal to p(t). Each of the four methods to compare uses a different mathematical expression of Hm(p(t)) and employs a number of parameters which must be identified by an error minimization method (see Section 5). The parameter identification process will identify the set of parameters that configure each of the hysteresis models Hm(p(t)) which are compared in this study. Due to the fact that they are phenomenological models, a set of M samples of measured output values of the tactile sensor {y(t1),…, y(tj),…y(tM)} are used to computationally derive the vector of parameters X minimizing the following error function [39].

Hysteresis Models
where ym are the samples of the output of each hysteresis model Hm(p(t)) with 1 jM  . In the next sections the models for Hm(p(t)) are called HGPI(p(t)), HMPI(p(t)), HPOLY(p(t)) and HELAM(p(t)), and their output values are named  

Generalized Prandtl-Islinskii Model (GPI)
The generalized model of Prandtl-Ishlinskii (GPI) [28] provides a model for H(p(t)) in the scheme of Figure 4, which we call HGPI(p(t)). This model is derived from a weighted superposition of a set of generalized play operators with different threshold values r with . Since the hysteresis curve of the sensor (see Figure 2a) is asymmetrical, it is necessary to use these generalized operators (see Figure 5a) instead of classical play operators (see Figure 5b). Note that an increase of input corresponds to an increase of output along the curve l  , while a decrease of input corresponds to a decrease of output along the curve r  . Both curves can be different, so that asymmetrical loops have to be modeled. The only conditions are that these so called envelope functions, l  and r  , be monotone and continuous functions with lr rr     . If a finite number n of operators is used [28], the output of the generalized Prandtl-Ishlinskii model can be approximated, by the expression.
where p(r) is the density function and H is the generalized play operator in Figure 5a defined from the envelope functions l  and r  as with y as the output value of the generalized play operator for the previous value of x.

26177
Therefore, the generalized Prandtl-Ishlinskii model will be obtained from the set of parameters X that define the envelope functions l  and r  , the density function p(r) and the threshold values of the play operators r. The GPI model output   GPI yt is used in Equation (1) to obtain the parameters X (see Section 5).

Modified Prandtl-Ishlinskii Model (MPI)
The modified Prandtl-Ishlinskii model (MPI) uses the classical play operators for modeling asymmetrical hysteresis such as the one presented by our tactile sensor (see Figure 2a). Although the classical Prandtl-Ishlinskii model (CPI) can only model symmetrical hysteresis curves, this MPI model proposes replacing the linear input function of the CPI model by a generalized input function. Thus, the asymmetric hysteresis nonlinearities can be determined, not only by the weighted superposition of a set of classical play operators but also by the generalized input function [30]. The main advantages of this model relative to the GPI are, as it continues using the classical play operators, that the mathematical description is simpler, the number of parameters to be identified is smaller and its inverse can be calculated analytically from the inverse of CPI model.
The classical Prandtl-Ishlinskii model (CPI) [27] is based on the combination of classical play operators (see Figure 5b) with different thresholds r. The modification of the play operator by a one-side operator (OSP) is proposed in [30] when the sensor or actuator only works with positive excitation signals. Then, the play operator (OSP) (see Figure 5c) is expressed as Using a finite number n of OSP operators the output of the MPI model is is a density function and the function g(x(t)) replaces the linear input function 0 () p x t  in the CPI model to achieve the modeling of asymmetric hysteresis [30].

Polynomial Based Model (POLY)
Methods have been proposed to compensate the hysteresis using mathematical structures that are not built with play operators. This is the case of the model described in [31], which is based on the construction of the model of hysteresis of a piezoelectric actuator from the external curves of the hysteresis data. This external loop consists of an ascending dominant curve when the input values increase, and a descending dominant curve when the input values decrease. All ascending curves converge at one point called upper target point, while all descending curves converge on the same lower converging point (see Figure 5d). The rest of the hysteresis curves can adopt their shape from these dominant curves, which can be expressed as two monotonically continuous functions, fra(x) and frd(x), respectively. Third-order polynomials to implement these dominant functions are proposed in [31], so in this paper, we call this model POLY.
The ascending curve function for any ascending trajectory starting from point 11 ( , ) xy and ending at xy is the lower converging point and ( , ) uu xy is the upper converging point of the dominant curves (see Figure 5d). Similarly, the descending curve for any trajectory starting point 11 ( , ) xy and ending at point 22 where The output of the model can be expressed as

External Loop Adaptation Model (ELAM)
The analysis of the GPI, MPI and POLY methods exposed in the latter sections, shows the presence of play operators, exponential functions, logarithmic, hyperbolic tangent, exponentiation and high degree polynomials, which anticipate a complicated implementation in devices such as FPGAs. Therefore, development of a method based on simple mathematical operations is essential to achieve fast and efficient real time tactile sensors. Furthermore, if one considers the possibility of working with arrays with a large number of taxels (for example, 16 × 16 = 256), each of which must be compensated in real time with its own model of hysteresis, depending on the application, compensation should be performed as fast as possible or with minimal resource consumption and high accuracy.

Direct External Loop Adaptation Model
The newly proposed approach is based on an adaptation of the external loop of the hysteresis curves to build all the inner hysteretic cycles. Linear interpolation of experimental data from this external loop (see Figure 2) is done to obtain the so called pattern curves, one for the ascending external trajectory and another for the descending external trajectory, so we call them a P and d P respectively. The adaptation is made by piecewise linear mapping of the pattern curves into the input interval of the internal target curves. Figure 6 illustrates the linear mapping of a generic pattern curve P between ( , ) iT fT xx , the corresponding value P x in [ , ] iP fP xx obtained by linear mapping is: iP fP yy is mapped into T y in [ , ] iT fT yy as y y y y y T x Y y y y y y y yy A key aspect of the ELAM model is that the target curve is split into pieces and the mapping in Equations (9) and (10) is done in each piece. It has been split into two segments in this work that are defined by the split point ( , ) ss xy with iT s fT x x x  and iT s fT y y y  . As a first simple approach, we propose the following linear expressions to find the split point: where 0 1   and 0 1   . The location of this point is determined by an error minimization algorithm as explained below, and the parameters  and  are the same for all target curves. This strategy allows different mappings to be performed at the left and right of the split point to achieve a more accurate fitting of the curve.
The coordinates of the split point ( , ) sd sd xy can be expressed from Equations (11) and (12) (13), the linear mapping in Figure 6 is carried out by replacing the initial and final values of the pattern curve P and of the target curve T in the Equations (9) and (10) (see Figure 6), by the initial and final values corresponding to the pattern curve and the target curve used in the construction of each of the segments. To do this, another split point has to be established in the pattern curve. From inspection of the experimental curves, it is set to ( , ( ))  (9) and (10) and  (9) and (10) to obtain (see Figure 7a) and Note that the operation to the left of the split point is simply the scaling of the pattern curve.
The output values in Equations (17) and (19) (14) and (15). An error minimization algorithm can fit the experimental data of the internal cycles in the values provided by Equations (17) and (19) (14) and (15) for any other possible descending curve in an internal cycle and the whole descending trajectory is given by Equations (17) and (19).
In this case, the adaptation of the curve a P from the external loop to the desired trajectory a T is performed. The split point ( , )  Figure 7b illustrates. Therefore, both split points do not now share the same x coordinate as in Figure 7a, but the same y coordinate.
For sa Tr u x x x  , the pattern curve a P is mapped into the target curve _ ar T (see Figure 7b) from the Equations (9) and (10) respectively, so we obtain  x x x , the pattern curve a P is mapped into the target curve _ al T (see Figure 7b) from the Equations (9) and (10) (24) and (26) (21) and (22).
Again, the error minimization algorithm can fit the experimental data of the internal cycles in the values provided by Equations (24) and (26) (21) and (22) for any other possible ascending curve in an internal cycle and the whole ascending trajectory is given by Equations (24) (17), (19), (24) and (26), respectively, and replacing ym(t) by yELAM(t) in Equation (1), the HELAM(p(t)) model for Figure 4 is:

Inverse External Loop Adaptation Model
The purpose of building a model of the hysteresis of the sensor is that the inverse model can be applied in cascade with the sensor output to compensate its hysteretic behavior. It is, therefore, necessary to invert the ELAM model proposed in this paper in such a way that for any output value of the tactile sensor () out T y T x  we can calculate the input value T x . The inverse model construction is performed similarly to the direct model described in the previous section. The following procedure is actually the same that is described in the previous section but with a swap of the coordinate axes. Since the split point is the same, the inverted curves are obtained just inverting the mapping in each piece. The inverted mapping is readily performed from Equations (9) and (10) The Equations (29) and (30) define the linear mapping for the inverse model ELAM, so the equations of the ascending or descending curves in the inverse model can be obtained from them in a similar way as that followed to build the direct model.

Parameter Identification
The four methods explained in the previous section are used to model hysteresis loops (see Figure 2a,b) of the tactile sensor. As a step prior to the construction, the validation and the comparison of the proposed models, it is necessary to carry out the identification of the parameters to adapt the models to the experimental data as accurately as possible. Many identification algorithms have been proposed for this purpose such as the least squares method, genetic algorithms and the particle swarm optimization method [40][41][42]. In this paper, we use genetic algorithms because they implement a parallel procedure able to simultaneously explore a wide range of solutions using probabilistic operators [43]. This feature allows them to discard local minima that do not correspond to an optimal solution. In the following, the different methods to obtain the parameters for the results in Section 6 are described.

GPI
The GPI model obtained from Equation (2) is completely defined by the number of generalized play operators used, the thresholds of the operators, the density function and the envelope functions.
Regarding of the number of operators n, from a theoretical point of view, the selection of a larger number of operators should obtain a more accurate approximation of the hysteresis loops. However, in real applications, it is found that further increase of the number of operators does not improve the fitting accuracy significantly. Since the complexity of the model is increased with the number of operators employed, it is advisable to use the smallest number of operators. According to the experiments reported in [44] for the same tactile sensor of this paper, it is enough to use a number of n = 4 play operators to obtain a good approximation.
For the results of this paper, the thresholds of the operators and the density function are obtained from [28] and their equations are i ri   (40) where i = 0,1,2,…,n and α is a positive real constant, and where ρ > 0 and τ are real constants.
Note that the density function p(r) vanishes for high values of the thresholds r and that there is no general criterion for its selection. Generally, it is completely selected by the designer. Once the structure of the density function p(r) is fixed, the parameters involved in the density function shall be determined by identification from experimental data.
The choice of the envelope functions used in the Equation (3) is decisive for a good fitting. Based on the results reported by other studies [44], which analyze and compare the use of different envelope functions, these functions have been chosen according to the following expressions in this paper:

MPI
The MPI model is also based on play operators, but of classical type in this case. Therefore, as with the model GPI, it will be necessary to select a number of operators n to use, thresholds for such operators, and a density function. The selection criteria for these items in the GPI model are also valid for the MPI model. In this article, we chose these elements according to what is proposed in [30], so a number of 10 classical play operators is selected. The expressions of the threshold operators ri and the density function p(r) are with i = 1,2,…,n and ( ) 1 vt   in the normalized case, and where ρ > 0 and τ are real constants.
To achieve the asymmetry observed in the hysteresis loop of the tactile sensor, it is necessary to define the generalized function g(x(t)) of the Equation (5). According to the proposal option in [30], this generalized function is chosen as a third-order polynomial 32 The selection of the generalized input function is not unique and other forms can also be chosen, but the use of a third-order polynomial seems a good choice to approximate curves of different shapes. Furthermore, third-degree polynomials are generally recognized as an effective choice to describe the hysteresis loops [45,46]. Actually, this polynomial is the one that achieves the best fitting to the hysteresis curve of the tactile sensor in this paper.

POLY
Regarding the POLY model, the only necessary parameter identification corresponds to the choice of coefficients of the polynomials fra and frd. Although third-order polynomials are proposed in [31] to be used in the model, sixth-order polynomials have been selected for the results of this paper. Therefore, the polynomials are 6 6 0 () Better results were achieved with tenth-order polynomials, but the choice of sixth-order polynomials seems most appropriate taking into account the balance between precision and the number of parameters required.

ELAM
The identification process for the ELAM model consists of determining the αd0, βd0, αa1, αa0 and βa0 parameters that allow the split points to be obtained in the Equations (14), (15), (21) and (22), respectively. The number of parameters to be identified is much lower than in the other models, so this identification process has a lower computational cost and the time employed is significantly shorter in this model.
It is remarkable that the GPI, MPI and POLY models require a prior step to choose the suitable envelope functions or polynomials to be used in the construction of the hysteresis loops. Moreover, in the case of models based on the Prandtl-Ishlinskii method, it is necessary to select the number of play operators and the structure of the density function. All this work is performed by the designer based on the results obtained in a selection process in order to find the best fit to the experimental data. Depending on the complexity of the hysteresis curve of the sensor, this process may be more or less costly. The proposed ELAM model does not require this previous step, because the pattern curves are obtained by linear interpolation of experimental data. In addition, the model provides the output from piecewise linear mapping of the pattern curves. Therefore, if the number of pieces to perform linear mapping or interpolation grows, the error always decreases and the model is robust against overfitting [47].
The parameters identified for each model that allow the best adaptation of the models to hysteresis data obtained experimentally for the descending and ascending curves are shown in Tables 2 and 3, respectively (see Figure 2a,b).

Results for the Output of a Single Average Taxel
The hysteresis loop in Figure 2a,b is really challenging to model because it is clearly asymmetrical, with high nonlinearity, and quite different ascending and descending external curves. Furthermore, the ascending curves inside the hysteresis loop do not have a similar shape to the external ascending curves, as can be seen in Figure 2b. Figures 9a-d and 10a-d show the models of the descending and ascending hysteresis curves, respectively obtained with the four methods implemented in this work. In addition, Tables 4 and 5 show the values of the average and maximum errors for each model with respect to the experimental data, both in absolute value and in percentage of the full scale. The error is evaluated according to the following expression.
where ym_i are the model samples, yi are the experimental samples and N is the number of samples. The column labeled Best refers to the minimum error value obtained after J(X) in Equation (1) is minimized during the parameter identification process.  The first thing that can be observed in Figures 9 and 10 is that the worse fit to the experimental data was obtained with the MPI model. This model, as explained in Section 3.2, consists of a set of symmetric classical play operators and a generalized function to model the asymmetry of the hysteresis loop. It is very difficult to find a single generalized function that fit to both external curves efficiently. Tests have been done with polynomial of different degrees, and the best result was achieved with the third-order polynomial in Equation (46). Since the results from this model are not good for the hysteresis curves in Figure 2, probably because classical play operators were primarily intended for adjustment of symmetrical loops, no further comment will be made about it.  The GPI model shows a very good performance for the descending curves, although some significant deviations are observed and circled in Figure 9a. The model is saturated for high values of the input signal, so that the adjustment in this area is not good. This behavior is repeated at the start of each inner descending curve. Regarding the ascending curves, the model shows a poor performance in fitting them, as Figure 10a depicts. The main reason behind these deviations is the difficulty in finding envelope functions that fit precisely to the external curves in the complete input range.
With respect to the POLY model, it shows a great difficulty in adjusting the descending curves (see Figure 9c). It is noted that the descending curves of the model are shifted to the left with respect to the experimental data. This is because it is not possible to find the polynomial in Equation (48) to fit the external descending curve of the tactile sensor. As for the ascending curves (see Figure 10c), the model cannot adjust the internal curves well for the same reason.
The ELAM models for the ascending and descending curves, are shown in Figures 9d and 10d respectively. It can be observed that the adjustment of the descending curves is better than with the other methods. The curves do not saturate as with the GPI model (see Figure 9a). The split point introduced in the ELAM model provides a larger flexibility than that of the POLY model to approximate highly nonlinear curves. Regarding the ascending curves, the ELAM model is also the one that achieves the best fitting, even in the internal ascending branches, where other models have large errors. Figure 11a shows the output voltage from samples measured and obtained with the ELAM method for different descending trajectories (see Figure 2a), as well as the error, whose zoom is displayed at Figure 11b. Figure 12 shows the same data for different ascending curves (see Figure 2b). The error is below 0.1% of the full scale sensor output in both cases. Therefore, the ELAM method provides the model that achieves the best fitting of the experimental data. The samples are the output voltage values obtained every two seconds. The number of measured samples is 211 for the descending curves and 451 for the ascending curves.  Furthermore, regarding the interest stated in the introduction in implementing the compensation algorithms in the local electronics based on an FPGA, the ELAM model is the most suitable because its complexity is similar to the POLY model but lower than that of the other two models based on Prandtl-Ishlinskii operators. Specifically, the GPI model shows a fit to the experimental curves close to that achieved with the ELAM model, but the mathematical operators involved are much more complex. The ELAM method is based on simple mathematical operations such as additions and multiplications, while the GPI and MPI models use exponential functions for the density function and for the envelope functions. These operations are more difficult to implement and require the use of more logical resources in devices such as FPGAs. The parameter identifying process in Section 5 is also simpler in the ELAM method. Finally, the GPI, MPI and POLY methods require of a prior selection step of the functions to be used in the construction of the models, while the ELAM method does not need it because it uses the experimental data for the construction of the hysteresis loops. Only the way to obtain the split point has to be determined, and simple linear (see Equations (11) and (12)) or quadratic (see Equation (21)) expressions provide good results. Therefore the ELAM model is inverted according to equations set out in Section 3.5 and it is used for compensating hysteresis nonlinearities presented by the experimental curves of Figure 2a,b. Figure 13 shows the output of the sensor for a set of loading-unloading cycles when it is compensated with the ELAM method. Figure 13a shows the direct ELAM model and the curves measured experimentally. Figure 13b displays the inverse ELAM model, wherein the pressure calculated by the model associated to the sensor output voltage is represented. Figure 13c displays the pressure calculated by the ELAM model versus the real exerted pressure, i.e., it shows the compensation of the descending cycles of the sensor hysteresis. The same data are shown in Figure 14 for the ascending curves. It is noteworthy that the performance of the ELAM model has been assessed from experimental data instead of using other simulated input that could provide artificially good results. In this respect, the errors observed in Figures 13c and 14c are due to noise in the flat areas of the experimental curves (note that the descending curves are saturated for high input pressure values). These errors limit the useful range of the sensor for a given resolution of the measurement.

Results with Tactile Sensor Matrix
This section shows the application of the ELAM method to compensate the hysteresis of a tactile sensor composed of 256 taxels distributed in 16 rows and 16 columns. Each taxel is an independent sensing unit that must be modeled individually, and its hysteresis nonlinearity has to be compensated with its ELAM model. Table 6 shows the frames obtained for different uniform pressures exerted on the sensor. A large mismatching between taxels is observed. Two frames measured for the same input pressure are compared in the same column of the table, one for an ascending sequence (top) and the other for a descending sequence (bottom). Table 7 shows the same frames once the compensation with the ELAM method is applied. It can be observed that not only the hysteresis nonlinearities but also the mismatching between different taxels are compensated. The value of the relative standard deviation with respect to the full scale output is shown for each frame to quantify the improvement regarding the mismatching after the compensation process. In addition, the difference between the frames in the ascending and descending paths is displayed, and its mean value with respect to the full scale output illustrates the improvement regarding the hysteresis.  In addition, Tables 8 and 9 show the sensor output when the force is exerted with a ring-shaped object. Table 8 shows the frames corresponding to the measured sensor output, while Table 9 shows the compensated frames with the ELAM method.

Conclusions
This paper presents a novel method to compensate for hysteresis nonlinearities observed in the response of a tactile sensor. The sensor shows a very pronounced hysteresis with high nonlinearity, and a large mismatching between the responses of different taxels. The proposed method builds a model that accurately fits the experimental data obtained in the characterization process. This so called ELAM method carries out a linear mapping of the external curves of the measured hysteresis loop to the curves of the inner cycles. Its main feature is the introduction of a split point in the curves to produce a different mapping to the left and right of this point, whose location is provided by an error minimization algorithm. The ELAM model is compared with the models obtained from three other approaches, the generalized Prandtl-Ishlinskii model (GPI), the modified Prandtl-Ishlisnkii model (MPI), and a model based on dominant curves built with polynomials (POLY). The results show that the ELAM model fits the measured data more accurately than the other three methods, especially in the ascending curves, where the other methods perform worse. Another very remarkable advantage of the ELAM method versus the other three methods is that the involved mathematical operations are simpler, so they can be implemented more easily in FPGA devices in order to cope with real-time applications. For instance, the ELAM method does uses neither play operators nor exponential functions that may require the use of additional computational resources. Moreover, the number of parameters to be identified by the error minimization algorithm is higher in the other methods, and they also require a prior selection of the appropriate functions to build the model. The performance of the proposed method is shown with data obtained from measurements of the sensor output when a uniform pressure is exerted on the entire matrix, and also when the force is exerted by objects with different shapes. The output of each taxel is compensated with its own ELAM hysteresis model and a significant reduction of the hysteresis, nonlinearity and mismatching errors, is observed. The ELAM method fits very well to complex and asymmetrical hysteresis curves, which cannot be characterized by mathematical functions in a direct way. This allows the application of the strategy followed by the ELAM method to other types of sensors or actuators. Moreover, it is a flexible method, since more split points can be added to divide the curves into a larger number of segments, so different mappings can be done for each segment and a good fitting can be achieved despite the complexity of the hysteresis loops.