A continuous switching contact model for virtual environment based teleoperation

Virtual environment (VE), as the proxy of slave contact environment, is the most promising technology to solve the time-delay problems in teleoperation. The accuracy of the predicted force depends not only on the reliability of the contact model but also on the estimation algorithm’s adaptability. A new contact model is proposed to be applicable in various materials, which includes both the Kelvin–Voigt model (KVM) and Hunt–Crossley model (HCM). An extra parameter is set in the model to express the capacity of continuous switching between KVM and HCM, whose rationality is proved based on the energy loss. The energy loss is proportional to a power of impact velocity, and the exponent is bounded at [2,3], which exactly lies between KVM and HCM. Furthermore, to estimate the parameters with a single-stage method, the nonlinear model is linearized approximatively with logarithm function and polynomials. Then, the recursive least squares (RLS) algorithm combining forgetting factor and self-perturbing action is designed to identify the four parameters online. Finally, the model’s continuous switching is verified with ideal simulation, and the model parameters are continuously changed without jumpy switch error. In the experiment, sponge, foam, and human hand represent the complex contact materials of the slave environment where the predicted force is shown to follow the real contact force with enough accuracy. Therefore, the virtual model can be considered the substitution of slave contact environment so that the feedback force in master can be calculated in real-time.


Introduction
Teleoperation, as the primary way to extend human skills, has wide-spread applications in bomb disposal, search and rescue, robotic surgery, etc [1]. However, due to the inevitable time-delay in network communication, the message flow between master and slave site cannot be synchronous, usually leading to wrong hand-eye coordination for the operator and even the unexpected instability in haptic teleoperation [2,3]. In [4], much conventional technique for solving time-delay problems are reviewed including passivity theory [5], wave variable method [6], which all concludes that the stability and transparency are two conflicting objectives: one guaranteed and another lost. In [7][8][9], prediction-based methods are proposed, predicting future motion and force based on his- torical data series or model parameters so that the transmitted message in the other end is synchronous. Although the stability and transparency are improved theoretically by these methods, it is difficult to have a good prediction strategy for sudden change and unknown environment.
With the development of computer graphics and image processing, virtual reality(VR) technology is a hot topic for researchers. In 1997, Tsumaki and Uchiyama joined VR technology into teleoperation and analyzed the effects of model errors [10]. As shown in Fig. 1, a virtual environment (VE) is located in the master site to imitate a remote slave environment, which can have a real-time motion and haptic interaction with the operator. To ensure that VE is precise enough to represent an isolated environment, the real slave image is utilized to correct the virtual model using virtualreal fusion technology, and the local force model is updated through real slave force-position data. Due to the real-time interaction with VE in the master site, the VE-based teleoperation is a promising solution to tackle time-delay problems in teleoperation where precise motion and force model are two key difficulties [11,12]. In terms of the force model con-sidered here, an accurate force impedance expression needs to be established for both integrity and rationality. Besides, an adaptive parameter estimation algorithm is necessary for the changeable unknown environment.
Generally, two main force models for 1-D contacts are the linear Kelvin-Voigt model(KVM) and nonlinear Hunt-Crossley model (HCM). KVM is represented by the parallel of a linear spring and a viscous damper, suitable for stiff materials [13,14]. Viscous item in HCM is dependent on penetration depth, which is more consistent with power exchange during contact and more suitable for soft materials [15]. Some standard estimation technology has been used to identify the model parameters, such as Recursive Least Squares(RLS), steepest descent method, and adaptive methods. In [16], a novel 2-stage method is introduced for online estimation for HCM where model parameters are identified linearly in two branches. In [17], the nonlinear HCM is linearized by logarithmic function, and then the parameters are estimated through the RLS family. Furthermore, in [18], corrected KVM and HCM are combined with fitting switch strategy, suitable for complex environments with stiff and soft materials.
Because of the above force model and its estimation method, a novel contact force model is proposed in this paper with enough parameters to describe complex tissue, including stiff and soft materials. Moreover, it is proved that the alternation between KVM and HCM is continuous instead of the jumpy switch according to a threshold in [18]. As it has one more parameter than HCM, this model is much nonlinear and cannot be simplified into a linear equation like HCM. In this paper, a linear combination of polynomials is used here to replace the nonlinear item so that the model is equivalent to a linear one. Meanwhile, a linear RLS-like estimation algorithm is designed to identify the model's parameters with both forgetting factors and self-perturbing. Finally, a VEbased teleoperation experiment is established, which proves that the proposed model can produce accurate real-time contact force using a position in the master site. This is a vital advantage, especially for complex contact environments like peg-in-hole, minimally invasive telesurgery, etc., because of the model's integrity and self-adaption to various materials.
The rest of the paper is organized as follows: in "A new dynamical contact model", a new dynamical contact model is proposed with its continuous switching between KVM and HCM, the rationality of the model is proved with power exchange during impacting. In "Identifying nonlinear parameters", the nonlinear model is linearized with logarithm function and polynomials first, and an RLS-like algorithm is designed to identify the linear parameters. In "Simulation and experiments", simulation and experiments are intended to verify the continuous switching feature and the feasibility of the model. In the final section, we conclude the whole work of the paper.

Model proposed
In VE-teleoperation, the contact model is a proxy of slave environment whose accuracy directly influences the transparency and stability of the system. However, the slave environment is previously unknown and multifarious. Therefore, it is impossible to describe the contact environment as a single case. Take robotic telesurgery, for example, where the instrument is inserted into the human's body to perform cutting, sewing, etc. Encountering with soft tissue is considered soft contact, and stiff regions like bone are rigid contact, and others may lie between the upper two. To describe such a contact environment, the desired model needs to satisfy the following requirements: -Physical consistency Model needs to generalize dynamical process during contact and satisfy the physical rule for both power exchange and momentum reaction; -Unified formula Single equation expression is needed to describe soft contact, stiff contact, and any other contact environment between the two; -Environment adaption Parameters in the model can be identified online to respond rapidly to the unknown and changeable slave environment.
The simplest dynamical contact model is the linear Kelvin-Voigt model, which has been confirmed to have the right solution in stiff contact. Its expression is such pure superposition with linear spring and linear damper and is described in 1-D space as where x(t) is the penetration depth into contact material,ẋ(t) is the penetration velocity, and K , B are the coefficients of spring and damper, respectively. From Eq. (1), we can see that the damper item is proportional to penetration velocity without any dependence on contact depth. That's to say, force extended by damper follows the change of speed whenever an impact is, which is usually inconsistent with applicable physical rule due to absolute velocity at beginning or end. According to the KVM calculation, there is also some contact force at this time, but the actual contact phenomenon does not produce the contact force at these two moments.
When simulating contact with a mass m on the material surface using KVM, the hysteresis loop and exchange power during contact are shown in Fig. 2.
In Fig. 2a, the compression phase corresponds to the upper arc, connecting point A to the maximum penetration x M , and  [16] H 1 represents energy that is extracted from the mass in this stage, which is equal to the area enclosed by the upper arc, the abscissa axis, the ordinate axis and the line x = x M . While the lower arc represents the restitution phase, and H 2 is the energy flowing from the touched material to the mass, hence the minus sign. The absolute value of H 2 is equal to the area enclosed by the lower arc, the abscissa axis and the line x = x M . H 3 represents energy that is extracted from the mass, even in the restitution phase, whereẋ(t) < 0, and this behavior is in contrast both with physical intuition and experimental results. These physical inconsistencies happen because when the penetration is small, the force F in Eq. (1) mainly depends on the damping term, which is assumed to be independent on x [16].
These energies are also plotted in Fig. 2b, showing the power flow obtained by multiplying the reaction force F(t) by the penetration velocityẋ(t)(P(t) = F(t)ẋ(t)). Positive areas represent energy flowing from the mass to the touched material and vice versa. We can see that all of the H i represent the areas enclosed by the axis and the arc, which correspond to Fig. 2a. A bizarre phenomenon appears that the power is extracted from contacting mass at the end of the restitution period because of the unexpected opposing force. Moreover, the KVM has terrible results in applying soft materials, which is widespread in many areas like operative procedures.
Because of shortcomings in KVM, Hunt and Crossley proposed a nonlinear dynamical model named HCM based on restitution coefficient in 1975 [15]. The physi-  [16] cal inconsistencies are avoided after using the exponential term of displacement in the damper item, and the nonlinear spring item is the same as Hertz's model. The expression is described as where k, λ are constant coefficients of nonlinear spring and damper, respectively, and n is exponential index which is always below 2.
As the nonlinear damper item is dependent on penetration depth, the shock force disappears at both beginning and end so that no power is extracted from contact mass at the end of the restitution period (see Fig. 3). The hysteresis loop shapes like a whole ellipse with its energy loss proportional to the cube of impact velocity. KVM is like a half ellipse, and energy loss is proportional to the square of impact velocity. Although its physical consistency to real dynamical contact process, HCM is not suitable for stiff materials as many experiments showed. Besides, due to the homology of exponential index n in spring and damper item, the estimation error of n have similar effects on both things, which leads to larger predicted force error in VE-teleoperation.
Considering the adaption of KVM and HCM to different contact materials, Achhammer and Weber proposed a threshold switching contact model (CSCM) where a threshold is used to switch between KVM and HCM. When the pre-dicted stiff parameter exceeds the threshold, the used model switches from HCM to KVM. The switching operator S(·) is where k kv is the estimated stiff parameter with KVM, and k th is the threshold of switching strategy. From Eq. (3), it is evident that the change from KVM to HCM is jumpy like 0-1 switching value and vice versa, at which moment the predicted force is usually shocking. What's more, according to the nature of gradually changeable materials, the transformation of the model structure should also be continuous. Encouraged by this idea, a novel continuous switching contact model (CSCM) is proposed as where n, α are adjustment parameters decided by materials for alteration between KVM and HCM(s = n+α), and others are the same as HCM. Considering the empirical value from previous researchers, n and α are objective to The damper includes the nonlinear item x n (t) so that the shock force in KVM is also avoided as HCM. Power loss during impacting is also consistent with physical rules as the quasi ellipse in Fig. 3. Besides, the continuous alteration between KVM and HCM is realized by adjusting parameters of n and α, which will be proved based on power loss in the following section.

Continuous switch between KVM and HCM
When a mass m contacting with viscoelastic material, the power of whole system is divided into three components: kinetic energy of mass, stored energy in material and dissipated power results from heat, irreversible deformation etc. This analysis is confined to compact solid bodies (normally regarded as elastic) impacting one another at low velocities, and this physical model holds if the object has a certain initial velocity and there is no external disturbance. According to the theorem of kinetic energy, the power exchange equation is where v i and v are the relative speeds respectively "in" and "out", E store and E diss are caused by the spring and damper item in model. Due to low impacting velocity when contacting, the E diss is small enough to be negligible contrast to E total . So when maximum penetration depth (assumed x M ) is reached, we can get that For intermediate process(set depth as x), the dynamical velocity is obtained by Eq. (6): Thus, energy loss decided by damper item within contact loop including compression and restitution can be acquired as then Eq. (9) becomes where . We can simplify this by making two substitutions. Let t = y/y M and u = t 2 , we have where F(β, 1.5; 2.5; 1) = h(β) is Gaussian hypergeometric function only relative to parameter β, whose function variation curve is drawn as Fig. 4a. In whole span of β, function h(β) increases with the increase of β and is bounded at [1, 2.5). According to Eqs. (7) and (10), we can get that v i ∝ y M . If the contact environment is given, the contact model, representing the dynamical behaviors, is also identified with four known parameters. Then, h(β) is a constant value and Eq. (12) becomes where γ = 3 − 2β and K 0 is a scale coefficient.
As the Eq. (13) shows, the energy loss ( E) is proportional to a power of impact velocity and the exponent is decided by γ , which changes with the different model parameters. According to inequations of Eq. (5), the available range of α − n can be drawn as the coloured area of Fig. 4b and the value of γ is calculated by where Concluded from the Fig. 4b and Eq. (14), the value of γ can be adjusted by (α, n) combination and is bounded at [2,3]: γ = 2 is the same as energy loss of KVM and γ = 3 is HCM. Besides, the varying of γ is continuous over the whole valid area of α − n, which also means that CSCM model can switch between KVM and HCM continuously.

Linearized scheme for CSCM
The four parameters (k, λ, n, α) in the CSCM model are decided by contact materials, which are usually unknown in advance. In VE-teleoperation, the local contact model is located at master to calculate the virtual haptic force in real-time, before which the parameters must be identified precisely according to slave position-force (P-F) data series. To be applicative in VE-teleoperation, estimation technology for parameters must satisfy two conditions: (1) The parameters are estimated online to follow the change of unknown environment; (2) The calculation is accurate and real-time, which means that the parameters will update immediately every time a new P-F data arrives, and the calculation is small.
Based on previous analysis, the CSCM model expressed as Eq. (4), is nonlinear with two exponential items (k · x s and λ·x n ·ẋ), which is more complex than HCM. For estimation of nonlinear parameters, the conventional strategy is to linearize the equation first with Taylor expansion, taking logarithm or others, and then use traditional estimation methods like RLS to identify linear parameters. The steepest descent method (SDM) is also a universal method for nonlinear estimation, which needs a good knowledge of initial values for its local convergence. Some intelligent algorithms, such as Genetic Algorithm or Neural Network, may be useful alternations for nonlinear estimation because of their iterative learning ability.
For CSCM model, the two exponential items (k · x s and λ · x n ·ẋ) can be easily linearized by moving one to the left side of equation by using two-stage method [16]. But due to the interconnection between both stages, the error ofk,ŝ will have big effects on next stage(λ,n), causing larger errors tô k,ŝ in turn. Considering this inapplicability, a reasonable linearized scheme is proposed as follows and can be estimated with single-stage method. To linearize the CSCM model, the following requirements need to be satisfied: (1) ϕ = λ k ·ẋ is near zero; (2) the value of x is small; First, taking logarithm function on both sides of Eq. (4): where ϕ = λ k ·ẋ is considered to be near zero ( ϕ < 0.1) [17], especially for teleoperation. And our experiments and simulations have shown that 0.1 is a rather conservative value in many cases. With regard to the feasibility of condition ϕ < 0.1, in many practical applications, the speed of operation within contact is not high and condition is often met. However, it should be noted that the speed of operation is not fully controlled by the user as it depends on the characteristics of the desired task and the capabilities of the robot. Therefore, the fact that condition is not guaranteed can be considered a shortcoming of this method.
And function g(ϕ) = ln (x α + ϕ) can be approximated around ϕ = 0 with Taylor Series expansion method: According to [19], the exponential term (x n , n ∈ [1, 2]) can be approximated with linear combination of primary and quadratic item if the value of x is small. For contact facts, the penetration depth is usually small enough in case of destroying the contact materials. So, the 1/x α in Eq. (16) can be replaced as following the linear combination where a 1 , a 2 are polynomial coefficients, depending on the value of x and α.
In addition, after careful observation, the two limiting cases a 1 = 0 and a 2 = 0 represents the form of KVM (α = 1) and HCM (α = 0) respectively. In other words, the Eq. (17) is a linear combination of KVM and HCM, another description for the continuous alternation of model.
Summarizing from Eqs. (15), (16), (17), can be linearized approximated as where b 1 = a 1 · λ k and b 2 = a 2 · λ k are the linear coefficients and can be estimated gradually within the range of x.
The Eq. (18), transformed from nonlinear model (Eq. (15)), is a good linear equation whose parameters (ln k, s, b 1 , b 2 ) can be handily estimated through many available linear methods. Before putting it into practice, its validity need be verified to ensure the accuracy of predicted force. Here, we use Matlab software to calculate force ideally with Eqs. (15), (18) respectively. The desired trajectory is set by The model parameters, having key influence on the approximated accuracy of Eqs. (16) and (17), are divided into six groups as Table 1.
Based on given α and trajectory, the a 1 and a 2 in Eq. (17) are identified successively by RLS method, which are also listed in the Table 1. With these model parameters, contact force is calculated with both Eqs. (15) and (18), which is expressed with line and dot sign in Fig. 5a respectively. Although constant estimated values (â 1 ,â 2 ) are used in Eq. (18), the approximated force is also lying closely on the whole curve of Eq. (15). Moreover, the error between both equations, as drawn in Fig. 5b, is too small to be sensed for haptic rendering (most below 0.01 N) and some big errors are mostly introduced by constant estimation valueâ 1 ,â 2 , which can be compensated through adaptive RLS algorithm. Therefore, the nonlinear model can be replaced by a linear equation Eq.14 Eq.17

Fig. 5 Errors of linearized equation
Eq. (15) within the permissible range, whose parameters will be identified with the following linear algorithm.

Algorithm design
After linearizing from Eqs. (15) to (18), the model parameters can be identified using RLS family methods. Firstly, the equation is summarized as where T k is regressor vector at the time of k, relative to the input of x,ẋ, ε k can be considered as Gauss noise consisting the measure and approximating error; θ k is the linear parameters which will be updating iteratively based on last time. Each item is expressed as For changeable contact environment, the model parameters may change gradually as the penetration depth, and even sometimes alter abruptly with a big step value. In order to be adaptive in different cases, estimation algorithm need make responds to these changes rapidly. Here, we design a hybrid RLS method (named FPRLS), which combines both forgetting factor and self-perturbing action to estimate the slowly and steeply changeable parameters respectively. The FPRLS algorithm is iterated as the following equations: where P k is 4 × 4 covariance matrix, K k is gain value representing contribution of estimation error on parameters. And the forgetting factor is adaptive according to penetration depth which is When penetration depth changes slowly or remains unchanged, the forgetting factor μ is near 1.0 having no forgetting effect on old data but depth changing fast makes a forgetting effect. This is a big advantage on changeable gradually parameters, also compensating for variable coefficients of Eq. (17). On the other hand, the sudden changed parameters are solved through self-perturbing action, which is defined as k+1 ·θ k and δ is a very big value. The covariance matrix P k will return back to be big which is set initially when estimated error is big enough; that's to say, the iteration will no longer consider old data once the sudden change happens. This is especially necessary for sudden changeable environment, such as from soft tissue to stiff bone in telesurgery.
After every step of estimation, the estimated parameters are substituted into Eq. (18) to calculate real-time contact force. The accuracy between estimated and real force depends on not only the validity of real P-F data series but also the adaptive ability of FPRLS. To obtain good convergence, the penetration trajectory need contain enough kinds of sine waves. The convergence and estimated error are analysed in the following section with both simulation and experiment.

Simulation verifying
In this section, we will make use of numerical simulation (Matlab) to verify the validity of CSCM model and FPRLS algorithm, which are proposed in above two sections. The contact action are abstracted by a needle penetrating into the surface of different materials and impact force is the function of penetration depth and velocity with a ideal contact model. Considering the conditions of convergence and low velocity in practice, the desired trajectory is set as  Fig. 6). Then, the contact force is calculated according to the depth and velocity with presumptive model which is the substitute of ideal contact material and is divided into following three types.

Estimation for KVM and HCM
According to given penetration trajectory, the ideal contact force is calculated by pure Kelvin-Voigt and Hunt-Crossley model, as depicted in Eqs.  Fig. 7a, the estimated k, s converge to set values at 20th iteration. As can be seen from the upper two subfigures, b 1 is not at zero but random value within (−1,1) because of the modelling error and force noise; b 2 , the coefficient of v/x, gradually goes to y = 0.083 line. For HCM, the ideal model is realized by setting α = 0 in Eq. (4), resulting in b 2 = 0 in Eq. (18). The estimated result is shown in Fig. 7a as pinkish red line with k, s also converging to set values rapidly. However, different from KVM, the value b 2 converges to zero at 10th iteration due to no "v/x" item in HCM. Therefore, the proposed model is applicable in both

Continuous switching from KVM to HCM
The most important advantage of proposed model is its adaptivity to continuous switching between models, which ensures that the estimated force is smoothly changing without the jump error in [18]. This part will imitate the continuous alteration ideally from HCM ( 2 ) to KVM ( 1 ) as follows Then, the simulated PVF datas are utilized in both CSCM model and TSCM [18] to obtain estimation force error as where F andF are the ideal and estimated force calculated with ideal and estimated parameters respectively.
The result of both methods is shown in Fig. 8 with curve and dot sign. For TSCM, the model switch from HCM to KVM just at the 198th iteration with a stiff threshold of 4000. Because of the 0-1 switching, the shock force appears from 4.5N to 0N, and at other times force error is bigger than CSCM because of the inaccuracy of the model. At each time, the varied model can be expressed by CSCM with different values of α, but for TSCM, only two forms of models are used to represent all continuous models: under the threshold is HCM and above is KVM, which contraries with the natural features of constant materials. Besides, after using forgetting factor μ in the FPRLS algorithm, the estimated parameters alter gradually between the two limit values (as shown in Fig. 7b), following the changes of parameters. Therefore, when the contact model changes, the CSCM model can alter gradually with an extra adjusting parameter, which is more reasonable and precise than the switch operator as proposed in [18].

Suddenly changed parameters
Sometimes, the contact environment may change suddenly in practice, such as step switch from soft to stiff material at a certain depth. For such case, the contact model must adjust its parameters immediately, responding to the new materials. Here, we assume that the sudden change from HCM ( 2 ) to KVM ( 1 ) is divided into following three steps When estimation errorê k ≥ 1N , the covariance matrix return to 10 6 · I as set at the beginning, which also means that the old data are not considered anymore, and parameters will converge to new ones rapidly. The estimation result is shown in Fig. 7c with two sudden changes at 70th and 150th iteration. At the first change, the model transforms to media form between KVM and HCM with α = 0.5 and the approximated polynomial coefficients in Eq. (17) is a 1 = 3.44, a 2 = 0.067. This is verified in the upper two subfigures of Fig. 7c with b 1 ≈ 0.5, b 2 ≈ 0. Besides, k, s also converge to real step values rapidly as desired at each sudden change, and the estimated force error is small enough to be ignored.

VE-teleoperation testing
To validate the practical applicability of the CSCM model, the VE-teleoperation experiment is constructed as Fig. 9 with a master device, a slave robot, and both controlling PCs. The interaction message flows are controlled as Fig. 1 with position sent from master to slave and force conversely via local area network (LAN). The slave robot is an RCM (remote center of motion) mechanism with a linear motion at its end, which can imitate the penetration action as a probe. Moreover, an ATI force sensor is mounted centrally to measure the contact force with its z-axis toward the penetration direction. The penetration depth and velocity can be acquired with the motor's encoder and are recorded in slave PC. For simplicity, Fig. 9 Experiment of VE-teleoperation the end-rod of the slave robot is kept vertical to reduce other motions' disturbance.
During the experiment, the operator moves the end of the master device to control the slave robot to contact with the environment materials. To enlarge the time-delay in force sending, a data buffer is used in the slave PC to store a certain number of force data, and the time-delay is near 300 ms approximately. But the time-delay of position is only 80ms because of the flow velocity and the small time-delay in LAN. In slave, the CSCM model parameters are updated with the FPRLS algorithm every time the measured PVF data are obtained. Then the parameters and force in a slave are sent together to the master via LAN. In master, the received parameters are used to predict the contact force according to the current position and velocity of the master device and then send to the operator for haptic rendering. The contact materials are set as foam, sponge, and the muscle of the human hand, whose results are illustrated as follows:

Sponge
As a soft and cellular material, the sponge has wide applications in shock absorption, cleaning, etc. Due to its delinquent recovery, elastic stiffness in restitution is usually smaller than in compression. In the experiment, the operator controls the end of the slave robot to penetrate the sponge with a certain depth and then repeatedly operates compression and restitution. The estimated parameters is recorded in real-time as shown in Fig. 10a. The parameters are stabilized at the 50th iteration, after which the varieties of parameters shape like a sawtooth wave. The sawtooth shape corresponds with the sponge's feature: the stiffness in compression is big and small in restitution. The predicted force in the master is shown in Fig. 11a in contrast to delayed force. It can be seen that the predicted force in master almost lies on real slave force at the same time but delayed force has a big offset. In other words, the rendering haptic with predicted force is much more synchronous than delayed force.

Foam
Unlike the sponge, the foam has many micropores and is stiffer when the penetration depth is bigger. Consequently, foam's contact property changes from soft to stiff in compression and inversely in restitution, which is especially suitable for the continuous testing of the CSCM model. The estimated parameters and predicted force are shown in Fig. 10b, 11b, respectively. The parameters change dramatically and can easily mutate. The predicted force in the master is also corresponding with the slave force within acceptable ranges, without incredible haptic sense on the operator.

Human hand
To prove the CSCM model's feasibility in a complex environment, a human hand is used here to represent the multivariate contact materials. In the experiment, a human hand is put statically with its thenar muscle just below the end of RCM. The slave robot is controlled to contact with the hand slowly with its force under 15N, in case of any harm on human. The result of this experiment is shown in Figs.10c, 11c. The estimated parameters have three convex peaks in each subfigure, where the penetration depth is up to the biggest in its cycle. At the second peak, the stiffness (k) is up to 150,000 N/m, which may happen when the end is contacting with the stiff bone. The other two peaks are smaller and can also be explained as the contact with the stiff muscle inside the hand. All these self-adaptiveness to the suddenly changed environment are owing to the self-perturbing of the FPRLS algorithm. As seen in Fig. 11c, the predicted force in the master site is accurate enough to be able to represent the factual slave contact force, either in the epidermal fat, stiff muscle tissue, or even the stiff bone. In other words, the proposed model is so integrated that the complex hand muscle can also be expressed, which is realized by adjusting its four parameters. Therefore, the master PC can calculate the real-time contact force based on the current parameters and kinestate of the master device, instead of the delayed force from the slave.

Conclusion
In this paper, a new contact force model is proposed to express the continuous changeable materials, avoiding the shock effect from sudden switching between KVM and HCM. The model is aimed to imitate the factual slave environment of VE-teleoperation with its parameters updating gradually. Firstly, the model is summarized based on the respective advantages of KVM and HCM in stiff and soft materials. By analyzing the energy during contact, the model has been proved to have a changeable energy loss between KVM and HCM. In other words, the model can switch continuously between the two models with its parameters. To estimate the four model parameters, the original nonlinear model is transformed to a linear equation firstly, whose approximated error is verified to be small enough to be negligible. Moreover, an FPRLS algorithm is adopted with forgetting factor and self-perturbing action to identify continuous or sudden changed parameters. To validate the feasibility of CSCM model and FPRLS algorithm, a simulation scheme is implemented in Matlab with an ideal exciting penetration trajectory. The perfect material models are set as pure KVM or HCM, continuous switching, and sudden change successively. It is confirmed that the model is integrated enough to express various contact models and has a unique advantage to avoid the shock force at switch moment. Finally, a VE-teleoperation is established with master and slave robots communicating through the LAN and an extra force data buffer. Sponge, foam, and human hand are used to represent the complex contact materials of the slave environment.
In each experiment, the estimated parameters have shown consistency with the feature of the material. The predicted force, calculated by the local contact model, is accurate for either changeable or sudden changed materials. It can be concluded that the proposed model is precise enough to represent the slave environment, which is an excellent advantage to improve the transparency and stability of the teleoperation system.