# Inverse Modeling of Nonlinear Artificial Muscle Using Polynomial Parameterization and Particle Swarm Optimization

Authors: Mohd Azuwan Mat Dzahir; Shin-ichiroh Yamamoto
Journal: Advances in Materials Science and Engineering
 (2020)

---

## Abstract

The properties of pneumatic artificial muscle (PAM) with excellent power-to-weight ratio and natural compliance made it useful for healthcare engineering applications. However, it has undesirable hysteresis effect in controlling a robotic manipulator. This behavior is quasistatic and quasirate dependent which changed with excitation frequency and external force. Apart from this, it also inherits frictional presliding behavior with nonlocal memory effect. These nonlinearities need to be compensated to achieve optimal performance of the control system. Even though an inverse modeling of PAM has limited application, it is important on certain control system implementation that requires the solution to the inverse problem. In this paper, the inverse modeling of PAM in the form of activation pressure was proposed. This activation pressure model was derived according to static pressure and extracted hysteresis components from pressure/length hysteresis. The derivation of the static pressure model follows the phenomenological-based model of third-order polynomial. It is capable of characterizing the nonlinear region of PAM at low and high pressure. The derivation of extracted hysteresis model follows the mechanism of dynamic friction. In this principle, the activation pressure model was dependent on regression coefficient of the static pressure model and dynamic friction coefficients of the extracted hysteresis model. The regression constants of these coefficients were characterized from the hysteresis dataset by using model parameter identification and the particle swarm optimization (PSO) method. The result from model simulation shows the root mean square error (RMSE) value of less than 10% error was evaluated at various excitation frequencies and external forces. This inverse modeling of PAM implemented a simple approach, but it should be useful in control design applications such as rehabilitation robotics, biomedical system, and humanoid robots.


---

## Body

## 1. Introduction

The configuration of PAM enables its pneumatic power to be converted into mechanical power with high power-to-weight ratio compared to rigid actuators. During inflation, the radial direction of PAM is enlarged, which caused its longitudinal direction to shorten due to contraction of braid muscle [1, 2]. This retraction strength is dependent on the total strength of individual fibers of braid muscle which produces unidirectional tensile force during the contraction of PAM. In exchange for its excellent power generation and natural compliance, it is known to have major drawbacks such as hysteresis, time variance, and pressure dynamics [3–5]. These drawbacks can cause major problem when designing a high-performance control system. The hysteresis of PAM increases the level of nonlinearities and leads to the complexity of the associated control system [6–8]. Many efforts have been delivered for modeling PAM, and several approaches have been implemented to compensate for the nonlinear behaviors [9–14].According to the literature, there are two common approaches widely used in modeling the nonlinear behaviors of PAM. The first is a physics-based model, and it is established based on the principles of physics such as geometric analysis, virtual work theorem, energy conservation, and empirical principle. In this approach, three models had received worldwide recognition in modeling of PAM. The first model was derived based on the empirical model introduced by Inoue et al. in 1988 [15]. In this model, the pressure force of PAM is separated into lateral and axial pressure forces. This model was modified by Tondu et al. in 1994 with shape irregularity qP coefficient which is dependent on the internal pressure of PAM to account for shape irregularity of the muscle [16, 17]. Later, the extended model was proposed by Kerscher et al. in 2006 by introducing shape slenderness με coefficient to the existing model [18]. The introduction of this coefficient was made to consider the nonlinear region of PAM at low- and high-pressure activation. The second model was derived by Chou and Hannaford in 1994–1996 using the geometrical analysis model [19, 20]. According to this model, it is stated that the pressure force of PAM can also be expressed by applying elementary volume variation. This model was equivalent to empirical modeling of PAM. The third model was derived by Caldwell et al. in 1995 which introduced the PAM model based on theorem of virtual work [21, 22]. The improved version of this model was made by Tsagarakis and Caldwell in 2000 by considering the effects of end-cap size/diameter on the PAM model [23]. Later, Davis and Caldwell in 2006 proposed another modification by considering braid effects on contractile range and account for its effect on friction modeling [24]. This model is identical to the model obtained by Chou and Hannaford through geometrical analysis [19]. However, they do not consider the end-cap surface area effect on modeling of PAM. The result with minimum error was reported in the literature which implemented the law of physics, and its performance can be further improved with an addition of the reliable intelligent control system.The second is phenomenological-based model, and it is not directly derived from first principles of physics. This model is a scientific model that describes the empirical relationship in a way which is consistent with fundamental theory and follows the phenomena of PAM. There are different type of phenomenological-based models such as operator based, differential equation based, data driven, and polynomial parameterization. Recently, the research has been shifted to the polynomial-based model in modeling of PAM. It has received many attentions from research experts in this field [25–29]. It was first intrigued by Petrovic et al. in 2002 which considered the implementation of the second-order polynomial model to describe the behavior of the static force model [30]. It was chosen because of the simplicity of the second order adjustment. This model was modified by Arrese et al. in 2010 with introduction of fourth-order correction term φq of the force element. This correction force element is dependent on displacement/contraction of PAM [31, 32]. Later, Tang et al. in 2012 proposed a similar model, but the coefficient of force element was reduced to second-order polynomial function [33]. It accounts for error caused by the thickness of the shell and the bladder. It is reported that the static model of second-order polynomial with correction force element provides accurate force estimation within the linear region of PAM. However, it was not able to characterize the asymmetricity of PAM at the nonlinear region. In response to this issue, Minh et al. in 2011 proposed a piecewise linearization method of static force model [34] and static pressure model [35] using the third-order polynomial model. They also included coefficient of force element and pressure element with third-order polynomial function. The result shows significant improvement at the nonlinear region of PAM, where it was able to characterize the asymmetric behavior of the PAM. Later, Schindele and Aschemann in 2014 proposed further improvement on the similar static force model by increasing the coefficient of force element to fourth-order polynomial function [36–38]. However, this increment was not significant as its accuracy was slightly improved.According to the literature, most of the studies were focusing on modeling of the contracting force model of the PAM. Not many research studies can be found which characterized the activation pressure of the PAM as reported in [11, 34, 39, 40]. In this study, an inverse modeling of PAM is described in the form of the activation pressure model. It utilizes a phenomenological-based model using the model parameter and system identification method. The model is theoretically derived according to activation pressure of the pressure/length hysteresis dataset. This model is dependent on both the static pressure and extracted hysteresis models. The inverse modeling is composed of three major sections including mathematical modeling, experimental dataset, and model simulation. In the mathematical modeling section, the proposed activation model is described using mathematical formulation according to polynomial regression and dynamic friction mechanism. In the experimental dataset section, the static pressure and extracted hysteresis of the hysteresis dataset are characterized using the analytical method. The regression coefficient of the static pressure model and dynamic friction coefficients of the extracted hysteresis model are identified using model parameter identification and the PSO method. These parameters were stored as data profile in the form of regression constants. In model simulation section, the proposed activation pressure model is evaluated using random wave input signal and then compared with the actual hysteresis dataset from the vertical load system with single PAM. This modeling approach should be useful in control design application that requires inverse model of PAM.


## 2. Materials and Methods

The inverse modeling of PAM is composed of mathematical modeling, experimental dataset, and model simulation sections. The mathematical modeling section comprised two parts. The first part involved derivation of the mathematical model for static pressure component using a phenomenological-based model of third-order polynomial function. The second part involved derivation of the mathematical model for extracted hysteresis components using the mechanism of dynamic friction. The experimental dataset section comprised four parts. The first part involved collection of the hysteresis dataset using the developed vertical load system with single PAM. The second part involved characterization of primary data (i.e., global minimum, global maximum, inflection, and tangent parameters) from the hysteresis dataset using the analytical method. The third part involved model parameter identification of the static pressure model using data enrichment. It was performed by characterizing its regression coefficient from primary data. The fourth part involved system identification of the extracted hysteresis model using the PSO method. It was performed by characterizing the dynamic friction coefficient based on dynamic response estimation. The code was programmed using MATLAB software (MATLAB R2019b, Math Works Inc., academic use). The model simulation section comprised two parts. The first part involved evaluation of the proposed activation pressure model based on quasistatic, quasirate, and nonlocal memory behaviors. It was performed using the random wave input signal of the independent variable obtained from the hysteresis dataset at various excitation frequencies and external forces. The second part involved evaluation of performance using correlation coefficient and RMSE values. It was performed based on the comparison of the proposed activation pressure model and actual pressure/length hysteresis dataset.


## 3. Mathematical Modeling

This section contains mathematical formulation of the proposed activation pressure model using the theoretical method. It can be observed from the literature that most of the existing hysteresis models are initially developed to describe a particular type of hysteretic system. However, its mathematical formulation could be implemented for multidisciplinary extension [41]. In general, there is no fundamental theory that allows for a general mathematical framework in modeling the hysteresis behavior. This is because the origins of this phenomenon are often unclear and complicated to understand [5].Figure1(a) shows an illustration of measured hysteresis of PAM according to the isotonic test experimental setup. This figure represents the hysteresis dependency of the measured activation pressure when the PAM experienced cyclic contraction and expansion. This activation pressure can be decomposed into virgin curve, major hysteresis, and minor hysteresis components. The virgin curve describes the quasistatic behavior of PAM which is dependent on external force (F) of the system, while the major hysteresis describes the quasirate behavior of PAM which is dependent on the excitation frequency (ω) of the system. Meanwhile, the minor hysteresis describes the nonlocal memory behavior of the PAM which is dependent on the previous values of contraction, force, and pressure. The proposed activation pressure model is derived from the constrained setup and is expressed as equation (1). It is an inverse model of PAM in the form of activation pressure (pi) which is dependent on both static pressure (yi) and extracted hysteresis (zi) models.(1)pi=yi+zi.Figure 1
Hysteresis of PAM: (a) measured activation pressure and (b) activation pressure model.
(a)(b)The activation pressure (pi) describes the pressure dynamics of PAM which is influenced by multiple factors such as compressibility, time variance, pressure build up, and friction. The static pressure (yi) describes the interpolated pressure of a virgin curve which is also known as constraint model. The mathematical formulation of the proposed static pressure model uses a polynomial regression of third order, as shown in equation (2), where yi is the dependent variable vector, εi is the independent variable vector, and βn=0,1,2,3 is the parameter vector. It has been reported in the literature that the second-order polynomial is enough to characterize the nonlinear region of PAM; however, to consider the nonlinear region at low- and high-pressure activation, third-order polynomial function is required.(2)yi=β0+β1εi+β2εi2+β3εi3.In general, the static pressure model can also be written in the form of matrix-based notation, as shown in equation (3), where the matrix X is the polynomial regression of third order and vector β⟶ is the regression coefficient of the static pressure model (y⟶). According to the literature, the pressure of PAM depends on both contraction and external force given to the system. Therefore, the regression coefficient of the static pressure model is proposed as a function of external force (F).(3)y⟶=Xβ⟶.To estimate the regression coefficientβ⟶ of the static pressure model, it requires four critical points ({εj, yj}, for j=1,2,3,4) known as primary data. It is obtained from observation of the hysteresis dataset. The primary data are composed of global minimum (ε1,y1), global maximum (ε2,y2), inflection (ε3,y3), and tangent (ε4,y4) parameters, as shown in Figure 1(b). The system of linear equations of this primary data is shown in equation (4) with single covariate and four observations. Based on these observations, the estimator of regression coefficient β⟶^ for the static pressure model can be obtained using the general theory of the least squares method as described in equation (5).(4)y1y2y3y4=1ε1ε12ε121ε2ε22ε231ε3ε32ε331ε4ε42ε43β0β1β2β3,(5)β⟶^=XℛTXℛ−1XℛTy⟶ℛ.The regressor for the polynomial regressionXℛ can be estimated based on regression of contraction parameters from the primary data as shown in equation (6), while the regressor for the dependent variable y⟶ℛ can be estimated based on regression of pressure parameters from the primary data as shown in equation (7).(6)ε1ε2ε3ε4=a1b1c1d1a2b2c2d2a3b3c3d3a4b4c4d4F3F2F1,(7)y1y2y3y4=a5b5c5d5a6b6c6d6a7b7c7d7a8b8c8d8F3F2F1.The regressorε⟶ℛ of the independent variable εj=1,2,3,4 and regressor y⟶ℛ of the dependent variable yj=1,2,3,4 can be written in simpler form, as shown in equations (8) and (9), where U is the polynomial regression which constitutes external force (F) and coefficient matrices α and δ are the regression constants. These regression constants characterize the quasistatic behavior of the PAM which is dependent on external force (F).(8)ε⟶ℛ=αU,(9)y⟶ℛ=δU.The extracted hysteresis (zi) describes the compensation of pressure during cyclic contraction and expansion of PAM. Its value can be measured by subtracting static pressure from activation pressure of the pressure/length hysteresis dataset. The mathematical formulation of the proposed extracted hysteresis model is described based on a mechanism of dynamic friction, as shown in equation (10). In this mechanism, the viscous friction (ν) during cyclic contraction and expansion of PAM is considered.(10)zi=νε˙i.Referring to this mechanism, extended form of the dynamic friction model, as shown in equation (11), is proposed. In this extended form, the extracted hysteresis model (z⟶) depends on two dynamic friction coefficients (i.e., γ and φ). According to the literature, the hysteresis of PAM system inherits quasistatic and quasirate dependencies. Therefore, the dynamic friction coefficients of the extracted hysteresis model are proposed as a function of excitation frequency (ω) and external force (F). The contraction velocity (ε˙⟶) can be obtained by using the first derivative of contraction.(11)z⟶=γφε˙⟶.To estimate the dynamic friction coefficients of the extracted hysteresis model, it requires dynamic response identification in between rate-dependent coefficientγ and static-dependent coefficient φ based on shape estimation of the major hysteresis dataset. The estimator of rate-dependent coefficient γ^ and static-dependent coefficient φ^ is described as a simple polynomial regression, as shown in equations (12) and (13).(12)γ^=a9b9c9d9ω3ω2ω1,(13)φ^=a10b10c10d10F3F2FF.These dynamic friction coefficients of the extracted hysteresis model can be written in simple form, as shown in equations (14) and (15). The coefficient matrices μ and ρ are the regression constants obtained from dynamic response identification, while the matrix W is the polynomial regression which constitutes excitation frequency of the system. The identification of these regression constants will be difficult without a sophisticated identification method. This is because the extracted hysteresis inherits multiple dependencies (i.e., quasistatic, quasirate, and nonlocal memory). Therefore, it is not an easy task to determine the dynamic response in between these two dynamic friction coefficients without implementing the metaheuristic method. In this identification, the PSO method is proposed due to its capability in finding the best parameters for the extracted hysteresis model (z⟶).(14)γ^=μW,(15)φ^=ρU.In addition, it has been reported that the effect of quasistatic dependency is often neglected due to its insignificant error produced on the performance of the control system. In this study, the rate-dependent coefficientγ^ uses only a simple linear regression. Therefore, the regression constants of both a9 and b9 were neglected.(16)p⟶=y⟶+z⟶,(17)p⟶=Xβ⟶^+γ^φ^ε˙⟶.The proposed activation pressure model (p⟶) is dependent on both static pressure (y⟶) and extracted hysteresis (z⟶) models as can be seen in equation (16), where its final form of the mathematical model is shown in equation (17). The regression coefficient β⟶^ of static pressure model and dynamic friction coefficients γ^φ^ of the extracted hysteresis model can be characterized from the experimental data. This experimental data contain the hysteresis dataset of single PAM at various excitation frequencies (ω = 0.5, 1.0, 1.5, and 2.0 Hz) and external forces (F = 0, 50, 100, 150, 200, 250, and 300 N). The dynamic friction coefficient γ^ is a function dependent on the excitation frequency (ω), while both the regression coefficient β⟶^ and dynamic friction coefficient φ^ are a function dependent on external force (F). These coefficients can be estimated based on the regression constants (i.e., α, δ, μ, and ρ) obtained from model parameter identification and the PSO method.


## 4. Experimental Dataset

This section contains model parameters and system identification of the proposed activation pressure model using the analytical method. An experimental test was conducted using a developed vertical load system with single PAM, as shown in Figure2. This vertical load system consists of PAM, laser displacement sensor, pressure sensor, load sensor, accelerometer, weights, digital regulator, air compressor, processing unit, and data acquisition system. Its operating system was developed and programmed using LabVIEW software (Academic Software-Spring 2015, National Instruments Corporation, Ireland). The PAM used in this experiment is self-developed McKibben muscle that consists of rubber tube, braided mesh sleeve, copper-ring, O-ring, and end-connectors. It has an initial diameter of 10 mm and an original length of 530 mm.Figure 2
Experimental setup for the vertical load system with single PAM.


### 4.1. Data Collection

Data collection is divided into three parts. The first part was collecting data from the experimental setup with quasistatic dependency. It was conducted by applying steady increment of vertical load from 0 N to 300 N. The second part was collecting data from the experimental setup with quasirate dependency. It was conducted by varying the excitation frequency of input signal from 0.5 Hz to 2.0 Hz. The third part was collecting data from the experimental setup with nonlocal memory behavior. It was conducted using random wave input signal which generates minor hysteresis loops. These data are composed of the contraction and pressure parameters at various excitation frequencies and external forces. The structural materials of the PAM inherently lead to hysteresis phenomena during cyclic contraction and expansion, as can be seen in Figure3. This nonlinearity was addressed as difficult error to handle especially for precise control system implementation.Figure 3
Sample of hysteresis data at various frequencies and external force conditions.


### 4.2. Data Characterization

Data characterization is divided into three parts. The first part involved characterization of local extremum of the hysteresis dataset using third-order polynomial regression. It was performed by characterizing the minimum and maximum values consisting of contraction-pressure parameters from the hysteresis dataset at various frequencies and external forces, as shown in Figure4. This characterization describes the quasistatic dependency of PAM. The major hysteresis dataset was implemented for this characterization method.Figure 4
Identification of local maximum and local minimum of hysteresis data.The second part involved characterization of the global extremum from the identified local extremum at various conditions. These are the first two critical points (primary data) from the hysteresis dataset: the first is global minimum (ε1, y1) and the second is global maximum (ε2, y2) parameters. These primary data allow for generalization of the virgin curves at various excitation frequencies (i.e., 0.5, 1.0, 1.5, and 2.0 Hz) using the single polynomial model, as can be seen in Figure 5. In this figure, the local extremum was distributed along with the identified single polynomial model. This characterization revealed that hysteresis dataset of various excitation frequencies was dependent on the quasistatic load.Figure 5
Identification of global maximum and global minimum of hysteresis data.The third part involved characterization of remaining two critical points (primary data) based on the obtained polynomial model solution: the first is inflection (ε3, y3) and the second is tangent (ε4, y4) parameters. These primary data can be obtained by solving the second derivative of the identified polynomial model and then finding its intersection from the origin. The unoptimized static pressure model based on characterization of the hysteresis dataset is shown in Figure 6. In this stage, the polynomial model has not been estimated based on the identified parameters (primary data). However, it was generated based on the curve fitting technique using the integral of error value.Figure 6
Unoptimized static pressure model using polynomial parameterization of third order.


### 4.3. Model Parameter Identification

The identified parameters from characterization of the hysteresis dataset were recognised as primary data for model parameter identification of the static pressure model. It composed of two stages of data enrichment; the first enrichment was performed based on the regressor of independent variableε⟶ℛ and regressor of dependent variable y⟶ℛ from identified eight parameters (primary data). These regressors were dependent on the external force (F) of the system. Its regression constants (i.e., α and δ) were organized as the first part of the data profile. This first part of data profile that characterized the static pressure component of the proposed activation pressure model is described in equation (1). The polynomial regression (X) of the static pressure model (y⟶)is shown in equations (2) or (3). It consisted of identification of regression coefficient (β⟶) which is characterized by the behavior of quasistatic dependency of PAM. The mathematical formulation of the regression coefficient β⟶^ was estimated according to equation (5). This regression coefficient was dependent on primary data, as shown in equation (4). Figure 7 shows regression constants for contraction (α) and pressure (δ) which controls the behavior of the proposed static pressure model. In this stage, global extremum (i.e., minimum and maximum) of the identified polynomial model was first enriched using the primary global minimum (ε1, y1) and global maximum (ε2, y2) parameters. It reconstructs the regression constants based on the optimized regression value of the global extremum as described in equations (6) and (8). This data enrichment reduces the deviation of quasistatic error and brings the polynomial model closer to primary global extremum.Figure 7
Regression constants of the static pressure model: (a) contraction and (b) pressure.
(a)(b)The second enrichment was performed based on the nonlinear region of PAM. As reported in the literature, PAM behaves highly nonlinear at low- and high-pressure activation. In the proposed static pressure model, the inflection parameter reflects nonlinearity boundary of low-pressure activation while the tangent parameter reflects nonlinearity boundary of high-pressure activation. This second data enrichment brings both inflection and tangent parameters of the polynomial model closer to the actual nonlinear region of PAM at low- and high-pressure activation as described in equations (7) and (9). Figure 8 shows an optimized static pressure model based on two stages of data enrichment.Figure 8
Optimized static pressure model using polynomial parameterization of third order.The measurement of quasistatic error during the two stages of data enrichment process is shown in Figure9(a). It was determined based on the shortest distance equation between the static pressure model and primary global extremum. The source of this error comes from the estimation used during the curve fitting technique. During the estimation, the model was not generated based on the identified primary data from the hysteresis dataset. However, it was performed based on the integration of minimal error in between the polynomial model and major hysteresis dataset. Figure 9(b) shows the estimated regression coefficient β⟶^ of the static pressure model after optimization, where its regression constants (i.e., α and δ) control the quasistatic dependency of the activation pressure model. The identified parameters of static pressure model are shown in Table 1.Figure 9
Model optimization: (a) quasistatic error and (b) regression coefficient of the static pressure model.
(a)(b)Table 1
Identified parameters of the static pressure model.
iaibicidi1−4.33E−093.58E−06−1.11E−031.33E−012−1.99E−098.70E−07−4.50E−042.04E−013−4.78E−093.74E−06−1.01E−035.48E−0246.98E−09−3.03E−06−3.61E−043.03E−015−2.85E−091.19E−063.86E−046.48E−026−4.89E−09−1.41E−071.30E−039.46E−0273.96E−10−1.84E−072.53E−052.77E−0287.40E−117.25E−08−2.91E−053.91E−01


### 4.4. System Identification

Apart from quasistatic and quasirate dependencies, the dynamic model of PAM also exhibits nonlocal memory behavior. It is historic dependent which changes drastically according to presliding domain of friction. This behavior is evident and can be observed as the shape of the hysteresis dataset drastically changed. Theoretically, this follows the mechanism of dynamic friction of presliding regime as described in [42, 43]. In response to this mechanism, extended form of the dynamic friction model was proposed for inverse modeling of PAM. Its mathematical formulation is described in equations (10) and (11). In this extended form, the proposed extracted hysteresis model (z⟶) was dependent on two dynamic friction coefficients (i.e., γ and φ). It is difficult to determine the dynamic response in between two dynamic friction coefficients without the use of metaheuristic methods such as genetic algorithm and PSO. There is no significant difference between these two methods in small scale; however, differences are seen in medium and large scale where genetic algorithms can only produce feasible solutions that are nearly optimal, while the PSO algorithm has ease of implementation and also has high calculation accuracy [44]. Therefore, system identification using the PSO method was proposed for model parameter identification of the extracted hysteresis model.The dynamic modeling of PAM needs to be properly established as it affects the accuracy of control system performance. However, it is almost impossible to derive scientific solution of mathematical problem due to these behaviors. Therefore, the dynamic response was estimated based on shape estimation of hysteresis dataset, as shown in Figure10. It was performed at various excitation frequencies (i.e., 0.5, 1.0, 1.5, and 2.0 Hz) and external forces (i.e., 0, 50, 100, 150, 200, 250, and 300 N). Even though the PSO method uses a simple algorithm, it is a powerful identification method that is suitable for this application [5, 41].Figure 10
Shape estimation at various frequencies: (a) external force 100 N and (b) external force 150 N.
(a)(b)To obtain the best solution for dynamic friction coefficients (i.e.,γ^ and φ^) of the extracted hysteresis model, the control parameters of the PSO method were utilized. It was initialized with a population size (Ns) of 50 particles and number of generation (G) of 100 iterations. The constant parameters of both cognitive learning (p1) and social learning (p2) gains were adjusted to 1.5 gain, while the inertia weight (Iw) which exhibits close relationship with balance of the PSO algorithm was set to 0.4 of the typical value for excellent result. Figure 11 shows the dynamic response of two dynamic friction coefficients to the change of excitation frequency (ω) and external force (F) given to the PAM system. Its regression constants (i.e., μ and ρ) were characterized and then organized as the second part of the data profile as described in equations (12) or (14) for rate-dependent coefficient and equations (13) or (15) for stati-dependent coefficient. This second part of data profile characterized the extracted hysteresis component of activation pressure model. The identified parameters of extracted hysteresis model are shown in Table 2.Figure 11
Regression of dynamic friction coefficients.Table 2
Identified parameters of the extracted hysteresis model.
iaibicidi90.00E+000.00E+002.67E−016.67E−0110−2.19E−069.55E−04−9.42E−033.51E+01


## 4.1. Data Collection

Data collection is divided into three parts. The first part was collecting data from the experimental setup with quasistatic dependency. It was conducted by applying steady increment of vertical load from 0 N to 300 N. The second part was collecting data from the experimental setup with quasirate dependency. It was conducted by varying the excitation frequency of input signal from 0.5 Hz to 2.0 Hz. The third part was collecting data from the experimental setup with nonlocal memory behavior. It was conducted using random wave input signal which generates minor hysteresis loops. These data are composed of the contraction and pressure parameters at various excitation frequencies and external forces. The structural materials of the PAM inherently lead to hysteresis phenomena during cyclic contraction and expansion, as can be seen in Figure3. This nonlinearity was addressed as difficult error to handle especially for precise control system implementation.Figure 3
Sample of hysteresis data at various frequencies and external force conditions.


## 4.2. Data Characterization

Data characterization is divided into three parts. The first part involved characterization of local extremum of the hysteresis dataset using third-order polynomial regression. It was performed by characterizing the minimum and maximum values consisting of contraction-pressure parameters from the hysteresis dataset at various frequencies and external forces, as shown in Figure4. This characterization describes the quasistatic dependency of PAM. The major hysteresis dataset was implemented for this characterization method.Figure 4
Identification of local maximum and local minimum of hysteresis data.The second part involved characterization of the global extremum from the identified local extremum at various conditions. These are the first two critical points (primary data) from the hysteresis dataset: the first is global minimum (ε1, y1) and the second is global maximum (ε2, y2) parameters. These primary data allow for generalization of the virgin curves at various excitation frequencies (i.e., 0.5, 1.0, 1.5, and 2.0 Hz) using the single polynomial model, as can be seen in Figure 5. In this figure, the local extremum was distributed along with the identified single polynomial model. This characterization revealed that hysteresis dataset of various excitation frequencies was dependent on the quasistatic load.Figure 5
Identification of global maximum and global minimum of hysteresis data.The third part involved characterization of remaining two critical points (primary data) based on the obtained polynomial model solution: the first is inflection (ε3, y3) and the second is tangent (ε4, y4) parameters. These primary data can be obtained by solving the second derivative of the identified polynomial model and then finding its intersection from the origin. The unoptimized static pressure model based on characterization of the hysteresis dataset is shown in Figure 6. In this stage, the polynomial model has not been estimated based on the identified parameters (primary data). However, it was generated based on the curve fitting technique using the integral of error value.Figure 6
Unoptimized static pressure model using polynomial parameterization of third order.


## 4.3. Model Parameter Identification

The identified parameters from characterization of the hysteresis dataset were recognised as primary data for model parameter identification of the static pressure model. It composed of two stages of data enrichment; the first enrichment was performed based on the regressor of independent variableε⟶ℛ and regressor of dependent variable y⟶ℛ from identified eight parameters (primary data). These regressors were dependent on the external force (F) of the system. Its regression constants (i.e., α and δ) were organized as the first part of the data profile. This first part of data profile that characterized the static pressure component of the proposed activation pressure model is described in equation (1). The polynomial regression (X) of the static pressure model (y⟶)is shown in equations (2) or (3). It consisted of identification of regression coefficient (β⟶) which is characterized by the behavior of quasistatic dependency of PAM. The mathematical formulation of the regression coefficient β⟶^ was estimated according to equation (5). This regression coefficient was dependent on primary data, as shown in equation (4). Figure 7 shows regression constants for contraction (α) and pressure (δ) which controls the behavior of the proposed static pressure model. In this stage, global extremum (i.e., minimum and maximum) of the identified polynomial model was first enriched using the primary global minimum (ε1, y1) and global maximum (ε2, y2) parameters. It reconstructs the regression constants based on the optimized regression value of the global extremum as described in equations (6) and (8). This data enrichment reduces the deviation of quasistatic error and brings the polynomial model closer to primary global extremum.Figure 7
Regression constants of the static pressure model: (a) contraction and (b) pressure.
(a)(b)The second enrichment was performed based on the nonlinear region of PAM. As reported in the literature, PAM behaves highly nonlinear at low- and high-pressure activation. In the proposed static pressure model, the inflection parameter reflects nonlinearity boundary of low-pressure activation while the tangent parameter reflects nonlinearity boundary of high-pressure activation. This second data enrichment brings both inflection and tangent parameters of the polynomial model closer to the actual nonlinear region of PAM at low- and high-pressure activation as described in equations (7) and (9). Figure 8 shows an optimized static pressure model based on two stages of data enrichment.Figure 8
Optimized static pressure model using polynomial parameterization of third order.The measurement of quasistatic error during the two stages of data enrichment process is shown in Figure9(a). It was determined based on the shortest distance equation between the static pressure model and primary global extremum. The source of this error comes from the estimation used during the curve fitting technique. During the estimation, the model was not generated based on the identified primary data from the hysteresis dataset. However, it was performed based on the integration of minimal error in between the polynomial model and major hysteresis dataset. Figure 9(b) shows the estimated regression coefficient β⟶^ of the static pressure model after optimization, where its regression constants (i.e., α and δ) control the quasistatic dependency of the activation pressure model. The identified parameters of static pressure model are shown in Table 1.Figure 9
Model optimization: (a) quasistatic error and (b) regression coefficient of the static pressure model.
(a)(b)Table 1
Identified parameters of the static pressure model.
iaibicidi1−4.33E−093.58E−06−1.11E−031.33E−012−1.99E−098.70E−07−4.50E−042.04E−013−4.78E−093.74E−06−1.01E−035.48E−0246.98E−09−3.03E−06−3.61E−043.03E−015−2.85E−091.19E−063.86E−046.48E−026−4.89E−09−1.41E−071.30E−039.46E−0273.96E−10−1.84E−072.53E−052.77E−0287.40E−117.25E−08−2.91E−053.91E−01


## 4.4. System Identification

Apart from quasistatic and quasirate dependencies, the dynamic model of PAM also exhibits nonlocal memory behavior. It is historic dependent which changes drastically according to presliding domain of friction. This behavior is evident and can be observed as the shape of the hysteresis dataset drastically changed. Theoretically, this follows the mechanism of dynamic friction of presliding regime as described in [42, 43]. In response to this mechanism, extended form of the dynamic friction model was proposed for inverse modeling of PAM. Its mathematical formulation is described in equations (10) and (11). In this extended form, the proposed extracted hysteresis model (z⟶) was dependent on two dynamic friction coefficients (i.e., γ and φ). It is difficult to determine the dynamic response in between two dynamic friction coefficients without the use of metaheuristic methods such as genetic algorithm and PSO. There is no significant difference between these two methods in small scale; however, differences are seen in medium and large scale where genetic algorithms can only produce feasible solutions that are nearly optimal, while the PSO algorithm has ease of implementation and also has high calculation accuracy [44]. Therefore, system identification using the PSO method was proposed for model parameter identification of the extracted hysteresis model.The dynamic modeling of PAM needs to be properly established as it affects the accuracy of control system performance. However, it is almost impossible to derive scientific solution of mathematical problem due to these behaviors. Therefore, the dynamic response was estimated based on shape estimation of hysteresis dataset, as shown in Figure10. It was performed at various excitation frequencies (i.e., 0.5, 1.0, 1.5, and 2.0 Hz) and external forces (i.e., 0, 50, 100, 150, 200, 250, and 300 N). Even though the PSO method uses a simple algorithm, it is a powerful identification method that is suitable for this application [5, 41].Figure 10
Shape estimation at various frequencies: (a) external force 100 N and (b) external force 150 N.
(a)(b)To obtain the best solution for dynamic friction coefficients (i.e.,γ^ and φ^) of the extracted hysteresis model, the control parameters of the PSO method were utilized. It was initialized with a population size (Ns) of 50 particles and number of generation (G) of 100 iterations. The constant parameters of both cognitive learning (p1) and social learning (p2) gains were adjusted to 1.5 gain, while the inertia weight (Iw) which exhibits close relationship with balance of the PSO algorithm was set to 0.4 of the typical value for excellent result. Figure 11 shows the dynamic response of two dynamic friction coefficients to the change of excitation frequency (ω) and external force (F) given to the PAM system. Its regression constants (i.e., μ and ρ) were characterized and then organized as the second part of the data profile as described in equations (12) or (14) for rate-dependent coefficient and equations (13) or (15) for stati-dependent coefficient. This second part of data profile characterized the extracted hysteresis component of activation pressure model. The identified parameters of extracted hysteresis model are shown in Table 2.Figure 11
Regression of dynamic friction coefficients.Table 2
Identified parameters of the extracted hysteresis model.
iaibicidi90.00E+000.00E+002.67E−016.67E−0110−2.19E−069.55E−04−9.42E−033.51E+01


## 5. Model Simulation

This paper proposed a polynomial parameterization of third order in characterizing static and dynamic behaviors of PAM. The simulation was performed based on evaluated data profile of static pressure and the extracted hysteresis dataset obtained from the model parameter and system identification method. The random wave signal of the independent variable obtained from the hysteresis dataset at various excitation frequencies and external forces was used for the model simulation of the proposed activation pressure model. This random wave signal generates various contracting pattern of cyclic contraction and expansion. It created different patterns of major and minor hysteresis loops with an effect of nonlocal memory behavior. To evaluate the reliability of the proposed activation pressure model, comparison between model simulation and actual hysteresis data was performed. The result obtained from model simulation of various excitation frequencies and external forces is shown in Figure12. This result shows that the proposed activation pressure model was able to generate comparable hysteresis of PAM during cyclic contraction and expansion. The statistical analysis of the data provides high significant correlation values from 0.875 to 0.988 at all evaluated conditions.Figure 12
Sample of hysteresis simulation at evaluated frequencies and external forces.Figure13 shows a sample of random wave simulation at various excitation frequencies and external forces. This result shows that the proposed activation pressure model was able to characterize the nonlinear behaviors of PAM when induced with quasistatic, quasirate, and nonlocal memory-dependent random wave input signal. The performance of the proposed activation pressure model was evaluated based on the RMSE value, as shown in Figure 14. The obtained RMSE value measures the reliability of the proposed activation pressure model in characterizing the nonlinear behaviors of PAM with less than 10% of error at all evaluated conditions. In providing a better solution in modeling dynamic behaviors of PAM, model parameter identification and the PSO method were presented to enhance the proposed activation pressure model. However, there was still some error remaining due to high nonlinearity of PAM that cannot be captured. This error cannot be avoided because the nonlinearities inherited by the PAM were varied and complicated. This remaining error could easily be compensated with an addition of an intelligent control strategy.Figure 13
Sample of random wave signal simulation at evaluated frequencies and external forces.Figure 14
Root mean square error (RMSE) at all evaluated frequencies and external forces.


## 6. Conclusions

In conclusion, the inverse modeling of PAM in the form of activation pressure model was presented. The mathematical formulation was derived according to the static pressure and extracted hysteresis components. This derivation follows the polynomial regression and dynamic friction mechanism. In this principle, the activation pressure modelp⟶ was dependent on polynomial regression β⟶ of the static pressure model and dynamic friction coefficients γφ of the extracted hysteresis model. These coefficients were estimated according to data enrichment of primary data obtained from characterization of the hysteresis dataset, where its regression constants (i.e., α, δ, μ, and ρ) were determined by using model parameter identification and the PSO method. Based on regression constants, the estimator of regression coefficient β⟶^ was obtained by the least squares method, while the dynamic friction coefficients γ^φ^ were obtained by dynamic response estimation. The final form of the activation pressure model (p⟶=Xβ⟶^+γ^φ^ε˙⟶) was designed to characterize the quasistatic, quasirate, and nonlocal memory behaviors of PAM. It was evaluated at various excitation frequencies (ω = 0.5, 1.0, 1.5, and 2.0 Hz) and external forces (F = 0, 50, 100, 150, 200, 250, and 300 N) using random wave input signal. The result from model simulation shows the RMSE value of less than 10% error evaluated at all conditions. This inverse modeling approach should be useful in solving some part of control system implementation of healthcare engineering applications such as rehabilitation robotics, biomedical system, and humanoid robots.



---

*Source: 8189157-2020-12-09.xml*