Control of Haptic Systems Based on Input-to-State Stability

Discretization of signals often generates additional energy in the haptic systems, and makes them unstable. The currently popular controls based on the passivity stabilize the systems by limiting the rendered force of the haptic systems so that the generated energy from the system stays below zero. The passivity-based methods are, however, often conservative and sacrifice the performance of haptic rendering. This paper proposes a control method to adjust the change rate of the stiffness of virtual environment connected to the haptic systems to satisfy the input-to-state stability (ISS) for better stability and transparency. The ISS conditions for the systems are derived by modeling the system as a linear time-varying system. The systems become dissipative if they satisfy the ISS conditions. The generated surplus energy of the dissipative system is less than a positive finite value while maintaining the stability. Since the passive system is a special case of the dissipative system, the proposed ISS-based control is less conservative than the passivity-based approaches. Performance of the proposed control is experimentally compared with the virtual coupling and the force-bounding method that are widely used passivity-based approaches. The energy generated from the system using the proposed method is positive finite whereas the generated energy using the passivity-based approaches is less than zero. This means that more energy is transferred to the operator. Increased stiffness corresponding to this additionally transferred energy can be rendered. Root-mean-square and maximum errors of the rendered forces are, therefore, reduced by at least 86 % and 50 %, respectively.


I. INTRODUCTION
Virtual environment connected to haptic systems may behave active due to signal discretization by the sampler, zero-order hold (ZOH), and velocity estimation [1], [2]. This active behavior incurs additional energy in the system and makes it unstable. The discretization effect becomes worse as stiffness of the virtual environment and sampling time increase. Physics-based simulations involving deformable objects have a slow force update rates of 300 Hz to 500 Hz due to large amounts of computation [3], which increases the active behavior of the haptic systems. It is important to enhance transparency while maintaining stability in systems that require high-fidelity haptic feedback, such as medical simulations.
Time-domain passivity approach (TDPA) stabilizes the systems by observing and controlling the energy additionally generated by the systems. Hannaford and Ryu [4] proposed a passivity observer and controller to preserve the system passivity. The passivity observer monitors the system energy in real-time, and excessive energy generated by the system is dissipated by the artificial damping force called the passivity controller. There are studies to improve the transparency of TDPA by predicting the system energy and reducing the observed energy [5]. The operator, however, may feel discrete feedback force because the controller works discontinuously depending on the observer state. The discontinuous feedback force sometimes cause contact oscillations especially when the haptic device is moving at a low velocity [6]. There is a possibility that the observer may not detect unstable motion even if the system diverges due to the accumulated energy even in stable environment and free space. This is called a memory effect, and it becomes worse in the environment where soft and stiff objects are mixed [7].
Output-limiting methods such as adjusting output-limiter [8] and force-bounding approach (FBA) [9]- [11] have been reported to limit the maximum environment force to guarantee the passivity of the system. These methods are free of discontinuous feedback force unlike TDPA since they limit the rendered force. The haptic system is modeled as one or two-port network, and the maximum environment force is computed based on the TDPA. The environment force is transmitted faithfully if it is less than the computed maximum force. The memory effect, however, occurs as in the TDPA since the system energy is observed in the time domain. A more conservative condition from scaling down the device damping [9] is used to handle the memory effect, thereby, deteriorating the transparency more than necessary.
There have been many researches to stabilize the haptic systems by limiting the environment impedance. Colgate et al. [12] proposed the virtual coupling that connects the device and the virtual tool with a spring and a damper. The haptic systems are modeled as a two-port network, and the spring and damper of the virtual coupling are designed to satisfy the passivity of the systems [13], [14]. The virtual coupling limits the environment impedance, and keeps the system passive. The virtual coupling is used in a wide range of applications including real-time simulation with deformable objects [15]- [17]. There is a disadvantage that the transparency of haptic rendering is degraded by the virtual coupling. The operator always feels a softer environment than the original environment since the spring of the virtual coupling is connected to the environment in series.
Methods to improve the transparency have been proposed. The input-output coupling method [18]- [20] keeps the system dissipative by controlling the coupled terms to satisfy the input-to-state stability (ISS) condition. These methods analyze the ISS considering hysteresis of the input-output [21], and control the input-output by coupling derivatives of the inputoutput pairs. The coupled pairs are controlled to satisfy the constraint expressed as the function of the derivatives pairs and the maximum slope of the hysteresis. The energy generated in the system is bounded to a finite value, and the system remains dissipative. Maximum renderable stiffness is increased compared to the previous control methods using the passivity. Haptic fidelity, however, deteriorates because bounded vibration may be included in the rendered force. Additional damping force, manually tuned by the operator, is required to prevent the vibration. This damping must be readjusted experimentally for new environment or operator. Model matching framework has been proposed to render the stiff environment [6]. The controller is designed based on the ∞ optimization problem to minimize the error between the rendered force and the desired force. Accurate knowledge of the physical properties of the environment is necessary to design the controller. Singh et al. [22], [23] proposed successive force augmentation (SFA) approach to improve transparency by adding an offset force to the force generated by the virtual coupling. The offset force is updated successively according to the pressing path and the target stiffness of the environment. This method depends on accurate knowledge of the physical properties of the environment.
Previous control methods based on the passivity condition, that limit the rendered force or impedance of the environment or simulation, achieve the goal of stabilization by keeping the energy generated from the environment less than zero. These approaches often lead to conservative control, and sacrifice the performance of haptic rendering. The previous methods developed to enhance the transparency of the haptic rendering, such as the input-output coupling method, model matching framework, and SFA approach, require accurate model and knowledge of the environment physics.
This paper proposes a control method which adjusts the change rate of the stiffness of the virtual environment to satisfy the input-to-state stability (ISS) conditions to enhance the transparency without a priori accurate model of knowledge of the environment physics. The nonlinear environment is modeled as a linear time-varying spring . The characteristics of the nonlinear environment are included since is calculated at each time step. The ISS conditions for the haptic systems are derived by modeling the system as a linear time-varying system. The change rate of is monitored at each time step, and it is adjusted to satisfy the ISS conditions. The haptic systems become dissipative if they satisfy the ISS condition [24]. The generated surplus energy of the dissipative system is less than a positive finite value [18]. The passive system is a special case of the dissipative system, where the generated energy is strictly less than zero. The proposed ISSbased control keeps the system dissipative, and the rendered force is controlled so that the generated energy has a positive finite value as long as stability is maintained. More energy of the dissipative system, compared to the passive system, can be transmitted to the operator. Hence increased stiffness of the environment corresponding to the surplus energy can be rendered to the operator, increasing the fidelity of the haptic system. The proposed control uses only the current magnitudes of the position of the haptic device and the environment force. The control does not require any exact model of the environment. Experiment results to compare performance of the haptic rendering of the proposed method with the previous virtual coupling and FBA methods show improved fidelity of the rendered force while maintaining the stability.
Modelling of the haptic systems as the linear time-varying systems is explained in Section II. The ISS conditions for the haptic systems are derived in Section III. Control design based on the ISS conditions to stabilize the systems is presented in Section IV-A. Control design to enhance the rendering performance is presented in Section IV-B. Experimental comparison with the passivity-based approaches is explained in Section V.

II. SYSTEM MODEL
The displayed stiffness ( ) felt by the operator at ℎ time step is computed as the environment force ( ) divided by the penetration depth ℎ ( ) of the virtual tool from the environment surface as shown in (1)  (1) Fig. 2 shows a typical model of the haptic systems including the operator, device and virtual environment with discretized and quantized signals. Muscle force ( ) is generated by the operator's intention. Contact force ℎ ( ) generated by the operator and position ℎ ( ) are transmitted to the haptic device. ℎ ( ) is bounded. ℎ ( ) is discretized into ℎ ( ) by the sampler. The device is modeled as mass and damping . ( ) = ( ) ℎ ( ) is generated by the virtual environment.
( ) turns into analog signal ( ) = ( ) which is the force rendered to the operator. Physics-based simulations involving deformable objects have slow force update rates of 300 Hz to 500 Hz due to large amounts of computation [3]. Sampling time is assumed to be 2 ms in this paper for the sake of analysis, which is considered typical and plausible in haptic systems. The proposed method can be also applied to systems with sampling times of less than 2 ms for stiffer envirionments. Physical damping of the passive human operator improves stability [25], [26]. Stability condition without the operator's impedance is, therefore, more conservative than ensuring the stability with human operator, which is adapted in this paper to show merits of the proposed method. Fig. 3 (a) shows the haptic systems excluding the operator impedance. Fig. 3 (b) shows the equivalent system in the continuous-time domain. The sampler and ZOH is given in [27]. The − is replaced with Taylor series as in (3). It is known that the human operator applies a position input of up to 20 Hz in typical haptic rendering [7]. The terms of third or higher order in the Taylor series (3) are assumed negligible considering the maximum frequency of 20 Hz is small compared to the sampling frequency, 500 Hz. Hence, the sampler and ZOH are modeled as in (4).  (4). Hence ( ) is expressed as a function of ( ), ℎ ( ) and . ), respectively. The haptic system is expressed as a linear time-varying system where ( ) and ( ) vary according to ( ) and .

III. INPUT-TO-STATE STABILITY CONDITIONS
Conditions for the system to satisfy the ISS are derived in this section. The haptic system is dissipative if it is ISS [24]. Equation (7) shows the energy stored in the system. The generated surplus energy of the dissipative system is less than a positive finite value α as shown in (7) [18]. The passive system is a special case of the dissipative system, where the generated energy is strictly less than zero, i.e., in (7) is zero. The proposed ISS-based control in section IV keeps the system dissipative, and the rendered force is controlled so that the generated energy has a positive finite value. More energy generated in the dissipative system compared to the passive system is transmitted to the operator. Increased stiffness of the environment corresponding to the surplus energy can be rendered to the operator while maintaining the stability of the system. The ISS is, therefore, less conservative than the passivity, and the fidelity of haptic rendering is improved using the ISS-based control.
Theorem 1 [24]: The system is dissipative if and only if the system is ISS.
The ISS conditions is applied to the linear time-varying system (6). Four conditions to satisfy ISS are derived. Consider the following system, where ∈ ℝ and F ∈ ℝ are the state and input, respectively. : ℝ × ℝ → ℝ is assumed as locally Lipschitz function with (0,0) = 0.

Theorem 2 [24]: The system in (8) is ISS if and only if there exists a positive definite and proper ISS-Lyapunov function
: ℝ → ℝ such that the following inequality is satisfied.
where and belong to a class ∞ function. |•| denotes Euclidean norm.
An ISS-Lyapunov function candidate is designed as in (10). (10) should be positive definite for to be positive definite. The other terms are positive definite. Condition I for the system to be ISS is as follows.
Condition I: Time derivative of the ISS-Lyapunov function candidate is as follows.
The , , and in (13) represent the coefficients of − 2 2 , 2 , and − 1 2 , respectively. The should be positive definite to satisfy the ISS condition (9) since it is a coefficient of − 2 2 . This is the condition II as in (14). 2 (13) is expressed as inequalities in (15) to (19).
Coefficients of − 2 2 , − 1 2 , and 2 should be positive definite for the inequality (19) to satisfy the ISS conditions (9). Two additional conditions are required so that the coefficients of − 1 2 and 2 become positive definite.
Condition IV: Table I summarizes the four conditions for ISS. The conditions are expressed as the function of the physical properties of the device ( , ), , ̇, and .

A. CONTROL DESIGN TO SATISFY ISS CONDITIONS
An ISS region that satisfies the four conditions in Table 1 is computed, and a controller for stabilizing the system is proposed in this section. The ISS region of the system (6) is shown in Fig. 4. The x and y axes represent and ̇, respectively. Condition ) corresponding to the condition II decreases as increases. There is a constraint on the maximum by the condition II. Equation (22) represent the function of that satisfy the condition II. The −(1 + 2 )/ (2 ) term in the condition III becomes a large negative value as the approaches +0. An asymptotic line of condition III is the same as the boundary line of condition II as shown in Fig.  4 (a). Hence, condition III is more conservative than condition II. Equation (23) represent the function of and ̇ that satisfy the condition III. Equation (24) is a boundary curve ( ) of the condition III (23), and it is represented in Fig.  4 (a). The (1 + 2 )/(2 ) term in the condition IV is positive definite if the is positive definite. The condition IV is satisfied if the condition II is satisfied.
Condition III: ̇ and cannot exceed the boundary curve ( ) to satisfy the ISS conditions. Stability is guaranteed by limiting ̇ and to the ( ). Value of ̇ is less than or equal to 0 N/(m•s) when the limit of maximum occurs. Equation (25) A controller 1 is proposed to adjust the change rate of to satisfy the ISS conditions in Fig. 4.
( ) is calculated as ( ) divided by ℎ ( ) in each time step when the contact occurs. ̇( ) is calculated using backward Euler as in (26). The controller 1 shown in Fig. 5 operates when contact occurs.
( ) is controlled to satisfy the conditions of ( ) and . ( ) denotes the controlled ( ) . The controlled force ( ) = ( ) ℎ ( ) is transmitted to the operator.
The ISS conditions may not be satisfied when increases rapidly or is larger than . The controller adjust the increase rate of to satisfy the ISS conditions, thereby ( ) is limited as ( ) . There are two ( ) candidates to satisfy the ISS conditions. ( ) candidate I and II are the conditions to satisfy ( ) and α , respectively. Equation (27) shows ( ) candidate I determined by ( ) . △ ( ) in (27) is restricted as � ( − 1)� if the ISS conditions are not satisfied due to the large increase rate of . Equation (28) shows ( ) candidate II determined as .
( ) is limited as (28) if it is larger than .
( ) shown in (29) is selected as the value that satisfies both candidates (27), (28). Fig. 4 (a) indicates a case in which Fig. 4 (b) shows a case in which ( ) is determined as since is less than ( ) candidate I: ( ) candidates II:

B. CONTROL DESIGN TO ENHANCE PERFORMANCE
The controller 1 cannot render the stiffness of over . This section presents an additional compensator for rendering of increased stiffness. The unintended damping, − 2 ( ) and spring, − 2̇( ) occur due to the discretization effect caused by ZOH and sampler as in (6). The compensator is designed to reduce the discretization effect. The stability of the system with the compensator is analyzed based on the ISS conditions. A controller 2 is proposed to satisfy the ISS conditions. The controller 2 bounds the change rate of ( ), and there is no maximum stiffness constraint that exists in the controller 1 .
( ) ( ) in (2) is approximated using the firstorder Pad é approximation as shown in (31), and ( ) ( ) is calculated using continuous-todiscrete-time conversion. ( ) is obtained by taking the inverse of ( ) ( ). Table II and Fig. 6 show ( ) ( ) using the continuous-to-discrete-time conversion methods when is 2 ms. First-order hold (FOH), Tustin approximation, and impulse-invariant mapping are suitable among the conversion methods for ( ) to become a causal compensator. The more magnitude and phase of ( ) ( ) is similar to that of ( ) ( ) , the more discretization effect is effectively compensated.
( ) ( ) using impulse-invariant mapping has larger magnitude than ( ) ( ) and cannot express the phase lag of ( ) ( ). There is a limitation to improving the stability of the system since it cannot express the phase lag.
( ) ( ) using Tustin approximation accurately expresses the phase lag of ( ) ( ) , but the magnitude becomes very small near the Nyquist frequency. Hence ( ) using Tustin approximation amplifies highfrequency noise.
( ) ( ) using FOH has an error in the phase near the Nyquist frequency, but the phase is expressed accurately in the low-frequency region. The magnitude error near Nyquist frequency is smaller than that of ( ) ( ) using Tustin approximation. ( ) ( ) using FOH, therefore, is used for ( ) as in (32). (32) Fig. 7 (a) shows the haptic systems with ( ). Fig. 7 (b) shows the equivalence in continuous time-domain.
( ) ( ) ( ) is defined as ( ) , and ( ) is approximated to the first order term using Taylor series as in (33). It is known that the operator applies a position input of up to 20 Hz in haptic rendering [7]. The terms of the second or higher order in the Taylor series (33) are considered negligible considering the maximum frequency of 20 Hz is small compared to the sampling frequency, 500 Hz.  inverse Laplace transform of ( ) ( ) 1− − and approximated as the inverse Laplace transform of ( )(1 + ′ (0) ) using (33). Hence ( ) is expressed as a function of ( ), ℎ ( ) and . An ideal ( ) is same as ( ) ℎ ( ). ( ) under the ( ) becomes ( ) ℎ ( ) with additional terms, ′ (0) ( ) ḣ ( ) and ′ (0)̇( ) ℎ ( ) as in (34), resulting in unintended force of damping and springs.
Differential equation of the system is shown in (35). The unintended damping, ′ (0) ( ) and spring, ′ (0)̇( ) occur due to ( ). Hence, the effective damping ( ) and spring ( ) of the system become + ′ (0) ( ) and ( ) + ′ (0)̇( ), respectively. The haptic system is expressed as the linear time-varying system where ( ) and ( ) vary according to ( ) and . Table III shows the ( ) and ( ) of the system with and without ( ). The effect of negative damping and spring caused by ZOH and sampler is reduced by adding ( ).
where ( ) , ( ) and ′ (0) denote + ′ (0) ( ) , ( ) + ′ (0)̇( ) and 5 × 10 −3 , respectively. (35) The ISS of the system with ( ) is analyzed using the ISS-Lyapunov function candidate (10). The four conditions in Table I should be satisfied for the system (35) to be ISS. Fig.  9 shows the region that satisfies the four conditions according to and ̇. Condition I ( = + = (1 + 5 × 10 −3 ) + ) is always positive definite since and are positive. The value of (= − = + 5 × 10 −3 − ) corresponding to the condition II increases as increases. Hence, there is no constraint on the maximum displayed stiffness. The system without ( ) has the constraint on the value of as shown in (22). The system with ( ) always satisfies the condition II by reducing the discretization effect, which improves the stable region. The (1 + 2 )/(2 ) term in the condition IV is positive definite if the is positive definite. The condition IV is satisfied if the condition II is satisfied. The system is, therefore, ISS if the condition III is satisfied. Equation (36) shows the function of the condition III. ̇2 term in (36) is neglected since the coefficient of ̇2 is very small as compared to other terms. The boundary line of (36) is defined as ℎ( ) and expressed in (37).
Condition III: The stiffness change rate is constrained as ℎ( ) to satisfy the conditions. It is confirmed that the ISS region with ( ) is expanded compared to the ISS region without ( ) shown in Fig. 4. The maximum displayed stiffness of the system without ( ) is limited to as shown in Fig. 4, whereas there is no constraint on the maximum displayed stiffness of the system with ( ).
A controller 2 that satisfies the ISS conditions is proposed. ( ) is calculated as ( ) divided by ℎ ( ) in each time step when the contact occurs. ̇( ) is calculated using backward Euler as in (26). The controller 2 shown in Fig. 10 operates when contact occurs.
( ) is controlled to satisfy the conditions of ℎ( ) . ( ) denotes the controlled ( ) . △ ( ) in (38) is limited as ℎ� ( − 1)� if the conditions are not satisfied due to rapidly increasing of ( ) . Fig. 9 shows how the controller operates. ( ) is bound to ( ) as in (39) if ( ) is larger than ( ) . In other cases, the measured ( ) is used as ( ) . The controlled force ( ) = ( ) ℎ ( ) is rendered to the operator. The controller 2 does not have the constraint on the maximum renderable stiffness (28) that the controller 1 has. Hence, the stiff environment can be rendered.

V. EXPERIMENTS
The ISS region is computed for an example using a commercial haptic device Omega 7 as shown in Fig. 11. The feasibility of the proposed method is confirmed by analyzing the ISS region of the system using Omega 7 in the 1-DOF of the z-direction, which is the direction of gravity. The proposed method can be extended straightforwardly to multi-DOF haptic rendering by defining in x, y, and z-directions, respectively.
The physical properties of the device with its local controller are measured to compute the ISS region. Force input of chirp signal that varies from 0.1 Hz to 30 Hz with an amplitude of 1.5 N is applied to the z-direction. The position of the end effector is measured from the encoder. The velocity and acceleration are calculated using a zero-phase filter. The properties in (40), mainly used for the model of the haptic device with its local controller [28], are optimized using the Levenberg-Marquardt algorithm. The device mass , device damping , and Coulomb friction are 0.285 kg, 0.739 Ns/m, and 0.243 N, respectively. The relative percent difference error between the input force and the estimated force using the optimized properties is 13.91%.
The ISS region of the system without ( ) in Section IV-A is computed. The four conditions in Table I are expressed as  the function of device properties the haptic device such as  mass , damping , stiffness , change rate of stiffness ( ) to satisfy ( ). ̇, and sampling time . Function value of each condition according to and ̇ is calculated by substituting the device properties, , ̇, and into the function of the four conditions. Fig. 12 shows the region that satisfies each condition according to and ̇ when is 2 ms. Colors of the graphs represent the function value of conditions that satisfy the ISS. The function value of each condition should be larger than zero to be ISS. Hence the value less than zero is not shown in the graph. Ranges of and ̇ are set to 0 ~ 600 N/m and -2000 ~ 2000 N/(m•s), respectively. X and y-axis represent and ̇, respectively.
Condition I is always positive definite regardless of and ̇. The black-dashed line shown in Fig. 12 (b) represents the boundary conditions of (22). An asymptotic line of condition III in Fig. 12 (c) is the same as the boundary line of condition II. The boundary curve ( ) indicated by red dots in Fig.  12 (c) represent the boundary conditions of (23).The (1 + 2 )/(2 ) term in the condition IV is positive definite if the is positive definite. The condition IV is satisfied if the condition II is satisfied.
The ISS region that satisfies the all conditions is shown in Fig. 12 (c). Equation (41) represents the boundary curve ( ). The maximum displayed stiffness (25) is 453 N/m. The ( ) and are used for the controller 1 as in (27)- (30). The ISS region in Fig. 12 (c) shows that the ISS cannot be satisfied if the ̇ increases rapidly beyond the ( ) or exceeds the maximum . The ISS region of the system with ( ) in Section IV-B is computed. The four conditions in Table I should be satisfied for the system (35) to be ISS. Fig. 13 shows the region that satisfies each stability condition according to and ̇ when is 2 ms. Omega 7 can render a stiffness of up to 14.5 kN/m [29], hence the range is set taking this into account. Ranges of and ̇ are set to 0 ~ 15 kN/m and -40 ~ 40 kN/(m • s), respectively. Colors of the graphs represent the function values of condition that satisfy the ISS. The function value of each condition should be larger than zero to be ISS. Hence the value less than zero is not shown in the graph. Condition I, II, IV are always positive definite regardless of and ̇. The system is ISS if the condition III is satisfied. The system without ( ) has the constraint on the value of (= − ) corresponding to condition II as shown in Fig. 12 (b). The system with ( ) always satisfies the condition II by reducing the discretization effect, which improves the stable region. Equation (42) The ISS region of the system with ( ) that satisfies all conditions is shown in Fig. 13 (c). It shows that the condition is always satisfied when the stiffness change rate is negative. The stiffness increase rate is constrained as ℎ( ) to satisfy the condition. It is confirmed that the ISS region with ( ) is expanded compared to the ISS region without ( ) shown in Fig. 12 (c). The ℎ( ) is used for the controller 2 as in (38), (39).
Performance of the proposed controller 1 and 2 is experimentally compared with the virtual coupling and the FBA method, which are widely used passivity-based approaches. Fig. 11 shows the experiment setup using Omega 7. The experiments are conducted for 1-DOF in the z-direction. There is the virtual environment at a position less than 0 mm in the z-axis direction. The operator moves the device in the zdirection and contacts the environment. Experiments are performed in which the operator maintains a force of 10 N while looking at an indicator of the environment force. is set to 2 ms, and is 0.739 Ns/m as in (40). A spring of the virtual  coupling is designed to satisfy Colgate's passivity condition [13] as in (43). A damping of the virtual coupling is set as 0 Ns/m to use the largest . The maximum stiffness that satisfy (43) is 739 N/m. Hence is set as 735 N/m. The FBA method considering the memory effect [9] is used in the experiments.
The virtual environment of the experiment is composed of the Hunt-Crossly model as in (44) that is widely used for nonlinear contact model [30]. , , and are set as 5000, 100, and 1.5, respectively. The of the Hunt-Crossly model increases as the insertion depth increases. Soft properties are rendered at the beginning of the contact, and stiff properties are rendered as the insertion depth increases.
Fig. 14 shows the experiment results using the virtual coupling. Smaller force is rendered compared to the environment force. The operator always feels stiffness softer than correct one when using the virtual coupling, since (= 735 / ) is connected serially between the virtual tool and the haptic device. As an example, the original at 2.5 s and 6.4 s is 300 N/m and 600 N/m, respectively. The using the virtual coupling at this time is 213.3 N/m (=1/(1/300 + 1/ )) and 330.3 N/m (=1/(1/600 + 1/ )), respectively. The original increases as the insertion depth increase in the Hunt-crossly model. The error of the rendered force increases as the original increases. Energy stored in the system using (7) after the experiment is 9.59 × 10 −5 N • m. Hence the generated energy is less than zero since the virtual coupling is designed to satisfy the passivity. Table IV shows the error between the rendered force and the environment force. A root-mean-square error (RMSE) is 3.26 N, and a maximum error is 4.88 N. Fig. 15 shows the experiment result using the FBA method considering the memory effect [9]. The computed environment force exceeds the passivity condition as the original increases. Force bounding does not occur at the initial contact, and soft properties are transmitted. Force bounding occurs when the insertion depth increases, and stiffer properties are transmitted. The rendered force is bounded at 4.608 s ~ 17.096 s. The magnitude of the bounded force is less than 5.832 N. RMSE is 3.02 N, and a maximum error is 4.79 N. The environment force is transmitted faithfully if it is less than the computed maximum force. Energy stored in the system using (7) after the experiment is 1.09 × 10 −4 N • m. Hence the generated energy is less than zero since the FBA method limits the maximum force to make the system passive. Fig. 16 shows the experiment result using the proposed controller 1 .
is controlled to gradually increase to satisfy the condition of ( )  95 N and 3.03 N, respectively. RMSE and maximum errors using the controller 1 are reduced by at least 35 % compared to using the virtual coupling and FBA method. Energy stored in the system using (7) after the experiment is −6.72 × 10 −3 N • m. It means that the energy generated from the haptic system is a positive finite value. More energy is transferred to the operator compared to the passivity-based methods. The increased stiffness corresponding to the additionally transferred energy is rendered compared to the passivity-based method. Performance of haptic rendering is, therefore, improved. Soft tissue models used in medical simulations generally have a stiffness up to 400 N/m [31], [32]. The proposed controller 1 is sufficient to render the haptic feedback of soft tissue since α is 453 N/m when T is 2 ms. There is, however, a limit to provide high fidelity haptic feedback when rendering the stiff environment with stiffness of more than 453 N/m. Fig. 17 shows the results of the proposed controller 2 . ̇ is adjusted when increase rapidly in 1.064 s ~ 2.464 s that is indicated by the dotted ellipse on the stiffness graph.
is controlled to gradually increase to satisfy the condition of ℎ( ). An error of the rendered force occurs in this region. The environment force is faithfully rendered to the operator after reaches the original . The RMSE and maximum errors are 0.43 N and 2.40 N, respectively. The error of the rendered force is summarized in Table IV. The RMSE and maximum error are reduced by at least 85.76 % and 49.90 %, respectively, compared to using virtual coupling and FBA method. The RMSE and maximum error are reduced by 77.95 % and 20.79 %, respectively, compared to using the controller 1 . Stored Energy in the system using (7) after the experiment is −8.20 × 10 −3 N • m. The energy generated in the system is increased compared to the system with controller 1 and passivity-based methods. This means that more energy is transferred to the operator, and the increased stiffness corresponding to the additionally transferred energy is rendered. Hence, the performance of haptic rendering is improved.

VI. CONCLUSION
This paper presents the control method to adjust the change rate of the environment stiffness based on the ISS conditions without the accurate model of the environment. The haptic systems becomes dissipative if they satisfy the ISS conditions. The generated surplus energy of the dissipative system may have a positive finite value. The proposed ISS-based control keeps the system dissipative, and the rendered force is controlled so that the generated energy has a positive finite value. More surplus energy of the dissipative system, compared to the passive system, can be transmitted to the operator. Hence increased stiffness of the environment corresponding to the surplus energy can be rendered to the operator, increasing the fidelity of the haptic system. The environment is modeled as the linear time-varying spring , and effective mass, spring and damper of the system are analyzed in continuous-time domain. The ISS conditions are analyzed, and the controller 1 is proposed to satisfy the conditions. The controller 2 is designed to enhance the performance using the compensator and ISS conditions. The compensator reduces negative damping and spring caused by the discretization effect, resulting in the wider ISS region. The is controlled to increase gradually to satisfy the ISS conditions if increases rapidly and deviates from the conditions. The environment force is faithfully rendered to the operator when the controlled reaches the original . The proposed controller uses only the current magnitudes of the device position and the environment force. The controller does not require any exact model of the environment. Hence, the proposed method guarantees the system stability with unknown nonlinear environment. Experiment results show that the energy generated from the haptic system using the   proposed method is a positive finite whereas the generated energy using the passivity-based approaches is less than zero. This means that more energy is transferred to the operator compared to the passivity-based method. The proposed controller reduces the RMSE and maximum error by at least 85.76 % and 49.90 %, respectively, compared to using the virtual coupling or FBA methods.