Modeling the Viscoelastic Hysteresis of Dielectric Elastomer Actuators with a Modified Rate-Dependent Prandtl–Ishlinskii Model

Dielectric elastomer actuators (DEAs) are known as a type of electric-driven artificial muscle that have shown promising potential in the field of soft robotics. However, the inherent viscoelastic nonlinearity makes the modeling and control of DEAs challenging. In this paper, we propose a new phenomenological modeling approach with the Prandtl–Ishlinskii (P–I) model to characterize the viscoelastic hysteresis nonlinearity of DEAs. Differently from the commonly used physics-based models, the developed phenomenological model, called the modified rate-dependent P–I model (MRPIM), produces behavior similar to that of physics-based models but without necessarily considering physical insight into the modeling problem. In this way, the developed MRPIM can characterize the asymmetric and rate-dependent viscoelastic hysteresis with a relative simple mathematical format using only the experimental data. To validate the development, experimental tests were conducted with seven different frequencies; four were selected to identify the model parameters and the rest of the data were used to further verify the model. Comparisons between the model prediction and experimental data demonstrate that the MRPIM can precisely describe the viscoelastic hysteresis effect of DEAs with a maximum prediction error of 9.03% and root-mean-square prediction error of 4.50%.


Introduction
Dielectric elastomer actuators (DEAs) are a promising actuation technology for soft robotics, owing to the advantages of large strain (>380%), high energy density (>3.4 MJ/m 3 ), and fast response (in the order of milliseconds) [1,2]. The working principle of DEAs is simple: under an applied electric field, the Maxwell stress squeezes the elastomer membrane in thickness and leads to expansion in area while keeping a constant volume [3]. On the basis of this principle, many DEA configurations have been invented to generate different degree-of-freedom (DOF) motions, such as rectilinear motion [4], rotational motion [5], and bending motion [6]. It is also interesting to find that these motions can mimic natural muscle, which makes them promising as artificial muscles for bioinspired soft robots, for instance, an arm wrestling robot [7], two locomotive robots [8,9], and a swimming fish [10].
Despite the diverse achievements in mechanism design, there still remains a challenge to accurately characterize the response of DEAs because of their inherent viscoelasticity. Figure 1a illustrates an experimental response of a DEA under a sinusoidal exciting input (Figure 1b). From the results, we can observe that the response can be divided into two regions: a transition region and To address this challenge, many efforts have been made in the literature. In general, the viscoelastic hysteresis can be considered as one of several dissipative phenomena, and some physicsbased models [11,[24][25][26] have been established to explain it. A main advantage of such physics-based models depends on the fact that they are built on the basis of understanding the physical behaviors of DEAs. However, they are generally very complicated in order to describe the experimental phenomena. For example, the creep effect was explained by a Gent model with 8 material parameters [27], but 20 parameters should be used to describe the creep, hysteresis, and the shift in the displacement's peak when the DEA is subjected to a cycle loading with a constant slope [11]. In addition, when the slope is variable, the model will be more complicated.
Rather than using physics-based models, some phenomenological models have been reported to describe the viscoelastic hysteresis effect of DEAs by using only experimental data without considering their physical behaviors. For example, by linearizing the input voltage, the hysteresis loops of a tubular DEA became symmetric and ellipse-like, and an elliptical modeling method was used to describe the hysteresis effect [28]. In addition, a polynomial-based approach was developed by using two polynomials to model the increasing and decreasing hysteresis loops [29].
However, the drawbacks of these existing phenomenological models are clear: they are either rate-independent or it is necessary to change the model parameters for a rate-dependent hysteresis description. It is known that the P-I model is a popular phenomenological model to describe the To address this challenge, many efforts have been made in the literature. In general, the viscoelastic hysteresis can be considered as one of several dissipative phenomena, and some physics-based models [11,[24][25][26] have been established to explain it. A main advantage of such physics-based models depends on the fact that they are built on the basis of understanding the physical behaviors of DEAs. However, they are generally very complicated in order to describe the experimental phenomena. For example, the creep effect was explained by a Gent model with 8 material parameters [27], but 20 parameters should be used to describe the creep, hysteresis, and the shift in the displacement's peak when the DEA is subjected to a cycle loading with a constant slope [11]. In addition, when the slope is variable, the model will be more complicated.
Rather than using physics-based models, some phenomenological models have been reported to describe the viscoelastic hysteresis effect of DEAs by using only experimental data without considering their physical behaviors. For example, by linearizing the input voltage, the hysteresis loops of a tubular DEA became symmetric and ellipse-like, and an elliptical modeling method was used to describe the hysteresis effect [28]. In addition, a polynomial-based approach was developed by using two polynomials to model the increasing and decreasing hysteresis loops [29].
However, the drawbacks of these existing phenomenological models are clear: they are either rate-independent or it is necessary to change the model parameters for a rate-dependent hysteresis description. It is known that the P-I model is a popular phenomenological model to describe the complex hysteresis nonlinearity of smart material actuators [16][17][18][19]22], but less attention has been paid to employ the P-I model for the viscoelastic hysteresis description of DEAs.
As a first attempt, this paper presents a modified rate-dependent P-I model (MRPIM) to describe the asymmetric and rate-dependent viscoelastic hysteresis effect of a DEA. The developed MRPIM consists of three parts: (i) a fourth-order polynomial is introduced to describe the asymmetric behavior; (ii) fixed weights and thresholds of the play operators with dynamical envelope functions are used to characterize the rate-dependent behavior; (iii) the second-order derivative of the input voltage is introduced to the fourth-order polynomial for describing the peak-to-peak displacements of the hysteresis loops that highly depend on the frequency of the input voltage. The advantages of our model lie in the facts that (1) it can precisely describe the viscoelastic hysteresis of DEAs with only experimental data; (2) compared with physics-based models, there are relatively fewer parameters to describe both asymmetric and rate-dependent hysteresis behavior of DEAs; and (3) differently from existing phenomenological models [28,29], our model is rate-dependent and only needs one step to identify all the parameters without data pre-processing. In order to verify the effectiveness of the MRPIM, a planar DEA was fabricated, and experiments were conducted. Experimental results demonstrate that the MRPIM can accurately predict the hysteresis effect of the DEA with different frequencies (in the range of 0.02 to 1 Hz). The maximum prediction error and root-mean-square prediction error were under 9.03% and 4.50%, respectively. In addition, it should be noted that our MRPIM can be used not only to describe the hysteresis effect of the planer DEA, but also can be applied for different DEAs (e.g., [28,29]) that have asymmetric and rate-dependent viscoelastic hysteresis behavior. The only difference depends on the model parameters, which need to be identified on the basis of experimental data of different DEAs.
The remainder of this paper is arranged as follows. Section 2 describes our fabricated planar DEA and the experimental setup; the viscoelastic hysteresis effect of DEAs is also analyzed in this section. Section 3 introduces the development, identification, and verification of the MRPIM. Finally, Section 4 concludes this paper.

Fabrication of the Planar DEA
In this work, a planar DEA as shown in Figure 2a was fabricated for proof-of-concept examining. The fabrication processes can be summarized as follows: (i) A prestrained DE membrane (3M VHB 4905, L 1 × L 2 = 100 mm × 40 mm, λ 1 × λ 2 = 6 × 4) is supported by a 3D printed frame. (ii) A platform (P 1 × P 2 = 80 mm × 5 mm) is fixed on the middle of the DE membrane to separate the DE membrane into passive and active regions. (iii) Carbon grease is used as electrodes covering both sides of the active region (the right part of the DE membrane). It should be noted that in order to avoid electric break down, there is an interval of 1 mm between the electrodes and the edges of other parts, such that each electrode is a rectangular area with a length of 98 mm and a width of 15.5 mm. (iv) The electrodes are connected to a high-voltage source.
The planar DEA can provide a linear actuation with one DOF. As shown in Figure 2b, its working principle can be described as follows: when a voltage is applied to the electrodes, the active region will expand along the Y-direction and the passive region works as a nonlinear spring [30,31], such that the platform will generate an output displacement. Figure 2c shows a photo of the fabricated planar DEA.

Experimental Setup
In order to capture and characterize the viscoelastic hysteresis effect of the planar DEA, we built an experimental setup as shown in Figure 3a. A high-voltage amplifier (TREK 10/10-HS, TREK INC., New York, NY, USA) with a fixed gain of 60 dB was used to apply an exciting voltage for the planar DEA, and a laser sensor (Micro-Epsilon ILD2300-100, range of 0-100 mm with an analogue output of 0-10 V, Micro-Epsilon, Ortenburg, Germany) was utilized to measure the output displacement of the planar DEA. A dSPACE-DS1103 (dsPACE, Paderborn, Germany) board equipped with 16-bit analogue-to-digital converters (ADCs) was adopted to generate a control signal for the high-voltage amplifier and to capture the output signal of the laser sensor. In this work, the sampling time was set to be 1 ms. Figure 3b gives a block diagram of the experimental setup. It should be noted that when the amplitude of the exciting voltage was higher than 3 kV, we could observe several failure phenomena, such as instability and wrinkling [2]. On the other hand, the output displacement of the planar DEA was too small to be observed when the exciting voltage was lower than 0.5 kV. Therefore, we kept the exciting voltage in the range of 0.5-2.5 kV in the study.

Experimental Phenomena
It is known that the viscoelastic hysteresis effect of DEAs is asymmetric and rate-dependent. In order to investigate the behavior of the viscoelastic nonlinearity, we firstly tested the dynamical response of the planar DEA. In this sense, the exciting voltages were selected to be same amplitude with different frequencies, which can be formed as where f represents the frequency of the exciting voltage. It should be noted that a phase delay of 0.5π in Equation (1) was used to reduce the response of the DEA because of the initial input voltage at t = 0. Of course, other phase delays could also be selected. Without loss of generality, we selected 0.5 φ π = . In this work, the relationships between the input voltage and displacement were measured under seven randomly selected different frequencies (0.02, 0.05, 0.1, 0.2, 0.4, 0.8, and 1 Hz). Remark: It is well known that the bandwidth of Very High Bond (VHB) based DEAs is usually limited to 3 Hz. For the specific planar DEA used in this work, the experimental results demonstrate that when the frequency of the exciting voltage increased from 0.02 to 1 Hz, the amplitude of the output displacement monotonously decreased by 27.26%, which means that the working frequency

Experimental Setup
In order to capture and characterize the viscoelastic hysteresis effect of the planar DEA, we built an experimental setup as shown in Figure 3a. A high-voltage amplifier (TREK 10/10-HS, TREK INC., New York, NY, USA) with a fixed gain of 60 dB was used to apply an exciting voltage for the planar DEA, and a laser sensor (Micro-Epsilon ILD2300-100, range of 0-100 mm with an analogue output of 0-10 V, Micro-Epsilon, Ortenburg, Germany) was utilized to measure the output displacement of the planar DEA. A dSPACE-DS1103 (dsPACE, Paderborn, Germany) board equipped with 16-bit analogue-to-digital converters (ADCs) was adopted to generate a control signal for the high-voltage amplifier and to capture the output signal of the laser sensor. In this work, the sampling time was set to be 1 ms. Figure 3b gives a block diagram of the experimental setup. It should be noted that when the amplitude of the exciting voltage was higher than 3 kV, we could observe several failure phenomena, such as instability and wrinkling [2]. On the other hand, the output displacement of the planar DEA was too small to be observed when the exciting voltage was lower than 0.5 kV. Therefore, we kept the exciting voltage in the range of 0.5-2.5 kV in the study.

Experimental Phenomena
It is known that the viscoelastic hysteresis effect of DEAs is asymmetric and rate-dependent. In order to investigate the behavior of the viscoelastic nonlinearity, we firstly tested the dynamical response of the planar DEA. In this sense, the exciting voltages were selected to be same amplitude with different frequencies, which can be formed as where f represents the frequency of the exciting voltage. It should be noted that a phase delay of 0.5π in Equation (1) was used to reduce the response of the DEA because of the initial input voltage at t = 0. Of course, other phase delays could also be selected. Without loss of generality, we selected φ = 0.5π. In this work, the relationships between the input voltage and displacement were measured under seven randomly selected different frequencies (0.02, 0.05, 0.1, 0.2, 0.4, 0.8, and 1 Hz). Remark: It is well known that the bandwidth of Very High Bond (VHB) based DEAs is usually limited to 3 Hz. For the specific planar DEA used in this work, the experimental results demonstrate that when the frequency of the exciting voltage increased from 0.02 to 1 Hz, the amplitude of the output displacement monotonously decreased by 27.26%, which means that the working frequency of the planar DEA was limited to about 1 Hz. Therefore, although the maximum effective frequency of the MRPIM was limited to 1 Hz, it could still satisfy the practical application of the planar DEA.
Polymers 2018, 10, x FOR PEER REVIEW 5 of 12 of the planar DEA was limited to about 1 Hz. Therefore, although the maximum effective frequency of the MRPIM was limited to 1 Hz, it could still satisfy the practical application of the planar DEA.  Figure 4 shows the experimental results when the frequency of the exciting voltage equaled 0.1 Hz. As mentioned in Figure 1, the displacement could be separated into a transition region and stable region. During the transition region, a creep effect (slow drift of the displacement [14]) was observed and lasted about 400 s before reaching the stable region. It should be noted that although such creep is caused by viscoelasticity of the DEA, it has been explained by many well-defined models [12][13][14]. Differently from the transition region, the viscoelastic hysteresis effect dominates the dynamical response in the stable region. However, how to model the viscoelastic hysteresis effect of DEAs is still challenging. Therefore, we were motivated to establish a MRPIM to describe the viscoelastic hysteresis effect of DEAs. To avoid the influence of the creep effect, we only took the experimental data of the stable region into consideration for developing the MRPIM. To further analyze the features of the viscoelastic hysteresis effect of the planar DEA, we compared the viscoelastic hysteresis loops with different frequencies (0.02, 0.1, 0.4, and 1 Hz). On the basis of the experimental results shown in Figure 5, we could observe the following:   Figure 4 shows the experimental results when the frequency of the exciting voltage equaled 0.1 Hz. As mentioned in Figure 1, the displacement could be separated into a transition region and stable region. During the transition region, a creep effect (slow drift of the displacement [14]) was observed and lasted about 400 s before reaching the stable region. It should be noted that although such creep is caused by viscoelasticity of the DEA, it has been explained by many well-defined models [12][13][14]. Differently from the transition region, the viscoelastic hysteresis effect dominates the dynamical response in the stable region. However, how to model the viscoelastic hysteresis effect of DEAs is still challenging. Therefore, we were motivated to establish a MRPIM to describe the viscoelastic hysteresis effect of DEAs. To avoid the influence of the creep effect, we only took the experimental data of the stable region into consideration for developing the MRPIM. To further analyze the features of the viscoelastic hysteresis effect of the planar DEA, we compared the viscoelastic hysteresis loops with different frequencies (0.02, 0.1, 0.4, and 1 Hz). On the basis of the experimental results shown in Figure 5, we could observe the following: (iii) The peak-to-peak displacements of the viscoelastic hysteresis loops are rate-dependent. When the frequency of the exciting voltage increased from 0.02 to 1 Hz, the peak-to-peak displacement of the viscoelastic hysteresis loop decreased by 27.26%.
Polymers 2018, 10, x FOR PEER REVIEW 6 of 12 (iii) The peak-to-peak displacements of the viscoelastic hysteresis loops are rate-dependent. When the frequency of the exciting voltage increased from 0.02 to 1 Hz, the peak-to-peak displacement of the viscoelastic hysteresis loop decreased by 27.26%. The above observations demonstrate that the viscoelastic hysteresis loops of the planar DEA are both asymmetric and rate-dependent. As discussed above, physical-based models may be too complicated to describe such viscoelastic hysteresis effects. Therefore, a phenomenological modeling approach was proposed to solve this challenge. At the first step, a MRPIM was developed to characterize the asymmetric and rate-dependent viscoelastic hysteresis effect. Finally, the MRPIM was identified and verified.

Viscoelastic Hysteresis Model
Classical P-I models have been widely used to describe hysteresis effects mathematically [16,17,22] because of their analytical inversion. The drawback of classical P-I models is that they are generally effective only for symmetric and rate-independent hysteresis effect description. In order to (iii) The peak-to-peak displacements of the viscoelastic hysteresis loops are rate-dependent. When the frequency of the exciting voltage increased from 0.02 to 1 Hz, the peak-to-peak displacement of the viscoelastic hysteresis loop decreased by 27.26%. The above observations demonstrate that the viscoelastic hysteresis loops of the planar DEA are both asymmetric and rate-dependent. As discussed above, physical-based models may be too complicated to describe such viscoelastic hysteresis effects. Therefore, a phenomenological modeling approach was proposed to solve this challenge. At the first step, a MRPIM was developed to characterize the asymmetric and rate-dependent viscoelastic hysteresis effect. Finally, the MRPIM was identified and verified.

Viscoelastic Hysteresis Model
Classical P-I models have been widely used to describe hysteresis effects mathematically [16,17,22] because of their analytical inversion. The drawback of classical P-I models is that they are generally effective only for symmetric and rate-independent hysteresis effect description. In order to The above observations demonstrate that the viscoelastic hysteresis loops of the planar DEA are both asymmetric and rate-dependent. As discussed above, physical-based models may be too complicated to describe such viscoelastic hysteresis effects. Therefore, a phenomenological modeling approach was proposed to solve this challenge. At the first step, a MRPIM was developed to characterize the asymmetric and rate-dependent viscoelastic hysteresis effect. Finally, the MRPIM was identified and verified.

Viscoelastic Hysteresis Model
Classical P-I models have been widely used to describe hysteresis effects mathematically [16,17,22] because of their analytical inversion. The drawback of classical P-I models is that they are generally effective only for symmetric and rate-independent hysteresis effect description. In order to describe asymmetric and rate-dependent viscoelastic hysteresis effects of the DEA, we developed a MRPIM in this work.

The Classical P-I Model
The classical P-I model can be expressed as follows [22]: where the input signal V(t) is any continuous function on the interval [0, t E ], y(t) represents output displacement, p is a positive constant, a(r) is a density function, and F r [V](t) with a threshold r ≥ 0 is a one-side play operator that can be written as follows: where 0 = t 1 < t 2 < · · · < t N = t E is a partition of [0, t E ], such that the function V(t) is a monotone on each of the subintervals [t i , t i+1 ]. Because the polynomial and play operators are rate-independent and symmetric, the classical P-I model cannot describe the asymmetric and rate-dependent viscoelastic hysteresis effect shown in Figure 5. Therefore, it is necessary to develop a new rate-dependent P-I model to describe the viscoelastic hysteresis effect of DEAs.

The MRPIM
The MRPIM can be expressed as where y(t) represents output displacement; V(t) and ..
V(t) are the input signal and its acceleration, respectively; q is a constant ratio; and g[V(t)] = p 1 V(t) + p 2 V(t) 2 + p 3 V(t) 3 + p 4 V(t) 4 is a polynomial input function with four constants p 1 , p 2 , p 3 , p 4 . The fourth-order polynomial is used to describe the asymmetric behavior. In this work, the order number 4 was selected on the basis of a trial-and-error method, as is generally used in the literature [22,32]. It should be noted that further increasing the order of the polynomial (higher than 4) contributes a little improvement in accuracy but makes the computation costlier. In the following development, we selected the order of the polynomial as 4 for a case study; a r represents a density function, and F h or [V](t) is the rate-dependent play operator with a constant threshold r that is defined as follows: where h l V(t), . V(t) and h r V(t), . V(t) are the dynamic envelope functions that depend on the input voltage V(t) and its derivative . V(t). On the basis of experimental observation, the dynamic envelope functions can be expressed as follows: where γ and λ are two constants. For the convenience of calculation, ]/T, and T is the sampling time.
As usual, the upper limit of the integration is infinite; it is not convenient to identify the parameters of the model, and thus we used a discrete form of the model to replace the integration. The discrete form of the MRPIM can be written as where N is the number of play operators, a i is the weighting factor of the ith play operator with a constant threshold value r i , and r i is equal to

Identification of the MRPIM
To experimentally verify the effectiveness of our development, the first step was to identify the parameters of the MRPIM, which were obtained by a modified particle swarm optimization (MPSO) algorithm [33]. In this work, the fitness function of the MPSO algorithm was selected as where m represents the number of the frequency; n m and k m are the values of the input signal and the weight ratio of ith frequency, respectively; y mi and y a mi are the predicted results and experimental data at the ith sampling time, respectively; and x = {γ, λ, p 1 , p 2 , p 3 , p 4 , q 1 , a 1 , a 2 , · · · , a 13 } is a vector that contains all identified parameters. Table 1 lists the parameters' values of the fitness function. MATLAB software was used to perform the identification process, and the identification results are shown in Table 2.   Figure 6a shows the comparisons between the experimental data (black lines) and the identified results (red lines). The identification errors are plotted in Figure 6b. We can see from Figure 6 that the identified MRPIM could accurately describe the hysteresis loops under different frequencies.
To further quantify the performance of the MRPIM, the maximum error e m and the root-mean-square error e rms are defined, respectively, as follows: where y(i) and y a (i) represent predicted results and experimental data, respectively. According to the definitions, the maximum identification error and root-mean-square identification error are listed in Table 3. We found that the maximum identification error and root-mean-square identification error were 9.03% and 4.50%, respectively.
Polymers 2018, 10, x FOR PEER REVIEW 9 of 12 Figure 6a shows the comparisons between the experimental data (black lines) and the identified results (red lines). The identification errors are plotted in Figure 6b. We can see from Figure 6 that the identified MRPIM could accurately describe the hysteresis loops under different frequencies.
To further quantify the performance of the MRPIM, the maximum error m e and the root-meansquare error rms e are defined, respectively, as follows:  Table 3. We found that the maximum identification error and root-mean-square identification error were 9.03% and 4.50%, respectively.  To further verify the effectiveness of the MRPIM, we used it to predict the viscoelastic hysteresis effect of the planar DEA with three different frequencies (0.05, 0.2, and 0.8 Hz) and compared the predicted results with experimental data. Figure 7 and Table 4 demonstrate that the MRPIM agreed well with the experimental data. The maximum prediction error and the root-mean-square prediction error were 6.95% and 3.25%, respectively. Therefore, the MRPIM can precisely describe the viscoelastic hysteresis effect of the planar DEA when the frequency of the exciting voltage is in the range of 0.02-1 Hz.  To further verify the effectiveness of the MRPIM, we used it to predict the viscoelastic hysteresis effect of the planar DEA with three different frequencies (0.05, 0.2, and 0.8 Hz) and compared the predicted results with experimental data. Figure 7 and Table 4 demonstrate that the MRPIM agreed well with the experimental data. The maximum prediction error and the root-mean-square prediction error were 6.95% and 3.25%, respectively. Therefore, the MRPIM can precisely describe the viscoelastic hysteresis effect of the planar DEA when the frequency of the exciting voltage is in the range of 0.02-1 Hz.

Conclusions
In this work, we propose a phenomenological modeling approach to describe the viscoelastic hysteresis effect of a planar DEA with a MRPIM. Firstly, we investigated the viscoelastic responses of the planar DEA on the basis of different experiments. The experimental results show that (i) the viscoelastic hysteresis loops are asymmetric; (ii) when the frequency of the exciting voltage increased from 0.02 to 1 Hz, the width of the viscoelastic hysteresis loop increased by 21.80%, while the peak-

Conclusions
In this work, we propose a phenomenological modeling approach to describe the viscoelastic hysteresis effect of a planar DEA with a MRPIM. Firstly, we investigated the viscoelastic responses of the planar DEA on the basis of different experiments. The experimental results show that (i) the viscoelastic hysteresis loops are asymmetric; (ii) when the frequency of the exciting voltage increased from 0.02 to 1 Hz, the width of the viscoelastic hysteresis loop increased by 21.80%, while the peak-to-peak displacement decreased by 27.26%. Thus, the MRPIM was developed, and the model parameters were identified by a MPSO algorithm. Finally, comparisons between experimental results and MRPIM simulation were performed to verify the effectiveness of the development.