Year: 2020

Abstract:
Although the Bayesian paradigm provides a rigorous framework to estimate the full probability distribution over unknown parameters, its online implementation can be challenging due to heavy computational costs. This paper proposes Adaptive Recursive Markov Chain Monte Carlo (ARMCMC) which estimates full probability density of model parameters while alleviating shortcomings of conventional online approaches. These shortcomings include: being solely able to account for Gaussian noise, being applicable to systems with linear in the parameters (LIP) constraint, or having requirements on persistence excitation (PE). In ARMCMC, we propose a variable jump distribution, which depends on a temporal forgetting factor. This allows one to adjust the trade-off between exploitation and exploration, depending on whether there is an abrupt change to the parameter being estimated. We prove that ARMCMC requires fewer samples to achieve the same precision and reliability compared to conventional MCMC approaches. We demonstrate our approach on two challenging benchmark: the estimation of parameters in a soft bending actuator and the Hunt-Crossley dynamic model. Our method shows at-least 70% improvement in parameter point estimation accuracy and approximately 55% reduction in tracking error of the value of interest compared to recursive least squares and conventional MCMC.

heading: 1 INTRODUCTION
text: Bayesian methods are powerful tools to not only obtain a numerical estimate of a parameter but also to give a measure of confidence (Kuśmierczyk et al., 2019; Bishop, 2006; Joho et al., 2013). In particular, Bayesian inferences calculate the probability distribution of parameters rather than a point estimate, which is prevalent in frequentist paradigms (Tobar, 2018). One of the main advantages of probabilistic frameworks is that they enable decision making under uncertainty (NoormohammadiAsl & Taghirad, 2019). In addition, knowledge fusion is significantly facilitated in probabilistic frameworks; different sources of data or observations can be combined according to their level of certainty in a principled manner (Agand & Shoorehdeli, 2019). Nonetheless, Bayesian inferences require high computational effort for obtaining the whole probability distribution and prior general knowledge on noise distribution before estimation.
One of the most effective methods for Bayesian inferences is the Markov Chain Monte Carlo (MCMC) methods. In the field of system identification, MCMC variants such as the one recently proposed by Green (2015) are mostly focused on offline system identification. This is partly due to computational challenges which prevent real-time use (Kuindersma et al., 2012). The standard MCMC algorithm is not suitable for model variation since different candidates do not share the same parameter set. Green (1995) first introduced reversible jump Markov chain Monte Carlo (RJMCMC) as a method to address the model selection problem. In this method, an extra pseudo random variable is defined to handle dimension mismatch. There are further extensions of MCMC in the literature, however, an online implication of it has yet to be reported.
Motion filtering and force prediction of robotic manipulators are important fields of study with interesting challenges suitable for Bayesian inferences to address (Saar et al., 2018). Here, measurements are inherently noisy, which is not desirable for control purposes. Likewise, inaccuracy, inaccessibility, and costs are typical challenges that make force measurement not ideal for practical use (Agand et al., 2016). Different environmental identification methods have been proposed in the literature
for linear and Gaussian noise (Wang et al., 2018); however, in cases of nonlinear models like HuntCrossley that does not have Gaussian noise (e.g. impulsive disorder), there is no optimal solution for the identification problem. Diolaiti et al. (2005) proposed a double-stage bootstrapped method for online identification of the Hunt-Crossley model, which is sensitive to parameter initial conditions. Carvalho & Martins (2019) proposed a method to determine the damping term in the Hunt-Crossley model. A neural network-based approach was introduced to control the contact/non-contact HuntCrossley model in (Bhasin et al., 2008)
This paper proposes a new technique, Adaptive Recursive Markov Chain Monte Carlo (ARMCMC), to address address certain weaknesses of traditional online identification methods, such as only being appllicable to systems Linear in Parameter (LIP), having Persistent Excitation (PE) requirements, and assuming Gaussian noise. ARMCMC is an online method that takes advantage of the previous posterior distribution, given that there is no sudden change in the parameter distribution. To achieve this, we define a new variable jump distribution that accounts for the degree of model mismatch using a temporal forgetting factor. The temporal forgetting factor is computed from a model mismatch index and determines whether ARMCMC employs modification or reinforcement to either restart or refine parameter distribution. As this factor is a function of the observed data rather than a simple user-defined constant, it can effectively adapt to the underlying dynamics of the system. We demonstrate our method using two different examples: soft bending actuator and Hunt-Crossley model and show favorable performance compared to state-of-the-art baselines.
The rest of this paper is organized as follows: In Sec. 2, introductory context about the Bayesian approach and MCMC is presented. Sec 3 is devoted to presenting the proposed ARMCMC approach in a step-by-step algorithm. Simulation results on a soft bending actuator with empirical results on a reality-based model of a soft contact environment capturing a Hunt-Crossley dynamic are presented in Sec. 4. Lastly, the final remarks and future directions are concluded in Sec 5.

heading: 2 PRELIMINARIES

heading: 2.1 PROBLEM STATEMENT
text: In the Bayesian paradigm, estimates of parameters are given in the form of the posterior probability density function (pdf); this pdf can be continuously updated as new data points are received. Consider the following general model:
Y = F (X, θ) + ν, (1)
where Y , X , θ, and ν are concurrent output, input, model parameters and noise vector, respectively. To calculate the posterior probability, the observed data along with a prior distribution are combined via Bayes’ rule (Khatibisepehr et al., 2013). The data includes input/output data pairs (X,Y ). We will be applying updates to the posterior pdf using batches of data points; hence, it will be convenient to partition the data as follows:
Dt = {(X,Y )tm , (X,Y )tm+2, · · · , (X,Y )tm+Ns+1}, (2) where Ns = Ts/T is the number of data points in each data pack with T, Ts being the data and algorithm sampling times, respectively. This partitioning is convenient for online applications, as Dt−1 should have been previously collected so that the algorithm can be executed from tm to tm + Ns+1 or algorithm time step t. Ultimately, inferences are completed at tm+Ns+2. Fig. 1 illustrates the timeline for the data and the algorithm. It is worth mentioning that the computation can be done in parallel by rendering the task of the adjacent algorithm step (e.g. phase A of algorithm t, phase B of algorithm t− 1 and phase C of algorithm t− 2 can all be done simultaneously) According to Bayes’ rule and assuming data points are independent and identically distributed ( i.i.d) in equation 1, we have
P (θt|[Dt−1, Dt]) = P ( Dt|θt, Dt−1 ) P (θt|Dt−1)∫
P ( D1|θt, Dt−1 ) P (θt|Dt−1)dθt , (3)
where θt denotes the parameters at current time step. P (θt|Dt−1) is the prior distribution over parameters, which is also the posterior distribution at the previous algorithm sampling time. P ( Dt|θt, Dt−1 ) is the likelihood function which is obtained by the one-step-ahead prediction:
Ŷ t|t−1 = F (Dt−1, θt), (4)
where Ŷ t|t−1 is the prediction of the output in (1). If the model in (4) is valid, then the difference between the real output and predicted should be measurement noise, (i.e., Y t|t−1 − Ŷ t|t−1 = ν). Therefore, the model parameter may be updated as follows:
P ( Dt|θt, Dt−1 ) = tm+Ns+1∏ tm+1 Pν ( Y t|t−1 − Ŷ t|t−1 ) , (5)
where Pν is the probability distribution of noise. Note that there is no restriction on the type of noise probability distribution.
Remark 1: As it was mentioned before, there is no need to know the exact probability distribution of noise. This probability distribution can be simply substituted with a Gaussian distribution, if one has minimal knowledge of the mean and variance of the data which can be easily obtained with preprocessing (Bishop, 2006).

heading: 2.2 MARKOV CHAIN MONTE CARLO
text: MCMC is often employed to compute the posterior pdf numerically. The multidimensional integral in (3) is approximated by samples drawn from the posterior pdf. The samples are first drawn from a different distribution called proposal distribution, denoted q(.), which can be sampled easier compared to the posterior. Brooks et al. (2011) discuss different types of MCMC implementations which may employ various proposal distributions and corresponding acceptance criteria. The main steps of the Metropolis-Hastings algorithm are listed as follows (Ninness & Henriksen, 2010):
1. Set initial guess θ0 while P (θ0|Y ) > 0 for iteration k = 1, 2. Draw candidate parameter θcnd, at iteration k, from the proposal distribution, q(θcnd|θk−1) 3. Compute the acceptance probability,
α(θcnd|θk−1) = min {
1, P (θcnd|D)q(θk−1|θcnd) P (θk−1|D)q(θcnd|θk−1)
} , (6)
4. Generate a uniform random number γ in [0, 1], 5. ‘Accept’ candidate if γ ≤ α and ‘ignore’ it if γ > α, 6. Set iteration to k + 1 and go to step 2.

heading: 2.3 PRECISION AND RELIABILITY
text: Two important notions in probabilistic framework to compare results are precision ( ) and reliability (δ). The former represents the proximity of a sample to the ground truth, and the latter represents the probability that an accepted sample lies within of the ground truth.
Lemma: Let Pk be k samples from MCMC, and E(Pk) denote its expected value. According to Chernoff bound (Tempo et al., 2012), given , δ ∈ [0, 1], if the number of samples (k) satisfies
k ≥ 1 2 2 log( 2 1− δ ), (7)
then Pr { {Pk − E(Pk)} ≤ } ≥ δ.
Algorithm 1 ARMCMC Assumptions: 1) roughly noise mean (µν) 2) roughly noise variance (σν) 3) desired precision and reliability ( 0, δ0) 4) desired threshold for model mismatch (ζth) Goal: Online calculation of parameters posterior distribution given the consecutive t-th pack of data (P (θt|Dt)) Initialization: Prior knowledge for θ01 , n=0 Consider desire precision and reliability ( , δ) repeat
Put t0 = n ∗Ns + 1 from (2), n++ Add new data pack to dataset Dt Model mismatch index: ζt from (10) if ζt < ζth then
Reinforcement: set prior knowledge equal to the latest posterior of previous pack Temporal forgetting factor: λt from (9)
else Modification: set prior knowledge θn1 Temporal forgetting factor: λt = 0 end if Set minimum iteration kmin from (12) for k = 1 to kmax do
Proposal distribution: • draw λk ∼ U(0, 1) • Variable jump distribution: qtk(.) from (8) Draw θt∗k ∼ qtk(.) Acceptance rate: α(.) from (6) Draw γ ∼ U(0, 1) if γ ≤ α then
‘Accept’ the proposal end if
end for Wait to build Dtm+Ns+1t0 (algorithm sample time)
until No data is obtained

heading: 3 ARMCMC ALGORITHM
text: At each time interval, ARMCMC recursively estimates the posterior distribution by drawing samples. The number of samples drawn is constrained by the desired precision and reliability, and the real-time requirement. On the other hand, the maximum number of data points in each data pack, Ns, is limited by the frequency of model variation, and the minimum is confined by the shortest required time such that the algorithm is real-time. We propose a variable jump distribution that enables both enriching and exploring. This will necessitate the definition of the temporal forgetting factor as a measure to reflect current underlying dynamics of the data. In other words, this parameter will show the validity of the previous model for the current data. We also prove that ARMCMC achieves the same precision and reliability with fewer samples compared to the traditional MCMC. Algorithm 1 summarizes ARMCMC.

heading: 3.1 VARIABLE JUMP DISTRIBUTION
text: We propose a variable jump distribution (also known as a proposal distribution) to achieve faster convergence, thereby enabling real-time parameter estimation:
qtk(θ t|θtk−1) = { P (θt−1|Dt−1) λk ≤ λt N(µD, σν) λk > λ t , (8)
where θtk−1 is the (k − 1)-th parameter sample which is given by the t-th data pack throughout the MCMC evaluation. Averaging the second half of this quantity will construct θt. P (θt−1|Dt−1) is the posterior distribution of the parameters at the previous algorithm time step, and N(µD, σν) is a Gaussian distribution with µD, σν computed using sample-based mean and variance of Dt−1.
The hyperparameter λt (temporal forgetting factor), is an adaptive threshold for the t-th pack that takes inspiration from in the classical system identification terminology; it regulates how previous knowledge affects the posterior distribution. Smaller values of λt intuitively means that there may be a large sudden change in θ, and thus more exploration is needed. Conversely, larger values of λt is appropriate when θ is changing slowly, and thus previous knowledge should be exploited. As more data is obtained, better precision and reliability can be achieved.

heading: 3.2 TEMPORAL FORGETTING FACTOR
text: Depending on whether the distribution of the parameter θ has changed significantly, a new sample can be drawn according to the modification or the reinforcement mode. Reinforcement is employed to make the identified probability distribution more precise when it is not undergoing sudden change. Modification is employed otherwise to re-identify the distribution ”from scratch”. Therefore, we define a model mismatch index, denoted ζt, such that when it surpasses a predefined threshold (ζt > ζth), modification is applied. In other respects, if ζt ≤ ζth, then ζt is used to determine λt as follows:
λt = e−|µν−ζ t|, (9)
where µν is an estimation of the noise mean, by calculation of the expected value in relation (1). Note that employing modification is equivalent to setting λt = 0. The model mismatch index ζt itself is calculated by averaging the errors of the previous model given the current data:
ζt = 1/Ns Ns∑ n=1 ( ytn − E θ∈θt−1 (F (Dt(n), θ)) ) , ζ0 =∞ (10)
Remark 2: The model mismatch index accounts for all sources of uncertainty in the system. To calculate ζth, one needs to precalculate the persisting error between the designated modeled and measured data. In other words, ζth is basically an upper bound for the unmodeled dynamics, disturbances, noises, and any other source of uncertainty in the system.
Remark 3: To avoid numerical issues, the summation of probability logarithms are calculated. In addition, each pair in the algorithm time sample is weighted based on its temporal distance to concurrent time. Therefore Eq. (5) is modified as
log ( P (·) ) = tm+Ns+1∑ tm+1 logPν(e t),
et = ( Y tn − F t(Dt−1(n), θt) ) e−ρ(Ns−n),
(11)
where ρ ∈ [0, 1] is a design parameter that reflects the volatility of the model parameters, and et = [et1, ..., e t n, ..., e t Ns
]. For systems with fast-paced parameters, ρ should take larger values. We are trying to solve the Bayesian optimization problem that gives us the pdf of the model parameters (θ), when given the data pairs in the presence of uncertainty.

heading: 3.3 MINIMUM REQUIRED EVALUATION
text: Theorem 1. Let and δ be the desired precision and reliability. . Furthermore, it can be assumed that the initial sample has enough number of evaluations (as in (7)). To satisfy the inequality in Eq. (7), the minimum number of samples k in ARMCMC is calculated using this implicit relation:
kmin = 1
2 2 log(
2
λt(1− δ) + 2(1− λt)e−2 2(1−λt)kmin ). (12)
Proof. Samples from previous pdf: According to the variable jump distribution in (8), given k samples, the expected number of samples drawn from the previous posterior probability (P (θ|Dt)) is λtk. By assumption, the algorithm has already drawn at least k samples in the previous algorithm time-step. Consequently, by (7), the expected number of samples with distances less than from E(Pk) drawn from a previous distribution is at least λkδ.
Samples from Gaussian: By (8), there are k0 = (1−λt)k samples drawn in expectation. According to (13), we have Pr { {Pk − E(Pk)} ≤ } ≥ δ0, where δ0 is given by rearranging (7):
δ0 = 1− 2e−2 2k0 . (13)
Thus, the expected number of samples with distances less than from E(Pk) are at least δ0(1−λt)k. Overall reliability: The total expected number of samples with distances less than from E(Pk) is the summation of the two parts mentioned above. Hence it is obtained through dividing by k:
δ1 = (λtkδ) + (δ0(1− λt)k)
k (14)
Given the new obtained reliability, which is greater than the desired one, helps us decrease the number of evaluations. For the sake of illustration, Fig. 2 presents the minimum required number of evaluations with respect to λ for different precisions and reliabilities. As it can be seen, the MCMC is equal to ARMCM if λ is always set to one. The number of evaluations in ARMCMC mitigates as the validity of the previous model increases.

heading: 4 RESULTS
text: In this section, we demonstrate the priority of the proposed approach given two different examples. First, we employ the proposed method to identify the soft bending actuator model and compare the results with a Recursive Least Squares (RLS). In the second example, we evaluate it on the Hunt-Crossley model given reality-based data and compare it with a simple MCMC and RLS.

heading: 4.1 SIMULATION RESULTS
text: For this part, we consider the dynamic model of a fluid soft bending actuator. The dynamic is given by the following relation: (Wang et al., 2019)
α̈ = q1(p− patm)− q2α̇− q3α uc sign(ps − p) √ |ps − p| =q4ṗ+ q5ṗp, ud = 0
ud sign(p− patm) √ |p− patm| = q6ṗ+ q7ṗp, uc = 0,
(15)
where α is the angle of the actuator, uc, ud are the control inputs, and p, ps, patm are the current, compressor and atmosphere pressure respectively. For this example, we assume q1 = 1408.50, q2 = 132.28, q3 = 3319.40 are known and patm = 101.3kpa, ps = 800kpa. We are trying to identify the four other parameters (q4, ..., q7). To this end, we assume the hybrid model below:
u sign(ps − p) √ |ps − p| = θ1ṗ+ θ2ṗp, u = {uc, ud} (16)
As the range of these parameters are small, we scale the input vector by the factor of 107 for RLS. Given the input (uc, ud) and the output (p, ṗ), we want to identify the parameter and estimate the current angle of actuator knowing its initial position at the origin. The data sample time is T = 1 ms and each data pack includes 100 samples which results in an algorithm sample time equal to Ts = 0.1 sec. The point estimation obtained by considering the mode at the modification phase and the medium during the reinforcement phase in ARMCMC is denoted as AR-MAPS. The point estimate results for the parameter estimation are shown in Fig. 3a. The true parameters are q5 = −2.14 × 10−4, q6 = 6.12 × 10−9, q7 = −9.76 × 10−5, q8 = −1.90 × 10−9. The second norm of the parameters’ errors are 0.0235, 6.0053 × 10−7 for θ1, θ2 in RLS and 0.0089, 1.1840 × 10−7 in AR-MAPS, respectively. Moreover, the estimation of the angle is plotted in Fig. 3b.

heading: 4.2 EMPIRICAL RESULTS
text: In this section, we demonstrate ARMCMC by identifying parameters of the Hunt-Crossley model, which represents an environment involving a needle contacting soft material. The needle is mounted as an end-effector on a high-precision translational robot, which switches between two modes: free motion and contact. Due to abrupt changes in the model parameters when the contact is established or lost, online estimation of the force is extremely challenging.

heading: 4.2.1 CONTACT DYNAMIC MODEL
text: Consider the dynamic of contact as described by the Hunt-Crossley model which is more consistent with the physics of contact than classical linear models such as Kelvin-Voigt (Haddadi & HashtrudiZaad, 2012). In order to overcome the shortcomings of linear models, Hunt & Crossley (1975) proposed the following hybrid/nonlinear model :
fe(x(t)) =
{ Kex p(t) +Bex p(t)ẋ(t) x(t) ≥ 0
0 x(t) < 0 , (17)
in which Ke, Bexp denote the nonlinear elastic and viscous force coefficients, respectively. The parameter p is typically between 1 and 2, depending on the material and the geometric properties of contact. Also x(t), ẋ(t), fe are the current position, velocity (as inputs X) and contact force (as output (Y in (1)) of a needle near or inside the soft material, with x ≥ 0 representing the needle being inside. This needle can move freely in open space or penetrate the soft material; the forces on this needle are modeled using the Hunt-Crossley model. The practical problem we consider is to estimate the force at the tip of the needle by identifying the model parameters. Ke, Be, p are three unknown parameters (θ in Eq. (1)) that needs to be estimated. An online estimate of environment force plays a pivotal role in stable interaction between robotic manipulators and unknown environments.
To make the parameters ready for the regression problem we have
log(fe) = log(Kex p s +Beẋsx p s),
log(fe) =p log(xs) + log(Ke +Beẋs). (18)
We will use the RLS method proposed in (Haddadi & Hashtrudi-Zaad, 2012) as a baseline of comparison where it is needed to make the assumption that Be/Keẋs << 1. It should be noticed that the vector of parameters (θ) in the following relation are not independent, which may lead to divergence. With this assumption, we have
log(1 +Be/Keẋs) ≈ Be/Keẋs, log(fe) = p log(xs)+ log(Ke) +Be/Keẋs.
(19)
φ =[1, ẋs, ln(xs)],
θ =[log(Ke), Be/Ke, p] T .
(20)

heading: 4.2.2 SETUP
text: The data structure is same as previous simulation. Prior knowledge of all three parameters (Ke, Be, p) are initialized to N(1, 0.1) (a normal distribution with µ = 1 and σ = 0.1) Moreover, more data is collected, the spread of the posterior pdf decreases. A bit after 5 seconds, the needle goes outside of the soft material, and experiences zero force; this is equivalent to all parameters being set to zero. The color-based visualization of probability distribution over time is used for the three parameters in Fig. 4a. During the period of time that the whole space is blue (zero probability density), there is no contact and the parameter values are equal to zero.
Since we are taking a Bayesian approach, we are able to estimate the entire posterior pdf. However, for the sake of illustration, the point estimates are computed by using AR-MAPS method which is shown in Fig. 4b for the time-varying parameters θ1 = Ke, θ2 = Be, θ3 = p. During the times that RLS results are chattering due to the use of saturation (if not, the results would have diverged), the needle is transitioning from being inside the soft material to the outside or vice versa. In addition, due to the assumption(19), performance can deteriorate even when there is no mode transition. Furthermore, in the RLS approach, estimated parameters suddenly diverged during free motion, since the regression vectors are linear dependent. In contrast, with the Bayesian approach, this issue can be easily resolved. The result of ARMCMC is presented in Fig. 5, which shows the force estimation with two different identification approaches. This probability of interest can be easily obtained by deriving the parameter density at one’s disposal.

heading: 4.2.3 QUANTITATIVE COMPARISON
text: Quantitative details of comparing a naive point estimate of the ARMCMC algorithm by averaging the particles (AR-APS) and the RLS method are listed in Table 1. This reveals more than a 70% improvement in the precision of all model parameters throughout the study by using the Mean Absolute Error (MAE) criteria and also more than a 55% improvements in the second norm of force estimation error. Among parameters, the viscose (Be) has the largest error in the RLS method since it is underestimated due to the restrictive assumption in Eq. (19). The AR-MAPS approach uplifts the performance of the parameter identification and the force estimation.
We also compare ARMCMC to MCMC. For the algorithm to run in real-time, MCMC requires more time to converge. For this example, with λ = 0.7, the value of Kmin is 15000 for MCMC but only 6000 for ARMCMC (more than two times faster) with = 0.01, δ = 0.9. Two approaches that can be used to fix this drawback are to reduce the number of samples, which results in worse precision and reliability compared to ARMCMC, or to increase the algorithm sample time which would cause more delay in the estimation result and slower responses to changes in the parameter.

heading: 5 CONCLUSIONS
text: This paper presented an algorithm for an online identification of full probability distribution of model parameters in a Bayesian paradigm using an adaptive recursive MCMC. Due to the abrupt
change of model parameters such as contact with a soft environment when it is established/lost, conventional approaches suffer from low performance. Empirical results on the Hunt-Crossley model as a nonlinear hybrid dynamic model was compared with a well-known conventional identification process and revealed proficiency of the proposed algorithm. The proposed method provides a systematic strategy for handling abrupt changes which relaxes the pre-requirement conditions in the parameters. As future work, deploying a fully probabilistic framework from identification to control and a decision-making stage is considered to exploit the full potentials of the Bayesian optimization. Additionally, employing a method to compensate the delay will be taken into consideration.

heading: A APPENDIX
text: According to Abolhassani et al. (2007), a nonlinear hybrid model based on a reality-based soft environment is considered as follows:
fe = fst(x, t, tp) + ffr(x, ẋ) + fct(x, t, tp), (21) where x is the needle tip position and tp is the latest time of puncture. Initial position of the environment is assumed to be at the origin. The stiffness of the force (fst) belongs to a pre-puncture and the friction (ffr) and cutting forces (fct) belong to a post-puncture. The stiffness force is modeled using a nonlinear Hunt-Crossley model:
fst(x, t, tp) =  0 x < 0 Kex p(t) 0 ≤ x ≤ x1, t < tp
0 x > x2, t ≥ tp (22)
where Ke, p are the same parameters defined in (17). The maximum depth that the soft environment yields before the puncture and its position after it is denoted by x1, x2, respectively (0 < x2 < x1). In this study, the needle can insert up to 16.65, 10.21 mm before and after penetration. A friction model is inspired from modified Karnopp model.
ffr(x, ẋ) =  Cnsgn(ẋ) +Bex pẋ ẋ ≤ −∆v/2 max(Dn, Fa) −∆v/2 < ẋ ≤ 0 max(Dp, Fa) 0 < ẋ < ∆v/2
Cpsgn(ẋ) +Bex pẋ ẋ ≥ ∆v/2
(23)
where Cn = −11.96 × 10−3 and Cp = 10.57 × 10−3 are negative and positive values of dynamic friction, Dn = −0.01823 and Dp = 0.01845 are negative and positive values of static friction, and Be, p are same as Eq. (17). The relative velocity between the needle and tissue is denoted by ẋ, and ∆v/2 = 0.005 is the value below where the velocity is considered to be zero, and Fa is the sum of non-frictional applied forces. The cutting force is considered as a static force constant for the tissue and the needle geometry if soft tissues are not inhomogeneous and anisotropic Cowan et al. (2011)
fct(x, t, tp) = { 0 x ≤ x1, t < tp 0.94 x > x2, t ≥ tp . (24)
According to the previous relations, the system is considered as a hybrid model while providing both free motion and in-contact environment. The manipulator is a translational mechanism with a friction, slip, and hysteresis loop for the actuator. To present the superiority of the proposed algorithm, the results are compared with the RLS method presented in (Haddadi & Hashtrudi-Zaad, 2012). To prevent the results of RLS from divergence in model mismatch sequences, saturation is applied in the outputs of the identifier.

Review summaries:
The paper introduces a new Markov chain Monte-Carlo (MCMC) algorithm to obtain and track the posterior distribution over unknown parameters in a non-linear system. Despite its simple elegance, i.e., the introduction of a data-driven _temporal forgetting factor_ into the usual Metropolis-Hastings algorithm, the approach is, to my knowledge, novel. Its discovery seems to be the result of the intersection between fields: system identification and Bayesian sampling techniques, leading to new bridges.