A Joint-space Position Control-based Approach to Haptic Rendering of Stiff Objects using Gain Scheduling

Haptic rendering often deals with interactions between stiff objects. A traditional way of force computing models the interaction using a spring-damper system, which suffers from stability issues when the desired stiffness is high. Instead of computing a force, this paper continues to explore shifting the focus to rendering an interaction with no penetration, which can be accomplished by using a position controller in the joint space using the encoders as feedback directly. In order to make this approach easily adaptable to any device, an alternative way to model the dynamics of the device is also presented, which is to linearize a detailed simulation model. As a family of linearized models is used to approximate the full dynamic model of the system, it is important to have a smooth transition between multiple sets of controller gains generated based on these models. Gain scheduling is introduced to improve the performance in certain areas and a comparison among three controllers is conducted in a simulation setup.


Introduction
Haptic devices are mechanisms that try to reconstruct the often missing haptic feedback in virtual reality (VR) interactions and reflect it to the user to make the whole interaction more realistic. Haptic feedback, though sounded unfamiliar, is one of the five main senses we rely on daily. Haptic technologies are everywhere from our mobile phones and laptops to highly specialized devices like dental training simulators [16,32]. Yang  Analogous to graphic/visual rendering, haptic rendering is an additional human-computer interface for the sense of touch so that the user can perceive virtual objects through haptic channels like skin, hands, and other receptors. In this paper, we focus more on the haptic rendering of interactions between a tool and a stiff object using desktop impedancetype force feedback devices. For such applications, it is commonly accepted that the refresh rate of the haptic rendering loop should be 1000 Hz and within each iteration there are 5 sequential steps [26]: 1. sense the configuration of the device in joint space; 2. map that into workspace using forward kinematics; 3. determine whether there is a collision and where it is; 4. compute the correct interaction force/torque (F/T) to render in workspace; 5. map that into joint space F/T as reference signals for low level control; Although interaction with soft objects can be one of the new frontiers of haptic research, stiff interaction is still a large part of the targeted rendering scenarios and quite challenging to render stably. Keeping a high refresh rate for the haptic rendering loop is one of the ways towards more stable rendering, which can be achieved by making computation more efficient. Different approaches are developed for different modeling methods. Zilles and Salisbury [36] develop the concept of god-object to determine where the interaction is. A similar concept is called "virtual proxy" [24]. This distinguishes where the tool should be, which is represented by the god-object, from where the tool actually is. The position of the god-object is solved by active constraints and it is suitable for the systems that use polyhedral representations. Ortega et al. [21] later extend the concept into a 6-DOF god-object for more complex haptic applications. As the computation gets more complex, they use asynchronous processes to compute the god-object motion and the rendering force, respectively, to meet the time requirement for each process. Otaduy and Lin [22] give another take at rendering complex objects with high polygon count. They linearize the contact model in complex contact scenarios to maximize the haptic refresh rate. Yu et al. [35] propose an extended version of the configuration-based optimization approach [31] to target sharp geometric features on objects. The approach uses a coarse-to-fine sphere tree to represent an object for faster collision detection and contact constraint formulation. It is designed for dental simulators and claimed better for rendering narrow space constraints and fine features in general.
Normally, for the whole haptic rendering loop to work, models of the device kinematics and dynamics are evaluated in real-time, which is another component to slow down the haptic rendering loop. Shortening this computation time is beneficial for stability and the methods are not limited within haptic applications. Park and Lee [23] introduce a simplified model to estimate the full model. To compensate the nonlinear dynamics and modeling errors, adaptive and robust control method are used to dynamically estimate and compensate these terms. Cheng et al. [6] experiment on a 2-DOF redundantly actuated planar parallel mechanism (PM) to compare the performance of 4 different control schemes. The device is also chosen for its relatively simple structure and the model of its dynamics can be obtained in closed-form. The comparison shows that the proposed control algorithm has smaller error in position/velocity tracking by moving a part of the computation offline to have a shorter execution time. Corbel et al. [11] propose a simplified dynamic model that is motivated by the high speed high acceleration application. The simplification is to ignore the inertia of some linkages and split their mass so that the calculation associated with these linkages can be avoided. Abdellatif et al. [2] propose a method to find a reduced set of dynamic parameters and use that to reduce computation. The optimal set of parameters is reduced to 10 rigid-body and 14 friction-related parameters for a 6-DOF PM. In a similar fashion, Diaz-Rodriguez et al. [12] propose an model identification method and demonstrate it on a 3-DOF PM. The dynamic model is first established with base parameters and then the parameter with the largest standard deviation is iteratively removed until the inertial matrix is positive definite. The results show that the controller based on the reduced model has less error and computational load, which enables the real-time control of the device.
Another instability factor is extra energy. Though the virtual environment is passive, extra energy can still be generated from interacting due to the sampling nature of the haptic device. Gillespie and Cutkosky [15] term the phenomenon "energy leak" and they propose passivitybased approaches to solve this. Hannaford and Ryu [17] propose time-domain passivity control with an observercontroller structure. An adaptive damping is introduced to precisely dissipate the extra energy. Ryu and Yoon [25] extend this into a memory-based passivity approach. In their case, a dedicated hardware FPGA is fast enough to log all the data needed for position versus force trajectory plot and the controller uses that to guarantee that the net energy coming from the passive virtual environment is less than 0. Note that most passivity-based controllers are considered a bit conservative as passivity is a sufficient condition for stability. Colgate et al. [10] propose "virtual coupling" as a filter to achieve stable behavior. This approach is quite popular though some studies also point out that it might lead to some "force artifacts" in force rendering [35]. Adams and Hannaford [3] later extend this approach to work for both impedance or admittance models of haptic interaction but the approach cannot handle a dynamically changing environment.
In summary, there are different ways to speed up the haptic rendering loop and cut down computations. For step 4 of force computing, stiff objects are mostly modeled as a spring-damper system and this can lead to unstable behavior if not dealt with carefully. To avoid the shortcomings of using a spring-damper system, we propose a positioncontrol based approach as a solution to achieve a stable interaction with stiff objects [33]. A new user scenario is introduced in this paper. In this scenario, we move the tool along a straight path on a stiff flat wall while we are applying a force perpendicular to the wall. This "pressing while sliding" scenario resembles the task of polishing in dental procedures. The following research questions (RQs) will be further explored: RQ 1: What are the benefits of using a family of piecewise linear models to develop a gain scheduling controller with smooth transitions? RQ 2: What performance can be obtained for the user scenario "pressing while sliding"?
The rest of this paper is structured as follows: the concept of the approach is revisited in Section 2. In Section 3, the control design process is laid out with the use of piecewise linear model and gain scheduling. Section 4 shows the simulation methodology on how to setup tests to evaluate the system. Results of the simulation are shown in Section 5, followed by a discussion in Section 6. Finally, conclusions and future work are presented in Section 7.

An Approach using Position Controller
In this section, we first introduce the test case device for the paper and then, the proposed approach is recapped as it originates from previous work [33,34].

Description of Test Case
The test case for this paper is TAU, a desktop 6-DOF PM haptic device developed at KTH as a dental training simulator l [4,5,29]. The design of the device is optimized in the previous work [4] to have a large enough singularityfree workspace (a 50×50×50 mm 3 cube) while maintaining a uniform performance within the whole workspace and a relatively small footprint. The whole system can be fit into a 250 × 250 × 300 mm 3 box. The design parameters can be found in Table 1 and other details can be seen in Fig. 1.
What makes TAU different from other well-studied PMs such as Delta [7] or Stewart [20] is that it has asymmetry in its kinematic structure. This is one of the reasons why we choose TAU as the test device and we will further explain this asymmetry and the challenges it opposes in Section 3.2

Recap of the Approach
For an impedance-type haptic device to render a stiff object, the traditional way is to model the virtual environment (VE) as a linear spring with high stiffness. With a high enough stiffness, the rendered spring behaves like a stiff object in real life because of the imperceptible deformation. When there is a collision and a penetration into the object, a penalty force following the Hooke's Law will push the device back out of the object. This is one of the common ways to do step 4 in the haptic rendering loop since the idea is introduced in [9], illustrated in Fig. 2a. It computes reference force signals in the Cartesian space and further converts them to joint torques signals through inverse dynamics. Here, we use the naming convention from [31],  Table 1 where we call the actual pose of the tool that is tracked by the encoders the haptic tool and the pose that satisfies no penetration the graphic tool. The haptic tool and the graphic tool are colored green and blue, respectively, in Fig. 2.
Our approach introduces a new way to compute the rendering force [34]. A stiff object behaves similarly to a position controller in the sense that it should keep the tool on the surface of the object with no penetration. In this way, given the collision position and the current position of the tool, a position controller can compute a penalty force to eliminate the difference between the two, bringing the tool back to the surface of the object. It is a closed-loop position control in joint space.
Compared to the traditional force rendering method, the newly proposed controller has the same input/output (I/O) signals as other controllers; hence the traditional one can be easily replaced by the new controller in haptic devices, which makes it easy to use the new controller in an old system. Because now the control goal is in joint space instead of the Cartesian space, some nonlinear and complex computation is avoided, which is good for maintaining the high haptic rendering loop rate. proposed. The blue ball symbolize the graphic tool while the green ball the haptic tool

Control Design
As recapped in Section 2.2, haptic rendering problem can be defined as generating the correct interaction force to create the illusion of directly interacting with virtual objects [19]. To solve the haptic rendering problem as a position control problem, the problem is first reformulated in Section 3.1. Then, for the remaining part of Section 3, we revisit the plant modeling and the discrete time control design. A linear-quadratic-regulator-based (LQR-based) full-state feedback structure is chosen because it can be optimized to prevent large control input. Given this choice, a linear time-invariant (LTI) state-space control plant is sufficient for the control design. In this paper, only key results are highlighted. A detailed version of this can be found in [34].

Equivalent Problem
With the discussion from the previous section, the haptic rendering of stiff interaction can be realized this way: given the current pose of tool handle in the Cartesian space, design a position controller that computes a motor torque reference so that the tool has minimal penetration into the stiff object when it is disturbed by the user. In this way, the user will feel like interacting with a real stiff object. Without losing generality or consistency with previous work [5], the pose in Cartesian space can be represented by a 6-by-1 vector where (p x , p y , p z ) denotes the position coordinates and (α, β, γ ) the orientation angles of the tool along x, y, and z axis, respectively. The pose in Cartesian space corresponds to joint angles of the active joints and the mapping from Cartesian space to joint space representation is the inverse kinematics. In the joint space, the pose of the device is also defined by a 6-by-1 vector where q i denotes the angular position of joint i. These are the active joints that govern the motion of the device while the rest of the joints are passive. Having passive joints is also one of the reasons why the forward kinematics of a PM is hard to derive.
Accordingly, the force to be rendered here is a 6-DOF force acting on the tool in Cartesian space, defined by a 6-by-1 vector and the joint torque to achieve that is defined as where τ i denotes the joint torque on joint i.
These are general definitions in kinematics and dynamics, and they can be used on both the graphic tool and the haptic tool. With these definitions, the problem becomes: given the current pose of the graphic tool X graphic determined by collision detection as reference signal and the current pose of the haptic tool X haptic estimated based on the q haptic measured by joint encoders, design a position controller that determines the joint torques τ to minimized the difference between X haptic and X graphic .

Plant Modeling
A model of the device dynamics is needed for this phase and traditionally, the analytical form of it is obtained through complex mathematical derivation. This is nontrivial, especially in our case, as the device is a PM with asymmetric chains. The asymmetry in the mechanism poses difficulties in modeling it in the sense that each chain has to be modeled differently. This not only is tedious but also further hinders the simplification of the equations which relies heavily on the symmetry within the mechanism. To tackle this problem and to have an approach that can be generalized to other devices, an alternative way to model the mechanism is needed. Note that unlike D-H parameterization for serial mechanisms forward kinematics [27], there are few mature methodologies to derive the dynamics for any PMs in analytical form. It mostly works on a case-by-case basis with the help of special structures or symmetry within each mechanism. Often times, the model is then simplified to some level so it can work with the chosen control design methodologies.
To speed up this phase and to develop a more general approach, we build the model of the device using Simscape, a module within the MATLAB/Simulink environment, and then obtain the targeted LTI state-space model by applying the MATLAB linearization tool. Simscape is chosen for its support for PM modeling and its ease of linearization. The approach has good generalization in the sense that it will also work on other devices with different structures.
The modeling process in Simscape is intuitive and close to how one assembles a mechanism in reality. As TAU (see Fig. 1) is a PM, the modeling starts from the base, the Icolumn, and continues to the end of each kinematic chains. Finally, the chains are connected by the moving platform, which forms kinematic loops. The tool handle is mounted onto the platform and the tip of the tool is aligned with the center of the moving platform. The I/O of the model is defined to match the physical machine and to have a good interface for the next step. The input u to the device is defined as τ and the output y = [q T , q T ] T . The next step is linearization. With the I/O defined in the previous step, the linearized model is obtained by simply applying the Linear Analysis Tool (LAT) from MATLAB. An operating point (OP) is chosen q OP = 0 6×1 ,q OP = 0 6×1 and the system is linearized around that point. Here, we use 0 m×n as a shorthand notation for an m × n matrix/vector that is full of zeros. The τ OP is obtained by running a static simulation of the mechanism at the OP. For any given OP, the tool should be static at the OP under the corresponding τ OP . The result of the LAT is an LTI state-space model where A, B, and C matrices are the constant coefficients in this form. Here, the state vector is δx instead of x because it is linearized at an OP and the result should be with respect to that OP where x is the state vector and x OP is the value the state vector should be at the OP. This is also true for δu, δy, and a few other variables below. Note that currently, the state vector is a 12 × 1 vector decided by LAT. It involves states related to things like passive spherical joints in the mechanism, which can be unreliable for state feedback and impossible to measure with the current hardware. Thus, a linear transformation is needed for the model to switch to another set of desired states that are only related to the active joints, monitored by the encoders. Because C is invertible, and I n×n is the n × n identity matrix.
Then, the state vector needs to be augmented with some integration of error states to eliminate the steady-state error.
These added integration of error states are and r is the reference signal in joint space. The augmented model becomes with

Discrete Time Control Design
The control plant in (2) is a continuous time model. It will be transformed into discrete time by the Zero-Order Hold (ZOH) method and then applied for discrete time control design.
where T s is the sampling time and it is set to be 1 ms (1000 Hz refresh rate) [1,14]. The model from (2) is transformed into The e here is the natural exponential constant, also known as Euler's number. The result is a discrete time LTI state-space model control plant, denoted by a subscript of auDT. So we formulate the control problem as an optimization problem to find the optimal gain K auDT such that the state-feedback control law minimizes the quadratic cost function with a LQR-based formulation The Q and R matrices are tuned to be The actual control command sent to the motors is

Gain Scheduling with Multiple OPs
So far, we revisit the approach introduced in [34] and the results with a single OP (SOP) in that paper are quite satisfying. However, we also noticed that the performance drops when the tool moves away from the SOP, which is the origin of the workspace. After a bit of tuning, we realized our performance is limited by the hardware capacity and the accuracy of the linearized model. Gain scheduling is a common nonlinear control technique [18] to design subcontrollers based on linearized models of the nonlinear plant at different OPs. This technique can be dated back to the 1960s and it also has some recent development. Over the years, it is sometimes combined with the study of linear parameter varying (LPV) systems.
One of the studies related to haptics is done by Walsh and Forbes [30]. In this paper, the authors introduce a very strictly passive gain-scheduled controller considering the passive properties of the controller. However, as most passivity-based control, it may behave too conservatively when rendering a stiff wall. In our study, as a first step to introduce gain scheduling in our control, we follow the classic gain scheduling method [18]. Therefore, a family of linearized models at different OPs is introduced to be the piecewise linear model to describe the full dynamics over the whole workspace. To achieve good performance even on the edge of the workspace, 27 ( where m, n, l are the OP indices in x, y, z directions, respectively. The number 25 is the half of the side length of the 50 × 50 × 50 mm 3 workspace and the unit is mm. Accordingly, the plant models and the controllers associated with these OPs are generated with the same procedures and Q, R matrices. The resulting family of linear controllers serves as a single controller whose parameters are changed based on the scheduling variables. The scheduling variables are all the variables in X haptic , which is estimated based on the joint encoder readings q haptic and forward kinematics. A naive switching method among these local linear controllers is based on the closest OP to the current X haptic . The switching mechanism is illustrated in Fig. 3a as a two-dimensional (2D) case. The backgound with different shades of gray indicates different sets of control gains based on the different OPs (red dots). Note that the size of the colored dots is not relevant and is only for illustration purposes. To have some consistency with Fig. 2, the light blue ball still represents where the graphic tool is. Note that in this view, the graphic and the haptic tool coincide because this view is looking at the wall and there are no constraints on the xy-plane. The area circled by black dashed line is a zoom-in illustration of what could be the estimated trajectory when a user move the tool right along the switching boundary. The zig-zag gray solid line is a possible trajectory for the tool on a millimeter scale due to user hand tremor, sensor noise, and numeric error, prior to accounting the effect of control. This trajectory will lead to rapid switching which is not an ideal behavior for gain scheduling controllers in general. The switching is more severe when the tool is working near the intersecting point of two switching boundaries because potentially it can switch among four sets of control gains. The discontinuity in the control signals based on this switching mechanism leads to oscillations in motion, which can be sensed by the user. To eliminate switching boundaries, a smooth transition can be achieved by using e.g. linear interpolation based on the distances to the closest set of OPs. Although we are using a multivariate linear interpolation, the mathematics is still based on the 1D case that where v is the value associated with x, and q v is the value to be interpolated based on q x. b x 1 , b v 1 and b x 2 , b v 2 are the value pairs of the break points, indicated by the superscript b. q v can be calculated as With higher dimensions, this 1D case is performed once in each direction. Because this is a 1D case, b x is only a scalar, whereas for nD case, b x is a nD vector The associated value b v is still a scalar and the value pair for a break point becomes b x, b v . We use f (. . .) to denote the associated value v with the corresponding set of x in different dimensions. A pseudo code version of the nD linear interpolation algorithm is shown in Algorithm 1.
Note that the order of the directions in which the interpolation is performed does not matter and the number of break points grows exponentially with the dimension. In the 2D case shown in Fig. 3b, the set contains 4 OPs and in 3D it will contain 8 (2 3 ) OPs. As also stated in the algorithm input data, for the j th dimension, i j can only take the value of either 1 or 2. This means that for all the OPs, the value for their j th coordinate can only be either b x 1 j or b x 2 j . Given there are n dimensions and for each dimension the value of the coordinate only has 2 choices, we can conclude that there are 2 n different OPs in the set of closest OPs for this Algorithm 1: nD linear interpolation.
Input: 2 n OPs with the notation Input: a QP with the notation  Because the value is interpolated, it is depicted in Fig. 3b with a gradual changing color of gray without clear switching boundaries as opposed to Fig. 3a. The interpolated values include the controller gain matrices K auDT as well as offsets at the OP such as q OP and u OP .
The new controller is named multi-OP position controller (MOPPC) to highlight the usage of multiple operating points (MOPs). As a comparison, we include the proposed controller in [34] and rename it to SOPPC to emphasize that it is generated from a control plant using a single OP. The baseline controller in that paper is also included as a comparison and renamed as open-loop impedance-type force controller (OLIFC).

Simulation Methodology
Simulations are conducted to investigate the performance of different controllers mentioned above. One of the tests frequently used for evaluating force rendering capability is to do a 1D pressing against the wall, as shown in [13,17], where other moving directions are constrained and therefore ignored. The device only moves in the direction perpendicular to the wall and the penetration in this direction is recorded and reported. Although this test demonstrates the performance of the device at one point in space in a certain direction with a near quasistatic setup, it does not represent the performance over the whole workspace or consider the effect of motions in nonpressing directions. If put in an application sense, this test resembles more to the working conditions during tasks like drilling but gives little insight into tasks like milling or surface polishing, where a major lateral motion is on-going while a pressing force perpendicular to the motion is also present.
To mitigate this, we test the performance of TAU under "pressing while sliding" conditions in this paper. This is also found in papers like [31]. It provides another layer of verification on the performance of the device in addition to "pressing while stationary". The object we are rendering is still a unilateral stiff wall. "Unilateral" means that the wall only constrains the motion of the tool in one direction [13]. Specifically, it prevents the tool from entering the wall. However, the tool should be free to leave the wall. A set of fixed trajectories is chosen in order to investigate how the performance varies along each trajectory. The trajectory is fed to the controller in real-time, mimicking the result from the last haptic rendering step, collision detection. Two different disturbance scenarios are also designed to represent typical disturbance in working conditions. For each scenario, 20 Monte-Carlo simulations are conducted on the full dynamic model of the device with randomly varying disturbance level and white noise. These specifications will be further explained in the following subsections.

Query Trajectory
Eight straight line segments are chosen as a set of representative paths to cover the rendered wall as shown in Fig. 4. Then, these 8 paths are turned into 16 query trajectories (QTs) given that one path can be traversed in two different directions. The performance along all QTs should give a more comprehensive view of the system instead of the performance at one point in the workspace. Each trajectory is parameterized with two waypoints, the start and stop pose of the tool, and the traveling time. The traveling time, regardless of the length of the path, is set to 2 s. The trajectory is also guaranteed to be smooth and second order differentiable. The 16 QTs are labelled in Fig. 4 with numbers close to its starting point.

Disturbance Scenario
Two scenarios with different disturbance models are introduced to test the performance under different kinds of workload.

Scenario 1: Force Model
Scenario 1 resembles the step test for general control systems. A pure force signal is applied onto the tool, Fig. 4 Illustration of all 16 QTs on the wall, labelled with number. Note that the two QTs on the same path are separated apart only for clarification. The red dots are OPs on the wall emulating the pressing force from the user. The force signal is a step signal overlaid with a white noise as shown below The step level (F c ) is randomly generated for each simulation uniformly distributed within the range of 8 -15 N. The step signal holds the whole duration of the simulation. The white noise follows a normal distribution parameterized by a mean of μ and a standard deviation of σ . The notation of N (μ, σ ) is adopted to represent that. The amplitude of the white noise is roughly 10% of the step level. As stated, the same simulation will be run along each QT with a randomization in the step level and added white noise as a Monte-Carlo simulation in order to demonstrate the robustness of the performance.

Scenario 2: Hand Model
The disturbance of Scenario 2 is applied through a hand model. As suggested by [28], the behavior of the hand can be modeled using a spring of 40 N/m and a damper of 3.6 Ns/m, connected in parallel. Unlike Scenario 1, where a pure force signal is generated, a hand motion is designed and applied to the hand. The relative motion and the springdamper system modeling the hand generate a disturbance on the tool, imitating how the user applies a disturbance.
Contrary to the force model, where the disturbance has a sharp edge to jump from 0 to F c , the hand model is aimed to apply a disturbance with a smooth transition. Other than that, the idea is still the same: the disturbance should stay constant for a while to confirm if there is a steady state error. This translates into having a hand motion that moves towards the wall and then stay stationary to maintain the level of pressing force, assuming the tool also stays stationary on the surface of the wall. The motion is designed to follow a constant acceleration, constant speed (zero acceleration), constant deceleration profile so that the acceleration follows The time period of acceleration/deceleration is chosen so that the generated disturbance has a smooth transition and little overshoot. With these constraints, a c is calculated based on F c and how much the hand needs to move. By designing the hand motion X hand , the disturbance is generated at run-time as a function of X hand and X haptic .

Post-processing and Performance Metrics
Stiffness related metrics are commonly used in evaluating the performance of haptic devices [8] since it considers the level of the disturbance as well as the penetration. However, it is also a calculated value based on a simple yet important assumption that there is indeed a penetration. In the calculation, the penetration is always in the denominator. If the penetration is close to zero and positive, the calculated stiffness will be huge. It does not necessarily mean that the stiffness rendered by the device can reach that level constantly. If the penetration is exactly zero, which is the ideal case for rendering a stiff wall, the calculated stiffness is infinity. This is a bit problematic if the next step is to take the average of the values. If the penetration is negative, which by definition means that there is still a distance between the tool and the wall, the stiffness calculation loses its physical meaning because there is no collision. Due to factors like numeric error, digitization and noise in sensors, oscillations, and so on, these three cases can all occur within one simulation. For this study specifically, a stiffness higher than 1.6 × 10 5 N/m is considered to be cut off, as the disturbance can reach no larger than 16 N and a penetration less than 0.1 mm is considered unidentifiable by the encoders. If this value is reached, we consider the device is rendering infinite stiffness and these data points will be processed separately. This is later referred to as the cut-off criteria. Note that the criteria is met mostly because there is too little penetration. So, when the criteria is met, it is also considered to have no penetration, which corresponds to the infinite stiffness.
To maintain consistency with the previous findings in [34], minimal stiffness K min along the QT is chosen for evaluating Scenario 1 and average stiffness over time K avg for Scenario 2 for each simulation.
K min is calculated by , ∀j ∈ {1, 2, ..., N} and it is chosen to report the worst-case scenario for Scenario 1. D z (t) is the penetration into the wall and N is the number of time steps for the simulation. Corresponding to the sharp step change in the disturbance, there is a relatively large initial penetration, as shown in Fig. 5. As the disturbance holds its constant level after the step, the penetration should be less noticeable. Because the initial penetration exists only for a short period of time compared to the whole simulation time, the average stiffness over time will filter out this low point. Because the disturbance in Scenario 2 has a smooth transition, as depicted in Fig. 6, the relatively large initial penetration found in Fig. 5 is not presented in Scenario 2. The slow change in disturbance value also gives controllers enough time to response, so the performance during the whole simulation time is quite consistent. Large penetrations are more likely as a result of the noise. Hence, K avg can better describe the performance along QTs.
As stated previously, no penetration is actually the aim of the controller and we consider stiffness that meets the cut-off criteria to have reached infinite stiffness. Infinite stiffness ratio (r IS ) is the ratio of the time period that the calculated stiffness meets the cut-off criteria over the total simulation time. Because the sampling frequency in discrete time is constant, the ratio is equivalent to Again, N is the total number of time steps for the whole simulation and N cut-off is the number of time steps that meet the cut-off criteria. r IS should be used as an auxiliary metrics when the results from the stiffness metrics are unclear. When r IS = 1, it means that the calculated stiffness can be considered infinite throughout the whole simulation and it should be considered no penetration during the simulation. When r IS = 0, it means that there exists a noticeable penetration throughout the whole simulation given the definition of the cut-off criteria. When r IS is between 0 and 1, it just means that the controller manages to render infinite stiffness sometimes.
r IS is accompanied with another auxiliary metrics called maximal off-wall distance (D mo ). When r IS is larger than 0, D mo is examined to make sure that the tool never gets too far off the wall, because that is an unwanted behavior.
These metrics will be used in the coming result and discussion sections.

Simulation Results
To demonstrate the value level as well as the distribution, (half) violin plots are chosen (see Figs. 7 and 8) to report the results corresponding to different scenarios, which can be considered as an enhanced box plot. Each violin plot includes a mini box plot with the circle as median, a rotated kernel density plot on the side, and a horizontal bar indicating the mean.

Scenario 1 with Force Disturbance
For Scenario 1, K min is chosen as a worst-case scenario estimation. As shown in Fig. 7, the results from three different controllers are shown in two different subplots. The value level of K min for OLIFC is quite consistent, at a value level lower than 3000 N/m. Thus, it is plotted in a separate subplot to show some details and variations in the distribution of the data. The results for both MOPPC and SOPPC are plotted in the top subplot since they are in a similar range. That being said, the results spread roughly between 3000 and 30000 according to the figure.
The results for MOPPC, indicated with the color green, can be further separated by the green dashed line at 20000. Similar separation can also be found for the results for SOPPC by the blue dashed line at 9000. The separation suggests that within the results for each controller the results for the first 8 QTs are on a different level from that for the last 8 QTs. The separation will be discussed in details in the discussion section. The red dashed line at 3000 gives a sense of where the results for OLIFC might be as 3000 is loosely the upper bound of the values. The results for QT 6 are omitted for all controllers because the results for MOPPC and SOPPC meet the cut-off criteria.

Scenario 2 with Hand Disturbance
K avg is the main metrics for evaluating Scenario 2. Note that based on how the post-processing is done, the data points that meet the cut-off criteria are not used for calculating K avg . In Fig. 8, the results for OLIFC still show a consistent stiffness at the level of 6000 N/m, which is the stiffness set for the wall for this controller. The results for OLIFC are still plotted separately from those for SOPPC and MOPPC due to the difference in the order of magnitude. A red dashed line of 6200 is added to the top plot just to illustrate where the numbers for OLIFC will be if the results for all three controllers are plotted together. For SOPPC and MOPPC, the results for the first 8 QTs are omitted because most of them meet the cut-off criteria and should be considered infinite stiffness. These results are also accompanied with an r IS , the infinite stiffness ratio, of 1 or close to 1. With this kind of r IS , D mo (maximal off-wall distance) ensures that the tool actually stays on the surface of the wall instead of breaking the collision and going off the wall back to free space. The D mo for all simulations along all 16 QTs is about 0.2 mm for MOPPC and 0.7 mm for SOPPC, which is an acceptable distance to be considered as on the surface of the wall. Meanwhile, the maximal penetration is around 0.3 mm and 0.9 mm for MOPPC and SOPPC, respectively.

Discussion
From the results for both scenarios, the performance is quite consistent for OLIFC. With a slow changing disturbance, it seems to reach its predefined stiffness of 6000 N/m. However, as shown in Scenario 1, with a sudden change in the disturbance, the worst-case stiffness is only about half of the targeted value.
The choice of the value 6000 N/m for OLIFC is not arbitrary in this study. An initial guess is based on guidelines in [1]. With some initial trial-and-errors and some margins, 6000 N/m is set for a balance between less oscillation and higher stiffness. Judging from the results for the other two controllers, higher stiffness cannot be reached due to instability instead of saturation. This might explain why the performance of the OLIFC is so uniform. The targeted 6000 N/m is well within the saturation limits at any point on the wall so that it can be reached everywhere easily.
Because OLIFC treats the wall as a linear spring with a constant stiffness, there is an unavoidable steady-state error in position, hence a penetration. With the maximal disturbance of 15 N and the predefined stiffness of the wall, a steady penetration of 2.5 mm can be expected. Similar results can be found from the simulations. Because the disturbance level uniformly varies from 8 to 15 N as a part of the Monte-Carlo simulation and the penetration varies with the disturbance level, the average number is slightly better, at around 1.9 mm for Scenario 1 and 1.5 for Scenario 2. The maximal penetration is reached in Scenario 1 with a value of 5.3 mm.
For MOPPC and SOPPC, many things are in common between the two when compared to OLIFC. They are both ambitious in trying to render the wall with the highest stiffness achievable for the physical device. The results are quite similar as shown in Figs. 7 and 8 with a magnitude of 10 4 in N/m. In comparison with each other, MOPPC seems to have a better performance across all scenarios and QTs, except for QT 10 in Scenario 2. The improvement can be interpreted as the benefit of having a more accurate model. With gain scheduling, the controller is using a set of more appropriate gain that is specifically designed for the current neighborhood. The exception on QT 10 can be caused by the cut-off criteria. The r IS for SOPPC in this scenario ranges from 0.56 to 0.74 while that for MOPPC ranges from 0.80 to 0.88. This means that more data points with high values are cut off and taken out in the calculation for Fig. 8. If they are added back, MOPPC can outperform SOPPC along QT 10.
In the result section, we raise the notion that the performance for the first 8 QTs seems better than that of the last 8 QTs for both MOPPC and SOPPC. This is quite clear with the dashed lines in Fig. 7. The results for first 8 QTs are omitted in Fig. 8 because they are much higher than that of the last 8 QTs. This difference indicates the QTs can be divided into 2 groups. Our hypothesis is that this division can be based on the position of starting point of the QTs, more specifically, whether the starting point is on the upper or lower half of the wall. Based on this criterion, the QTs can be regrouped into 3 groups: 1. QT 1 -4 and 6 -8 (upper half); 2. QT 9 -12 and 14 -16 (lower half); 3. QT 5 and 13 (boundary cases) e starting point is on the upper or The boundary of upper and lower half of the wall is the x axis (p y = 0), where QT 5 and 13 are on. This division is close to the division discovered in the figures. The reason behind this division can be that a good initialization is important for the performance along the whole QT. In the previous study [34], a tendency of the upper half of the wall showing a better performance was identified. Combining these findings, we can assume that the performance along a QT tends to be similar to the performance at the starting point. It maintains the performance that way even when the tool later moves to a point that is queried to show a different performance.
According to this hypothesis, most of the first 8 QTs are from the first group so their performance metrics is higher than the last 8 QTs which are mostly from the second group. The performance for QT 5 and 13 can be considered slightly different because they actually belong to the third group. The results for QT 6 and 14 are higher than the rest of their group, respectively. That is probably enhanced by the mechanical structure of the device as Chain 3 has a revolute joint instead of a universal joint for a passive joint and yzplane is considered a symmetric plane for the device. Note that the hypothesis can be over-fitting the data because the number of QTs in different groups are limited.

Conclusion and Future Work
In this paper, we have revisited our proposed approach to solve the problem of rendering stiff interaction as a position controller design problem and extended this approach with gain scheduling so that the performance drop in the lower half of the wall can be improved. The haptic rendering case is still the interaction with a stiff wall. The performance improvement is investigated through Monte-Carlo simulations under new "pressing and sliding" conditions. From the results and discussions, we can conclude that the proposed method is able to render a stiffness that is on the magnitude of 10 4 N/m on TAU device stably. This is a huge improvement against the traditional OLIFC as both guidelines and tests suggest that the stable stiffness it can render is on the level of 10 3 N/m. Note that similar conclusions are also reported in [13] that they also achieve high stiffness beyond stability limits suggested by papers like [1,10]. Compared to the previous study [34], the new scenario of sliding investigate performance under a different condition which resembles more to the task of polishing in dental procedures instead of drilling. The performance is consistent with previous findings in magnitude but it proves hard to estimate the performance along the trajectory only with the knowledge of the performance at each point on the trajectory. Instead, the performance at the starting point of the trajectory gives a good estimate of the performance along the whole trajectory. With the newly proposed MOPPC, the transition between areas using different control gains proves to be smooth. Compared to SOPPC, MOPPC shows significant improvement in both scenarios by adopting classic gain scheduling.
Currently, the proposed approach can be used as an alternative solution to deliver high stiffness in rendering stiff objects instead of the traditional spring-damper controller. Nonetheless, it is still at its early phase of development and further research is needed to explore the full potential of the approach being an all-purpose haptic rendering approach, suitable for deformable objects. Furthermore, we expect that more complex scenarios, which require more accurate modeling and faster execution, will reveal more advantages of our proposed approach.

Conflict of Interests
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.