Dynamic parameters identification of a haptic interface for a helicopter flight simulator

The haptic interface force feedback is one of the key factors for a reliable flight simulation. This paper addresses the design and control implementation of a simple joystick-like haptic interface to be used for a helicopter flight simulator. The expression of the haptic interface force is obtained by dynamic analysis of the haptic interface operation. This paper proposes a new strategy aiming at avoiding the use of an expansive and complex force/torque sensor. Accordingly, specific dynamic model is implemented by including Stribeck friction to describe the friction moment. Experimental data are processed as based on a genetic algorithm for identifying the dynamic parameters in the Stribeck friction model. This allows to obtain the friction moment parameters of the haptic interface, as well as the torque distribution due to gravity and the rotational inertia parameters of the haptic interface for the calculation of the haptic interface force. Experimental tests are carried out and results are used to validate the proposed dynamic model and dynamic parameter identification method and demonstrate the effectiveness of the proposed force feedback while using a cheap photoelectric sensor instead of an expansive force/torque sensor.

1 Introduction

Due to the high risks and costs of training on a real helicopter, pilots often need to complete a preliminary training on a helicopter flight simulator. Real flight training is carried out only after a successful training on a simulator. In helicopter flight simulators, the joystick is a key element, which is used to achieving a reliable helicopter driving while providing a force feedback feeling, which needs to properly mimic the force feedback feelings of a real helicopter flying (Prachyabrued and Robert, 2018). In the real flight control process, the pilot's hand will feel the feedback force and make judgments, then adjust the flight status through the feedback force. In the whole process, the current force of the haptic interface needs to be used as the simulation calculation of the feedback force. Thus, its measurement accuracy is particularly important. Most of the operating forces are measured by using torque sensors. The papers Cavallo et al. (2004), Rinaldi and Beckham (1983), Hutton et al. (1994), and Fergani et al. (2016) installed a pull-up/pressure sensor at the end of the rod to read the applied force at each moment. However, torque sensor occupies large space in both mechanical installation and wiring, which affects the normal operation of the pilots. This paper proposes a different strategy aiming at avoiding the use of an expansive and complex force/torque sensor. Namely, the proposed approach is based on establishing a reliable dynamic model of the haptic interface to replace the role of the force/torque sensor. Accordingly, a proper dynamic model and a parameter identification strategy have been proposed.

In a servo system, there are complex rolling friction and sliding friction phenomena between gear train, bearing, input/output shaft and seal. Therefore, the modelling of friction torque is very important, which directly affects the dynamic accuracy of the haptic interface. A variety of empirical models have been proposed for friction calculation of servo system, which can be divided into two phases: static friction model and dynamic friction model (Bona and Indri, 2005). The most widely used static friction model is the Stribeck model. Several papers analyzed the friction of a servo system by establishing its Stribeck friction model (Lichun and Pavelescu, 1982; Xu et al., 2011; Chen et al., 2016; Márton and Lantos, 2009; Broel-Plater et al., 2018). Considering the advantages and drawbacks of the existing models in literature, it has been decided to implement a Lugre model to analyze the friction of the servo system, as proposed in Canudas de Wit et al. (1995), Wang et al. (2019), Azizi and Yazdizadeh (2019), Freidovich et al. (2010), and Ishikawa et al. (2010). Main advantage of Lugre model is its relatively high accuracy, while its implementation is relatively complex, since it introduces an unknown quantity which can be difficult to be measured directly and requires a high precision servo system. The Stribeck model is currently used in many fields to model the friction moment. Since there are unknown parameters in the haptic interface dynamics model, such as gravity moment, Stribeck friction moment model parameters and rotational inertia, these parameters need to be identified. The main available parameter identification methods are tracing methods, least square methods, and genetic algorithms. The genetic algorithms are very effective evolutionary random search methods, which are proven in literature to be very effective in managing non-linear problems with a high level of computational robustness and calculation parallelization. Accordingly, genetic algorithms can be seen as the most convenient approach for addressing this identification problem. In this paper, a dynamic analysis is carried out on the a 2-dof flight haptic interface for a helicopter simulator. A nonlinear Stribeck model is established for the friction torque of the servo system. A curve fitting and genetic algorithm method are used for processing the data to identify the gravity torque, four parameters of the Stribeck model and the rotational inertia. Finally, the available parameters identified are brought into the dynamic model to calculate the haptic interface force in real time and an experimental validation is carried out to proof the effectiveness of the proposed approach.

2 Haptic interface operating principle and dynamics analysis

2.1 The proposed haptic interface and its operating principles

This paper adopts the helicopter simulator at laboratory in Yanshan University which is shown in Fig. 1. The helicopter flight simulator is mainly composed of four modules: a helicopter simulation cockpit system, a virtual scene system, a dynamic motion system and an audio system. The simulator's haptic interface is shown in Fig. 2. The haptic interface is driven by two servomotors, and two disc speed reducers are used for speed regulation. The motor adopts Kollmorge CKM04 AC permanent magnet synchronous servo motors, and its corresponding AKD series servo driver supports Modubus TCP communication. The parameters of the servo motors and the reducers are shown in Tables 1 and 2. The hardware system composition of the haptic interface is shown in Fig. 3.

Figure 1 The helicopter simulator in Yanshan University.

Figure 2 Force feedback test bench of haptic interface. (a) Control handle; (b) servo motor and rotating mechanism.

Figure 3 Hardware system composition of haptic interface.

First, the computer is connected to the switch, and the switch is connected with two motor drivers respectively through Enternet port. Then, C$++$ programmed code has been written to establish communication between servo motor and computer through Modbus TCP protocol. When the driver operates the haptic interface, the computer can obtain the current electricity value from the servo motor at each time. Then, the driving torque of the servo motor is obtained. At the same time, the photoelectric encoder of the servo motor can measure the stick's position and speed at the current time, which can be used for the dynamic calculation of the haptic interface.

Table 1 Parameter of the Kollmorge CKM04 servo motor.

Table 2 Parameter of the reducers.

2.2 Dynamic analysis of the haptic interface

First, coordinate system for the haptic interface is established for the haptic interface, as shown in Fig. 4. The intersection point of the two motor axes is defined as origin O. The two motor axes are set as x axis and y axis respectively. z axis is perpendicular to the xy plane.

Figure 4 Haptic interface coordinate system.

The free-body diagram is shown in Fig. 5. The angle θ, θ_x, θ_y formed between the haptic interface and the z axis is the output angle of reducer 1 and reducer 2 respectively. the relationship of the three parameters can be expressed as follows:

$$\begin{array}{}\text{(1)}& {\dot{\mathit{\theta }}}^{\mathrm{2}}={\dot{\mathit{\theta }}}_x^{\mathrm{2}}+{\dot{\mathit{\theta }}}_y^{\mathrm{2}}.\end{array}$$

When the haptic interface moves, the control handle, the haptic interface and the base of the haptic interface move at an angle θ, the motor 1, the reducer 1 and the rotating mechanism move at an angle θ_y.

Figure 5 The proposed free-body diagram.

In Fig. 5, m is the mass of the haptic interface system, l_G is the distance between the center of mass of the haptic interface and the origin, and l is the distance between the top of the stick and the origin. τ_a1 and τ_a2 are the output torques by motor 1 and 2 through the disc reducer respectively; τ_f1 and τ_f2 are the friction torques on the x axis and y axis of the servo system respectively. The manual force of the pilot in the direction of x axis, y axis and z axis defined as F_x, F_y, F_z.

The Lagrange method is adopted to conduct dynamic analysis of the system. The Lagrange function can be formulated as

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\partial E_{\mathrm{k}}}{\partial {\dot{q}}_i}}\right)-{\displaystyle \frac{\partial E_{\mathrm{k}}}{\partial q_i}}+{\displaystyle \frac{\partial E_{\mathrm{p}}}{\partial q_i}}=F_i(i=\mathrm{1},\mathrm{2}).\end{array}$$

In Eq. (2), E_k is the kinetic energy of the haptic interface system; E_p is the potential energy of the haptic interface system; ${\dot{q}}_{\mathrm{1}}$, ${\dot{q}}_{\mathrm{2}}$ in the haptic interface system represent the output angles θ_x, θ_y of the two reducers respectively.

The total kinetic energy E_k for the haptic interface system can be expressed as:

$$\begin{array}{}\text{(3)}& E_{\mathrm{k}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_{\mathrm{1}}{\dot{\mathit{\theta }}}_y^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_{\mathrm{2}}{\dot{\mathit{\theta }}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_{\mathrm{2}}{\dot{\mathit{\theta }}}_x^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(J_{\mathrm{1}}+J_{\mathrm{2}}){\dot{\mathit{\theta }}}_y^{\mathrm{2}}.\end{array}$$

In Eq. (3), J₁ is the sum of the moment of inertia of motor 1, reducer 1 and rotating mechanism; J₂ is the sum of the moment of inertia of the control handle, haptic interface and haptic interface base.

The potential energy of the haptic interface E_p can be expressed as:

$$\begin{array}{}\text{(4)}& E_{\mathrm{p}}=\mathrm{mg}{l}_{\mathrm{G}}\cos \mathit{\theta }.\end{array}$$

In combination with Fig. 5 and Eq. (1), the following equation can be obtained:

$$\begin{array}{}\text{(5)}& \cos \mathit{\theta }=\mathrm{1}/\sqrt{\mathrm{1}+t{g}^{\mathrm{2}}{\mathit{\theta }}_x+t{g}^{\mathrm{2}}{\mathit{\theta }}_y}.\end{array}$$

According to the comprehensive Eqs. (2)–(5), can be obtained:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}{J}_{\mathrm{2}}{\ddot{\mathit{\theta }}}_x-\frac{\mathrm{mg}{l}_{\mathrm{G}}tg{\mathit{\theta }}_x{\sec }^{\mathrm{2}}{\mathit{\theta }}_x}{(\mathrm{1}+t{g}^{\mathrm{2}}{\mathit{\theta }}_x+t{g}^{\mathrm{2}}{\mathit{\theta }}_y{)}^{\mathrm{3}/\mathrm{2}}}={\mathit{\tau }}_{\mathrm{1}}+{\mathit{\tau }}_{\mathrm{a}\mathrm{1}}-{\mathit{\tau }}_{\mathrm{f}\mathrm{1}}\\ (J_{\mathrm{1}}+J_{\mathrm{2}}){\ddot{\mathit{\theta }}}_y-\frac{\mathrm{mg}{l}_{\mathrm{G}}tg{\mathit{\theta }}_y{\sec }^{\mathrm{2}}{\mathit{\theta }}_y}{(\mathrm{1}+t{g}^{\mathrm{2}}{\mathit{\theta }}_x+t{g}^{\mathrm{2}}{\mathit{\theta }}_y{)}^{\mathrm{3}/\mathrm{2}}}={\mathit{\tau }}_{\mathrm{2}}+{\mathit{\tau }}_{\mathrm{a}\mathrm{2}}-{\mathit{\tau }}_{\mathrm{f}\mathrm{2}}\end{array}\right..\end{array}$$

In Eq. (6), τ₁ and τ₂ are the torques exerted by the pilot on the x axis and y axis during the control process respectively.

Especially, the torque output τ_a by the motor through the disc reducer is related to the reduction ratio i of the reducer and the output torque T_a of the motor, as:

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_{\mathrm{a}}=i{T}_{\mathrm{a}},\end{array}$$

T_a can be deduced from Eq. (7):

$$\begin{array}{}\text{(8)}& T_{\mathrm{a}}=\mathrm{1.5}{P}_n[{\mathit{\psi }}_{\mathrm{r}}{I}_q+(L_d-L_q\left)I_d{I}_q\right].\end{array}$$

In Eq. (8), P_n is polar logarithm; ψ_r is the flux chain of rotor magnetic pole; L_d, L_q is the armature inductance of axis d and axis q; I_d, I_q is the current of d axis and q axis.

Since the armature inductance of d axis and q axis of the motor used same value, namely:

$$\begin{array}{}\text{(9)}& L_d=L_q.\end{array}$$

From Eqs. (7), (8) and (9) can be obtained as follows:

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{a}}=\mathrm{1.5}i{P}_n{\mathit{\psi }}_{\mathrm{r}}{I}_q.\end{array}$$

At this time, the torque and the q axis current have a linear relationship, which can be expressed as:

$$\begin{array}{}\text{(11)}& {\mathit{\tau }}_{\mathrm{a}}=i{K}_t{I}_q.\end{array}$$

Combined with Eqs. (6) and (11), and using mgl_G replaced by the gravity moment T_G, the following equation can be obtained:

$$\begin{array}{}\text{(12)}& \left\{\begin{array}{l}{\mathit{\tau }}_{\mathrm{1}}=J_{\mathrm{2}}{\ddot{\mathit{\theta }}}_x-\frac{T_{\mathrm{G}}tg{\mathit{\theta }}_x{\sec }^{\mathrm{2}}{\mathit{\theta }}_x}{(\mathrm{1}+t{g}^{\mathrm{2}}{\mathit{\theta }}_x+t{g}^{\mathrm{2}}{\mathit{\theta }}_y{)}^{\mathrm{3}/\mathrm{2}}}-i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}+{\mathit{\tau }}_{\mathrm{f}\mathrm{1}}\\ {\mathit{\tau }}_{\mathrm{2}}=(J_{\mathrm{1}}+J_{\mathrm{2}}){\ddot{\mathit{\theta }}}_y-\frac{T_{\mathrm{G}}tg{\mathit{\theta }}_y{\sec }^{\mathrm{2}}{\mathit{\theta }}_y}{(\mathrm{1}+t{g}^{\mathrm{2}}{\mathit{\theta }}_x+t{g}^{\mathrm{2}}{\mathit{\theta }}_y{)}^{\mathrm{3}/\mathrm{2}}}-i_{\mathrm{2}}{K}_t{I}_{q\mathrm{2}}+{\mathit{\tau }}_{\mathrm{f}\mathrm{2}}\end{array}\right..\end{array}$$

Thus, the manual force of the pilot can be expressed as:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}F& =\left[\begin{array}{l}{F}_x\\ F_y\\ F_z\end{array}\right]={\displaystyle \frac{\mathrm{1}}{l}}\left[\begin{array}{ll}\mathrm{0}& \cos {\mathit{\theta }}_y\\ \cos {\mathit{\theta }}_x& \mathrm{0}\\ -\sin {\mathit{\theta }}_x& -\sin {\mathit{\theta }}_y\end{array}\right]\left[\begin{array}{l}{\mathit{\tau }}_{\mathrm{1}}\\ {\mathit{\tau }}_{\mathrm{2}}\end{array}\right]\\ & ={\displaystyle \frac{\mathrm{1}}{l}}\left[\begin{array}{ll}\mathrm{0}& \cos {\mathit{\theta }}_y\\ \cos {\mathit{\theta }}_x& \mathrm{0}\\ -\sin {\mathit{\theta }}_x& -\sin {\mathit{\theta }}_y\end{array}\right]\\ & \cdot \left[\begin{array}{l}{J}_{\mathrm{2}}{\ddot{\mathit{\theta }}}_x-\frac{T_{\mathrm{G}}tg{\mathit{\theta }}_x{\sec }^{\mathrm{2}}{\mathit{\theta }}_x}{(\mathrm{1}+t{g}^{\mathrm{2}}{\mathit{\theta }}_x+t{g}^{\mathrm{2}}{\mathit{\theta }}_y{)}^{\mathrm{3}/\mathrm{2}}}-i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}+{\mathit{\tau }}_{\mathrm{f}\mathrm{1}}\\ (J_{\mathrm{1}}+J_{\mathrm{2}}){\ddot{\mathit{\theta }}}_y-\frac{T_{\mathrm{G}}tg{\mathit{\theta }}_y{\sec }^{\mathrm{2}}{\mathit{\theta }}_y}{(\mathrm{1}+t{g}^{\mathrm{2}}{\mathit{\theta }}_x+t{g}^{\mathrm{2}}{\mathit{\theta }}_y{)}^{\mathrm{3}/\mathrm{2}}}-i_{\mathrm{2}}{K}_t{I}_{q\mathrm{2}}+{\mathit{\tau }}_{\mathrm{f}\mathrm{2}}\end{array}\right].\end{array}\end{array}$$

2.3 Friction torque model of the haptic interface

In the haptic interface structure, there are complex frictions within the servo motor and between the reducer. The friction torque generated by these frictions has a great influence on the accurate modelling of the haptic interface. Therefore, the friction torque at each moment must be accurately obtained to ensure the accuracy of the calculated force and torque at each moment.

The classic Coulomb friction model friction indicated that friction force is related to the positive pressure acting on the object. However, due to the lubrication in the servo system, the influence of viscous friction must be considered.

Stribeck proposed to model the servo friction as a combination of Coulomb friction and viscous friction. He has given a quantitative formulation for friction with respect to velocity (Balogh and Krstic, 2004). Similarly, the Stribeck formula can be converted into the formula of friction torque regarding motor speed as proposed in (Iwasaki et al., 1999):

$$\begin{array}{}\text{(14)}& {\mathit{\tau }}_{\mathrm{f}}=[{\mathit{\tau }}_{\mathrm{c}}+({\mathit{\tau }}_{\mathrm{s}}-{\mathit{\tau }}_{\mathrm{c}}\left)e^{-(\frac{\mathit{\omega }}{{\mathit{\omega }}_{\mathrm{s}}}{)}^{\mathrm{2}}}\right]\mathrm{sgn}\left(\mathit{\omega }\right)+B_{\mathit{\omega }}\mathit{\omega },\end{array}$$

In Eq. (14), ω is the motor speed; τ_c is Coulomb friction torque; τ_s is the maximum static friction torque; ω_s is Stribeck speed; B_ω is the coefficient of viscous friction.

The Stribeck curve is related to the speed of the servo system. In the low-speed area, it is affected by Coulomb friction and viscous friction, and the friction torque decreases. This area is called boundary lubrication friction area. With the increase of speed, the influence of viscous friction is much higher than Coulomb friction. At this time, the friction torque is approximately proportional to the rotational speed. This operation condition is defined as liquid lubrication condition.

There are two problems in applying the Stribeck friction model to the control of haptic interface force control: on the one hand, it is difficult to describe the change of haptic interface friction torque when the velocity is zero; on the other hand, the change of velocity near the zero point will cause a sudden change of friction torque, which will lead to a jitter in the reversing process. Therefore, approximate treatment of the low speed phase of the Stribeck friction model is required. This paper proposes to use a sigmoid function to deal with Stribeck friction model (Ciliza and Tomizuka, 2007).

The Stribeck friction torque model improved by sigmoid function can be expressed as follows:

$$\begin{array}{}\text{(15)}& {\mathit{\tau }}_{\mathrm{f}}={\displaystyle \frac{\mathrm{1}-e^{-\mathit{\gamma }\mathit{\omega }}}{\mathrm{1}+e^{-\mathit{\gamma }\mathit{\omega }}}}\left[{\mathit{\tau }}_{\mathrm{c}}+({\mathit{\tau }}_{\mathrm{s}}-{\mathit{\tau }}_{\mathrm{c}})e^{-(\frac{\mathit{\omega }}{{\mathit{\omega }}_{\mathrm{s}}}{)}^{\mathrm{2}}}\right]\mathrm{sgn}\left(\mathit{\omega }\right)+B_{\mathit{\omega }}\mathit{\omega }.\end{array}$$

According to the requirements of the model and accuracy, the value of γ is defined as 500 in the proposed model.

3 Identification of dynamic parameters

3.1 Identification of the gravity torque model parameter

According to Eq. (13), the total force applied to the haptic interface is directly related to the magnitude of the gravity torque T_G. In order to obtain an accurate model, it is necessary to identify the gravity torque T_G. This is achieved with an experiment conducted on the motor 1 when its movement is on xz plane. Set the haptic interface to rotate at a constant speed, and then ${\ddot{\mathit{\theta }}}_x$ is zero. The haptic interface was not moved in the y direction, so θ_y is zero. The driver has no additional force to the haptic interface, accordingly also τ₂ is zero. Since the velocity is constant, τ_f1 is constant. In addition, sin θ_x is relatively small in practical application, accordingly θ_x can be used instead. The experimental process is shown in Fig. 6 and Eq. (6) can be converted into:

$$\begin{array}{}\text{(16)}& i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}+T_{\mathrm{G}}{\mathit{\theta }}_x-{\mathit{\tau }}_{\mathrm{f}\mathrm{1}}=\mathrm{0}.\end{array}$$

Figure 6 The experimental process.

At this time, θ_x maintain a liner relationship with I_q1 have a linear relationship. Discretizing Eq. (16) one can obtain:

$$\begin{array}{}\text{(17)}& {\mathit{\theta }}_x\left(k\right)={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{f}\mathrm{1}}}{T_{\mathrm{G}}}}-{\displaystyle \frac{i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}\left(k\right)}{T_{\mathrm{G}}}}.\end{array}$$

A following experiment consists in setting the motor speed at 40 rpm (Rotation per minute), make the haptic interface rotate at a constant speed, take the values of θ_x and I_q1 in multiple groups. A curve fitting of the results in such experiments is shown in Fig. 7. The identified gravity torque is 11.14 N m. The toolbox of Curve Fitting Tool in Matlab is used to get the result.

Figure 7 Identification process of heavy torque.

3.2 Identification of friction torque model parameters via genetic algorithm

When we choose Stribeck friction model, it contains four unknowns: τ_c, τ_s, ω_s and B_ω, these four parameters need to be identified. Similarly, the motor 1 is set as the speed mode and a specific motor speed is given to make the haptic interface rotate at a constant speed. At the same time the driver has no additional force to the haptic interface, in this case, θ_y, ${\ddot{\mathit{\theta }}}_x$and τ₁ are equal to zero. Therefore, according to Eq. (6), the following equation can be obtained:

$$\begin{array}{}\text{(18)}& i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}+T_{\mathrm{G}}\sin {\mathit{\theta }}_x-{\mathit{\tau }}_{\mathrm{f}\mathrm{1}}=\mathrm{0}.\end{array}$$

According to Sect. 2.1, the current value, rotation angle and rotation speed of the motor can be read at each moment. It is can be concluded from Eq. (18), the friction torque τ_f1 at the corresponding rotational speed ${\dot{\mathit{\theta }}}_x$ can be obtained. The value of friction torque τ_f1 can be obtained by the different output angle ${\dot{\mathit{\theta }}}_x$.

Figure 8 Flow chart of genetic algorithm.

Figure 9 Parameter identification process of friction model. (a) The objective function changes with generations; (b) the elite individual changes with generations.

Genetic algorithm is an algorithm designed according to Darwinian evolution theory: through continuous natural selection, survival of the fittest, so as to obtain the most adaptable results (Liu, 2006). The process of parameter identification mainly includes generation of initial population, selection of race parameters, setting of fitness function, selection, crossover, variation and termination.

The search scope of parameters τ_c, τ_s, ω_s and B_ω can obtained by tracing method. The iteration interval is [4–5], [4–5], [0–0.1] and [20–40] respectively. The population size was set as 100, and the crossover probability and mutation probability were set as 0.4 and 0.1 respectively. The selection of GA parameters has been carried out by considering the peculiarities of the specific identification problem as well as several preliminary identification tests. In particular, after preliminary identification tests we selected 0.4 as crossover parameter and 0.1 as variation parameters. The selected crossover parameter gives a preference to parental individuals on filial individuals. The selected variation value is motivated by a limited influence of variation on this specific parameter identification. Numerical simulations have been carried out also with other values of the crossover parameter (e.g. 0.4 and 0.8). Results led to similar outcomes of the identification procedure.

The multi-group friction torques obtained through Eq. (18) can be expressed as:

$$\begin{array}{}\text{(19)}& {\widehat{\mathit{\tau }}}_{\mathrm{f}\mathrm{1}}=\left[{\widehat{\mathit{\tau }}}_{\mathrm{c}}+({\widehat{\mathit{\tau }}}_{\mathrm{s}}-{\widehat{\mathit{\tau }}}_{\mathrm{c}})e^{-(\frac{{\dot{\mathit{\theta }}}_x}{{\widehat{\mathit{\omega }}}_{\mathrm{s}}}{)}^{\mathrm{2}}}\right]\mathrm{sgn}\left({\dot{\mathit{\theta }}}_x\right)+{\widehat{B}}_{\mathit{\omega }}{\dot{\mathit{\theta }}}_x.\end{array}$$

In Eq. (8), ${\widehat{\mathit{\tau }}}_{\mathrm{c}}$, ${\widehat{\mathit{\tau }}}_{\mathrm{s}}$, ${\widehat{\mathit{\omega }}}_{\mathrm{s}}$ and ${\widehat{B}}_{\mathit{\omega }}$ are the four parameters to be identified respectively, and ${\widehat{\mathit{\tau }}}_{\mathrm{f}\mathrm{1}}$ is the theoretical friction torque calculated by the motor speed. Take the target function:

$$\begin{array}{}\text{(20)}& J=\sum _{i=\mathrm{1}}^n({\mathit{\tau }}_{\mathrm{f}\mathrm{1}}-{\widehat{\mathit{\tau }}}_{\mathrm{f}\mathrm{1}}{)}^{\mathrm{2}}.\end{array}$$

The actual identification task of genetic algorithm is to solve the minimum value of objective function.

The minimum target value is defined as 10⁻⁶ here, and the maximum evolutionary algebra is 10 000 generations. In each generation, the population is sorted in ascending fitness order, and the optimal solution of this generation is recorded as the elite. Parental individuals were selected in the fitness sequence for crossover, and the offspring were mutated. Finally, when the fitness of a certain generation of elite meets the requirements ($J\le {\mathrm{10}}^{-\mathrm{6}}$) or the evolutionary algebra is more than 10 000 generations, the whole selection process of genetic algorithm is terminated.

The process of genetic algorithm is shown in Fig. 8, and its identification process is shown in Fig. 9. Figure 10 shows the comparison between the theoretical friction moment calculated by the identified parameters and the actual friction moment. The measurement result of the friction torque at the moment of ${\dot{\mathit{\theta }}}_x<\mathrm{0}$ is the same as the moment of ${\dot{\mathit{\theta }}}_x>\mathrm{0}$, but the symbol is opposite, so it was omitted. Here we only figure out the absolute value of the friction torque. The identified parameters at the end of program operation are shown in Table 1. The identification of friction model parameters in the y axis direction was the same as that in the x axis. The identification results are listed in Table 3. Several cases have been considered in the identification process for the friction torque model parameters to obtain a reliable mean value of all the parameters. Considering the characteristics of the proposed identification process we can estimate a confidence interval of 0.0001 % (the minimum target value is in Eq. 20), which can be considered suitable for the specific application, also in comparison with commonly available accuracy of commercial sensors.

Figure 10 Curve of friction torque and rotational speed.

Table 3 Parameter identification.

3.3 Identification of rotational inertia model parameters

According to Eq. (14), the Stribeck friction torque model can be simplified as:

$$\begin{array}{}\text{(21)}& {\mathit{\tau }}_{\mathrm{f}}={\mathit{\tau }}_{\mathrm{c}}\mathrm{sgn}\left(\mathit{\omega }\right)+B_{\mathit{\theta }}\mathit{\omega }.\end{array}$$

Firstly, the motor 1 is analyzed. When the pilot does not operate the haptic interface, Eq. (6) can be converted into:

$$\begin{array}{}\text{(22)}& i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}+T_{\mathrm{G}}{\mathit{\theta }}_x-{\mathit{\tau }}_{\mathrm{c}}\mathrm{sgn}\left({\dot{\mathit{\theta }}}_x\right)-B_{\mathit{\theta }}{\dot{\mathit{\theta }}}_x=J_{\mathrm{2}}{\ddot{\mathit{\theta }}}_x.\end{array}$$

Thus, the structure block diagram of Eq. (21) can be described and is shown as in Fig. 11.

Figure 11 Structural block diagram of motor driven haptic interface.

After sorting out Fig. 11, the simplified structure block diagram of motor driven haptic interface is shown in Fig. 12.

Figure 12 Simplified structure diagram.

According to Fig. 12, the continuous function can be expressed as:

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}\left(s\right)}{i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}\left(s\right)-{\mathit{\tau }}_{\mathrm{c}}\left(s\right)}}={\displaystyle \frac{\mathrm{1}}{J_{\mathrm{2}}{s}^{\mathrm{2}}+B_{\mathit{\theta }}s-T_{\mathrm{G}}}}.\end{array}$$

The zero order holder method is adopted for Z transformation in Eq. (22), and the following equation can be obtained:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& {\displaystyle \frac{{\mathit{\theta }}_x\left(z\right)}{i_{\mathrm{1}}{K}_t{I}_{q\mathrm{1}}\left(z\right)-{\mathit{\tau }}_{\mathrm{c}}\left(z\right)}}=Z\left[{\displaystyle \frac{\mathrm{1}-e^{-Ts}}{s}}G\left(s\right)\right]\\ & ={\displaystyle \frac{\mathrm{1}}{\sqrt{B_{\mathit{\theta }}^{\mathrm{2}}+\mathrm{4}{J}_{\mathrm{2}}{T}_{\mathrm{G}}}}}\left[{\displaystyle \frac{z-\mathrm{1}}{p_{\mathrm{2}}(z-e^{p_{\mathrm{2}}T})}}-{\displaystyle \frac{z-\mathrm{1}}{p_{\mathrm{1}}(z-e^{p_{\mathrm{1}}T})}}\right.\\ & \left.+\left({\displaystyle \frac{\mathrm{1}}{p_{\mathrm{1}}}}-{\displaystyle \frac{\mathrm{1}}{p_{\mathrm{2}}}}\right)\right].\end{array}\end{array}$$

In which,

$$\begin{array}{}\text{(25)}& \left\{\begin{array}{l}{p}_{\mathrm{1}}=\frac{-B_{\mathit{\theta }}-\sqrt{B_{\mathit{\theta }}^{\mathrm{2}}+\mathrm{4}{J}_{\mathrm{2}}{T}_{\mathrm{G}}}}{\mathrm{2}{J}_{\mathrm{2}}}\\ p_{\mathrm{2}}=\frac{-B_{\mathit{\theta }}+\sqrt{B_{\mathit{\theta }}^{\mathrm{2}}+\mathrm{4}{J}_{\mathrm{2}}{T}_{\mathrm{G}}}}{\mathrm{2}{J}_{\mathrm{2}}}\end{array}\right..\end{array}$$

By discretizing Eq. (24), the following equation can be obtained:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathit{\theta }}_x\left(k\right)& =(e^{p_{\mathrm{1}}T}+e^{p_{\mathrm{2}}T}){\mathit{\theta }}_x(k-\mathrm{1})-e^{(p_{\mathrm{1}}+p_{\mathrm{2}})T}{\mathit{\theta }}_x(k-\mathrm{2})\\ & +{\displaystyle \frac{i_{\mathrm{1}}{K}_t}{p_{\mathrm{1}}{p}_{\mathrm{2}}\sqrt{B_{\mathit{\theta }}^{\mathrm{2}}+\mathrm{4}{J}_{\mathrm{2}}{T}_{\mathrm{G}}}}}\left(I_{q\mathrm{1}}\right(k-\mathrm{1})-d(k-\mathrm{1}\left)\right)\\ & \cdot (p_{\mathrm{2}}-p_{\mathrm{1}}+p_{\mathrm{1}}{e}^{p_{\mathrm{2}}T}-p_{\mathrm{2}}{e}^{p_{\mathrm{1}}T})\\ & +{\displaystyle \frac{i_{\mathrm{1}}{K}_t}{p_{\mathrm{1}}{p}_{\mathrm{2}}\sqrt{B_{\mathit{\theta }}^{\mathrm{2}}+\mathrm{4}{J}_{\mathrm{2}}{T}_{\mathrm{G}}}}}\left(I_{q\mathrm{1}}\right(k-\mathrm{2})-d(k-\mathrm{2}\left)\right)\\ & \cdot [p_{\mathrm{1}}{e}^{p_{\mathrm{1}}T}-p_{\mathrm{2}}{e}^{p_{\mathrm{2}}T}+(p_{\mathrm{2}}-p_{\mathrm{1}}\left)e^{(p_{\mathrm{1}}+p_{\mathrm{2}})T}\right].\end{array}\end{array}$$

In which,

$$\begin{array}{}\text{(27)}& d\left(k\right)=\left\{\begin{array}{ll}{\mathit{\tau }}_{\mathrm{c}+}& {\dot{\mathit{\theta }}}_x>\mathrm{0}\\ {\mathit{\tau }}_{\mathrm{c}-}& {\dot{\mathit{\theta }}}_x<\mathrm{0}\\ \mathrm{0}& {\dot{\mathit{\theta }}}_x=\mathrm{0}\end{array}\right..\end{array}$$

In Eq. (27), τ_c+ represents positive Coulomb friction and τ_c− represents negative Coulomb friction.

Set the torque mode through the motor driver, set the specific current size, and drive the haptic interface to rotate. When the motor speed reaches a certain value, take multiple groups of current value I_q1 and rotation angle θ_x.

For the obtained data, the genetic algorithm in Sect. 3.2 can be used to identify the moment of inertia. The estimated J₂ value among search range [0.1–1.1]. The population size was set as 100, and the crossover probability and mutation probability were set as 0.4 and 0.1 respectively. The selection of GA parameters has been carried out by considering the peculiarities of the specific identification problem as well as several preliminary identification tests. In particular, the selected crossover parameter has been defined to give a preference to parental individuals on filial individuals. The selected low variation value is motivated by a limited influence of variation on this specific parameter identification.

The multiple groups of rotation angles ${\widehat{\mathit{\theta }}}_x\left(k\right)$ calculated by Eq. (26) can be expressed as:

$$\begin{array}{}\text{(28)}& J=\sum _{k=\mathrm{1}}^n\left({\mathit{\theta }}_x\right(k)-{\widehat{\mathit{\theta }}}_x(k){)}^{\mathrm{2}}.\end{array}$$

The minimum target value is given as 10⁻⁴, and the maximum evolution algebra is 10 000 generations. The evolution of its optimal solution is shown in Fig. 13. The J₂ of the recognized moment of inertia is 0.2887 kg m².

Figure 13 Identification process of moment of inertia. (a) The objective function changes with generations; (b) the elite individual changes with generations.

Performming the same operation on the motor 2 as above, and the identified moment of inertia J₁+J₂ is 0.2989 kg m². Therefore, J₁ is 0.0102 kg m².

4 Experimental validation

Based on the experimental platform in Fig. 2, the driver is connected to the switch through the X11 terminal, and the switch is connected to the computer. Use Visual Studio C$++$ for programming. The computer communicates with the driver via Modbus TCP. Specifically, the program adopts the dynamic mapping function of Modbus. Modbus allows the driver to map any fixed register address to a new register address, and read-write access to the remap parameters by reordering the sequence block. The computer gives instructions to read or change the motor current state by reading or modifying the value in the register.

In this paper, the sampling time of the motor is given as 15 ms, and the quadrature axis current I_q1 and I_q2 size are set as 0. After noise removal, the rotation speed ω₁, ω₂ and position θ_m1, θ_m2 of the two motors at each moment were measured. ω₁ and θ_m1 can be converted to the reducer 1. the reducer 1 output angular speed ${\dot{\mathit{\theta }}}_x$ and position θ_x are corresponding to them respectively; ω₂ and θ_m2 can be converted to the reducer 2. The reducer 2 output angular speed ${\dot{\mathit{\theta }}}_y$and position θ_y are corresponding to them respectively; by substituting ω₁ and ω₂ into Eq. (14), the friction torque τ_f1 and τ_f2 of two degrees of freedom on the haptic interface at the current moment can be obtained. The angular acceleration ${\ddot{\mathit{\theta }}}_x$, ${\ddot{\mathit{\theta }}}_y$ of two degrees of freedom of the haptic interface can be calculated by difference method, and they can be expressed as:

$$\begin{array}{}\text{(29)}& \left\{\begin{array}{l}{\ddot{\mathit{\theta }}}_x=\frac{{\dot{\mathit{\theta }}}_x\left(k\right)-{\dot{\mathit{\theta }}}_x(k-\mathrm{1})}{T}\\ {\ddot{\mathit{\theta }}}_y=\frac{{\dot{\mathit{\theta }}}_y\left(k\right)-{\dot{\mathit{\theta }}}_y(k-\mathrm{1})}{T}\end{array}\right..\end{array}$$

In Eq. (29), T is the sampling time of the motor.

At this time, the driver's force on the haptic interface can be obtained through Eq. (13).

In the experiment, the pilot manipulates the haptic interface diagonally and nearly sinusoidal, as shown in Fig. 14a, and records the control force after the kinetic solution. The calculation result of haptic interface force in the whole process is shown in Fig. 14b.

Figure 14 Control process of the haptic interface approaching sinusoidal motion. (a) Operating conditions of the haptic interface; (b) calculation results.

Then take the take-off stage and the helicopter hovering to the left for example.

For the control process in the take-off stage, after lifting the total pitch bar, the pilot slowly and uniformly pulls the haptic interface backward. The haptic interface needs to be operated negatively along the x axis to make the helicopter climb steadily, and holds the haptic interface steadily after pulling back to a certain degree. When the helicopter climbs to the target altitude, lay down the haptic interface. The haptic interface needs to be operated forward along the x axis to make the helicopter attitude. The specific control process is shown in Fig. 15. From Fig. 15b, it can be seen that the pilot in t₁ moment is pulling the haptic interface in a negative direction towards the x axis and F_x, F_z both are negative. In the t₂ stage, the pilot holds the haptic interface for a certain seconds. At this time, because the motor speed is zero, F_x and F_z only equal to gravity. F_x, F_z are all constant positive values. In the t₃ stage, the pilot pulls the haptic interface forward along the x axis, therefore, both F_x, F_z are positive. At the t₄ stage, the haptic interface returns to its original position, and the values of F_x, F_z back zero. In addition, because the driver did not apply force in the direction of the y axis, F_y was constant zero in the whole process.

Figure 15 Helicopter take-off stage control process. (a) Operating conditions of the haptic interface; (b) calculation results.

Figure 16 Helicopter hovering to the left. (a) Operating conditions of the haptic interface; (b) calculation results.

For the control process of helicopter hovering to the left, the pilot slowly and uniformly pulls the haptic interface. The haptic interface needs to be operated forward along the y axis so that the helicopter fuselage tilts to the left. After pulling the haptic interface to a certain extent, the pilot stabilizes the haptic interface. At this time, with the operation of the helicopter pedal, the helicopter can make circular circling. The specific control process is shown in Fig. 16. As can be seen from Fig. 16b, in the t₁ stage, the pilot pulls the haptic interface forward along the y axis, which F_y is positive and F_z negative. In the t₂ stage, the pilot holds the haptic interface firmly, which F_y is negative value and F_z positive value. F_x was constant zero in the whole process.

It can be seen from the experiments that the control process curve of take-off stage and left hover stage is consistent with the actual driving situation. It can be indicated that proposed identification method for the haptic interface model parameters are feasible and can be applied in the helicopter flight simulator.

5 Conclusion

In this paper we have addressed the dynamic modelling and parameter identification of a joystick-like haptic interface to be used for helicopter flight simulators. Careful attention has been addressed at the Stribeck friction model and its parameters for an accurate modelling of servomotors. A specific genetic algorithm has been implemented for an accurate identification of the Stribeck friction parameters. A full dynamic model has been developed and implemented as based on Lagrange formulation. Experimental tests show that the proposed identification procedure and haptic interface allow a control of both take-off and left hover stage with results being consistent with the actual helicopter driving. Accordingly, the proposed haptic interface can provide a reliable force feedback even without using expansive force sensors. As future work, we aim at improving the proposed model for a more accurate modelling of friction especially when the velocity is close to zero. Future activities will also include further experimental testing.