Forecast for Artificial Muscle Tremor Behavior Based on Dynamic Additional Grey Catastrophe Prediction

Recently, bio-inspired artificial muscles based on ionic polymers have shown a bright perspective in engineering and medical research, but the inherent tremor behavior can cause instability of output response. In this paper, dynamic additional grey catastrophe prediction (DAGCP) is proposed to forecast the occurrence time of tremor behavior, providing adequate preparation time for the suppression of the chitosan-based artificial muscles. DAGCP constructs various dimensions of time subsequence models under different starting points based on the threshold of tremor occurrence times and peak-to-peak values in unit time. Next, the appropriate subsequence is selected according to grey correlation degree and prediction accuracy, then it is updated with the newly generated values to achieve a real-time forecast of forthcoming tremor time. Compared with conventional grey catastrophe prediction (GCP), the proposed method has the following advantages: (1) the degradation of prediction accuracy caused by the immobilization of original parameters is prevented; (2) the dynamic input, real-time update and gradual forecast of time sequence are incorporated into the model. The experiment results show that the novel DAGCP can predict forthcoming tremor time earlier and more accurately than the conventional GCP. The generation mechanism of tremor behavior is illustrated as well.


Introduction
Artificial muscle, as one of the most promising smart materials, has attracted great attention for its potential applications in intelligent robots, artificial organs and biomedical devices in recent decades [1][2][3][4][5]. Due to the impressive characteristics of light weight, high flexibility, super agility and long durability [6][7][8][9][10], plentiful research results have been reported regarding the properties of artificial muscle. Li et al. [11] proposed a novel cantilever beam artificial muscle using single-walled carbon nanotubes, which showed superfast response and ultrahigh mechanical output power density. Lu et al. [12] put forward an artificial muscle based on multi-walled carbon nanotubes, ionic liquids and biopolymer chitosan, performing an excellent bio-compatibility. Jager et al. [13] studied the electrically-controllable characteristics of conjugated polymers artificial muscles and applied them to the drug injection therapy. Then, Escudero et al. [14] designed a magnetically sensitive actuator with a ball screw and rare-earth magnet coreless, which was compact enough to be used in practical artificial prosthesis, and the limitations of high-price and lower-power in conventional electrochemical polymer actuations were overcome. Recently, Kim and Kwon [15] developed a hybrid muscle powered

Materials and Fabrication of Artificial Muscle
Chitosan powder (deacetylation degree 85%) was purchased from Jingchun biochemical technology limited company (Shanghai, China). Dilute acid (HA) (Hengxing chemical reagents limited company, Tianjin, China), multi-walled carbon nanotubes aqueous dispersion (Boyu material technology limited company, Beijing, China), glycerin (Dongli chemical reagent factory of Tianjin University, Tianjin, China), etc. were applied as the ingredients.
The basic processes of muscle include the preparation of actuation film, the synthesis of electrode film and hot embossing course as shown in Figure 1. Chitosan powder was added to 2% dilute acid solution, and the mixture was stirred for 30 min at 60 • C water bath, then multi-walled carbon nanotubes (MCNTs) aqueous dispersion was gradually added into the above mixture for another 15 min. After ultrasonic oscillation, a homogeneous electrode solution was obtained and poured into a square mold (50 mm × 50 mm × 4 mm). The electrode film is formed in an oven at 80 • C for 6 h. Meanwhile, the actuation solution was prepared through appending chitosan powder and glycerin into the HA solution, and the mixture was stirred for 15 min at 60 • C. After ultrasonic oscillation, the solution was poured into the mold and placed in the oven at 80 • C for about 4 h. Finally, the electrode film and actuation film were fabricated by the hot-pressing method for 20 min at 30 • C.

Measurement Methods and Bending Mechanism
The experiment system and deflection displacement model were established for testing the deflection displacements of the prepared artificial muscle samples, as shown in Figures 2 and 3, respectively. Under the effect of direct current (DC) voltage, the muscle sample achieved cathode deflection by overcoming the inherent bending stress, and real-time measurements d were collected from a laser displacement sensor device. After enough data was collected, a curve of time dependent displacement was drawn.  As a polymer material with good chemical stability, its side chains distribute a large number of amino (-NH2) after taking off the acetyl from the 40% sodium hydroxide solution. In general, dilute acids (HA) with concentrations of 2% are used as solvents. Figure 4 presents the schematic diagram of artificial muscle motion mechanism. During hydrolysis, several free amino groups (-NH2) from the molecular side chains are combined with hydrogen ions (H + ) in solvent to form charged polyelectrolyte (-NH3 + ); free anions (-A − ) are in a scattered state. Available dilute anions (-A − ) include

Measurement Methods and Bending Mechanism
The experiment system and deflection displacement model were established for testing the deflection displacements of the prepared artificial muscle samples, as shown in Figures 2 and 3, respectively. Under the effect of direct current (DC) voltage, the muscle sample achieved cathode deflection by overcoming the inherent bending stress, and real-time measurements d were collected from a laser displacement sensor device. After enough data was collected, a curve of time dependent displacement was drawn.

Measurement Methods and Bending Mechanism
The experiment system and deflection displacement model were established for testing the deflection displacements of the prepared artificial muscle samples, as shown in Figures 2 and 3, respectively. Under the effect of direct current (DC) voltage, the muscle sample achieved cathode deflection by overcoming the inherent bending stress, and real-time measurements d were collected from a laser displacement sensor device. After enough data was collected, a curve of time dependent displacement was drawn.  As a polymer material with good chemical stability, its side chains distribute a large number of amino (-NH2) after taking off the acetyl from the 40% sodium hydroxide solution. In general, dilute acids (HA) with concentrations of 2% are used as solvents. Figure 4 presents the schematic diagram of artificial muscle motion mechanism. During hydrolysis, several free amino groups (-NH2) from the molecular side chains are combined with hydrogen ions (H + ) in solvent to form charged polyelectrolyte (-NH3 + ); free anions (-A − ) are in a scattered state. Available dilute anions (-A − ) include

Measurement Methods and Bending Mechanism
The experiment system and deflection displacement model were established for testing the deflection displacements of the prepared artificial muscle samples, as shown in Figures 2 and 3, respectively. Under the effect of direct current (DC) voltage, the muscle sample achieved cathode deflection by overcoming the inherent bending stress, and real-time measurements d were collected from a laser displacement sensor device. After enough data was collected, a curve of time dependent displacement was drawn.  As a polymer material with good chemical stability, its side chains distribute a large number of amino (-NH2) after taking off the acetyl from the 40% sodium hydroxide solution. In general, dilute acids (HA) with concentrations of 2% are used as solvents. Figure 4 presents the schematic diagram of artificial muscle motion mechanism. During hydrolysis, several free amino groups (-NH2) from the molecular side chains are combined with hydrogen ions (H + ) in solvent to form charged polyelectrolyte (-NH3 + ); free anions (-A − ) are in a scattered state. Available dilute anions (-A − ) include

Tremor Behavior Analysis
In this work, the tremor behavior of muscles refers to the rhythmic and involuntary contraction characterized by oscillations of the muscle. It is an urgent issue for improving the stability of output response.
The surface microscopic structures of actuation film samples are illustrated in Figure 5. It can be seen that the overall structure of the actuation film is complete and compact without cracks, which can ensure expected muscle deflection. However, there are many obvious wrinkled areas and spatial agglomerations after it was magnified 5000 times. It will affect the normal movement of ions and cause the uneven distribution of electrode ions, resulting in unstable output response, namely tremor behavior. Based on the above bending mechanism, time dependent displacement data of sample 1 (40 mm × 5 mm × 0.5 mm) was collected at 3 V voltage as shown in Figure 6. It is easy to discover that the changes in displacement are fluctuating, but the general trend is monotonically rising within the time period. Regarding the displacement, if it decreases suddenly and then reverts to its last increasing state and the maximum peak-to-peak value or duration value of this wave is no less than a specific

Tremor Behavior Analysis
In this work, the tremor behavior of muscles refers to the rhythmic and involuntary contraction characterized by oscillations of the muscle. It is an urgent issue for improving the stability of output response.
The surface microscopic structures of actuation film samples are illustrated in Figure 5. It can be seen that the overall structure of the actuation film is complete and compact without cracks, which can ensure expected muscle deflection. However, there are many obvious wrinkled areas and spatial agglomerations after it was magnified 5000 times. It will affect the normal movement of ions and cause the uneven distribution of electrode ions, resulting in unstable output response, namely tremor behavior.

Tremor Behavior Analysis
In this work, the tremor behavior of muscles refers to the rhythmic and involuntary contraction characterized by oscillations of the muscle. It is an urgent issue for improving the stability of output response.
The surface microscopic structures of actuation film samples are illustrated in Figure 5. It can be seen that the overall structure of the actuation film is complete and compact without cracks, which can ensure expected muscle deflection. However, there are many obvious wrinkled areas and spatial agglomerations after it was magnified 5000 times. It will affect the normal movement of ions and cause the uneven distribution of electrode ions, resulting in unstable output response, namely tremor behavior. Based on the above bending mechanism, time dependent displacement data of sample 1 (40 mm × 5 mm × 0.5 mm) was collected at 3 V voltage as shown in Figure 6. It is easy to discover that the changes in displacement are fluctuating, but the general trend is monotonically rising within the time period. Regarding the displacement, if it decreases suddenly and then reverts to its last increasing state and the maximum peak-to-peak value or duration value of this wave is no less than a specific Based on the above bending mechanism, time dependent displacement data of sample 1 (40 mm × 5 mm × 0.5 mm) was collected at 3 V voltage as shown in Figure 6. It is easy to discover that the changes in displacement are fluctuating, but the general trend is monotonically rising within the time period. Regarding the displacement, if it decreases suddenly and then reverts to its last increasing state and the maximum peak-to-peak value or duration value of this wave is no less than a specific value, the tremor of artificial muscle is considered to happen. In particular, tremor behaviors are marked in circles. In addition, the thresholds of peak-to-peak value and duration are set as 0.2 mm and 1.5 s, respectively.
Appl. Sci. 2018, 8, x 5 of 12 marked in circles. In addition, the thresholds of peak-to-peak value and duration are set as 0.2 mm and 1.5 s, respectively. The displacement data in the first 65 s is selected as the original sequence and the last 20 s is used to test the accuracy of prediction. Table 1 shows the original and testing sequence of tremor occurrence times and the maximum peak-to-peak values with a time unit of 5 s, respectively. Without loss of generality, if no tremor occurrs in a specific time unit, the peak-to-peak value is set as 0 mm. The general idea of grey prediction is to establish the first-order differential equation about generated time sequence named GM (1, 1) model and forecast the occurrence time of forthcoming abnormal values based on the inherent regularity [32]. In particular, the time sequence is divided into equal time intervals. The threshold for determining abnormal values mainly depends on the actual situation and experiment methods when implementing the catastrophe prediction.

Conventional GCP Method
By specifying a threshold, the time sequence data of tremor behavior from 0 s to 65 s is considered as the original time sequence , where Q (0) is a nonnegative sequence and n The displacement data in the first 65 s is selected as the original sequence and the last 20 s is used to test the accuracy of prediction. Table 1 shows the original and testing sequence of tremor occurrence times and the maximum peak-to-peak values with a time unit of 5 s, respectively. Without loss of generality, if no tremor occurrs in a specific time unit, the peak-to-peak value is set as 0 mm. The general idea of grey prediction is to establish the first-order differential equation about generated time sequence named GM (1, 1) model and forecast the occurrence time of forthcoming abnormal values based on the inherent regularity [32]. In particular, the time sequence is divided into equal time intervals. The threshold for determining abnormal values mainly depends on the actual situation and experiment methods when implementing the catastrophe prediction.

Conventional GCP Method
By specifying a threshold, the time sequence data of tremor behavior from 0 s to 65 s is considered as the original time sequence Q (0) = (q(1), q(2), · · · , q(n)), where Q (0) is a nonnegative sequence and n is the number of the outliers.
The first-order accumulative generation operation (1-AGO) is given by Then, the generated mean sequence Z (1) of Q (1) is defined as: The least square estimate of GM (1, 1) model is defined as: where a and b represent the development coefficient and the grey action of the prediction model, respectively. Finally, the estimated equation of the training sequence can be expressed aŝ

DAGCP Method
The whole process of dynamic additional prediction is completed as follows:

Establishment of different time subsequences
The GM (1, 1) model typically requires at least four points, so the maximum starting point should be set to the (n − 4)th observation to ensure that there are no less than four points in the shortest subsequence. The subsequence size of n − 4 can be established at the moment. Abnormal values time sequence and time subsequences are expressed as: Q (0) = (q(1), q(2), · · · , q(i), · · · , q(n)) where n is the number of the outliers as well as the dimension of original sequence and i is the ordinal number of calculation starting point. For example, i = 1 represents that the distances between the first and subsequent other points are computed to obtain a new data subsequence under the dimension of n − 1. Normally, Response equation of catastrophe time subsequences can be obtained:

Selection of time subsequence
The optimal time subsequence is selected based on the grey correlation degree and model precision. Essentially, grey correlation degree is applied to analyze the model's fitting degree by making a comparison between differential subsequences and original sequence in terms of the approximate level. As known, the fitting degree is improved with the increase of grey correlation value. The relevant equation about grey correlation degree is shown as follows: where ξ i (k) is grey correlation coefficient and it can be calculated by where ρ is the distinguishing coefficient of the correlation axiom and the value is set as 0.5 based on previous experience. If the grey correlation degree is greater than 0.6, the model is considered as a satisfaction. Accordingly, the appropriate subsequence is picked out according to the model precision and correlation degree.

Achievement of additional forecast
Conventional GCP adopts a fixed model to make a forecast for multiple times, which can continuously cause the degradation of prediction precision owing to the immobilization of inherent parameters. In view of the immediacy and contingency in tremor behavior, the newly generated predicted value is added into the original subsequence for future prediction purpose after rounding down decimal places. The decimal places are omitted so as to be consistent with the original nonnegative integer data sequence.

Examination of Prediction Accuracy
To analyze the reliability of DAGCP, the accuracy of models has to be tested. The data tested comes from the China's hydropower production from 2000 to 2015, and prediction results are given by Wang et al. [25] by seasonal autoregressive integrated moving average (SARIMA) method and GM (1, 1) are provided for a direct comparison (see Figure 7). It can be seen that the present method is in close agreement with the actual values and SARIMA. Thus, compared with conventional methods, DAGCP shows a fine prediction precision and reliability. The reason for excellent prediction performance is that various dimensions of time subsequences constructed based on different starting points weakens the effects of the initial value, and the selected subsequence with continuously updated values captures time-varying performance of the system. Next, this method will be used for the forecast of tremor occurrence time.

The Investigation of Tremor Times
During the response process of artificial muscle, tremor occurrence times in each time unit are extracted, which characterizes the occurrence frequency of the tremor, and can be used as an important index for the evaluation of tremor behavior.
By referencing Table 1, the lower limiting outlier of occurrence times covered by tremor behavior was assigned as 1. Thus, the corresponding tremor data sequence within the first 65 s was: (2,5,6,7,8,9,11,12,13) When i = 1, the distances between the first and subsequent other points were computed to obtain a new data subsequence: Q (0) (1) = (3,4,5,6,7,9,10,11)  And then i = 2, the next subsequence was given: (1,2,3,4,6,7,8) To analyze the reliability of DAGCP, the accuracy of models has to be tested. The data tested comes from the China's hydropower production from 2000 to 2015, and prediction results are given by Wang et al. [25] by seasonal autoregressive integrated moving average (SARIMA) method and GM (1, 1) are provided for a direct comparison (see Figure 7). It can be seen that the present method is in close agreement with the actual values and SARIMA. Thus, compared with conventional methods, DAGCP shows a fine prediction precision and reliability. The reason for excellent prediction performance is that various dimensions of time subsequences constructed based on different starting points weakens the effects of the initial value, and the selected subsequence with continuously updated values captures time-varying performance of the system. Next, this method will be used for the forecast of tremor occurrence time.

The Investigation of Tremor Times
During the response process of artificial muscle, tremor occurrence times in each time unit are extracted, which characterizes the occurrence frequency of the tremor, and can be used as an important index for the evaluation of tremor behavior.
By referencing Table 1, the lower limiting outlier of occurrence times covered by tremor behavior was assigned as 1. Thus, the corresponding tremor data sequence within the first 65 s was:  Figure 8a-e, respectively. Meanwhile, model response, model accuracy, grey correction degree and fitting degree were calculated to evaluate the availability of the candidates (Table 2).
Analogously, based on the dimension of the original sequence, five new data subsequences were constructed under dimensions from 4 to 8 as the candidate models. The estimated values and original values of each subsequence with respect to ordinal k were shown in Figure 8a-e, respectively. Meanwhile, model response, model accuracy, grey correction degree and fitting degree were calculated to evaluate the availability of the candidates (Table 2).  Table 2 and Figure 8, the results showed that the 4-D model had the highest prediction accuracy; however, it was not available for prediction because its grey correlation degree was less than 0.6 and the number of data samples was too small. Additionally, the lower correlation degree and accuracy of 5-D and 6-D models limited their applications in the prediction. Although both of them had high accuracy at the initial stage, compared with the 7-D subsequence, the model accuracy and grey correlation degree of 8-D were better as a whole, which showed a better fitting degree between the subsequence and original sequence. Therefore, the 8-D subsequence model about tremor occurrence times in unit time was used for prediction and the starting time point q(1) was set as 2. From Table 2 and Figure 8, the results showed that the 4-D model had the highest prediction accuracy; however, it was not available for prediction because its grey correlation degree was less than 0.6 and the number of data samples was too small. Additionally, the lower correlation degree and accuracy of 5-D and 6-D models limited their applications in the prediction. Although both of them had high accuracy at the initial stage, compared with the 7-D subsequence, the model accuracy and grey correlation degree of 8-D were better as a whole, which showed a better fitting degree between the subsequence and original sequence. Therefore, the 8-D subsequence model about tremor occurrence times in unit time was used for prediction and the starting time point q(1) was set as 2. The estimated values and relative errors of conventional GCP and DAGCP were shown in Table 3. Specifically, the conventional GCP was constructed by the original time sequence. It is obvious that the actual number of tremors only occurred once during 75~80 s and was affected by process parameters, experimental conditions and human factors, as well as existing contingency. Hence, it was not necessary to predict ahead in this circumstance.
It was obvious that the occurrence times in 65~70 s, 70~75 s and 80~85 s were 2, 3 and 2, respectively. As shown in Table 3, there was a little difference between GCP and DAGCP in the relative error during 65~70 s. However, the relative error of DAGCP was smaller than GCP after 70 s. Therefore, the prediction accuracy of DAGCP was superior to that of GCP as a whole.

The Investigation of Peak-to-Peak Values
The maximum peak-to-peak values (the difference between adjacent crests and troughs) in unit time are extracted, which characterizes the severe degree of tremor, and also can be used for the evaluation of tremor behavior. Here, the lower limiting outlier was assigned as 0.2 mm. In addition, the corresponding data sequence within 65 s was: (2,5,6,7,8,11,12) When i = 1, a new data subsequence was: Q (0) (1) = (3,4,5,6,9,10). Thus, three new data subsequences were constructed under the dimensions from 4 to 6.
The estimated values and original values were shown in Figure 9a-c, respectively. In addition, the model response, model accuracy, grey correction degree and fitting degree were displayed in Table 4.  Table 4 and Figure 9, although the grey correlation degree had achieved the basic requirement, the prediction precision of 4-D and 5-D models limited their applications in the prediction. Among them, 6-D model had the highest model accuracy and grey correlation degree, which showed a more accurate fitting degree. Thus, the 6-D model about tremor peak-to-peak value in unit time was used for the prediction and q(1) was also set as 2.
The estimated values and relative errors were shown in Table 5. The actual peak-to-peak value was less than 0.2 mm during 75~80 s, which was allowed to happen. Therefore, it was not necessary to predict ahead in this circumstance.  In the above analysis, the peak-to-peak values in 65~70 s, 70~75 s and 80~85 s were 0.6 mm, 0.5 mm and 0.2 mm, respectively. As shown in Table 5, the relative error of DAGCP was obviously less than that of GCP.
To validate the effectiveness and accuracy of the proposed method, the catastrophe prediction results of 10 groups of chitosan-based artificial muscles with different process parameters based on GCP and DAGCP were shown in Table 6. Here, the data in the first 70 s was selected as original sequence and the last 10 s was used to assess the prediction accuracy. The lower thresholds of occurrence times and peak-to-peak values were assigned as 2 and 0.2 mm, respectively. If either of the conditions was satisfied, the tremor is considered to happen.   In the above analysis, the peak-to-peak values in 65~70 s, 70~75 s and 80~85 s were 0.6 mm, 0.5 mm and 0.2 mm, respectively. As shown in Table 5, the relative error of DAGCP was obviously less than that of GCP.
To validate the effectiveness and accuracy of the proposed method, the catastrophe prediction results of 10 groups of chitosan-based artificial muscles with different process parameters based on GCP and DAGCP were shown in Table 6. Here, the data in the first 70 s was selected as original sequence and the last 10 s was used to assess the prediction accuracy. The lower thresholds of occurrence times and peak-to-peak values were assigned as 2 and 0.2 mm, respectively. If either of the conditions was satisfied, the tremor is considered to happen. From the comparative analysis, the proposed DAGCP can better identify the occurrence tendency of tremor behavior and individual variation of original data, and achieve smaller relative errors. Therefore, this novel method is viable and reliable for predicting the tremor behavior of artificial muscle.

Conclusions
In this paper, we illustrate the tremor mechanism of artificial muscle and forecast the tremor occurrence time using GCP based methods. In order to further improve the prediction accuracy, we propose a new method called DAGCP by integrating sequence optimization mechanism and dynamic additional strategy into the conventional GCP. Its excellent prediction performance indicates that various dimension of time subsequences constructed based on different starting points weakens the effects of the initial value, and the selected subsequence with continuously updated values captures time-varying performance of the system.
As shown in the case study of artificial muscle, the novel DAGCP can earlier and more accurately forecast the forthcoming tremor time than conventional GCP. Therefore, it can be selected to provide the preparation time for parameter adjustment of artificial muscles' online behaviors and realize the suppression of tremor behavior. In summary, this method has the following advantages: the degradation of prediction accuracy caused by the immobilization of original parameters is prevented; dynamic input, real-time update and gradual forecast of the model can be realized and widely used in the industry, agriculture and medicine, etc.