Stability Analysis and User Perception of Haptic Rendering Combining Virtual Elastic, Viscous and Inertial Effects

: Virtual Reality environments are being used on a mass scale in ﬁelds, such as Industry and Medicine. These virtual scenarios serve very different purposes such as prototyping, gaming and exercising. Interaction with the virtual environment is mainly achieved by senses of sight and hearing through devices, such as a mouse or VR glasses. To this end, haptic research started a few decades ago with the aim of improving this interaction through a sense of touch. A key element, hitherto not researched, is the effective combination of virtual elastic, viscous, and inertial effects in haptic feedback restored to the user and the safety implications of these feedback effects. It is of particular importance in neurological rehabilitation exercising, as interaction realism and safety are of great importance in therapy and for the patient. Therefore, this work addresses the stability analysis of the combination of three haptic effects—elastic, viscous, and inertial—and the subjective feeling on the part of users regarding different combinations of these effects. A theoretical analysis is presented with a view to establishing stable control principles, and a user-study was carried out in order to help understand the perception of users to different combinations of haptic effects.


Introduction
Neurological rehabilitation is the key to helping hemiparetic patients regain movement. Rehabilitation therapy usually consists of extensive task-specific and repetitive movements of the patient's limbs in order for them to regain lost mobility. Traditionally, these therapies are carried out by physiotherapists and nurses and they require intensive physical and occupational resources. However, the advent of robotic rehabilitation is enabling patients to perform these exercises with the sole assistance of a robotic device, with similar or improved effectiveness [1]. Moreover, IoT technology allows these therapies to be carried out, even at the patient's house or at nearby clinical centers, further developing the concept of tele-rehabilitation and showcasing its efficiency [2,3]. Robotic rehabilitation is usually undertaken together with VR environments in order to engage and push patients to keep training, as patients' motivation and cognitive involvement have a great impact on rehabilitation outcomes [4,5]. These platforms have also shown their ability for learning purposes [6].
Over the last few years, many robotic systems have been developed in order to perform tailored rehabilitation exercises at home [7]. Usually, these platforms are connected to a computer, where serious games (mazes, puzzles, etc.) are shown in two-dimensional (2D) and three-dimensional (3D) scenes. The user is able to interact physically with these games through the robotic systems. Thus, it is desirable to properly recreate the sensation of handling objects in a virtual world [8,9]. A haptic device or a collaborative robot with force feedback can trace the user's movements in a virtual world and restore appropriate reaction forces during contact with virtual objects. Furthermore, force feedback can also be used in order to help patients complete the exercise (i.e., assist-as-needed approaches [10]), or it can be established as resistive forces that change during the course of therapy to compel the user to exert a physical effort (in terms of strength level and duration), according to the patient's actual muscular capacity [11].
Rehabilitation games are relevant scenarios in which to incorporate complex haptic effects, because virtual tasks are easier to understand if the haptic feedback is realistic [12]. These effects have several benefits, among them: the possibility to enhance the virtual scene and the ability to facilitate the rehabilitation, such as adding or subtracting the weight of the device by modifying the virtual mass or even making it weightless [13]. These virtual effects add immersiveness and realism to virtual interactions.
A wide variety of realistic haptic effects can be simulated by combining the three main mechanical components: spring, damper, and mass. Or, in terms of their characteristic parameters, by combining stiffness, viscosity, and inertia. To use all of these effects, an elementary prerequisite is to ensure the stability of the human-robot interaction and, to this end, Section 2 provides a summary of stability conditions. Some are already known in scientific literature, but others are novel-mainly analytical stability conditions for virtual mass, virtual mass-damper, virtual spring-mass, and virtual spring-mass-damper systems. However, despite there being an awareness of stability limits, little or no research in haptics has focused on subjective discrimination between elastic, viscous, and inertial effects. In Section 3, we present a user study that was carried out to gain some insight in this area. Section 4 gathers the results and Section 5 discusses them before the conclusion.

Stability Analysis of Haptic Interaction
In order to improve the architecture of a robotic rehabilitation platform where the patient must interact with virtual environments, virtual effects must be programmed with realistic tactile or kinesthetic sensations. However, modifying the physical parameters of the haptic effect can lead to instabilities in the robotic device. Therefore, we shall now analyze the most significant parameters. Based on a closed-loop force control, let us consider a mechanical interface with mass m and viscous damping b. The impedance haptic loop ( Figure 1) implements a discrete-time force law C(z) using a position sensor, while the sampling period is T.
Zero-order hold Impedance The characteristic equation of this system is: When defining a dimensionless number δ, which only depends on physical parameters of the system, the characteristic equation of the system is: The roots of the characteristic equations will be the poles of the system, with latter being considered stable if-and only if-all of the poles are inside the unit circle of the z-plane, i.e., the magnitudes of all poles are smaller than 1. In the following subsections, the stability of all combinations of the three main haptic effects (spring, damper, and mass) are analyzed.

Virtual Spring
The simplest haptic effect involves the implementation of an elastic force with stiffness coefficient K, which is C(z) = K. In this case, the characteristic equation of the system will be a second-order polynomial in the complex variable z: The Routh-Hurwitz criterion is applied following a bilinear transformation in order to analyze the stability of the system [14]. Consequently, the range of stable values of K will be: This expression contains two conditions that stiffness K must satisfy as upper and lower limits. The dimensionless parameter δ is very low in haptic systems [15] and, thus, the following Taylor series can be used in order to obtain some approximations as follows: Using the Taylor series (7), a good approximation of the previous stability condition (6) is: where dimensionless virtual stiffness α is defined, as follows: The results that are obtained from this section are very well known in literature [14,16,17]. If physical damping b is low, then the upper limit in (8) will also be relatively low. In this case, there is little room to simulate rigid virtual springs without the addition of other virtual effects.

Virtual Damper
If the haptic effect is a viscous damping, with damping coefficient B, and velocity is estimated while using backwards difference without filtering, the characteristic equation of this system will be a second-order polynomial: When applying the Routh-Hurwitz criterion, the exact stability condition is: Using the Taylor series (7), a good approximation for the stability condition (13) is: where dimensionless virtual damping β is defined as follows: In contrast to the previous section, if the haptic system behaves ideally (as depicted in Figure 1), it is possible to implement and simulate virtual dampers while using large damping coefficients. Thus, it is possible to use this haptic effect to dissipate energy from the system.

Virtual Mass
If the haptic effect is an inertia with virtual mass M and acceleration is estimated while using double backwards difference, the characteristic equation of this system will now be a third-order polynomial: The exact analytical stability condition is: While using the Taylor series (7), an approximation of (19) is: where dimensionless virtual mass γ is defined, as follows: The analytical result shown in (19) is new in the literature, and it is consistent with [18], where the stability region is depicted by numerically solving the characteristic equation while using experimental values of m, b, and δ from different devices.
The lower stability limit of (20) could be used in order to eliminate inertia from the mechanism almost completely. Hence, the haptic effect with negative inertias could be used to enhance the simulation, increasing the transparency of the system (the mechanical interface is not perceived by the user). However, because parameter M is set closer to the stability limit, there would be more poorly damped oscillations, reducing the transparency of the haptic simulation. As for the upper limit in (20), critical virtual inertia is only twice the physical inertia.

Virtual Spring-Damper
If the haptic interaction consists of an elastic force with stiffness coefficient K plus a viscous effect with damping coefficient B, then the characteristic equation of this system will be a third-order polynomial: The condition that is required in order to preserve stability [17] is: Figure 2 shows the stability boundary (25) in the dimensionless (α, β)-plane for two values of δ. The shape and size of the stability region only changes with parameter δ. However, typical values of δ are so low [15] that the influence of this parameter on the shape of the complete stability region is barely noticeable. The abscissa axis of Figure 2, for α = 0, meets the previous condition (15), although the lower limit is so low as compared to the upper one that the stability boundary would seem to start at the origin (0, 0) instead of (−δ, 0). Similarly, the ordinate axis for β = 0 meets condition (9). Again, the stability boundary would seem to start at the origin (0, 0) instead of passing through the point (0, 2δ).

Virtual Mass-Damper
If the haptic interaction consists of a viscous effect with damping coefficient B plus an inertia with virtual mass M, then the characteristic equation of this system will be a third-order polynomial: When applying the Routh-Hurwitz criterion, the condition required to preserve stability is: This analytical result is new in literature. Figure 3 shows the exact stability boundary of the haptic system in the (β, γ)-plane. The abscissa axis for γ = 0 meets condition (15), while the ordinate axis for β = 0 meets condition (21).

Virtual Spring-Mass
If the haptic interaction consists of an elastic force with stiffness coefficient K plus an inertia with virtual mass M, then the characteristic equation of this system will be a fourth-order polynomial: When applying the Routh-Hurwitz criterion, two conditions must be met at the same time for the system to be stable. The first is simply K > 0, although the second is a very long expression. When defining the following parameters, this second stability condition is: This analytical result is new in literature. Figure 4 shows the exact stability boundary of the haptic system in the (α, γ)-plane. The abscissa axis for α = 0 meets condition (21), while the ordinate axis for γ = 0 meets condition (9). In [18], the poles of this characteristic equation are numerically solved by substituting experimental values m, b, and δ from different devices. Consequently, several experimental stability boundaries are obtained, although the exact analytical stability condition is not reported.

Virtual Spring-Mass-Damper
If the haptic interaction contains the three virtual effects described so far-stiffness, damping, and inertia-, it can be expected that the stability region of this system meets that of previous sub-regions. In fact, these three regions match perfectly in a three-dimensional plot ( Figure 5). Curiously, the complete stability region of the positive dimensionless parameters (α, β, γ) would seem to be confined within a cube with a side length equal to 2.  The characteristic equation of this system with virtual stiffness, virtual damping, and virtual inertia will be a fourth-order polynomial: When applying the Routh-Hurwitz criterion, two conditions must be met at the same time for the system to be stable. The first is simply K > 0, although the second is a very long expression. When defining the following parameters, this second stability condition is: This analytical result is new in the literature. Note that, by substituting B = 0, the coefficients c i exactly become the previous coefficients a i . Figure 6 shows the exact stability boundary of the haptic system in the (α, β, γ)-plane while using condition (43) with δ = 10 −2 .

Recommendations for Developing Safe Haptic Controllers
Before implementing any haptic effect within the force loop of the rehabilitation system, some useful recommendations can be put forward from the results of the previous theoretical background. Regarding the inertial effect, the range of permissible virtual inertias is quite small, and it is greatly conditioned by the physical effect. The virtual mass cannot be more than double the physical mass (20), and this limit is lowered even further when virtual damping is introduced. In addition, the neighborhoods of the critical limits should be avoided. At most, in our rehabilitation device, we were able to simulate a virtual mass equal to the physical mass (M = m) that was free of disturbing oscillations and, therefore, the overall inertia (physical plus virtual) was able to be doubled. Most desktop haptic interfaces have small physical inertia. For example, all of the devices analyzed in [19] have only 250 g or less, as double the inertia of all these mechanisms is still a small amount.
Conversely, the viscous effect can be chosen for up to relatively high values (14), with moderate damping values (around 10 N·s/m), giving rise to perceivable forces, even with low-speed movements of the mechanism.
The elastic effect does not pose any special problems either. While using our rehabilitation device, and even without additional virtual damping or virtual inertia, moderate or high stiffness coefficient values can be simulated (around 1000 N/m).

Materials and Methods
Stable values for each haptic parameter or combination of parameters were obtained in the previous section. The theoretical analysis allows for safe controllers to be implemented, although it is still necessary to select the appropriate values from the user perception standpoint. To this end, a user study was carried out with the aim of better understanding how users perceive and discriminate between sensations of virtual stiffness, damping, and inertial haptic effects.

Rehabilitation Device
The HomeRehab robotic tele-rehabilitation system [7] was used for the experiments (Figure 7), and Table 1 shows its main technical specifications. Although the weight of the mechanism is 8.2 kg, the inertia of the moving parts of the device that is felt by the users is m ≈ 3 kg. This is relatively high when compared to other mechanical interfaces [19], although the device is very easy to handle, because the motors compensate for the gravitational components of the mechanism [13]. HomeRehab has the option of working in two-dimensional (2D) and three-dimensional (3D) workspaces. When working in 2D, patients exercise with HomeRehab sitting in a chair and training 2D movements in a planar workspace. In this case, experiments were carried out in 2D and not all of the device workspace was used.

Participants
Twenty healthy subjects took part in the experiments-four women and sixteen men-from 22 to 47 years old. All were right-handed and reported normal tactile function, and all of the subjects had prior experience with haptic applications.

Experimental Setup
A virtual environment was developed with two independent areas, referred to as 1 and 2 ( Figure 8), with each area being of 200 × 200 mm in size and having two virtual walls in the upper and lower zone. The users are able to move freely within wall limits, and they can also move from one area to the other. Virtual force law is implemented as a different combination of the three haptic parameters subject to analysis: stiffness, damping, and inertia. Note that, while stiffness is only felt when users collide with the virtual walls, damping and inertia are also felt within the complete area.
Users feel a different combination of haptic effects in each area. Indeed, the aim of the experiment was to analyze whether they are able to discriminate between each haptic parameter and recognize differences. The values of the virtual stiffness coefficients that are tested in the experiments are 500 and 1000 N/m, the values of the virtual inertias are 0 and 1 kg, and the virtual viscous coefficients 0 and 10 N·s/m (almost equal to the physical effect).
At the beginning of the experiment, a brief report is given to all users describing the robotic device, configuration of the virtual environment, and the meaning of the three haptic parameters being tested. They also allowed for testing the virtual environment and trying out some haptic feedback effects. Afterwards, each of the twenty users carries out seven different trials, with haptic parameters of Area 1 being the same in all of the trials: virtual stiffness coefficient of the walls K = 500 N/m and no additional effects in free movement (M = 0 and B = 0). However, the haptic parameters of Area 2 are modified between trials (Table 2), with the values being hidden from users and presented in randomized order.   In each trial, the users are allowed to interact with the scene as many times as desired with no time restrictions, and it is explained to them that haptic parameters of Area 2 change (increase, decrease, or maintain), among trials. Whenever they feel confident, they are requested to answer whether stiffness, inertia, and viscosity of Area 2 are higher or lower than, or equal to, those of Area 1.

Results
In each of the sven trials, the users reported three answers (one per haptic parameter). This yields 21 responses per user and 420 responses for the whole study. Discriminating according to the haptic effect, we have 140 responses per each parameter. Table 3 presents the outcomes of the experiment, in which the results of the responses per haptic parameter have been grouped together in order to better evaluate the outcomes of the experiment. The stiffness evaluation group, for example, gathers all 140 responses regarding the discrimination of virtual stiffness between trials, and it groups together responses that were obtained from the trials where the virtual stiffness was equal in both areas (trials 2, 3, and 6) or different (trials 1, 4, 5, and 7). The overall success rate (all trials and users) in discriminating between haptic effects is 75.71%. This means that the haptic effects are quite distinguishable for users and they should be applied for haptic rendering in order to enhance force feedback realism (with the stability limits that are pointed out in Section 2). Note that virtual inertial effects, for example, are not usually applied in most haptic applications.

Discussion
Traditionally, haptic rendering has only relied on virtual stiffness in order to render objects rigid in virtual scenarios, and rarely also included viscosity. However, the realism of any virtual contact can be greatly improved if all three effects are properly combined. To this end, two relevant issues have to be taken into account: the stability implications of combining the three haptic effects and the perception of each effect by the user. Stability control principles have been established in Section 2) and they do not require further discussion, rather than extending them to future analysis of other sources of instability (e.g., vibration modes of the mechanical device). As to how users are able to discriminate between these effects when they are displayed virtually, the results that are gathered in Table 3 are very open to further discussion.
Inertia, viscosity, and stiffness are physical parameters that, once properly explained to users, would seem to be easily distinguishable from each other. However, when analyzing the results of Table 3 for each haptic parameter independently, it can be seen that virtual damping is the effect that is best discriminated between trials, with a 87.14% success rate. This is quite a surprising result as users in general reported that virtual stiffness was the parameter that they felt to be most recognizable and were able to discriminate easier. In our opinion, this may be because stiffness is the physical parameter that is the easiest to understand, but, surprisingly, the results that were obtained from the experiment showed better discrimination rates for viscosity. Inertia was confirmed to be the most difficult haptic effect to recognize, although it should be noted that it was the only parameter that did not double the value among trials; rather, it only changed from 3 kg physical to 3 kg physical, plus 1 kg virtual.
Analyzing the results more in-depth, the fact of whether it is easier to discriminate the equal or different values of the parameters between areas can be evaluated. According to the results, virtual stiffness was properly identified as remaining equal in 45% of the trials, while the success rate for identifying different stiffness values was 96.25%; for virtual damping, the successful identification of equal values in both areas was 75% and, for different values, it was 96.25%; lastly, for inertia, the results were 51.67% and 76.25% success rates for equal and different values between areas, respectively. Thus, most users were capable of detecting changes in the stiffness and viscosity with ease (both with 96.25% success rate). It is also remarkable that the least accurate discrimination was that of stiffness when it remains equal in both areas. We think that these results may be because variations in the other two parameters are mainly evaluated in free motion, while virtual stiffness is discriminated when colliding with virtual walls where the effect of all parameters is mixed.
In summary, while the results show that how users perceive and discriminate between the three haptic effects is still not fully understood and requires further experiments, the obtained discrimination rates would seem to be high enough to consider including all three effects when designing virtual haptic scenarios. In our own research using the HomeRehab robotic tele-rehabilitation system, we are applying virtual stiffness and damping in order to render walls that guide patients through the appropriate paths and movements, and we use inertial effects when users hold elements virtually.

Conclusions
This work presents the stable limits for rendering the haptic effects that combine virtual stiffness, damping, and inertial parameters, and it also provides some insight into how humans perceive and discriminate these effects when they are presented together. The proper implementation of haptic feedback that combines these parameters can provide enhanced realism and improved performance in many applications, with particular relevance in domains, such as robotic rehabilitation, where the haptic stimuli presented to the user may play a relevant role in therapy outputs. Nevertheless, further studies are necessary in order to better understand how humans perceive and discriminate all haptic effects that are taken into consideration.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: IoT Internet of Things  2D two dimensional 3D three dimensional