Modeling of Rate-Independent and Symmetric Hysteresis Based on Madelung’s Rules

Hysteresis is a kind of nonlinearity with memory, which is usually unwanted in practice. Many phenomenological models have been proposed to describe the observed hysteresis. For instance, the Prandtl-Ishlinskii (PI) model, which consists of several backlash operators, is the most widely used. On the other hand, the well-known Madelung’s rules are always used to validate hysteresis models. It is worth pointing out that the PI model obeys Madelung’s rules. In this paper, instead of considering these rules as criteria, we propose a modeling method for symmetric hysteresis by directly constructing the trajectory based on Madelung’s rules. In the proposed method, turning points are recorded and wiped out according to the input value. After the implementation of the recording and wiping-out mechanisms, the curve which the current trajectory moves along can be determined and then the trajectory can be described. Furthermore, the relationship between the proposed method and the PI model is also investigated. The effectiveness of the presented method is validated by simulation and experimental results.


Introduction
Hysteresis is the phenomenon that can be found in a wide variety of smart materials [1][2][3], such as piezoceramics [4][5][6], magnetostrictive materials [7], and shape memory alloys [8]. It brings considerable problems to the application of smart materials because the output of the system actuated by these materials cannot be predicted without the knowledge of the hysteresis behavior of the system. Generally, it is hard to give an accurate and general model through basic physical principles of a material. To simulate the observed hysteresis behavior, many phenomenological models have been developed in the literature, e.g., Bouc-Wen [9][10][11][12], Duhem, Preisach model [1], Maxwell [4], and Prandtl-Ishlinskii (PI) models [13].
The above phenomenological models can further be classified into two categories: differential-based models and operator-based models. The Bouc-Wen and Duhem models are differential-based models since they contain the derivative of the output in the model. Compared with the Duhem model, the Bouc-Wen model is more popular since it has only a few parameters to be identified. For this reason, it has been widely used in precision control systems actuated by piezoelectric ceramics [10,11]. However, it was shown in Reference [12] that the Bouc-Wen model is inferior to the Preisach model in the trajectory prediction of piezoelectric actuators (PEAs). The inverse of a hysteresis model is also important since it can be used as a feedforward controller to reduce the hysteretic effect. Unfortunately, the construction of the inverse for the Bouc-Wen model is a difficult task since it cannot be inversed analytically [3]. In Reference [9], a Bouc-Wen least square support vector machine was proposed for hysteresis compensation without the need to model the inverse. A multiplicative inverse structure was proposed in Reference [10] to compensate the hysteresis modeled by the Bouc-Wen model. The Preisach model, Maxwell model, and PI model all belong to operator-based models. The main idea of operator-based models is to use the weighted sum of several simple motions to describe complex movements. For example, the simple relay hysteron is the basic operator in the Preisach model. Generally, the number of adopted operators influences the accuracy of the Preisach model. The accuracy grows with the number of operators. However, the computation time also increases. It is not convenient to use the Preisach model in real-time applications since it cannot be inversed analytically. Some strategies were proposed to overcome these disadvantages of the Preisach model, such as the inverse multiplicative technique [14] and the recursive approach [15]. As for the Maxwell model, the basic operator is the elasto-slide element, which consists of a massless linear spring and a massless block that suffers from Coulomb friction. To give a physical insight into PEAs, the elasto-slide element was replaced with the resistive capacitance element in Reference [4]. Among the operator-based models, the PI model is the most widely used in practice [10,13]. It is a weighted superposition of several backlash operators. One of the advantages of the PI model lies in the fact that its inverse exists and is also of the PI type. The parameters of the inverse PI model can be derived from those of the PI model, and this property facilitates the design of feedforward controllers. It should be emphasized that the above mentioned phenomenological models are only good at describing rate-independent and symmetric hysteresis.
However, the observed hysteresis is always rate-dependent or asymmetric, or both of them. A rate-dependent circuit model was proposed in Reference [16], where the creep phenomenon was also considered. Note that this model has an equivalent hardware circuit. In References [17,18], the rate-dependent PI model was developed by dynamically changing its threshold values according to the input rate. In Reference [19], the Bouc-Wen model cascaded with a linear dynamic model was used to describe rate-dependent hysteresis. Based on the classical PI model, many modified forms have been presented to describe asymmetric hysteresis [13,20,21] and rate-dependent hysteresis [22][23][24]. In essence, the core of these models is still the classical PI model. Therefore, the investigation for rate-independent and symmetric hysteresis models is important, since they are the basis of describing rate-dependent and asymmetric hysteresis.
On the other hand, German physicist Madelung presented three rules based on his observation of the hysteresis phenomenon in the early 20th century. These rules are always used to validate hysteresis models [25], such as the Preisach model and the PI model [13]. In this paper, instead of considering Madelung's rules as criteria, we propose a modeling method for the rate-independent and symmetric hysteresis based on these rules. In our method, turning points are recorded as the record of movement history and then wiped out according to the input value. After the implementation of the recording and wiping-out mechanisms, a key turning point named as the current turning point (CTP) is determined. Then the trajectory of the symmetric hysteresis can be described by this CTP. The contributions of this paper are listed as follows:

•
We propose a modeling method to describe the symmetric hysteresis by directly constructing its trajectory based on Madelung's rules rather than considering these rules as criteria. Furthermore, this method is translated into an algorithm that can be run by digital processors.

•
The relationship between the proposed method and the PI model is investigated.
The remainder of this paper is organized as follows: In Section 2, the proposed hysteresis modeling method is explained in detail. Simulations and experiments are conducted in Section 3 to show the effectiveness of this method. At last, some discussions and conclusions are drawn in Section 4.

Madelung's Rules and Their Applications in Trajectory Construction
In this subsection, the definition of turning point is first given. A turning point is a point on the trajectory where the input changes its direction, for instance, point A 0 and point A 1 as shown in Figure 1a,b, respectively. For point A on the trajectory, we use x A and y A to represent its x− and y− coordinates, respectively (see Figure 1a). Before introducing our method, we list Madelung's rules as follows [13,25]:

1.
Any trajectory starting from a turning point is uniquely determined by the coordinates of this point. For example, the turning point A 1 in Figure 1b is the starting point of curve A 1 A 0 , as we will see in the next that the function of A 1 A 0 can be uniquely described by x A 1 and y A 1 .

2.
If any point A 2 on the trajectory, as shown in Figure 1c, becomes a new turning point, then the trajectory leads back to the previous turning point A 1 . In other words, any hysteresis loop is closed. 3.
If the trajectory moving along curve A 2 A 1 is continued beyond A 1 , then it coincides with the continuation of curve A 0 A 0 as if hysteresis loop A 1 A 2 − A 2 A 1 did not exist, as shown in Figure 1d.
In addition to the above three rules, a fourth rule can be given for symmetric hysteresis.

4.
The hysteresis loops of symmetric hysteresis are centrally symmetric.
Sensors 2018, 18, x 3 of 17 Figure 1a,b, respectively. For point A on the trajectory, we use A x and A y to represent its x − and y − coordinates, respectively (see Figure 1a). Before introducing our method, we list Madelung's rules as follows [13,25]: 1. Any trajectory starting from a turning point is uniquely determined by the coordinates of this point. For example, the turning point 1 A in Figure 1b is the starting point of curve 1 0 A A , as we will see in the next that the function of 1 0 A A can be uniquely described by A A A A did not exist, as shown in Figure 1d.
In addition to the above three rules, a fourth rule can be given for symmetric hysteresis.
4. The hysteresis loops of symmetric hysteresis are centrally symmetric.  In our method, it is assumed that the starting point of the hysteretic trajectory is point A 0 which corresponds to the minimum input value, as shown in Figure 1a. Starting from point A 0 , the trajectory moves along the major ascending curve A 0 A 0 . Note that the function of A 0 A 0 is f 0 (x), where x ∈ D f 0 = [A 0 , A 0 ] and D f 0 is the domain of f 0 (x). According to Rule 2, the trajectory leads back to A 0 along A 1 A 0 after it is reversed at point Figure 1b). We can see from Rule 4 that A 1 A 0 and A 0 A 1 are centrally symmetric as illustrated in Figure 1b, and, thus, point (x, f 1 (x)) on A 1 A 0 and its corresponding point (x t , f 0 (x t )) on A 0 A 1 has the following relationship: (1), the function of A 1 A 0 can be expressed by If point A 2 on A 1 A 0 becomes a new turning point, the new curve A 2 A 1 can be described by At this time, the trajectory moves along curve A 2 A 1 as illustrated in Figure 1c. If the trajectory is continued beyond the former turning point A 1 , the trajectory will then move onto curve A 0 A 0 as if hysteresis loop A 1 A 2 − A 2 A 1 did not exist at all (see Rule 3), which means A 1 and A 2 can be wiped out. Otherwise, the trajectory should be reversed at turning point A 3 which belongs to A 2 A 1 , and the function of the new curve A 3 A 2 can be written as Similarly, we can construct the following trajectory for any input signal. Without loss of generality, when turning points A 0 , A 1 , · · · , A m are generated in sequence, the function of curve A k A k−1 can be expressed by It is noted that Equation (5) is not suitable for real-time applications since it is an iterative equation. In fact, Equation (5) can be simplified as follows: Note that there are three terms in Equation (6) for each value of k, the first term is f k−1 (x A k ) and the other two terms can be computed by f 0 (x). This means that curve A k A k−1 can be considered to be determined by the coordinates of point A k and, thus, Rule 1 is validated. According to Rule 3, the trajectory is transferred from Until now, we have shown how to use Madelung's rules to describe the trajectory of symmetric hysteresis. Next, we will translate this method into a computer algorithm. Before we proceed, some observations are made as follows:

•
As the recorders of the movement history, all the turning points should be identified and recorded.

•
If the trajectory surpasses previous turning point A k−1 , it will be transferred to the previous curve A k−2 A k−3 . Point A k−1 and its previous point A k should be wiped out since they are of no use to describe the future trajectory.

•
The domain range of the curve described by Equation (6) is decreasing with k, which means that • The trajectory at any time instant must be on one of the curves described by Equation (6), namely the current curve. If the current curve is determined, the trajectory can then be described. In fact, this curve can be determined by finding the minimum domain that contains the input value x among the existing curves. Mathematically, we have where D c is the domain of the current curve, σ(D f ) represents the length of D f , and • If k is odd, the curve described by Equation (6) is a descending curve, and it is an ascending curve when k is even.

Algorithm and Complex Analysis
Based on the above analysis, it is obvious that the trajectory at any time instant must be on the current curve. The starting point of this curve is defined to be the current turning point (CTP). Once the CTP is determined, the trajectory can be described by Equation (6) according to Rule 1. The flowchart of the algorithm is shown in Figure 2, where each block is labeled to facilitate the presentation.  , , At the beginning, the initialization is conducted in 0 S . The function of the major ascending curve is 0 ( ) f x , and it is described by the following polynomial is the coefficient and n is the degree of the polynomial. In the algorithm, TPs are further classified into two types: left turning points (LTPs) and right turning points (RTPs). At the beginning, the initialization is conducted in S 0 . The function of the major ascending curve is f 0 (x), and it is described by the following polynomial f 0 (x) = c n x n + c n−1 x n−1 + · · · c 1 x + c 0 (9) where c i (i = 0, 1, · · · , n) is the coefficient and n is the degree of the polynomial. In the algorithm, TPs are further classified into two types: left turning points (LTPs) and right turning points (RTPs). For instance, A 0 , A 2 , A 4 , · · · are LTPs and A 1 , A 3 , A 5 , · · · are RTPs. Array X LTP and array Y LTP are used to record the x− and y− coordinates of LTPs, respectively. Similarly, X RTP and Y RTP are used to record the coordinates of RTPs. Auxiliary variables x 1 , y 1 , and x 2 are utilized to identify turning points.
The Boolean variable TP_ f lag is the indicator for the type of the latest turning point. For instance, TP_ f lag = 0 implies that the latest turning point is an LTP and the trajectory moves along some ascending curve. The integer variable index represents the position of the CTP in turning point arrays. Furthermore, temporary variables le f t, right, and mid are used in the binary search method. In S 1 , the new input value x is read into the program. The remainder of the algorithm can mainly be divided into three parts: The first part consists of S 2 , S 3 , S 4_1 , and S 4_2 . The first step aims to identify turning points. To identify whether the previous point (x 1 , y 1 ) is a turning point, the following inequality is utilized.
where x 1 and x 2 are two input values previous to x. If Inequality (10) holds, point ( Next, this point is restored in either S 4_1 or S 4_2 according to its type. The second part contains S 5 , S 6_1 , S 6_2 , and S 7 . The objective of this part is to determine the CTP and to wipe out the turning points after the CTP. As mentioned earlier, the current curve and the CTP can be determined by Equations (7) and (8). Herein, we will give simplified results. We can see from Equation (7) that X LTP and X RTP are naturally sorted in ascending order and descending order, respectively. To see the fact, the domains of the functions described by Equation (6) are plotted in Figure 3, where we have implies that the latest turning point is an LTP and the trajectory moves along some ascending curve. The integer variable index represents the position of the CTP in turning point arrays. Furthermore, temporary variables left , right , and mid are used in the binary search method. In 1 S , the new input value x is read into the program. The remainder of the algorithm can mainly be divided into three parts: The first part consists of 2 S , 3 S , 4 _1 S , and 4_ 2 S . The first step aims to identify turning points.
To identify whether the previous point 1 1 ( , ) x y is a turning point, the following inequality is utilized.
where 1 x and 2 x are two input values previous to x . If Inequality (10)  The second part contains 5 S , 6 _1 S , 6 _ 2 S , and 7 S . The objective of this part is to determine the CTP and to wipe out the turning points after the CTP. As mentioned earlier, the current curve and the CTP can be determined by Equations (7) and (8). Herein, we will give simplified results. We can see from Equation (7) that LTP X and RTP X are naturally sorted in ascending order and descending order, respectively. To see the fact, the domains of the functions described by Equation (6) are plotted in Figure 3, where we have On the contrary, if = _ 1 TP flag , the turning point m A is an RTP and Expression (8) can be rewritten as Figure 3. Domains of functions f k (x), k = 1, 2, · · · m.
If TP_ f lag = 0, the latest turning point A m is an LTP and the current curve is an ascending curve. Furthermore, we have x A 0 < x A 1 < · · · < x A m < x, and, thus, Expression (8) can be simplified as follows: On the contrary, if TP_ f lag = 1, the turning point A m is an RTP and Expression (8) can be rewritten as Apparently, the CTP can be determined by locating the position of x in X LTP or X RTP . As shown in S 6_1 and S 6_2 of Figure 2, the binary search method is used to locate the position of x. After the "while" operation, the value of right is exactly the target position. Furthermore, all the turning points after the CTP are wiped out by setting index = right in S 7 .
The third part comprises S 8 , S 9_1 , S 9_2 , and S 10 . In this part, the output y is computed using the CTP (see Equation (6)). In addition, auxiliary variables are updated only when x = x 1 , since successive identical input values will influence the identification of turning points.
Owing to the wiping-out mechanism, the number of turning points cannot be increased infinitely. The maximum number of turning points equals to N u , where N u is the number of different input values. For instance, a data acquisition system with a 12-bit digital-to-analog converter (DAC) has N u = 2 12 = 4096 different input values. On the other hand, since LTPs and RTPs are generated alternatively and eliminated in pairs. Thus, we can simply suppose that the number of LTPs is the same as that of RTPs. Note that the number of turning points should be no more than N u . Therefore, the length of each array in the algorithm can be set to N u /2.
Apparently, the binary search method is the main body of the algorithm, and its computational complexity is T(L) = O(log 2 L), where L is the length of the input array. In our algorithm, the maximum value of L is L = N u /2 = 2 n s −1 , where n s is the resolution of the used DAC. Then we can see that the computational complexity of the algorithm is O(log 2 2 n s −1 ) = O(n s ) and the maximum computation time (MCT) of the algorithm corresponds to the case when L = 2 n s −1 . In addition to computational complexity, the space complexity of the algorithm is O(N u ) = O(2 n s ).

Parameter Identification
The least mean square method is used to identify the parameters of f 0 (x) in Equation (9). The oscillating input signals, such as the sinusoidal and triangular waves, between the minimum and maximum input are used to get the input and output data. When the input is increased from the minimum value x A 0 to the maximum value x A 0 , the trajectory moves along the major ascending curve. However, note that when the input value decreases, the trajectory moves along the descending curve f 1 (x) rather than f 0 (x). In order not to lose the information of the descending curve f 1 (x), the point on the descending curve has to be translated to its corresponding point on the major ascending curve. Note that f 0 (x) and f 1 (x) are centrally symmetric, the point (x, f 1 (x)) can be translated to its corresponding point (x t , f 0 (x t )), as shown in Figure 4, where Based on Equation (14), the output-input data can be reproduced by the algorithm shown in Figure 5. In the algorithm, the output-input characteristics of f 0 (x) can be totally captured by the reproduced arrays X _data and Y _data. Then the least mean square method can be used to derive the coefficients of f 0 (x). Sensors 2018, 18, x 9 of 17  Simulation Example 1: We assume that 2 x ( ∈[0, 1] x ) is the function of the major ascending curve. After the injection of the sinusoidal wave with constant amplitude 1, the input and output data are collected in _ X data and _ Y data , respectively. The plot of _ Y data versus _ X data is shown in Figure 6, where the output data are affected by noises. Setting = 2 n in Equation (9)    Simulation Example 1: We assume that 2 x ( ∈[0, 1] x ) is the function of the major ascending curve. After the injection of the sinusoidal wave with constant amplitude 1, the input and output data are collected in _ X data and _ Y data , respectively. The plot of _ Y data versus _ X data is shown in Figure 6, where the output data are affected by noises. Setting = 2 n in Equation (9) and applying the identification algorithm, we have = − + 2 0 ( ) 1.0005 0.0008 0.0002 f x x x (see Figure 6), which is very close to the function 2 x . Simulation Example 1: We assume that x 2 (x ∈ [0, 1]) is the function of the major ascending curve. After the injection of the sinusoidal wave with constant amplitude 1, the input and output data are collected in X_data and Y_data, respectively. The plot of Y_data versus X_data is shown in Figure 6, where the output data are affected by noises. Setting n = 2 in Equation (9) and applying the identification algorithm, we have f 0 (x) = 1.0005x 2 − 0.0008x + 0.0002 (see Figure 6), which is very close to the function x 2 .

Relationship with the PI Model
Generally, the PI model consists of a set of backlash operators. Considering the fact that input signals for PEAs are positive, one-side backlash operators [20,21] are used in this paper. The PI model with one-side backlash operators can be written as where [ ]( ) i H x t is the one-side backlash operator, i w is the weight value, i r is the threshold value, and PI n is an integer which represents the number of operators. The function of the major ascending curve the PI model is x r x r r (17) It is assumed that = < < <  PI 1 2 0 n r r r . At this time, we can present a question as follows: If PI ( ) f x is used to replace 0 ( ) f x , is the trajectory construction method able to compute the output of the PI model? The answer is "yes" since the PI model is symmetric and obedient to Madelung's rules. In other words, the output of the trajectory construction method is equivalent to that of the PI model. This equivalence is not accidental since the trajectory construction method captures the symmetric property of the PI model and guarantees the equivalence in principle. Simulation Example 2: To simulate the hysteresis phenomenon in Example 1, we use ten operators to construct the PI model. The parameters of the PI model are listed in Table 1. Next, a series of different types of input signals are used to excite both models. The output signals of these two models are also plotted in Figure 7, where we can see that these two methods have the same output signals.

Relationship with the PI Model
Generally, the PI model consists of a set of backlash operators. Considering the fact that input signals for PEAs are positive, one-side backlash operators [20,21] are used in this paper. The PI model with one-side backlash operators can be written as where H i [x](t) is the one-side backlash operator, w i is the weight value, r i is the threshold value, and n PI is an integer which represents the number of operators. The function of the major ascending curve the PI model is where It is assumed that 0 = r 1 < r 2 < · · · < r n PI . At this time, we can present a question as follows: If f PI (x) is used to replace f 0 (x), is the trajectory construction method able to compute the output of the PI model? The answer is "yes" since the PI model is symmetric and obedient to Madelung's rules. In other words, the output of the trajectory construction method is equivalent to that of the PI model. This equivalence is not accidental since the trajectory construction method captures the symmetric property of the PI model and guarantees the equivalence in principle.
Simulation Example 2: To simulate the hysteresis phenomenon in Example 1, we use ten operators to construct the PI model. The parameters of the PI model are listed in Table 1. Next, a series of different types of input signals are used to excite both models. The output signals of these two models are also plotted in Figure 7, where we can see that these two methods have the same output signals.   Figure 7. Output signals of the PI model and the presented method.

HIL Simulations
Hardware-in-the-Loop (HIL) simulations are conducted in this subsection to investigate the computation time of the trajectory construction method with a DSP28335 development board. One pin of the DSP is programmed to denote the state of the algorithm. The pin turns to be a high level when the algorithm starts a new cycle and becomes low at the end of the cycle. The duration of the high level represents the computation time of the trajectory construction method.
It is obvious from Figure 2 that the computation time depends on the number p of numerical comparisons in the binary search. According to the property of the binary search, there is a logarithmic relation between p and m, with the function given by where m is the number of non-empty elements in array X LTP or array X RTP . It is noted that the computation time increases with the amount of turning points. Fortunately, thanks to the logarithmic relation of Equation (18), p grows slowly with the increase of m. On the other hand, owing to the wiping-out mechanism, m cannot increase infinitely. As mentioned above, the maximum value of m is N u /2, where N u represents the number of different input values in the data acquisition system. To cope with the following experiments, we assume N u = 2 n s = 2 12 = 4096 in the HIL simulations. Figure 8 depicts the relation between computation time and the number of comparisons in the binary search. The data in Figure 8 can be fitted by a straight line 0.19p + 0.57, which means that it takes about 0.19 µs to complete a comparison. According to Equation (18), the maximum value of p is p = log 2 (N u /2) + 1 = n s = 12 (19) and the computation time corresponding to this p value is 2.80 µs, which is also the maximum computation time (MCT). Since the computation time varies with the number of turning points, it is important to estimate the MCT of the trajectory construction method. Furthermore, we can see from Equation (19) that the MCT is related to n s , which is the resolution of the DAC.

Experiments
In this subsection, we will show how to predict the output of PEAs. The experimental system is set up as shown in Figure 9, where the micro positioner 20VS12 actuated by PEA is from Tomorrow Core Company. The positioner has an integrated resistance strain-gauge transducer. The high voltage amplifier XE-650.CA is also from the same company. The data acquisition system consists of

Experiments
In this subsection, we will show how to predict the output of PEAs. The experimental system is set up as shown in Figure 9, where the micro positioner 20VS12 actuated by PEA is from Tomorrow Core Company. The positioner has an integrated resistance strain-gauge transducer. The high voltage amplifier XE-650.CA is also from the same company. The data acquisition system consists of the above used DSP28335 development board and an AD/DA board with 12-bit analog-to-digital and digital-to-analog converters. The running frequency of the DSP28335 is set to 150 MHz. The algorithm is compiled on the host computer and then downloaded into the DSP28335.

Experiments
In this subsection, we will show how to predict the output of PEAs. The experimental system is set up as shown in Figure 9, where the micro positioner 20VS12 actuated by PEA is from Tomorrow Core Company. The positioner has an integrated resistance strain-gauge transducer. The high voltage amplifier XE-650.CA is also from the same company. The data acquisition system consists of the above used DSP28335 development board and an AD/DA board with 12-bit analog-to-digital and digital-to-analog converters. The running frequency of the DSP28335 is set to 150 MHz. The algorithm is compiled on the host computer and then downloaded into the DSP28335. First, a triangular wave with constant amplitude 40 V as shown in Figure 10 is used to excite the positioner system to identify the function 0 ( ) f x . The input-output data are collected and stored in _ X data and _ Y data , respectively. Figure 11 depicts the plot of _ Y data versus _ X data . First, a triangular wave with constant amplitude 40 V as shown in Figure 10 is used to excite the positioner system to identify the function f 0 (x). The input-output data are collected and stored in X_data and Y_data, respectively. Figure 11 depicts the plot of Y_data versus X_data.
The function 0 ( ) f x is also plotted in Figure 11, where we can see that the measured data can be well described by 0 ( ) f x . Next, a test signal is used to validate the effectiveness of the proposed method, as shown in Figure 12. When ∈[0, 13.4] t , the test signal is a triangular wave with a decreasing amplitude, which begins to increase at = 13.4 t s. The position prediction results are shown in Figure 13, where the maximum error is 0.12 µm. Setting n = 3 and applying the parameter identification algorithm, we derive the function f 0 (x) as follows: The function f 0 (x) is also plotted in Figure 11, where we can see that the measured data can be well described by f 0 (x). Next, a test signal is used to validate the effectiveness of the proposed method, as shown in Figure 12. When t ∈ [0, 13.4], the test signal is a triangular wave with a decreasing amplitude, which begins to increase at t = 13.4 s. The position prediction results are shown in Figure 13, where the maximum error is 0.12 µm.  Figure 14 depicts the number of turning points in the above hysteretic movement. The number of turning points increases one at a time and decreases in pairs (an LTP and an RTP, respectively). In the beginning, there is only one turning point, which is the starting point. When the input signal leads back to 0 V at = 3 t s, the number of turning points becomes to zero. At ∈[3, 13.4] t , the number of turning points is monotonically increased. When > 13.4 t s, the amplitude of the test signal increases with time. The number of turning points is decreased to zero finally when the test signal arrives at 0 V. As shown in Figure 14, the maximum number of turning points is 14, and more precisely, there are 7 LTPs and 7 RTPs. Thus, we have = 7 m and = 4 p . As shown in Figure 8, the computation time corresponding to = 4 p is 1.23 µs. As a matter of fact, the computation time in the whole movement is no greater than 1.23 µs, which is much smaller than the MCT.
Voltage (V) Figure 12. Test signal. Figure 14 depicts the number of turning points in the above hysteretic movement. The number of turning points increases one at a time and decreases in pairs (an LTP and an RTP, respectively). In the beginning, there is only one turning point, which is the starting point. When the input signal leads back to 0 V at t = 3 s, the number of turning points becomes to zero. At t ∈ [3, 13.4], the number of turning points is monotonically increased. When t > 13.4 s, the amplitude of the test signal increases with time. The number of turning points is decreased to zero finally when the test signal arrives at 0 V. As shown in Figure 14, the maximum number of turning points is 14, and more precisely, there are 7 LTPs and 7 RTPs. Thus, we have m = 7 and p = 4. As shown in Figure 8, the computation time corresponding to p = 4 is 1.23 µs. As a matter of fact, the computation time in the whole movement is no greater than 1.23 µs, which is much smaller than the MCT. number of turning points is monotonically increased. When > 13.4 t s, the amplitude of the test signal increases with time. The number of turning points is decreased to zero finally when the test signal arrives at 0 V. As shown in Figure 14, the maximum number of turning points is 14, and more precisely, there are 7 LTPs and 7 RTPs. Thus, we have = 7 m and = 4 p . As shown in Figure 8, the computation time corresponding to = 4 p is 1.23 µs. As a matter of fact, the computation time in the whole movement is no greater than 1.23 µs, which is much smaller than the MCT.

Discussions
As shown in Equation (9), there are 1 n + parameters if the degree of the polynomial function 0 ( ) f x is n . In the experiments, a third order polynomial function (see Equation (20)) is used to describe the major ascending curve. Thus, only four parameters need to be identified. For comparison, if the PI model with ten backlash operators is used, we have to identify twenty parameters. Therefore, the proposed method is superior to the PI model in parameter identification.
However, since the turning points have to be recorded in the proposed method, the worst case has to be considered. From the above analysis, we know that the maximum number of turning points is related to the number of different input values. In other words, a certain amount of memory is needed. However, it is worth pointing out that the demand for memory resources is affordable to most of industrial computers, as shown in Table 2. Different from traditional operator-based method, the trajectory construction method uses the observed rules for hysteresis (Madelung's rules) and a measured curve (the major ascending curve) to simulate the symmetric hysteresis. As we have known that hysteresis is a kind of nonlinearity with memory. The concept of turning point is important in the proposed method since turning point is the recorder of movement history. In other words, the turning point is the memory of the

Discussions
As shown in Equation (9), there are n + 1 parameters if the degree of the polynomial function f 0 (x) is n. In the experiments, a third order polynomial function (see Equation (20)) is used to describe the major ascending curve. Thus, only four parameters need to be identified. For comparison, if the PI model with ten backlash operators is used, we have to identify twenty parameters. Therefore, the proposed method is superior to the PI model in parameter identification.
However, since the turning points have to be recorded in the proposed method, the worst case has to be considered. From the above analysis, we know that the maximum number of turning points is related to the number of different input values. In other words, a certain amount of memory is needed. However, it is worth pointing out that the demand for memory resources is affordable to most of industrial computers, as shown in Table 2. Different from traditional operator-based method, the trajectory construction method uses the observed rules for hysteresis (Madelung's rules) and a measured curve (the major ascending curve) to simulate the symmetric hysteresis. As we have known that hysteresis is a kind of nonlinearity with memory. The concept of turning point is important in the proposed method since turning point is the recorder of movement history. In other words, the turning point is the memory of the hysteretic movement. The proposed method is an intuitive and easy-to-use hysteresis modeling method, and it is also beneficial to the in-depth understanding of the hysteresis phenomenon.

Conclusions
In this paper, a symmetric hysteresis modeling method has been proposed based on Madelung's rules. The turning points are recorded and eliminated according to the input value. The key point of the proposed method is to determine the CTP. Then the current can be determined, and the output can be computed by Equation (5). Furthermore, the relationship between our method and the PI model has also been investigated. It is found out that the outputs of both methods are equivalent to each other when f 0 (x) is replaced with f PI (x). It should be noted that the proposed method has a relatively high demand for memory resources. Fortunately, the memory consumption is affordable to most of industrial computers. Even embedded processors can be used with low resolution requirements. Simulation and experiment results have shown the effectiveness of the proposed algorithm. In the future work, we will focus on how to reduce the memory usage of the trajectory construction method.
Author Contributions: K.C. and R.L. proposed the theoretical analysis, and designed the experiments; K.C. conducted the simulations and wrote the paper.