A hybrid model for fast and efficient simulation of fluid power circuits with small volumes utilizing a recurrent neural network

The modeling and simulation of fluid power systems are essential parts of the real-time simulation of virtual prototypes of mobile working machines. In several cases in the dynamic simulation of such fluid power systems, a longer simulation time is required. This makes the traditional mathematical models inefficient for real-time simulations, particularly when simulating fluid power systems because of the small volumes in stiff differential equations of pressure. The main reason for such an issue in the simulation of fluid power systems is the existence of small oil volumes in stiff differential equations of pressure. To overcome this issue, a novel hybrid model is proposed for stiff fluid power systems simulation. The main feature of the model is the utilization of a recurrent neural network instead of stiff differential equations of pressure with small volume. At the same time, the dynamics of the rest system are traditionally presented with algebraic and differential equations. The testing results of the introduced hybrid model showed that the novel method can reduce the simulation time, which makes the model suitable for real-time applications. Moreover, the accuracy of the model remains at a fairly high level compared to traditional mathematical models.


I. INTRODUCTION
F LUID power systems are widely used in the real-life modeling of various mobile machines such as logging harvesters, cranes, excavators, and industrial robots. The recent trend in modeling digital twins [1] and virtual prototypes [2]- [4] has shown that mathematical modeling of fluid power systems plays a vital role in the development of industrial simulators of such mobile working machines. For this purpose, real-time and faster than real-time techniques [5] are used to get a fast response in the system. However, fluid power circuits can contain singularities that directly affect the computational speed of the whole system and make a real-time simulation of the system impossible.
In general, there are two main problems in fluid power circuits modeling and simulation. The first problem is related to the pressure drop approaching zero. This phenomenon is associated with difficulties in the use of the traditional turbulent flow orifice equation because of the infinite value of the flow rate derivative. To solve this problem, several combined orifice models were proposed by Ellman et al. at the end of the 1990s [6], [7]. Another computationally efficient solution was proposed by Aman et al. in 2008 [8] in which the polynomial relation between flow rate and pressure drop was derived for cases when the pressure drop approaches zero. The model was named the two-regime flow model in which the third-order polynomial is used for describing the laminar and transition flow areas, whereas the traditional square root turbulent orifice equation of flow is used for the turbulence regime.
Another problem in the dynamic simulation of fluid power circuits can be associated with the numerical stiffness of ordinary differential equations [9]. The numerical stiffness in fluid power systems can be explained by the fact that such circuits include volumes of different orders of magnitude.
This results in difficulties in the numerical integration of ordinary differential equations, and classical explicit integration solvers are not able to generate a stable dynamic response at a high integration time step, slowing down the simulation by implementing significantly small time steps. The numerical stiffness of the system can often be associated with small volumes in a hydraulic circuit. The solution for the systems with small volumes should be derived by using specific numerical solvers to avoid the problem of instability in a dynamic model of a fluid power circuit [10].
The problem of small volumes and the resulting numerical stiffness of fluid power circuits arises in various scientific papers, and as a result, special solvers are developed. This problem was first mentioned in 1990 by Bowns and Wang [11], who proposed an iterative technique to solve the problem of numerical stiffness due to the small volumes in hydraulic pipes. However, the technique was still computationally costly.
Another iterative technique was first presented by Aman and Handroos [12]- [14]. This implicit solver was named the pseudo-dynamic solver. This solver consists of two loops-an outer (or main) loop and an inner (or pseudo) loop. The main loop contains algebraic and differential equations of the whole fluid power system, except for the pressure and flows that are associated with the small volume. The integration of such pressure occurs in the inner loop of the solver with an artificially enlarged volume instead of the real small volume. The idea of the inner loop is to receive the steady-state response of the pressure by iterating such pressure with artificial volume until the convergence criterion is reached.
The advanced version of the pseudo-dynamic solver (Ad-vPDS) was proposed by Malysheva et al. in 2020 [15]. The main difference of this approach from the original pseudodynamic method is the implementation of adaptive criterion of convergence. The criterion depend on the pressure levels in the system, and it was experimentally proven that with small levels of pressure, the smaller criterion have more accurate results; and at the same time, the bigger criterion can accelerate the computational speed of the fluid power system. Ultimately, the method is a trade-off between accuracy and speed suitable for real-time or faster than real-time applications. The method has also been improved by Malysheva et al. [16], where it was investigated that different integrators can also affect the computational speed of AdvPDS.
One of the most reliable methods proven to be excellent in simulating fluid power systems in theory and practice is the method proposed by Kiani-Oshtorjani et al. [17]. This method is based on the singular perturbation theory. The main idea of this approach is to substitute stiff ordinary differential equations with small volumes for algebraic equations modified for fluid power systems simulation in accordance with this theory. The method was tested in simulating multibody systems [18], and an accurate and fast response of the system was achieved. However, the method has its own drawbacks related to the accumulative error in certain cases of use. To overcome the error, a special corrector factor for the model should be used. To solve this problem and allow the simulation without a corrector factor definition, a novel method based on the combined singular perturbation theory and Method of multiple scales (MMS) was proposed by Kiani-Oshtorjani et al. [19]. This robust method allows the elimination of the accumulative error and makes the simulation faster and more accurate than the original method based on the singular perturbation theory.
In several works [20]- [22], the components of fluid power systems were introduced as neural networks with basic inputs and outputs. This type of system modeling showed that the use of neural networks in fluid power circuits modeling provides a fast response in the system that allows a simulation of the system in real-time or faster than real-time applications. In addition, predictive neural network models of various dynamic systems with ordinary differential equations (ODE) were studied in [23] and [24]. Both papers are based on the studies of new methods for improving the performance of neural networks based on simple ODEs and complex nonlinear dynamic systems. Nevertheless, the most significant results were shown by using recurrent neural networks (RNN) [25]- [27] in the modeling and simulation of different systems.
However, the systems based on neural networks are simulated as a black box, and the dynamics of the system are totally neglected. Substituting a stiff differential equation with small volumes with a computationally efficient and accurate neural network can solve the problem of singularity arising in the simulation of such stiff fluid power systems. It means that only the stiff equation will be replaced with the neural network and the system dynamics of the entire system will be saved in the mathematical model. And at the same time, the simulation is supposed to be faster due to the absence of stiff differential equations in the areas with small volumes.
The main purpose of this research paper is to introduce a novel hybrid model for a fast simulation of fluid power systems with small volumes using the RNN only for the pressure continuity equation with small volume in calculations. Section II describes the fluid power systems under consideration in the research. Section III provides information about the RNN used in the development. Section IV describes the development in detail, from the collected data and the RNN training to the implementation of the neural network in the fluid power circuits from Section II. Finally, sections V and VI provide the results of the simulations, and a conclusion is made about the model developed.

II. FLUID POWER CIRCUIT DESCRIPTION
This section describes the systems studied in the current research. Overall, two fluid power circuits were investigated. Both systems have one thing in common-the content of one small volume inside. The first circuit is a simple fluid power system with multiple orifices located in series. The second fluid power circuit is a simple pressure-compensated system with a 2-chamber hydraulic cylinder controlled by a directional control valve.
The systems were modeled in accordance with the lumped fluid theory, which assumes that a fluid power system can be divided into sections with separate volumes where the pressure can be distributed. A differential equation is formed for each pressure in such a fluid power system where the derivative of pressure can be expressed with the general formula:ṗ where p i is the pressure in at i th section, B e is the effective bulk modulus, V i is the volume in the same section, Q i and Q i+1 are the inlet and outlet volumetric flows, and dVi dt is the time rate changes of volume V i . The detailed description of the fluid power circuits studied will be in the following subsections.

A. CIRCUIT 1: MULTIPLE ORIFICES
The first fluid power circuit is related to the system of three identical orifices in series. The system contains the constant pressure pump, which is assumed as an ideal pressure source and tank for recovery of the fluid. Fig. 1 shows the schematic representation of the circuit. The position of a small volume in such a system is between the orifice number 2 and 3. The total simplicity of this circuit allows the use of a large integrator time step for the dynamic simulation of the circuit even with a small volume inside, and it requires less computational time than more complex fluid power systems.
Each pressure in the circuit can be represented by the lumped fluid theory and a general formulation (1). Therefore, pressures in this circuit are integrated from the following formulations:ṗ where V 1 and V 2 are pipe volumes in two pressure sections; Q 1 , Q 2 , and Q 3 are volume flows through orifices 1, 2, and 3, respectively, and B e1 and B e2 are pressure-dependent effective bulk modulus for pressures p 1 and p 2 , respectively. The effective bulk modulus in such a system can be derived as follows: where E max is the maximum bulk modulus of the oil, p max is the maximum pressure in the system, p i denotes pressure corresponding to the bulk modulus, and a 1 , a 2 , and a 3 are the empirical constants [28]. Volume flows in the circuit can be calculated with the orifice equation as follows: where C d denotes the discharge coefficient, A i is the crosssectional area of orifice i (where i = 1, 2 or 3), ρ is the density of the hydraulic fluid, and p i and p i−1 are the outlet and inlet pressures, respectively. The system has a total of 3 volumetric flow variables corresponding to each orifice. The values of all parameters used in the modeling and simulation of the multiple orifice fluid power circuit are represented in Table 1.

B. CIRCUIT 2: HYDRAULIC CYLINDER CONTROLLED BY A DIRECTIONAL CONTROL VALVE
The second circuit under consideration is the system with a more complex architecture and which is more practical and applicable for real-time simulation. The circuit consists of a two-chamber double-acting hydraulic cylinder with a mass attached to the end of the horizontal cylinder's rod. The mass is not totally fixed, having 1 degree of freedom to slide in a piston movement direction. The cylinder is controlled by the 4/3-proportional directional control valve. The pressure in the system is supported by the constant pressure pump, which is assumed to be an ideal pressure source. A pressure compensator is also used in the circuit between the pump and the directional control valve. All the components in the system are connected with the hydraulic pipes of different volumes.
The circuit contains one extremely small pipe volume, located between the pressure compensator and the directional control valve. Fig. 2 depicts the whole circuit; the small volume is denoted as V v . All the parameters of this system are presented in Table 2. The mathematical model of the circuit is represented with differential and algebraic equations. The system is controlled and activated by the voltage signal U , obtained from the following equation: where K v is the valve gain, U s is the signal proportional to the valve spool displacement, ζ is the valve damping ratio and ω n is the natural angular frequency. The flow Q v through the throttle of the pressure compensator is calculated according to a semi-empirical approach, as in [27]: where K denotes the semi-empirical flow coefficient, C 1 , C 2 , C 3 , and C 5 denote empirical constants [29], p s and p v are the constant pump pressure and the pressure in small volume, respectively. The output pressure of the shuttle valve VOLUME 4, 2016  between way A and B (Fig. 2) p shuttle is dependent on the maximum value of pressures p 1 and p 2 and can be expressed as p shuttle = max(p 1 , p 2 ).
The volume flow rates Q 1 and Q 2 of the 4/3-proportional directional control valve are calculated with the use of the turbulent orifice model with a triangular groove cross-section, as follows: where C ν is the flow constant that accounts for crosssectional areas and the geometry of the valve orifices, U d is the insensitivity area for the applied signal, and p 1 , p 2 , p v , and p t are the pressures in two cylinder chambers, the pressure in small volume V v , and the pressure in the tank, respectively. Pressures in the system are integrated from the pressure continuity equations in accordance with the lumped fluid theory, as follows: where V v is the small volume between the pressure compensator and directional control valve, Q v is the volume flow through the throttle of pressure compensator,ẋ is sliding speed of the piston, A 1 and A 2 are cross-sectional areas of chambers of the cylinder, V 1 and V 2 are volumes of the pipelines and chamber in way A and way B, respectively. B e1 , B e2 , and B e3 are effective bulk moduli for each pressure in the system and are represented by (3). Volumes V 1 and V 2 are calculated as follows: where S c is the full stroke of the cylinder, x is the position of the piston, A 1 and A 2 are areas of the first and second chamber, respectively, and V dead is the dead volume, which represents the volume of pipelines. The outlet flow Q 3 from (8) depends on the valve position, which can be determined as follows: The equation of motion for the system with a total force acting on it can be represented by the following relation: where F µ denotes friction between the walls of the cylinder and piston, m is the mass of the load attached to the end of the rod, andẍ is the acceleration of the piston. The friction model used for the simulation is based on the LuGre friction model [30]- [32]. This model can be introduced with the following set of equations:ż where F µ is the total friction force; z is the non-measurable internal state, F c is the Coulomb friction force,ẋ is the sliding velocity of the piston, v s is the sliding speed coefficient, F s is the static friction force, k v is the viscous friction coefficient, and σ 0 and σ 1 are the flexibility and damping coefficients, respectively.

III. RECURRENT NEURAL NETWORK ARCHITECTURE DESCRIPTION
In this section, the features of the RNN selection are presented. In addition, the description of the selected RNN is provided in detail. The most common RNN architectures used in modeling dynamic systems are the nonlinear finite impulse response (NFIR), the nonlinear autoregressive network with exogenous inputs (NARX), and the nonlinear autoregressive moving average network with exogenous inputs (NARMAX). There are also more complex and advanced architectures used in dynamic systems modeling, which are Long Short-Term Memory (LSTM) and Gated Recurrent Unit (GRU) neural networks. The NFIR architecture, described in [33], is the simplest RNN architecture of all the architectures mentioned above. The operating process of this network occurs through feeding all values of past inputs to achieve the current output value. The defining equation for the network is formulated as follows: where y(t) is the RNN output vector at time t, Ψ H is mapping performed by a multilayer feedforward network, d is the past values of series x(t), and x(t) is the RNN input vector at time t. The main advantage of the NFIR architecture is its stability while all past inputs are fed to the network.
The NARX neural networks principle is related to the utilization of the outputs of the network for feeding the input with past states of the outputs and inputs while remembering the state of the system at every step of the network operation. The ordinary NARX RNN can be defined by the following equation [34]: where y(t) is the RNN output vector at time t. The main feature of the NARX RNN is an accurate approximation of output values. However, in certain cases, it can be inherently less stable due to operation in a closed loop using the past values of the output. Another architecture that can be presented as an advanced NARX structure is NARMAX RNN. The main difference with this RNN architecture is the ability to use error of previous values in the loop. Thus, the defining equation for NARMAX networks is the following [35]: where e(t) is an RNN error vector at time t − 1. In the NARMAX architecture, all elements defined by x and e are sometimes called "controlled" and "uncontrolled" inputs [25]. This means that NARMAX is the most prevalent architecture in cases with real-world data, including the system error generated by noise. The structure also uses the error as an input dataset, which makes the method more complex than the above-mentioned NFIR and NARX. In addition to the mentioned RNN architectures, more complex architectures are used in dynamic systems modeling. These architectures are more suitable for training and may easily store long-term dependencies. Such architectures, which are LSTM and GRU neural networks, were studied and compared [36]. Both architectures show significantly good performance in case of complex dynamics persisting in the system. In the case of LSTM networks, the accuracy of predictions is at a high level, but selecting numerous hyperparameters can affect the performance of the network. GRU networks are similar to LSTM due to their functionality, however, in several cases, it can be applied less time to train the network.
The features of the selected RNN architecture for fluid power system simulation were also described in [25]. According to the author, the data obtained from the simulation model of a similar fluid power system does not contain real system noise, and the use of a large number of RNN parameters can affect the simulation time and speed. At the same time, the maximal accuracy of the system with the embedded RNN is required. In such a case, the NARX RNN is the most suitable architecture of neural network for fluid power system simulation even in the case of the simulation of one stiff differential equation since it is more accurate than NFIR due to the use of output data as feedback. Ultimately, VOLUME 4, 2016 the NARX architecture provides a trade-off between speed and accuracy in simulating the final system, which is the main objective of the research. Fig. 3 illustrates the basic structure of the NARX network used in the modeling of the system. The RNN consists of one input layer with four input values, and one feedback value that can be used only during the training of the network, two hidden layers with 40 neurons, and one output layer with one output activated with a linear function. Because the system has to be modeled as a mathematical model of one differential equation of pressure, the number of layers and neurons were manually selected by trial and error during the training process to ensure the accuracy of the network. The sigmoid function σ(x) was selected as the activation function for hidden layers of the network; the function is defined as follows: where x is the argument of the function σ(x) and e is the Euler's number.

IV. HYBRID MODEL DEVELOPMENT
The whole development of the hybrid model was divided into two stages. The first stage of development included the collection of the training data and training of the NARX RNN. The second stage of the development is related to the implementation of the neural network to the fluid power systems from Section II. The modeling and simulation of the systems were performed in MATLAB R2020a software in a form of MATLAB code, and the formulation of the RNN was performed through the embedded MATLAB Deep Learning Toolbox.

A. DATA COLLECTING AND TRAINING OF THE RNN
The data collected for the training was based on the simulation results of the practical fluid power circuit 2 described in Section II. The original system was simulated for 3,000 seconds with an integrator time step of 1.0×10 -5 s and an input voltage that was supplied randomly in a range between -10 and 10 Volts. The integrator time step was selected empirically to ensure the correct response of the system in the presence of small volume. The data set of 300 million samples was created from several parameters of the system, where each sample displayed the data of parameters obtained at every time step of the simulation. The number of samples was reduced to 3 million by saving each 100th sample to reduce the computational load of the computer and provide relatively fast training of the neural network. The input data chosen for the training were data arrays of volumetric flows Q 3 and Q 4 , effective bulk modulus B e3 and the fixed small volume V v obtained from the simulation of the original system mentioned above. Pressure in small volume p v was also saved and utilized as the output data for the training of the neural network. This data was chosen for training, validation, and testing of the neural network. Since all the input data, except the small volume, is variable, the neural network based on such data will work with any system with similar variable parameters. In case of changes of the small volume, new training of the network might be required.
At the beginning of the training, the training data was distributed for training, validation, and testing sets in the proportion of 70/15/15 percent. The input and output data were normalized in a range between -1 and 1 to achieve effective training results. The NARX RNN was trained multiple times in order to find the appropriate number of neurons in hidden layers for the most accurate and effective simulation. The Levenberg-Marquardt algorithm was utilized in the training process as a main backpropagation-based training algorithm due to its relatively fast training of the network and accurate results [37]. The Early-Stopping technique was also utilized in cases when the generalization stopped improving. The network was trained a total of ten times, for a maximum of 1,000 epochs with a different number of neurons in each hidden layer. The results of the training are displayed in Table 3. The results show that the number of neurons affects the training time of the neural network, and at the same time, the most accurate validation performance was obtained at training 8 and 9 with 40 and 45 neurons in a hidden layer, respectively. Based on the obtained data, the number of neurons in the RNN was selected to achieve the most accurate result, and the most accurate network in terms of Mean-Square error (MSE) was selected for the Hybrid model simulation.
The most accurate network (see Fig. 4a and 4b) contains 40 neurons in hidden layers. The validation performance of the network was expressed in the form of MSE, equal to 0.000216 (normalized data). The training time of the selected network was 3 hours and 33 minutes.

B. IMPLEMENTATION OF THE RNN IN THE HYBRID MODEL
After the training, the most accurate and fastest network was implemented in the MATLAB code of the traditional or reference Runge-Kutta mathematical model of both circuits described in Section II. The network was added to the code in the form of a MATLAB function as a substitute for the numerically stiff equation with small volume. First, the traditional mathematical model was simulated to obtain the input dataset for the hybrid system simulation. Such inputs are effective bulk modulus B e3 , volumetric flows Q 3 and Q 4 , and volume V v .
Note that the simulation of each hybrid system should be performed after obtaining the inputs for the network, and this means that the simulation of the traditional mathematical model should be completed and input values should be saved. After that, the hybrid system can be simulated an unlimited number of times.
The simulation of the hybrid system includes the stage of preprocessing the RNN before the main model simulation with the use of the data obtained from the previous simulation. The whole system is automatically simulated using the pressure data obtained from the RNN preprocessing stage at the corresponding simulation time step.

V. RESULTS AND DISCUSSION
This section presents the results of the simulation of the fluid power circuits described in Section II. In the case of Circuit 1, the simulation was performed with two variations of the integrator time step to show the responses of the pressure of the reference and hybrid systems at different simulation speeds. Circuit 2 was simulated with a different input signal to ensure the ability of the hybrid system with the RNN was accurate and fast with any set of input signals. The results are presented in the form of a comparison of plots related to the responses of the traditional system modeling and the hybrid approach of system modeling. The results obtained in the simulation also provide an understanding of the advantages and features of the hybrid method compared to the classical mathematical modeling and simulation of fluid power circuits. of the small volume described in [10], [11]. In this case, the response of the pressure at the time step of 1.0 × 10 −4 is correct. In other words, it is impossible to use the model at large time steps due to its numerical instability.
The Table 4 represents the results of the simulation in the form of the simulation time and Relative Root-Mean-Square  The hybrid model was simulated at a time step of 1.0 × 10 −4 . The simulation shows an accurate pressure response with an RRMSE of 0.895% with respect to the reference system (time step of 1.0 × 10 −4 ), which is observed from the plot in Fig. 6a. The simulation time of the hybrid system is 23.56 s, which is shorter than the reference system (24.63 s) even at the same simulation time step.
The plot for the next simulation of the hybrid model performed at a time step of 1.0 × 10 −3 s is depicted in Fig.  6b. The plot shows the stable response of the hybrid model in relation to the reference model at the integrator time step of 1.0 × 10 −4 . The simulation time of the hybrid system is only 2.32 s which is 10 times shorter than the numerically stable response of the reference system, which was obtained at the 1.0 × 10 −4 s integrator time step (24.63 s). The RRMSE for the pressure response of the hybrid system at the time step of 1.0 × 10 −3 s is 7.333%.
According to the plots in Figures 6a and 6b, the only difference between the responses of the hybrid system at different time steps is the longer pressure rise time, which indicates the stability of the hybrid model. In other words, the difference between the responses of the hybrid model at time steps of 1.0 × 10 −4 and 1.0 × 10 −3 is neglectable compared to the reference model at the same time steps. At the same time, responses of the hybrid model are close to the reference at the time step of 1.0 × 10 −4 . It means that the hybrid model can be used without significant accuracy losses at larger simulation time steps, which makes this model a trade-off between accuracy and simulation speed.

B. CIRCUIT 2: RESULTS
Circuit 2 was simulated with a constant pressure supply of 1.4 × 10 7 Pa and two different sets of input signals ( Fig.   7a and 7b) using both traditional (reference) and hybrid models. First, the circuit was simulated for 50 seconds with a randomly distributed input signal in a range between -10 and +10 volts (Fig. 7a). The numerical stiffness of the reference system with a small volume between the pressure compensator and the directional control valve allowed the maximal simulation time step of 1.0 × 10 −5 s to obtain a numerically stable response. The stiffness of the hybrid model was generally reduced by using the RNN instead of the stiff differential equation related to the pressure in small volume, which allowed an increase in the time step of the simulation to 1.0 × 10 −4 s.  As the result, the simulation time for the real 50 seconds was 501.47 seconds for the reference model. Compared to the simulation time of the hybrid model (53.27 s), the reference was almost ten times slower. However, the simulation time difference between the two models is associated with the time step difference at the same level of accuracy, which can be VOLUME 4, 2016 observed among the simulation results in Table 5 (see Fig.  8).
The responses of pressure p v related to small volume as well as the cylinder piston position x s for the reference and hybrid models with a random input signal are illustrated by plots in Fig. 8. The accuracy of the model is well-observed in the plot. The RRMSE for the reference model defined in (17) and represented in the simulation results (Table 5) was calculated for the hybrid system and equals 0.089% for the cylinder piston position x s and 1.479% for pressure in small volume p v .
The second simulation of Circuit 2 was performed within 10 seconds utilizing the reference and hybrid models. The repeatable input signal in a range between -6 and +9 Volts with a time period of 2.5 seconds was set for both systems. The simulation time steps remained the same for the reference and hybrid models, which were 1.0 × 10 −5 s and 1.0 × 10 −4 s, respectively. The responses of cylinder piston position x s and pressure in small volume p v are plotted in Fig. 9. The plots clearly represent the accurate responses of the hybrid model with respect to the reference with a resultant RRMSE of 0.260% and 2.441% for cylinder piston position x s and pressure p v , respectively.
The simulation times of the second simulation for the reference and hybrid models are presented in Table 5. Hybrid model simulation was performed for only 9.89 seconds, while the simulation of the reference model required almost ten times as much time (102.96 seconds).
At the end of the simulation, it can be concluded that at both the random and repeatable sets of input signals, the performance of the hybrid model is stable and accurate despite the RNN utilized instead of the stiff differential pressure equation was trained on the random data obtained from the 3,000-second simulation. Moreover, the simulation time of the hybrid system is significantly shorter than with the traditional mathematical model. And the conversion of the model, for instance to C++ code, may allow simulating it in real-time, or even faster than real-time applications.

VI. CONCLUSION
This paper proposes a novel hybrid method to solve a fundamental problem in the dynamic simulation of fluid power circuits. The problem is associated with the existence of stiff differential equations in the mathematical model of the system in the presence of small volumes that restrict the simulation speed for real-time applications. The method is based on the use of the NARX recurrent neural network in the model, which substitutes the stiff pressure continuity equation. The dynamics of the remaining parts of the system are preserved and modeled with conventional differential and algebraic equations. The network was trained using data from the practical fluid power circuit simulation in order to receive an accurate response.
To demonstrate the features and advantages of the method, the classical mathematical model and the hybrid method were compared. The simple circuit of three orifices in series and the practical circuit with a two-chamber cylinder controlled by a directional control valve were modeled and simulated using the two above-mentioned approaches. A simple circuit was tested at different simulation time steps, whereas the practical circuit was tested at random and repeating inputs. In both case studies, the response of the hybrid method was several times faster than the conventional reference model due to the elimination of the numerical stiffness problem. At the same time, the accuracy of the hybrid method is relatively high, which allows a trade-off between accuracy and simulation speed. Finally, the hybrid method performance allows us to use it in real-time applications, and future research can be associated with the testing and comparison of different time series RNN architectures (e.g. LSTM, GRU architectures in comparison with the NARX architecture) to find the most efficient method for real-time simulation of fluid power circuits.   HEIKKI HANDROOS has been a Full Professor in Machine Automation with LUT University since 1992. He received M.Sc (Eng.) degree and D.Sc. (Tech.) degree from Tampere University, Finland, in 1985 and 1991 respectively. He is a member of IEEE and ASME. He has served as an Associate Editor for ASME Journal of Dynamic Systems, Measurement and Control from 2013 to 2020 and himself has published more than 300 journal and conference articles in the field of mechatronics, robotics and control. He has been the supervisor of about 30 doctoral students and he has served as a principal investigator of international and domestic projects, whose total funding exceeds C30M.