Comparison between Genetic Programming and Dynamic Models for Compact Electrohydraulic Actuators

Abstract

A compact electrohydraulic actuator (C-EHA) is an innovative hydraulic system with a wide range of applications, particularly in automation, robotics, and aerospace. The actuator provides the benefits of hydraulics without the expense and space requirements of full-sized hydraulic systems and in a much cleaner manner. However, this actuator is associated with some disadvantages, such as a high level of nonlinearity, uncertainty, and a lack of studies. The development of a robust controller requires a thorough understanding of the system behavior as well as an accurate dynamic model of the system; however, finding an accurate dynamic model of a system is not always straightforward, and it is considered a significant challenge for engineers, particularly for a C-EHA because the critical parameters inside cannot be accessed. Our research aims to evaluate and confirm the ability of genetic programming (GP) to model a nonlinear system for a C-EHA. In our paper, we present and develop a GP model for the C-EHA system. Furthermore, our study presents a dynamic model of the system for comparison with the GP model. As a result, by using this actuator in the 1-DOF arm system and conducting experiments, we confirmed that the GP model has a better performance with less positional error compared with the proposed dynamic model. The model can be used to conduct further studies, such as designing controllers or system simulations.

1. Introduction

Hydraulic actuators are among the most common in heavy industries. These actuation systems can provide a high amount of force while consuming a minimum amount of energy. However, most of these actuation systems are composed of a number of components that may impose limitations on the system. As a result, these types of actuators might be very difficult to use in applications where space, cleanliness, and weight are limited, such as in robotics and aerospace. A compact electrohydraulic actuator (C-EHA), which offers a number of advantages over conventional hydraulic actuators, has been developed to overcome this issue. Figure 1 shows a commercial type of this actuator. These actuators are significantly smaller and lighter compared with conventional hydraulic actuators, can run on battery power, are isolated, very clean, plug-and-play, and can be integrated into a wide range of systems. These characteristics distinguish this actuator and open up a broad variety of applications for it, particularly in positioning systems and robotics applications where high force is required in a limited space. Despite this product being well developed and already commercially available, engineers continue to improve the control system’s performance because these actuators are associated with a high level of non-linearity and uncertainty. Having an accurate model of the system could help improve tracking accuracy and enhance control performance, especially for model-based controllers [1].

The dynamic modeling approach is considered the most popular method of determining the system’s model. Modeling the system using physical rules offers its own set of benefits and provides engineers with a unique perspective on the system’s features. However, this strategy requires an extensive understanding of the system, which may not always be feasible when dealing with highly complex nonlinear systems, specifically hydraulic systems [2]. A variety of linear and nonlinear controllers have been well developed over the past few decades to more precisely control the conventional type of electrohydraulic actuators [3,4,5,6]. Most of these systems monitor and measure critical internal parameters, such as oil pressure and flow rate, during the actuator’s operation. These parameters can be utilized to considerably increase the controller’s precision and performance. Moreover, a standard controlled valve has been used to control most of these actuators’ cylinder rods. However, no such components to control the output of an C-EHA exist, which are also known as pump-controlled hydraulic or electrohydrostatic systems, or self-contained electrohydraulic cylinder drives (Figure 2). Because these actuators rely solely on the pump motor’s speed as input, it is difficult to determine an accurate dynamic model for the control system because oil pressure and flow rate are not accessible to measure in C-EHA systems.

This type of actuator has been the subject of several studies over the past decade with the aim of improving their characteristics, modeling, and identification, or developing higher-precision controllers. The following are some highlights from previous studies. Hagen et al. listed various types of internal structures for valveless types of actuators [7]. Ling et al. conducted a study on identifying the system and improving the accuracy of an EHA’s position-tracking system. The authors applied stimulus–response data and implemented system identification processes without applying physics laws for modeling, resulting in a linear model for the nonlinear system [8]. Nie et al. developed an extended-state-observer-based backstepping control system for underwater electrohydrostatic actuators [9]. Martin and Louis developed an electrohydrostatic actuator model; however, they only made simulations and did not compare the results with the actual actuator [10].

The some of research in the area of electrohydrostatic actuators is focused on controlled-valved actuators; valveless actuators have not received significant attention. Even though there have been some significant achievements in these studies, finding an accurate model of the system remains a challenge. As previously mentioned, measuring hydraulic actuator internal parameters such as oil pressure or flow rate can significantly increase the accuracy of the model and control. Due to the fact that these parameters are not accessible in this type of isolated system, the only parameter that can be measured is the pump speed. Therefore, developing an accurate model of the system is still considered a challenge for researchers.

Alternatively, artificial intelligence techniques can be utilized to facilitate the development of new modeling techniques without requiring the accurate dynamical models of various systems [11,12,13], specially an EHA [14,15]. Genetic programming (GP), an extension of genetic algorithms (GA), is a well-known solution for nonlinear problems. It is an optimization method that can be applied to optimize complex systems with nonlinear structures, as well as to present systems and their parameters [16]. GP provides a number of benefits, such as the possibility of utilizing system pre-knowledge to prepare a better initial candidate solution that could help to achieve better results [17]. GP has numerous applications in engineering, and it can be applied to the design and identification of various systems. For example, Santos utilized GP to identify a nonlinear model for an experimental ball-and-tube system [18]. Recently, Raibaudo and his team implemented GP on an experimental setup to optimize unsteady feedback controllers. The GP method also has additional applications in engineering and non-engineering fields (see [19,20,21,22]).

A primary objective of our research is to examine the effectiveness of using GP methods to develop an accurate model for a case study of a commercial type of C-EHA (Figure 1). We investigated a GP model that was proposed and reported in a previous paper, and we evaluated the accuracy and performance of this method [23]. Our study develops an appropriate dynamic mode for comparison. We applied the same input command to both dynamic and GP models, as well to a real actuator system. We performed a comparison of the position-tracking error data between the two models to check which has better accuracy.

An introduction to the C-EHA is provided in Section 2, followed by details on the development of a proper dynamic model. Section 3 briefly describes the GP method and the parameter settings developed in the previous paper, which we used in our study to model the compact electrohydraulic actuator. The experiment and final results for both the GP and dynamic models are presented and discussed in Section 4. The final section concludes our study.

2. Dynamic Modeling

2.1. Compact Electrohydraulic Actuator (C-EHA)

C-EHAs have a broad range of applications, particularly in automation and robotics, due to their unique features. In contrast to conventional actuators, this actuator does not use an actuated valve to control the amount of oil and, consequently, pressure on the piston sides. Instead, automatic valves, which control the flow of oil based on the rational direction of the pump, are integrated with the system. C-EHA technology simply connects two ports of a bi-directional gear pump to either side of a piston, along with a small oil reservoir as a means of generating an initial pressure. As a result of a motor driving the gear pump, the piston moves in either direction. Figure 3 shows that this actuator’s motor rotates the gear pump in one direction to pressurize the hydraulic fluid and extend the cylinder. Fluid is drawn from both the small reservoir inside the actuator and the cylinder’s head side as the auto-valves open, resulting in the actuator rod’s extension. The motor runs in the other direction during the retraction phase, reversing the operation and returning the fluid to the piston’s opposite side. The leakage is assumed to be negligible because of the system’s isolation and its inability to measure any internal leakage. The motion equation for this actuator can be derived using the continuity equation shown in (1) following the assumption that oil is not compressible [24,25].

$$\begin{array}{c}\hfill Q_1=\frac{d{V}_1}{dt}+\frac{V_1}{{\beta }_e}\frac{d{P}_1}{dt}\\ \hfill Q_2=-\frac{d{V}_2}{dt}-\frac{V_2}{{\beta }_e}\frac{d{P}_2}{dt}\end{array}$$

(1)

where $Q_i$, $P_i$, and $V_i$ are the flow rate, cylinder pressure, and cylinder volume, respectively. ${\beta }_e$ is defined as the effective bulk modulus of oil. The cylinder volume can be rewritten as in (2), where $V_0$ is the cylinder’s initial volume, $x_p$ is the piston position, and $A_i$ is the piston’s cross-section area.

$$\begin{array}{c}\hfill V_1=V_{01}+A_1{x}_p\\ \hfill V_2=V_{02}-A_2{x}_p\end{array}$$

(2)

Consequently, by substituting (2) into (1), the continuity equation can be rewritten as (3).

$$\begin{array}{c}\hfill \dot{P_1}=Q_1\frac{{\beta }_e}{V_1}-A_1\dot{x_p}\frac{{\beta }_e}{V_1}\\ \hfill \dot{P_2}=-Q_2\frac{{\beta }_e}{V_2}+A_2\dot{x_p}\frac{{\beta }_e}{V_2}\end{array}$$

(3)

The difference in pressure between the piston’s two sides generates external force. Equation (4) describes this external force in the presence of friction between the cylinder and piston.

$$\begin{array}{c}\hfill P_1{A}_1-P_2{A}_2=F_{ext}+F_{friction}\end{array}$$

(4)

As previously mentioned, the inside pressure value is unavailable for measurement. Thus, we can consider the derivative of (4) as follows.

$$\begin{array}{c}\hfill \dot{P_1}{A}_1-\dot{P_2}{A}_2=\dot{F_{ext}}+\dot{F_{friction}}\end{array}$$

(5)

The flow rates $Q_1$ and $Q_2$ are circulated by the internal gear pump, described as follows, where $c_p$ is the pump’s volumetric capacity and ${\omega }_p$ is the pump’s rotational speed.

$$\begin{array}{c}\hfill Q_1=c_p{\omega }_p,Q_2=-\frac{A_2}{A_1}{Q}_1\end{array}$$

(6)

The generated force’s final equation can be expressed as (7).

$$\begin{array}{cc}\hfill & {\beta }_e{c}_p\left(\frac{A_1^2}{V_1}+\frac{A_2^2}{V_2}\right)\dot{x_p}+{\beta }_e{c}_p\left(\frac{A_1}{V_1}-\frac{A_2^2}{V_1{V}_2}\right){\omega }_p={\dot{F}}_{ext}+{\dot{F}}_{friction}\hfill \end{array}$$

(7)

Stribeck effects are applied to determine the C-EHA’s friction model. Figure 4 shows the Stribeck effect for the viscus friction force as a function of velocity. Contact surfaces can have different friction forces depending on their velocities. There are four major regions; steady-state velocity, boundary lubrication, and partial fluid and full fluid lubrication. The friction model and its derivatives can be defined by (8) and (9), respectively [26].

$$\begin{array}{c}\hfill F_{friction}=\left(F_k+(F_s-F_k)e^{-{\left(\frac{{\dot{x}}_p}{{\dot{x}}_s}\right)}^2}\right)sgn\left({\dot{x}}_p\right)+F_v{\dot{x}}_p\end{array}$$

(8)

$$\begin{array}{c}\hfill {\dot{F}}_{friction}=-2(F_s-F_k)e^{-{\left(\frac{{\dot{x}}_p}{{\dot{x}}_s}\right)}^2}\frac{\dot{x}\ddot{x}}{{\dot{x}}_s^2}sgn\left({\dot{x}}_p\right)+F_v{\ddot{x}}_p\end{array}$$

(9)

where $F_k$, $F_s$, and $F_v$ are the coefficients of kinetic, static, and viscous frictions, respectively, and ${\dot{x}}_p$ is the Stribeck velocity. We can determine these unknown parameters by performing specific tests on the actuator. Section 4.1 provides a detailed explanation of this procedure.

2.2. Equation of Motion for Arm Manipulator

Figure 5a shows a one-DOF arm manipulator that we developed to provide the experimental datasets. This arm is equipped with the commercial-type C-EHA made by the Parker Hannifin Company with a maximum payload of $2669N$ (Figure 1).Both mounting sides of the C-EHA can be easily adjusted in this setup. A strong bolt is also used to attach the weight as an external load to the head of the arm. Weights of different sizes are provided, which can be easily changed or even installed together using longer bolts. Section 4 shows the results of a number of experiments we conducted to provide training and validation datasets for the GP model and find friction parameters for the dynamic models. We used these data to demonstrate the proposed models’ effectiveness and accuracy.

Here, we present the motion equation for the arm manipulator. We chose to apply the Lagrange method to this system to find the model. Based on the parameters defined in Figure 5b, (10) and (11) present the system’s kinematic and potential energy equations, respectively.

$$\begin{array}{c}\hfill T=\frac{1}{2}\left(m_1{L}^2+\frac{m_2{L}^2}{3}\right){\dot{\theta }}^2\end{array}$$

(10)

$$\begin{array}{c}\hfill U=\left(m_1+\frac{m_2}{2}\right)gLsin\theta \end{array}$$

(11)

Lagrange’s method is described as follows

$$\begin{array}{c}\hfill \frac{d}{dt}\left(\frac{\partial T}{\partial \dot{\theta }}\right)-\frac{\partial T}{\partial \theta }+\frac{\partial U}{\partial \theta }=\tau \end{array}$$

(12)

where $\tau$ represents the applied torque at the rotational joint. The relation between the external force and torque generated by the hydraulic actuator can be defined as follows:

$$\begin{array}{c}\hfill \tau =F_{ext}lsin\alpha \end{array}$$

(13)

By substituting Equations (10)–(13) into Lagrange’s method, Equation (12), the motion equation of the arm manipulator and its derivative, can be obtained, as expressed in Equation (14) and Equation (15), respectively.

$$\begin{array}{cc}\hfill & \left(m_1{L}^2+\frac{m_2{L}^2}{3}\right)\ddot{\theta }+\left(m_1+\frac{m_2}{2}\right)gLcos\theta =F_{ext}lsin\alpha \hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & \left(m_1{L}^2+\frac{m_2{L}^2}{3}\right){\theta }^{\left(3\right)}-\left(m_1+\frac{m_2}{2}\right)gLsin\theta ·\dot{\theta }={\dot{F}}_{ext}lsin\alpha \hfill \end{array}$$

(15)

By substituting Equation (7) in Equation (15), the final motion equation for the arm manipulator can be summarized as in Equation (16)

$$\begin{array}{cc}\hfill & {\theta }^{\left(3\right)}=\frac{H_2}{H_1}sin\theta ·\dot{\theta }-\frac{H_3}{H_1}{T}_1{\dot{x}}_p-\frac{H_3}{H_1}{\dot{F}}_{friction}+\frac{H_3}{H_1}{T}_2{\omega }_p\hfill \end{array}$$

(16)

where

$$\begin{array}{cc}\hfill & H_1=m_1{L}^2+\frac{m_2{L}^2}{3},H_2=\left(m_1+\frac{m_2}{2}\right)gL,H_3=lsin\alpha \hfill \\ \hfill & T_1={\beta }_e\left(\frac{A_1^2}{V_1}+\frac{A_2^2}{V_2}\right),T_2={\beta }_e{c}_p\left(\frac{A_1}{V_1}-\frac{A_2^2}{V_1{V}_2}\right)\hfill \end{array}$$

3. Genetic Programming and Its Application for Hydraulic Actuator Problem

Genetic programming is considered an extension of genetic algorithms. Generally, genetic programming refers to an automated method of generating computer programs that carefully solve specific problems by the principle of natural selection. GP was created using evolutionary algorithms and was first developed by Barricelli to simulate evolution [27]. Later, Rechenberg used GP in the early 1970s to solve more complex engineering problems through evolution strategies. Nicholas Cramer made the first “tree-based” genetic programming declaration [28]; furthermore, John Koza played a key role in the development of GP and its application to a number of different situations [29]. The flowchart in Figure 6a illustrates a general GP procedure. The initial populations are made based on pre-knowledge of the system and basic algebra operations. A fitness function is used to evaluate each candidate solution based on the available training-set data, which is obtained from experiments. As each candidate solution undergoes GP operations, their fitness values impact their chances of survival. Good candidate solutions can be improved with each successive generation and, ultimately, the best result can be achieved after many generations. The final model’s validation can be conducted using the validation-set data obtained from the experiments. In our study, we used previously proposed GP model structures and parameters, which are briefly discussed as follows [23].

3.1. Encoding and Initial Population

Figure 6 shows all the individuals are encoded as a tree structure in a mathematical model. Terminal (basic algebra operations and subfunctions) and root nodes (constant numbers and inputs) exist in every structure. GP individuals do not have a fixed length, and they may change from generation to generation. Randomly selecting initial candidates allows GP to generate a wide range of candidate solutions [16]. However, C-EHA requires only one input command to increase the GP model’s accuracy, with the motor pump’s current speed as an input also being considered. Additionally, the arm’s current position is also added to the model as an input to obtain the arm’s absolute position. Considering these three inputs and one output as the next arm position, the GP model is considered a multi-input–single-output model (MISO). Considering the available functions and conditions, we assume that each node can have a maximum of two branches and a minimum of none. The candidate’s final value can then be calculated by using the input values and substituting them into the roots.

3.2. Function Library

Generally, function libraries can be imagined as pools containing functions and terminal sets for making initial population or mutation branches. The terminal set consists of statements, operations, basic algebra, and sub-functions while the root set is composed of the system inputs and constant numbers. Any prior knowledge of the system should also be considered for the initial population. Therefore, in addition to the basic algebraic and sinusoidal functions, the hyperbolic tangent function is also selected as part of the function library because it has similar properties to the friction model.

3.3. Fitness Function

Fitness functions are a type of objective function that measures how well a solution candidate performs in comparison with the actual data available. Several well-known fitness functions exist for GA and GP depending on the conditions [30]. We define fitness function in our study as the sum of squared errors between the model output and experimental data.

3.4. Operations

Similar to the GA method, the GP method uses three major operations to create the main components of the GP: selection, crossover, and mutation [31]. The selection process, such as roulette-wheel, tournament, ranking, and elitist selection, is the process in which candidates from one generation are selected according to the selection method and fitness values to generate new candidates for the next generation. Crossover refers to the process of combining two individuals (parents) to create new individuals (offspring). A mutation is a genetic operation that changes one part of an individual by using a randomly generated sub-tree. In each generation, after the selection operation, both crossover and mutation operations are usually partially applied to the candidates during the evaluation based on their probability factors. These factors might vary depending on the optimization. In our proposed GP model, in accordance with elitist selection theory, 25 percent of the population with the highest fitness values are directly passed down to the next generation in order to avoid the effects of crossover and mutation. The remainder of the population is then subjected to crossover and mutation with rates of 80% and 3%, respectively. We used the roulette-wheel method to select candidates for crossover and mutation operations. These operations may result in a longer or shorter tree structures because of the random selection of nodes for crossover and mutation. However, long structures pose some challenges to GP. Consequently, we defined a 20-node limit for individuals to avoid creating long tree structures. An operation that produces offspring with a length greater than the limitation is not acceptable, and the selected node should be changed again until the new tree structure length is less than the limitation.

3.5. Termination

Using termination rules, a genetic algorithm determines whether to continue searching for a better answer or to stop searching. The termination criterion is checked after each generation to determine whether the answer meets the criterion. Based on the problem conditions, there are many ways to define the criteria for termination. The most well-known criteria are generation number, fitness threshold, and fitness convergence [32]. We used all three criteria simultaneously in our study to determine termination during the GP generations.

4. Experiment, Result, and Discussion

Our research requires various experiments to be designed and performed in order to produce data sets. Three major experiment setups are designed and conducted as follows.

4.1. Friction Parameter Identification

To find the unknown friction coefficients expressed in Equation (8) and check how much the friction term can affect the system response, it is essential to obtain the relation between the friction force and actuator speed. The first step in conducting such a test is to activate the manual release valve in the C-EHA to make a passive motion of the piston while the actuator is not in operation mode. Next, the system should be moved at a constant speed by applying controlled external force. We measured and recorded the external force and piston velocity for post-processing. The external load represents the amount of friction inside the actuator because it is not in operation. Figure 7 shows the results of a number of tests conducted at various speeds. We obtained the arm manipulator’s estimated friction parameters using Matlab curve-fitting toolbox and listed them in

Table 1. Parameters of C-EHA used in this study.

Parameter	Symbol	Unit	Value
Extra mass	$m_1$	kg	46.7
Mass of arm	$m_2$	kg	1.64
Length of arm up to weight	L	m	0.488
Length of arm up to actuator head pin	l	m	0.225
Cross-section area of piston with rod	$A_1$	m${}^2$	$6.319\times {10}^{-4}$
Cross-section area of piston without rod	$A_2$	m${}^2$	$8.042\times {10}^{-4}$
Initial length of hydraulic actuator	$x_0$	m	0.3055
Volumetric capacity of the pump	$c_p$	m${}^3$/rad	$1.050\times {10}^{-7}$
Effective bulk modulus of the oil	${\beta }_e$	Pa	$1.747\times {10}^9$
Kinetic friction coefficient	$F_k$	N	1.87
Static friction coefficient	$F_s$	N	13.15
Viscous friction coefficient	$F_v$	Ns/m	741.36
Stribeck velocity	$\dot{x_s}$	m/s	0.0028

along with other geometric parameters. Figure 7 shows the range of inner-friction force ($60N$) is relatively smaller than the actuator maximum payload when it is in operation ($2669N$). Friction terms can be neglected depending on the system and the model’s purpose. Nevertheless, we did not ignore them in our study.

4.2. GP Dataset

The GP method requires two sets of data for training and validation purposes. We conducted a large number of tests on the real system while measuring and recording inputs and outputs using a computer to create these datasets. The piston’s speed in this actuator is controlled exclusively by the pump speed, which is directly related to the motor input voltage. We conducted five sets of 80-s tests for GP training, as well as three sets for validation. A variety of frequencies and amplitudes are incorporated into input commands for creating datasets to ensure that all of the various actuator modes are stimulated (Figure 8).

4.3. Dynamic Model Validation

We implemented the same validation tests utilized for the GP model to validate our proposed dynamic system model. Furthermore, the results can be used to compare the performance of the two methods. However, the motor speed rather than motor voltage is required for dynamic models. We attached a proximity sensor to the motor side prior to conducting tests to determine motor speed (Figure 9). We conducted additional tests to provide a relationship between the motor voltage and speed. Figure 10 shows that this relation can be considered linear.

4.4. Results

We used the GP method to obtain a model of the system based on the training data set. The following Equation (17) represents the best answer found by the GP method:

$$\begin{array}{cc}\hfill & y=x_1+0.0002027\times x_3+0.002027\times tanh(cos\left(x_2\right))-0.002027\times tanh(x_2\times (x_1+8.97))\hfill \\ \hfill & -0.004659\times tanh(cos(cos\left(x_2\right)))+0.004659\times x_2\times cos(cos\left(x_3\right))+0.004659\times x_2\hfill \\ \hfill & +0.004659\times x_2\times tanh(x_1+8.97)+0.0001014\times x_1\times cos(2\times cos\left(x_3\right))+0.00214\hfill \end{array}$$

(17)

where, $x_1$, $x_2$, and $x_3$ represent the current motor position, the voltage input, and the motor velocity, respectively. In our study, each time step is 1 millisecond. This is followed by applying the three different validation input command sets to the dynamic and GP models, resulting in the arm angle and speed illustrated in Figure 11 and Figure 12, respectively. As can be seen in these figures, the GP model has slightly better accuracy in comparison to the dynamic model. Position accuracy is of greater concern in most C-EHA applications because C-EHA actuators operate at low speeds compared with other actuators. Therefore, for ease of comparision, the position error between these two models and experimental data presented in Figure 11 are calculated and presented in Figure 13. As can be seen in these figures, the GP model has slightly better accuracy in comparison to the dynamic model. This could be confirmed with the error bar presented in Figure 14. Moreover, the root mean square error (RMSE) of each dataset for each model is calculated by following equation;

$$\begin{array}{c}\hfill I_{error}=\frac{{\sum }_{i=1}^N\sqrt{e_i^2}}{N}\end{array}$$

(18)

where $N=1,2,3$. The results in

Table 2. Root mean square error of GP and dynamic models.

Set	Dynamic Model	GP Solution 1
Validation 1	1.355	0.735
Validation 2	1.048	0.556
Validation 3	1.217	0.715

confirm that the GP model has a slightly better performance in comparison with the dynamic method.

5. Conclusions

We employed two different modeling methods in our study to obtain a compact electrohydraulic actuator model and analyze its performance and accuracy. We developed a one-DOF arm with a payload to provide experimental data where all input and output can be measured and recorded by a computer. First, we created a dynamic actuator model along with a friction model for the system. Because the critical parameters, such as pressure or flow rate, are not accessible for measurements during the actuator’s operation, we need to be consider some assumptions that could reduce the dynamic model’s accuracy. We identified friction parameters through the specific experiment and curve fitting technique. Figure 8 shows that the friction-force contribution to the system response is considerably small regarding the payload rate; therefore, we consider it negligible. However, we decided to include it in the final model to keep the accuracy of the dynamic model as high as possible. We used the parameters of a previously developed GP model to obtain our model. We obtained various sets of training and validation data by utilizing an experimental setup. Finally, we employed both models to estimate the system’s response to the input command using a variety of validation datasets, as shown in Figure 11. The error of these two models is presented in Figure 12 and Figure 13 and Table 2. Our final results confirm that the GP model has a relatively effective performance, which could be used later for further simulations or improving control accuracy. Although it may be difficult to represent these GP models at a state–space level or to rewrite them as frequency responses, they can still be used for system simulations or in reference model controllers. More work can be performed in this field, and the exact practical applications of this model will be investigated in future works.