Stability Boundary for Haptic Rendering: Influence of Damping and Delay Jorge Juan Gil∗ Thomas Hulin Assistant Professor Senior Scientist ´ Emilio Sanchez Carsten Preusche Assistant Professor Senior Scientist Applied Mechanics Department, CEIT Gerd Hirzinger and TECNUN, University of Navarra Professor ´ E-20018 San Sebastian, Institute of Robotics and Mechatronics Spain German Aerospace Center (DLR) Email:
[email protected] Oberpfaffenhofen-Weßling, D-82234 Germany Email:
[email protected]
ABSTRACT The influence of viscous damping and delay on the stability of haptic systems is studied in this paper. The stability boundaries have been found by means of different approaches. Although the shape of these stability boundaries is quite complex, a new linear condition which summarizes the relation between virtual stiffness, viscous damping and delay is proposed under certain assumptions. These assumptions include a linear system, short delays, fast sampling frequency and relatively low physical and virtual damping. The theoretical results presented in this paper are supported by simulations and experimental data using the DLR Light-Weight Robot and the LHIfAM.
1 Introduction A haptic interface can be used to link a human operator to a virtual environment in such way that the user is able to perceive feedback from the environment though the sense of touch. The ability to perceive contact forces of virtual environments is essential in many industrial applications, such as virtual prototyping [1] and maintainability analysis [2, 3], surgery training [4, 5] and driving simulators [6]. Maintaining stability is an elementary prerequisite for all of these applications. This paper studies the stability conditions of impedance-type haptic systems and presents a relation between virtual stiffness, viscous damping and delay. Section 2 presents a mathematical model of haptic systems and Section 3 summarizes previous stability and passivity conditions for such systems. In Section 4 a new stability condition is presented. The validity of this condition is verified in sections 5-8 in four different ways. Section 9 discusses the validity of the linear stability condition and conclusions are stated in Section 10.
2 System Description From the control point of view, a haptic system is a sampled-data controlled mechatronic device. Fig. 1 shows the model of the haptic device colliding with a virtual wall with delay td . This time delay can be the sum of several effects: computations, communications, etc. The interface has a mass m, a physical damping b and a Coulomb friction c. A viscoelastic impedance model with stiffness K and damping B is used to compute the force of the environment. It is assumed that the system has no velocity sensor; therefore, the backwards difference is used to estimate velocity. The sampling period is T and the position sensor resolution q.
∗ Address
all correspondence to this author. This paper was presented in part at the IEEE International Conference on Robotics and Automation, Rome, Italy, April 10-14, 2007.
+DSWLF LQWHUIDFH
&RXORPEIULFWLRQ
0DVVDQG 9LVFRXVGDPSLQJ
B_ L
L
B
'HOD\
~
d EX L K +ROGHU
d X
0QX
Fig. 1.
Table 1.
d X
(QYLURQPHQW
4XDQWL]DWLRQ
R
F d
L4 xF x
xX
Model of a haptic system
Dimensionless parameters
Parameter
Variable
Dimensionless variable
Sampling period
T
-
Mass
m
-
Physical damping
b
δ=
Virtual stiffness
K
α=
Virtual damping
B
β=
Delay
td
d
bT m KT 2 m
BT m td = T
+DSWLFLQWHUIDFH
B_ L B
d EX1L KX 'HOD\
0Q X
Fig. 2.
+ROGHU
d X
xX
(QYLURQPHQW
L4
F d
xF x
Linear model of a haptic system
Some phenomena, such as possible saturation and quantization in the actuators or internal vibration modes of the interface, are not taken into account. Although user dynamics are also involved in the loop, their influence makes the system more stable [7, 8]. Therefore, these dynamics have not been included in the block diagram, and it can be considered as “a worst-case scenario”. A simplified model of the system is shown in Fig. 2. Since this model contains only linear phenomena, it is valid only if Coulomb friction and quantization are negligible. Moreover, in [8] and [9] it was found that Coulomb friction can dissipate the amount of energy introduced by quantization. Under these assumptions, the linear model can be used as good approximation to find stability conditions. In accordance with [10], the dimensionless parameters that will be used in this paper are shown in Table 1. Although some authors [8] have used the virtual stiffness to normalize the parameters, we prefer to use the mass (or the inertia for rotary degrees of freedom), which means that the values of the device do not change with the contact force law. Both real and dimensionless parameters can theoretically take any value (m > 0, b > 0, T > 0 and td ≥ 0). However, typical sampling rates in haptics are quite fast (≥ 1 kHz) and the relation between physical damping and mass cannot be b 1 s−1 for all investigated arbitrarily large. For example, some experimentally-acquired values given in [8] show that m haptic devices. Therefore, the dimensionless physical damping δ should be quite small in haptic systems. In this paper we will suppose that δ <10−3 . As shown later, some conclusions in this paper can only be drawn based on this assumption.
3 Previous Analyses This section discusses previous approaches for analyzing the system introduced above, based on stability and passivity. 3.1 Stability approaches Classical control tools have been applied to the linear system in order to obtain stability conditions. In [7] it was stated that, with no delay, d = 0, the stability condition of the linear system (using the dimensionless parameters of Table 1) is: α < δ(δ + β)
(1 − )(β + βδ − β + δ 2 ) , (1 − − δ)(β + βδ − β + δ 2 )
(1)
where is a dimensionless number: = e− m = e−δ . bT
(2)
Stability condition (1) is consistent with [11]. Substituting (2) with its Taylor approximation: 1 1 = 1 − δ + δ 2 − δ 3 + O(δ 4 ), 2 6
(3)
makes it possible to linearize (1) around the origin and to obtain the following more compact stability condition [7]: α < 2(δ + β).
(4)
The fact that the dimensionless physical damping δ is a considerably small number enforces the validity of this approximation. Stability condition (4) was also experimentally found in [12]. Further experimental studies confirm the conclusion that increasing both the physical viscous damping, either electrically [13, 14] or magnetically [15], and the virtual damping [16, 17] allows for larger stable stiffness coefficients. If the system contains a delay of one sampling period, d = 1, the stability condition that has been proposed in [18] using the Pad´e approximation is: α<
2 (δ + β). 3
(5)
3.2 Passivity approaches Another way to guarantee the stability of the system is by ensuring its passivity [19]. The passivity condition for the linear system without delay proposed by Colgate [20] is: α < 2(δ − |β|).
(6)
The influence of the virtual damping on the passivity condition differs from the stability condition (4) for β > 0. Since passivity is a more restrictive condition than stability, the passive region in the (α, β)-plane is a subregion inside the stable region. However, it is worth noting that the same condition for both passivity and stability can be obtained if no virtual damping is included. Passivity analysis has been successfully extended [8] to the non-linear system depicted in Fig. 1, but without including the effect of the virtual damping β. In [9], equivalent results were obtained without the influence of the delay. 4 Stability Condition In this paper, we propose a stability condition for the linear system, including the effect of both delay and virtual damping. Our stability condition may be seen as a generalization of previous conditions (4) and (5) for any delay, consistent with the study of the non-linear system [8], but including the effect of the virtual damping β: α<
2 (δ + β). 1 + 2d
(7)
Using the physical values of the parameters, the proposed stability condition is:
K<
2 (b + B). T + 2td
(8)
The validity of this formula will be checked in four different ways: 1) with theoretical analysis, 2) by solving the characteristic equation numerically and performing a graphical approach, 3) by running simulations and 4) with experimental results. Rearranging the stability equation, which relates the real parameters of the system, yields:
K<
b+B . + td
T 2
(9)
Taking into account that the effect of the sampling and hold in the control loop can be approximated by a delay of half the sampling period T2 , we can interpret our stability condition (7) with the following statement: Damping . Critical stiffness = Delay
(10)
The critical stiffness of a haptic system is equal to the overall damping of the mentioned system divided by the total delay of the loop. Therefore, a double viscous damping in the system, physical plus virtual, will allow for double stiffness while a double delay in the haptic loop, regardless of its nature, will cut the maximum stable stiffness by half.
5 Routh-Hurwitz Analysis In the first stage, the Routh-Hurwitz criterion is applied to the discrete-time model shown in Fig. 2. In the second stage, the same method is applied to an equivalent continuous-time model. 5.1 Exact discrete-time model The methodology followed in [7] can be used to obtain the analytical stability condition from the characteristic equation of the system. In the Z domain, this equation consists of a polynomial if the delay td is a multiple of the sampling period T , causing d to take on natural values: δ 2 (z − )(z − 1)z d+1 − (1 − − δ)(α + β)z 2 + [(1 − − δ)(α + β) + (1 − − δ)β] z −(1 − − δ)β = 0.
(11)
The order of this polynomial (11), d+3, increases with d and using the Routh-Hurwitz criterion becomes quite complex. In Section 3 the analytical solution of the stability boundary for d = 0 was presented (1). The following lines derive the stability condition for a delay that is equal to the sampling period (d = 1). In this case, the characteristic equation is: δ 2 z 4 + (1 − )δ 2 z 3 + δ 2 − (1 − − δ)(α + β) z 2 + [(1 − − δ)(α + β) + (1 − − δ)β] z −(1 − − δ)β = 0.
(12)
After the bilinear transformation used in [7], it is possible to apply the Routh-Hurwitz criterion and identify the most restrictive condition. This exact condition for stability is: a1 a2 a3 − a4 a21 − a3 a20 > 0,
(13)
0.7 d=0 0.6 Unstable 0.5 0.4 α 0.3 0.2 d=1 0.1 Stable 0
Fig. 3.
0
0.5
1 β
1.5
2
Stability boundaries using the analytical conditions, for small dimensionless physical damping (δ
and with a delay equal to the sampling period (d
< 10−3 ), without delay (d = 0)
= 1)
where: a0 =(1 − ) δα 2 ,
a1 =(1 − )(δβ + δ 2 − α) + δα, a2 =(1 − )(α − 2β) + δ(3δ − δ + 3β − β − α), a3 =(1 − )(α + 4β) + δ(3δ + δ − 3β − β − α), a4 =−(1 − )(α + 2β) + 2δ (1 + )(α + 2β + 2δ). It is not possible to isolate parameter α as was done in (1) because some powers of this parameter appear in (13). Fig. 3 shows the stability boundaries using analytical conditions (1) and (13), and for small dimensionless physical damping (δ < 10−3 ). Taking different values of δ within this range, it is quite difficult to appreciate any change in the d=0 d=1 ≈ 0.686 and αmax ≈ 0.144 can be taken overall shape of the stability boundaries. Therefore, their maximum values αmax as constants [21]. As expected, the stability region with delay is smaller than that without delay. Although it has been stated that the overall shape of the stability boundaries depicted in Fig. 3 barely changes with δ, it is possible to detect some minor but important differences for relatively small values of β. Fig. 4, which is a close-up of Fig. 3 close to the origin, summarizes this behavior. The stability boundaries start at a virtual damping equal to the negative physical damping (β = −δ), making it possible to introduce a negative virtual damping up to this value. Substituting the Taylor approximation (3) in stability condition (13), yields: α<
2 (δ + β), 3
(14)
where it was possible to isolate parameter α. This condition fits the beginning of the corresponding stability boundary in Fig. 4 and also confirms result (5), which was obtained with different approximations and confirmed with experimental results. 5.2 Equivalent continuous-time model The Routh-Hurwitz criterion is applied to an equivalent continuous-time model of the haptic system, which is presented in Fig. 5, to obtain stability condition (8) for any delay td . In this linear continuous model, the effect of the sampler and the holder has been replaced by a delay of T2 seconds and the backwards difference has been replaced by the exact derivative. The characteristic equation of the continuous system is: 1+
K + Bs −(td + T )s 2 e = 0. ms2 + bs
(15)
6 Unstable 5 4 d=0 α δ
3 d=1 2 Stable 1 0 −1
0
1
2
3
4
5
β δ
Fig. 4. Close-up of the stability boundaries using the analytical conditions, for small dimensionless physical damping (δ delay (d = 0) and with a delay equal to the sampling period (d = 1)
< 10−3 ), without
+DSWLFLQWHUIDFH
B_ L B 'HOD\
Fig. 5.
d EX1L KX
0QX
6DPSOHU DQGKROGHU (QYLURQPHQW
x 1X
L 4X
Linear continuous model of the haptic system
To apply the Routh-Hurwitz criterion, the following approximation of the delay [22] is used: e−tr s ≈
1 , 1 + tr s
(16)
which is valid only for low frequencies, and the characteristic equation is: mtr s3 + (btr + m)s2 + (b + B)s + K = 0, where tr = td +
T 2
(17)
. The system is stable if the following two inequalities hold: b + B > 0, m(b + B − Ktr ) + btr (b + B) > 0.
(18) (19)
By neglecting the second addend in condition (19), the value of which is small in haptic systems compared to the first one, that condition can be reduced to linear stability condition (8). The continuous model behaves similarly to the discrete-time model at low frequencies, far below the Nyquist frequency. This fact may justify the use of the continuous model to find the critical stiffness of the system. Consider the continuous and the discrete-time models for B = 0 and td = 0. Their critical stiffness values are (20) and (21) respectively, where Gm{.} represents gain margin of the transfer function within brackets. 1 ms2 + bs 1 1 − e−T s KCR =Gm Z s ms2 + bs
KCR =Gm
(20) (21)
20
Gain (dB)
0 −20 −40 −60
1 s2 +s
−80
Z
1 −e−0.002s 1 s s2 +s
−100 −90 −105 Phase (°)
−120 −135 −150 −165 −180 −195 −1 10
Fig. 6.
0
1
2
10 10 Frequency (rad/s)
10
Bode diagram of the continuous transfer function (solid) and the discrete-time transfer function (dashed)
Fig. 6 shows the Bode diagrams of these transfer functions, using some hypothetical values: m = 1 kg, b = 1 Ns/m and T = 2 ms. The continuous transfer function (solid) has an infinite gain margin, therefore it is stable for any value of virtual stiffness K. The discrete-time system (dashed) has 60 dB of gain margin at 31.6 rad/s, therefore it has a finite critical value for virtual stiffness KCR ≈ 1000 N/m. b Assuming m Tπ rad/s in haptic systems, the crossover frequency is placed far below the Nyquist frequency as well. Therefore, the effect of the sampler and the zero-order-holder can be substituted by a delay of half the sampling period to obtain the gain margin and critical stiffness. In the example, the crossover frequency is 31.6 rad/s, while the Nyquist frequency is 1570.8 rad/s. Using this approximation, the critical stiffness for the system with certain delay is: KCR = Gm
1 −(td + T2 )s . e ms2 + bs
(22)
Fig. 7 shows the Bode diagram of the continuous transfer function with an overall delay of 5 ms (dashed). It is straightforward to conclude that the presence of delay in the control loop reduces the critical stiffness of the system. For moderate T values of td , the overall delay e−(td + 2 )s introduces a small phase at the crossover frequency and therefore can be replaced by the first-order transfer function (16) to find the critical stiffness.
6 Graphical Analysis The validity of (7) is checked in this section using the graphs of the stability boundaries. The critical stiffness of the linear system with delay has been obtained and depicted in two different ways. The first one follows [7] and directly obtains the critical stiffness for different values of the virtual damping by evaluating: ⎧ ⎨ α < Gm
⎩
1 z d (z−1)(z−)δ 2 (−1+δ)z+1−−δ
⎫ ⎬ ⎭ + β z−1 z
.
(23)
The second method, used in [10] and [23], numerically solves the poles of the characteristic equation (11) and finds the stiffness coefficients, which place all the poles just within the unit circle. R only if the delay is Although both methods obtain the same results, the gain margin can be easily computed in Matlab a multiple of the sampling period T , while the other method allows for introducing fractional numbers for the delay. Fig. 8 shows the stability boundaries for different delays d and setting δ < 10−3 while Fig. 9 shows a close-up of Fig. 8 near the point of origin. The shown boundaries in Fig. 9 perfectly fit the linearized stability condition (7).
20
Gain (dB)
0 −20 −40 −60
1 s2 +s 1 e−0 .005 s s2 +s
−80 −100 −90 −105 Phase (°)
−120 −135 −150 −165 −180 −195 −1 10
Fig. 7.
0
1
2
10 10 Frequency (rad/s)
10
Bode diagram of the continuous transfer function without delay (solid) and with a delay of 5 ms (dashed)
0.7 d=0
Unstable
0.6 0.5 0.25
0.4 α 0.3
0.5 0.75
0.2
1 0.1 0
Fig. 8.
1.5 2 3 0
0.5
1 β
Stability boundaries for small dimensionless physical damping (δ
1.5
2
< 10−3 ) and delays d = [0,0.25,0.5,0.75,1,1.5,2,3]
The initial slope of the stability boundaries becomes smaller with the delay. Therefore, the critical stiffness without virtual damping β = 0 also decreases with the delay. This means that, using the physical parameters, the critical stiffness depends on both the physical damping and the delay. Notice that this result is a graphical analysis which solely shows the consistency of (7). It can only be considered as proof after showing that it fits the exact analytical conditions, like the ones presented in the previous section.
7 Simulation Results A number of simulations have been performed to check stability condition (8). However, since they all hold that condition, only a few of them have been included in this paper. Table 2 shows the critical virtual stiffness for several values of delay and virtual damping. The sampling period was set equal to 1 ms, the mass equal to 1 kg and the physical damping 0.1 Ns/m. Simulations show that the virtual damping increases the critical stiffness in the same way as the physical one does. All the values obtained during the simulations hold condition (8) well.
6
d=0
0.25
0.5
5 0.75
Unstable
α δ
4
1
3
1.5 2
2
3
1 Stable 0 −1
0
1
2
3
4
5
β δ
Fig. 9. Close-up of the stability boundaries near the point of origin for small dimensionless physical damping (δ d = [0,0.25,0.5,0.75,1,1.5,2,3] Table 2.
Critical stiffness KCR versus delay and virtual damping, with m
KCR (N/m)
Delay (ms)
< 10−3 ) and delays
= 1 kg, b = 0.1 Ns/m and T = 1 ms
Virtual damping (Ns/m) 0.1
0.2
0.3
0.4
0
399.89
599.78
799.75
999.22
0.25
266.61
399.89
533.17
666.26
0.5
199.97
299.93
399.89
499.73
0.75
159.98
239.94
319.91
399.79
1
133.32
199.95
266.58
333.16
1.5
99.99
149.97
199.93
249.89
2
79.99
119.97
159.93
199.89
3
57.14
85.69
114.23
142.76
4
44.44
66.64
88.83
111.03
8 Experimental Results Two different haptic interfaces have been used to perform experiments: the DLR Light-Weight Robot III [24] and the LHIfAM [25]. A bilateral virtual wall consisting of a virtual spring and damper was implemented using one joint of each interface. Limit-stable parameter values were obtained when sustained oscillations were observed increasing the stiffness. No user was involved in the experiments. 8.1 DLR Light-Weight Robot The DLR Light-Weight Robot III (Fig. 10) is a 7 DoF robot arm with carbon fiber grid structure links [24]. Though it weighs only 14 kg, it is able to handle payloads of 14 kg throughout the whole dynamic range. The electronics, including the power converters, is integrated into the robot arm. Every joint has an internal controller which compensates gravity and Coulomb friction. Since high-resolution position sensors are used to measure link orientation (quantization q ≈ 20”), non-linear effects can be neglected. The virtual wall was implemented in the third axis of the robot, indicated by the rotating angle φ in Fig. 10. The environment was implemented using a computer connected to the robot via Ethernet. The sampling rate was 1 kHz and the overall loop contained a delay of 5 ms. Fig. 11 shows the experimental results, introducing several fixed values for the virtual
G
Active joint
Fig. 10. Third generation of the DLR Light-Weight Robot arm
2000
Critical stiffness KCR (Nm/rad)
1800 td = 5 ms
1600
td = 6 ms
1400 1200 1000
td = 10 ms
800 600 400 200 0
0
5 10 15 Virtual damping B (Nm s/rad)
20
Fig. 11. Experimental stability boundaries for a delay td of 5, 6 and 10 ms (pluses and solid) and theoretical boundaries (dashed)
damping. A set of experiments was performed with only the system delay of 5 ms, while additional delays were artificially introduced into the loop to obtain an overall delay of 6 and 10 ms. The theoretical behavior is depicted with dotted lines. The experimental stability boundaries fit the linear condition remarkably well. A significantly long delay was also introduced into the system in order to obtain a curved stability boundary. Fig. 12 shows the experimental stability boundary for an overall delay of 55 ms. The beginning of the stability boundary for a delay of 10 ms is also shown in the same figure. The theoretical stability curve has been computed using the device’s inertia in the configuration selected for the experiments: 0.8 kg·m2 . 8.2 LHIfAM The LHIfAM (Fig. 13) is a haptic interface with a large workspace developed in CEIT [25]. The mechanism consists of a counterbalanced parallelogram moving on a linear guide 1.5 m in length. The virtual wall was implemented in the direction of the translational movement of the guide (x axis in Fig. 13). In this direction, both the inertia of the device and the sensor resolution are quite high: 5.4 kg and quantization q ≈ 3.14 μm, respectively. The Coulomb friction is compensated by the controller as well. The controller can acquire the information from the sensor, compute the force of the virtual wall and command the motor
200
Critical stiffness KCR (Nm/rad)
180 td = 10 ms
160 140 120
td = 55 ms
100 80 60 40 20 0
Fig. 12.
0
5 10 15 20 Virtual damping B (Nm s/rad)
25
Experimental stability boundaries for a delay td of 10 and 55 ms (pluses and solid) and theoretical boundaries (dashed)
?
Fig. 13. LHIfAM haptic interface
within the same sampling period, that is, theoretically without delay in the loop. Therefore, significant stiffness coefficients can be implemented with stable behavior. However, the motor is saturated with few mm of penetration in the virtual wall. In order not to saturate the actuator in the overall critical oscillation, artificial delays of 3, 6 and 12 ms have been introduced into the control loop. Fig. 14 shows the experimental stability boundaries. The theoretical stability boundaries have been computed using a physical damping of 4.6 Ns/m. Although the physical damping is quite high, since the sampling period was 1 ms, the dimensionless damping of the LHIfAM in the direction of x was δ = 0.85 × 10−3 and therefore still within the δ < 10−3 range. The angular frequency ωCR of the critical oscillation has been registered for several experiments. The angular frequency ωCR can be calculated for an ideal mass-spring-damper system by using the well-known resonance equation ωr =
ko b2 − o2 . mo 4mo
(24)
Table 3 shows the theoretical oscillation frequency, when assuming that the virtual spring corresponds to its physical
3000
Critical stiffness KCR (N/m)
2500 td = 3 ms 2000 td = 6 ms
1500
1000
td = 12 ms
500
0
0
2
4 6 8 10 Virtual damping B (Ns/m)
12
14
Fig. 14. Experimental stability boundaries for the LHIfAM and a delay td of 3, 6 and 12 ms (pluses and solid) and theoretical stability boundaries for a physical damping b of 4.6 Ns/m (dashed) and same delays
counterpart ko = KCR and the damper is the sum of the physical and discrete damping bo = b + B. For mass and physical damping, the device’s values are taken: m = mo = 5.4 kg and b = 4.6 Ns/m. Therefore, by knowing the virtual damping and the overall delay within the loop, it is possible to approximatively determine the inertia and the physical damping of the mechanical interface with this kind of experiment.
9 Valid Range of the Linear Condition The shape of the stability boundary can be divided into two different parts. The first one follows the linear condition (8) for relatively small values of virtual damping (Fig. 9). The second one is a curve (Fig. 8) which can only be obtained graphically, although we know the exact analytical conditions of two cases: without delay (1) and with a delay of one sampling period (13). In this section we discuss the boundary between these two parts. Fig. 15 shows the points in which the relative error of the linearization (8) to the exact boundaries is 2% and 5% for cases δ = 0 and δ = 0.01. Obviously, the linear condition (8) presented in this paper can be used even for high physical damping δ = 0.01. For this value, the 2% relative error is only violated for delay d ≥ 3 from the beginning. Note that the points of the fixed relative errors are located in the plane of the normalized parameters α and β. The range inside which the linear condition is held, in terms of non-normalized parameters, increases for systems with faster sampling rate T1 or greater mass m. Regarding the delay, several factors, such us computation of the collision detection algorithms for complex virtual environments, digital to analog conversion and amplifier dynamics, introduce a certain delay in the haptic system that is usually equal to or less than one sampling period: d ≤ 1 [8]. Therefore, the linear condition (8) is appropriate for haptic devices. In other kinds of systems, which usually involve longer delays, the linear stability condition should not be used. For example, it is quite common to suppose a delay equal to hundreds of ms in teleoperated systems.
10 Conclusions and Future Work This paper studies the influence of viscous damping and delay on the stability of haptic systems. The stability boundaries have been found by means of different approaches: 1) analytically using the Routh-Hurwitz criteria; 2) numerically computing the poles of the characteristic equation, 3) performing simulations and 4) experimentally. Although analytical expressions of the stability boundaries are quite complex, a linear condition relating stiffness, damping and system delay can be used for impedance haptic rendering under certain assumptions. Virtual and physical damping have the same influence on the system inside the valid range of the linear condition. Thus, in combination with the delay, the virtual spring causes energy gain [26] and not virtual damping. Comparing the linear stability condition to the exact boundary demonstrates that it is valid for a wide parameter range of the virtual environment. Yet, since the analyses presented in this paper assume the linearity of the system, its results can only be taken as an approximation if non-linear phenomena (like Coulomb friction and sensor resolution) are not negligible. Another limit is
Table 3.
Critical oscillations using the LHIfAM
Experimental
Theoretical
KCR (N/m)
B (Ns/m)
td (ms)
ωCR (rad/s)
ωCR (rad/s)
700
0
6
11.21
11.34
800
1
6
12.08
12.13
950
2
6
13.36
13.21
1300
4
6
15.32
15.45
1600
6
6
17.45
17.14
1900
8
6
18.47
18.67
2250
10
6
20.20
20.31
2600
12
6
21.66
21.83
350
0
12
7.75
8.01
520
2
12
9.59
9.77
650
4
12
10.83
10.91
800
6
12
12.17
12.10
1000
8
12
13.54
13.52
1200
10
12
14.68
14.80
1300
12
12
15.32
15.40
the required low-frequency of the system compared to the sampling rate, which may be violated if the haptic device collides with a real environment, for example. In terms of future research, the investigation of these nonlinear effects is necessary to be carried out. Also the robustness against uncertain physical parameters and external disturbances has to be examined.
Achnowledgements This work has been supported in part by the EU Government, Enactive Network of Excellence, project number IST2002-002114 and the Basque Government mobility program MV-2006-1-6.
References [1] Chen, E., 1999, “Six Degree-of-freedom Haptic System for Desktop Virtual Prototyping Applications,” Proc. First International Workshop on Virtual Reality and Prototyping, Laval, France, pp. 97–106. [2] Savall, J., Borro, D., Gil, J. J., and Matey, L., 2002, “Description of a Haptic System for Virtual Maintainability in Aeronautics,” Proc. IEEE/RSJ Int. Conf. Intell. Robot. Syst., Lausanne, Switzerland, pp. 2887–2892. [3] Hulin, T., Preusche, C., and Hirzinger, G., 2005, “Haptic Rendering for Virtual Assembly Verification,” Proc. World Haptics Conf., Pisa, Italy. [4] Madhani, A. J., Niemeyer, G., and Salisbury, J. K., 1998, “The Black Falcon: A Teleoperated Surgical Instrument for Minimally Invasive Surgery,” Proc. IEEE/RSJ Int. Conf. Intell. Robot. Syst., Victoria B.C., Canada, pp. 936–944. [5] Li, M., and Liu, Y.-H., 2006, “Haptic Modeling and Experimental Validation for Interactive Endodontic Simulation,” Proc. IEEE Int. Conf. Robot. Autom., pp. 3292–3297. [6] Lee, W.-S., Kim, J.-H., and Cho, J.-H., 1998, “A Driving Simulator as a Virtual Reality Tool,” Proc. IEEE Int. Conf. Robot. Autom., Leuven, Belgium, pp. 71–76. ´ and Fl´orez, J., 2004, “Stability Analysis of a 1 DOF Haptic Interface Using the [7] Gil, J. J., Avello, A., Rubio, A., Routh-Hurwitz Criterion,” IEEE Trans. Control Syst. Technol., 12(4), pp. 583–588.
0.1
δ=0
d=0 0.25
0.08 Unstable 0.06
0.5
0.04
0.75 1
α
2 3 6
0.02 0 −0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
β 0.12
δ = 0.01
d=0
0.1
0.25 Unstable
0.08 0.5 α 0.06
0.75 1
0.04 2 3 6
0.02 0 −0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
β Fig. 15. Exact stability boundaries (solid) for different delays and places where the relative error of the linear condition (8) (dashed) is equal to 2% (circles) and equal to 5% (squares) for δ = 0 and δ = 0.01
[8] Diolaiti, N., Niemeyer, G., Barbagli, F., and Salisbury, J. K., 2006, “Stability of Haptic Rendering: Discretization, Quantization, Time-Delay and Coulomb Effects,” IEEE Trans. Robot., 22(2), pp. 256–268. [9] Abbott, J. J., and Okamura, A. M., 2005, “Effects of Position Quantization and Sampling Rate on Virtual-Wall Passivity,” IEEE Trans. Robot., 21(5), pp. 952–964. [10] Hulin, T., Preusche, C., and Hirzinger, G., 2006, “Stability Boundary for Haptic Rendering: Influence of Physical Damping,” Proc. IEEE/RSJ Int. Conf. Intell. Robot. Syst., Beijing, China, pp. 1570–1575. [11] Gillespie, R. B., 1996, “Haptic Display of Systems with Changing Kinematic Constraints: The Virtual Piano Action,” Ph.D. thesis, Stanford University. [12] Minsky, M., Ouh-young, M., Steele, O., Brooks Jr., F., and Behensky, M., 1990, “Feeling and Sensing: Issues in Force Display,” Comput. Graph., 24(2), pp. 235–243. [13] Mehling, J. S., Colgate, J. E., and Peshkin, M. A., 2005, “Increasing the Impedance Range of a Haptic Display by Adding Electrical Damping,” Proc. First WorldHaptics Conf., Pisa, Italy, pp. 257–262. [14] Tognetti, L. J., and Book, W. J., 2006, “Effects of Increased Device Dissipation on Haptic Two-Port Network Performance,” Proc. IEEE Int. Conf. Robot. Autom., Orlando, Florida, USA, pp. 3304–3311. [15] Gosline, A. H., Campion, G., and Hayward, V., 2006, “On The Use of Eddy Current Brakes as Tunable, Fast Turn-On Viscous Dampers For Haptic Rendering,” Proc. Eurohaptics Conf., Paris, France, pp. 229–234. [16] Colgate, J. E., and Brown, J. M., 1994, “Factors Affecting the Z-Width of a Haptic Display,” Proc. IEEE Int. Conf. Robot. Autom., San Diego, California, USA, 4, pp. 3205–3210. [17] Janabi-Sharifi, F., Hayward, V., and Chen, C.-S. J., 2000, “Discrete-Time Adaptive Windowing for Velocity Estimation,” IEEE Trans. Control Syst. Technol., 8(6), pp. 1003–1009. [18] Bonneton, B., and Hayward, V., 1994, “Implementation of a Virtual Wall,” Technical Report, McGill University. [19] Adams, R. J., and Hannaford, B., 1999, “Stable Haptic Interaction with Virtual Environments,” IEEE Trans. Robot. Autom., 15(3), pp. 465–474. [20] Colgate, J. E., and Schenkel, G., 1997, “Passivity of a Class of Sampled-data Systems: Application to Haptic Interfaces.” J. Robot. Syst., 14(1), pp. 37–47. [21] Hulin, T., Preusche, C., and Hirzinger, G., 2006, “Stability Boundary and Design Criteria for Haptic Rendering of
Virtual Walls,” Proc. 8th International IFAC Symposium on Robot Control, Bologna, Italy. [22] Ogata, K., 2002, Modern Control Engineering, 4th ed., Prentice-Hall Inc., Upper Saddle River, New Jersey. [23] Salcudean, S. E., and Vlaar, T. D., 1997, “On the Emulation of Stiff Walls and Static Friction with a Magnetically Levitated Input/Output Device,” J. Dyn. Syst. Meas. Control-Trans. ASME, 119, pp. 127–132. [24] Hirzinger, G., Sporer, N., Albu-Sch¨affer, A., H¨ahnle, M., Krenn, R., Pascucci, A., and Schedl, M., 2002, “DLR’s Torque-controlled Light Weight Robot III - Are We Reaching the Technological Limits Now?,” Proc. IEEE Int. Conf. Robot. Autom., Washington D.C., USA, pp. 1710–1716. [25] Borro, D., Savall, J., Amundarain, A., Gil, J. J., Garc´ıa-Alonso, A., and Matey, L., 2004, “A Large Haptic Device for Aircraft Engine Maintainability,” IEEE Comput. Graph. Appl., 24(6), pp. 70–74. [26] Basdogan, C., and Srinivasan, M. A., 2002, “Haptic Rendering In Virtual Environments,” Virtual Environments HandBook, K. M. Stanney, ed. Lawrence Erlbaum Associates, pp. 117–134.
List of Figures 1 Model of a haptic system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Linear model of a haptic system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 3 Stability boundaries using the analytical conditions, for small dimensionless physical damping (δ < 10−3 ), without delay (d = 0) and with a delay equal to the sampling period (d = 1) . . . . . . . . . . . . . . . . . 5 4 Close-up of the stability boundaries using the analytical conditions, for small dimensionless physical damp6 ing (δ < 10−3 ), without delay (d = 0) and with a delay equal to the sampling period (d = 1) . . . . . . . . 5 Linear continuous model of the haptic system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 6 Bode diagram of the continuous transfer function (solid) and the discrete-time transfer function (dashed) . . 7 7 Bode diagram of the continuous transfer function without delay (solid) and with a delay of 5 ms (dashed) . 8 8 Stability boundaries for small dimensionless physical damping (δ < 10−3 ) and delays d = [0,0.25,0.5,0.75,1,1.5,2,3] 9 Close-up of the stability boundaries near the point of origin for small dimensionless physical damping 9 (δ < 10−3 ) and delays d = [0,0.25,0.5,0.75,1,1.5,2,3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Third generation of the DLR Light-Weight Robot arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 11 Experimental stability boundaries for a delay td of 5, 6 and 10 ms (pluses and solid) and theoretical boundaries (dashed) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 12 Experimental stability boundaries for a delay td of 10 and 55 ms (pluses and solid) and theoretical boundaries (dashed) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 13 LHIfAM haptic interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 14 Experimental stability boundaries for the LHIfAM and a delay td of 3, 6 and 12 ms (pluses and solid) and theoretical stability boundaries for a physical damping b of 4.6 Ns/m (dashed) and same delays . . . . . . . 12 15 Exact stability boundaries (solid) for different delays and places where the relative error of the linear condition (8) (dashed) is equal to 2% (circles) and equal to 5% (squares) for δ = 0 and δ = 0.01 . . . . . . . . . 14
List of Tables 1 Dimensionless parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Critical stiffness KCR versus delay and virtual damping, with m = 1 kg, b = 0.1 Ns/m and T = 1 ms . . . 3 Critical oscillations using the LHIfAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 9 13
8
