Classical robots perturbed by Lévy processes: analysis and Lévy disturbance rejection methods

The stability and convergence of state, disturbance and parametric estimates of a robot have been analyzed using the Lyapunov method in the existing literature. In this paper, we analyze the problem of stochastic stability and also prove some results regarding behavior of statistically averaged Lyapunov energy function in the presence of jerk noise modeled as the sum of independent random variables hitting the robot at Poisson times. This type of noise is also called jerk noise in contrast to white Gaussian noise. Jerk noise is a Lévy process, i.e., a process with stationary independent increments, and is the natural non-Gaussian generalization of white Gaussian noise. Jerk noise can successfully be used to model hand tremor occurring at sporadic time intervals. We also study here the problem of time delay estimation in Lévy noise.


Introduction
A lot of literature exists on analyzing a d-link robot when the external torque applied to it through motors has a noise/disturbance component. Various kinds of disturbance observers have been designed to estimate the disturbance based on measurements of the angular positions and velocities of the links. After designing the disturbance observer using such real-time algorithms, we subtract it from the torque input, thereby reducing the disturbance. Apart from motor noise, disturbance can come from hand tremor or jerks if we are dealing with only a master robot or it can come from jerky environment if we are dealing only with slave robots. Modeling the disturbance as white Gaussian noise processes and estimating the parameters of the robot based on such a Gaussian statistical model have been discussed in the existing literature [1,2].
Some work [1] also exists in proving stability theorems for robots with Gaussian noise using stochastic Lyapunov theory. This analysis is based on calculating the rate of change on the ensemble averaged Lyapunov energy function using Ito's formula for Brownian motion. However, in many instances, it is not possible to model the disturbance as a Gaussian process. This is especially true if the robot links experience sudden jerks of random amplitudes at random times. If the successive jerks are independent and the inter-jerk time is an exponential random variable, then the accumulated effect of the jerks can be successfully modeled as a compound Poisson process or more generally as a Lévy process, i.e., a stochastic process with independent increments. Examples include Gaussian distributed, uniformly distributed and exponentially distributed jump amplitudes.
The Lévy process is Markovian and has independent increments. However, there is another property [3] which states that any Lévy process is a limit of finite superposition of independent Poisson processes and Poisson processes are pure jump processes. When a human operator handles a robot, jerk tremors will be present (especially if the hand has something like Parkinson's disease) and these jerks are pure jumps of varying magnitude which can be modeled using Poisson processes.
Another situation involves modeling the motion of a slave robot in a rough environment which generates jerk noise acting on the slave robot.
In this paper, we study the effect of Lévy processes on a two-link robot manipulator. We simulate Lévy processes and study the solution of the robot differential equation subject to Lévy noise using simulations.
We also calculate the average Lyapunov energy for the trajectory tracking error and also the statistical fluctuations in the Lyapunov energy. These calculations are based on Ito's formula for the Poisson process. In the case of white Gaussian noise, differentials of the tracking error amplitude have negligible third and higher moments owing to Ito's formula, not so for general Lévy processes. Differentials of the error of degree three and higher play a role in the computation of the fluctuations in the Lyapunov energy, and we therefore get a higher value of the energy fluctuations in the Lévy case than in the Gaussian case. This effect can also be explained by noting that an accumulated sequence of finite independent random variables has a large fluctuation as compared to an accumulated sequence of small random variables owing to the applicability of the central limit theorem in the latter case.
In this paper, we also study the maximum likelihood method of parameter estimation of a two-link robot based on measurements of the total number of jerks suffered by the robot angular variables in a fixed time duration. Specifically we show how the times at which the jerks take place and the size of the angular jumps can be estimated. Our formula for the maximum likelihood estimate should be compared to the case of white Gaussian noise.
In "Appendix," we discuss algorithms for estimating parameters in SDEs driven by Lévy processes and general theory for statistical analysis of a Lévy processdriven stochastic differential equation (SDE) including an infinite series formula for the rate of change in the average of any Lyapunov function of a process satisfying a SDE driven by a Lévy process that is representable as a sum of a Poisson number of independent random vectors.
In summary, this paper analyzes statistical properties of classical robots perturbed by Lévy noise and derives parameter estimation algorithms.
As a final application, we address the important problem of rejecting the Lévy disturbance from the system dynamics using techniques of nonlinear filtering theory, i.e., estimate the current state based on past noisy measurements, and then use this state estimate to estimate the current Levy disturbance which is subtracted from the dynamics at the next time instant.

Review of literature
In [4], the authors consider an SDE in R driven by a Lévy process characterized by a Brownian motion process and a Poisson field. The process considered has a special form, namely the coefficients of the Brownian motion and Poisson fields are linear functions of the state. By applying Ito's formula for Brownian and Poisson processes, they arrive at a formula for the rate of change in a Lyapunov function of the error between the state process of the Lévy-driven SDE and a stationary solution of the same. Using this, they are able to characterize the Lyapunov of exponent of the error. Our robot model is a bit more complicated in that the coefficients of the Lévy noise are vector-valued nonlinear functions of the state process. Nevertheless, the results of [4] may be extended and applied to our model by considering a linearized approximation. We have in our paper also made a small attempt to bound the Lyapunov function by an exponentially increasing function of time, and this result may be compared to [4] in the context of Lyapunov exponents.
In [5], an ordinary delay differential equation with time-varying delays has been considered. This dynamics has been stabilized by adding a Lévy noise component with state-dependent coefficient. For this, the authors apply Ito's formula for the Brownian and Poisson fields to compute the average Lyapunov function. They prove that under certain assumptions on the infinitesimal generator of the Lévy process applied to the Lyapunov function, its rate of increase is lower and upper bounded by certain nice functions and then using this result, they deduce asymptotic stability of the stochastic Lévy modified system. To establish this result, they make use of Markovian stop time theory like Doob's optional sampling theorem. This theory can definitely be generalized to multivariate state variable differential systems with delay, and we can hope therefore to apply it to master and slave robots communicating via time-varying teleoperation delays. This will be the subject of a future paper by us.
In [6], stochastic exponential stability of differential equations perturbed by Lévy noise that can be represented as the sum of Brownian motion and a superposition of Poisson processes is investigated. Under local Lipschitz conditions on the coefficients of the SDE, they prove zero crossing theorems as well as Martingale convergence theorems which are interpreted as strong laws of large numbers. They then perturb the system by large jump Poisson processes and prove that the system is exponentially stable, i.e., almost sure asymptotic boundedness of log |x(t)| t . We can apply the results of [6] to our problem, since [6] deals with vector-valued processes driven by vector-valued Poisson processes. In order to do so, we must verify appropriate Lipschitz conditions. Applebaum and Siakalli [6] presents a beautiful computation of the Lyapunov exponent lim sup t→∞ log |x(t)| t , and we need to apply it to the twolink system. Moreover in [6,7], a special form of the coefficient matrices of the Lévy noise has been chosen to prove exponential stability. We need to generalize this for general coefficients appearing in robotics.
In [8], existence and uniqueness theory for general stochastic differential equations driven by a continuous superposition of Poisson fields with the coefficients in the superposition being causal functionals of the state vector which takes values in a separable Hilbert space has been developed. This theory can be applied to the robot problem when the dynamical equations have the form where F(t, u, q, q ) for any u ∈ E is a functional of {q(s), q (s), s ≤ t}, i.e., F is a causal functional of the state vector. Such a model is valid when the intensity of the jerk noise is state dependent and has memory, i.e., the intensity depends on the past orientation and angular velocities of the robot. This may happen when a feedback is used to control the motion of the robot.
In [9], the authors study linear SDE's driven by bivariate Lévy process as generalized Ornstein-Uhlenbeck (OU) processes associated with Brownian motion. If we consider a linearized robot around a given trajectory with Lévy jerks, then the "error process" as discussed in the body of our paper satisfies a similar equation. More generally consider a SDE let dξ 0 (t) = f (ξ 0 (t))dt be the "equilibrium trajectory" and e(t) = ξ(t)−ξ 0 (t) the error process. Then we have approximately, where V t is Lévy process. This can be expressed in the form [9] studies a generalization of this linear SDE in which the differential process is replaced by the differential G t dW t where W t is a Lévy process different from V t . We can thus apply the results of [9] to our linearized robot dynamics. In [10], the heat equation with Lévy noise is studied. Existence results are established generalization including the study of partial differential equations driven by Lévy process, i.e., field u(t, x) satisfying where L 1 and L 2 are partial differential operators in the variable x ∈ R n and w(t, x) is a white noise field, not necessarily Gaussian. Thus is a Lévy noise field, i.e., its increments L(t k+1 − t k . . .), k = 1, 2, . . . are independent fields for t 1 , t 2 , . . .. This idea can be generalized to the situation of a distribution of robots or flexible robots whose angular configuration is specified by q(t, x). There may be an interaction potential energy between two robots at locations x, y or equivalently an energy between the configurations at x and y. The Lagrangian then has the form the result of applying the Euler-Lagrange equation to this action functional is In addition if a Lévy noise field w(t, x) is present, the right-hand side of (10) acquires a random torque term of the general form and the resulting equation can be analyzed for the statistics of fluctuations.
In [11], the authors construct the ML estimator for Gaussian signals in the presence of a random medium which varies impulsively causing the received signal field to have a Lévy distribution of the α-stable type. In effect, the parameters in a signal model having the representation as a superposition of Gaussian and Poisson processes are estimated by the ML method. The model considered is the direction of arrival (DOA) model in which the source signal statistics and the DOA's are to be estimated. We can use the results of [11] in our robot problem by considering α-stable Lévy noise and constructing the ML estimator for the robot parameters, namely the link masses and lengths from measurements on the angular positions. This would be a block processing approach. We could go a step further and consider estimating the robot parameters on a real-time bases from noisy measurements of the robot angular positions or end effector positions using an extended Kalman filter (EKF) developed for α-stable Lévy process white noise and measurement noise.
Asymptotic equivalence of Gaussian white noise processes [12,13] and Lévy processes has been studied in [14]. However, in general the behavior of Lévy process over a short time duration is very different from that of Gaussian noise.
If we assume that the jerks are statistically independent of each other, then these jerks must be modeled as the differentials of a process having independent increments (i.e., non-Gaussian white noise). If we further assume that this jerk noise is a stationary stochastic process, then the process increments must be stationary and we know from the basic Kolmogorov-Levy-Khintchine theorem that any such process can be represented as the limit in probability of a sequence of compound Poisson processes. A delightful proof of this fact is contained in [15]. Thus as an approximation, it is worth modeling a jerk process as a superposition of compound Poisson processes. Some of the useful references containing these ideas are [16][17][18].
The main contributions of our paper are as follows: (a) derivation of a linearized SDE for the tracking error process and parametric estimation error process-driven by a compound Poisson process. (b) Solution to the compound Poisson process-driven SDE using a state variable approach. (c) Theoretical expressions for the average rate of increase/decrease in a quadratic Lyapunov energy function of the tracking and parametric estimation errors with simulation plots. (d) Theoretical expressions for average rate of increase/decrease in arbitrary functions of the tracking and parametric estimation errors that explicitly display the difference between Ito's calculus for Brownian motion and Poisson processes. (e) Path integral expressions for the likelihood function for varying time delay and robot parametric uncertainties. (f) General approximate expressions for the likelihood function of parametric uncertainties and fixed time delay. (g) Verification of the Jerk theory developed using hardware Phantom Omni Bundle as well as MATLAB-based simulational studies.

Problem formulation
We first formulate the SDE for a robot with the perturbation to the input torque being a Lévy process modeled as a compound Poisson process, i.e., as a sum of i.i.d. random variables taken N (t) times where N (.) is a Poisson process. The input torque is the computed torque for a desired trajectory q d (t) with proportional derivative feedback control included for tracking purpose. We also assume parametric uncertainty for the robot. By linearizing around the true parameter values, and assuming the trajectory tracking error e(t) = q d (t) − q(t) to be small, we derive a secondorder linear differential equation for the error process e(t) which contains the parametric estimation error δθ(t). The model for the robot dynamics is whereθ(t) is the estimate of the robot parameters θ 0 based on observation collected up to time t.
is a Poisson process. The above model for the robot dynamics is derived from the Lagrangian with a damping torque N (q,q) and noise torque included in τ . The Euler-Lagrange equations [19] are given by d dt This results in the given differential equation, where M(q) is the mass-moment of inertia matrix and 1 2q T M(q)q is the kinetic energy of the robot. V (q) is the gravitational potential energy. K p and K d are, respectively, position and velocity error gain coefficients. Thus, w(t) is a special kind of Lévy noise differential. Writingθ(t) − θ 0 = δθ(t), we get approximately, from (12) or equivalently, An explicit computation for W is given in [1]. Further, we assume a parameter update equation

Lyapunov energy function
As in the existing literature [1] we define a quadratic Lyapunov function of e(t) and δθ(t) which decreases monotonically to zero as t → ∞ in the noiseless case provided that we assume an appropriate parametric update equation. Like the robot dynamics, we also assume that the parametric update equation also contains a Lévy noise component and then obtain a expression for the average rate of change in the Lyapunov function using Ito's formula for Poisson process. In passing, in order to illustrate the difference between the effects of white Gaussian noise and compound Poisson noise on a robot, we give a sample computation of the average rate of change in an arbitrary function of e(t) and δθ(t). Unlike the white Gaussian noise case where this average rate of change involves only up to two partial derivative of the Lyapunov function, in the compound Poisson case, all partial derivatives come into the picture. Consider now a Lyapunov energy function Equation (16b) can be approximately expressed as , .
This approximation involves neglecting terms of order (q d − q) ⊗ δθ and higher.
Remark 1 If this linearized approximation is not made, then we would get an equation of the form i.e., the coefficients of the Lévy-driven SDE would depend on the stateẽ(t).
To solve this, we would have to write where δ is a small perturbation parameter. We can then expand the solutioñ substituting this expression into the SDE and equating equal powers of δ, gives This gives us an iterative stochastic integral algorithm for calculating the processesẽ k (t). The expressions will be in term of multiple stochastic integrals with respect to the processes The successive approximation e k (t) will have mean varying as the kth moment of the Y k s and variance as the 2kth moment of the Y k s. Hence, our linearized approximation (21) will be good enough if the Y k s have small variances as compared to the initial error energy ||ẽ(0)|| 2 .
Likewise Eq. (18) can be expressed as Then, Ito's formula gives We assume that We note more generally that if V (ẽ, δθ) is any function ofẽ and δθ, then and hence, . Equation (32) shows the explicit dependence of the rate of change in the average Lyapunov function on the Poisson rates λ,λ as well as on the noise amplitude correlations matrices R Y , R Z .

Remark 2
We can use (34) to arrive at an upper bound on EV (ẽ(t), δθ (t)) which is a step toward proving stability.
Suppose that the Lyapunov function V (ẽ, δθ) satisfies the inequalities Suppose in addition, the partial derivatives of V (ẽ, δθ) satisfy the inequalities where it is assumed that So letting We get giving So we can put a upper bound on the rate of growth of V (ẽ(t), δθ (t)). Stability may be hard to prove but we can assert that log E[V (ẽ(t), δθ (t))] grows at most linearly with t.

Parameter estimation
Here, the problem involves computation of the probability density functional of the trajectory error process as a function of the parametric uncertainty. This expression may be used to obtain a maximum likelihood estimate of δθ in the block processing format. The case of block processing parameter estimation based on discrete time measurements is also considered. We measureẽ(t) (i.e., q d (t),q d (t), q(t),q(t) or equivalently e(t),ė(t),ë(t)) and let f δt (w) denote the probability density of w(t)dt, i.e., of Y N (t)+1 dN (t). Then from Eq. (16b), the approximate likelihood function of δθ (assuming that it is a constant uncertainty) is given by where , δ being the time discretization step size and by 0≤t≤T ψ(t) we actually mean Thus, we have an alternate description of the Likelihood function as a path integral: Alternately, if we measureë(t k ), e(t k ),ė(t k ), k = 1, 2, . . . , N , then we can write the likelihood function for δθ in the form Feynman path integrals are used extensively in quantum physics to calculate probabilities of events by solving the Schrodinger equation. In classical probability, they appear in the form of Kac formulas as the solution to the Heat equation with potentials. The maximization of this function of δθ is easily implementable.
N ote : Equations (45) and (46) are meant to express the probability law of the error process in path space aa a path integral. For this, we need the characteristic functional of the compound Poisson process We have dξ(t) = Y N (t)+1 dN (t) and so From which we declare that By infinite-dimensional Fourier inversion of the characteristic functional, we obtain the probability law as a path integral.

Time delay estimation in Lévy noise
Here, the problem considered is approximate time delay estimation from trajectory error measurements. The desired trajectory is read with some delay and the maximum likelihood estimation for the joint time delay and parameter uncertainty is constructed. q d is known except for a time delay, i.e., the reference signal is q d (t − T ) where T = T 0 + δT with T 0 known but δT unknown and small. We write q d0 (t) in place of Then, it follows from Eq. (17) that Sö If the noise were white Gaussian then the maximum likelihood estimator of δT, δθ would be We write δT δθ = δϕ and then setting the derivative of the above w.r.t δθ and δϕ to zero gives In the compound Poisson case of Eq. (53), i.e., when the noise is which is white but non-Gaussian, we have approximately where so that the probability density function of w(t)Δt is given by The probability density functional of {e(t), 0 ≤ t ≤ T } is then approximately given by This is evaluated using a finite pulse approximate for the δ-function and maximized with respect to δT and δθ to obtain their maximum likelihood estimates.
Remark 3 In Sects. 3.2 and 3.3, we have pointed out how parameter uncertainties and time delays can be accounted for and even estimated in this jerk noise model. For example, we could consider an SDE where θ is the parametric uncertainty and τ is the time delay. N k (.), k = 1, 2, . . . , p are independent Poisson processes. X (.) is no longer a Markov process but as in [1], we can transform this SDE in to an infinitedimensional Markovian SDE by introducing a vector wherẽ and likewise forg k . Then, ξ(t) is clearly a Markov process and the parameter θ can be estimated by constructing the likelihood function of ξ(.) given θ . Markovian systems are more easily real-time implementable and the infinitesimal generator of a Markov process can directly be used to construct its probability distribution by solving the Chapman-Kolmogorov forward equation [20][21][22][23].

Simulation and numerical results
The differential equation that describes the tracking error evolution around a fixed point q(t) = q 0 = q d0 ∀t and henceq(t) = 0,q(t) = 0∀t along with the parameter update equation has been carried out in our experiments. We assumed that Y k are constant vectors in R 2 and Z k are constant vectors in R 4 . In other words, the robot process noise is assumed to be a Poisson process and so is the noise process in the parameter update equation. Both the Poisson processes are assumed to be independent of each other. We could also assume that the Y k s are independent and identically distributed (i.i.d.) random vectors and also the Z k s. In that case, our process noise and parameter update noise processes are compound Poisson processes which is a special case of a Lévy process. The state and parameter evolution Eqs. (21) and (66) can be expressed as where The solution to the above state equation is expressed as a Poisson integral: Plots of ξ(t) based on approximating this integral by a discrete stochastic sum are shown in Figs. 1, 2, 3, 4, 5, 6, 7 and 8. The plots show sudden spikes in the angular position and velocity tracking error. These spikes represent abrupt deviation of the error about the stable point q 0 arising from the Poisson jerk noise. Figure 9 shows This energy plot is randomly fluctuating with spikes but the value around which these oscillations take place are roughly the same, namely the expected value of

E(V (t)) = T r(QE(ξ(t)ξ(t) T ))
where It is easily seen that since the eigenvalues of P all have negative real part, V (t) converges as t → ∞ to a constant and this constant represents the limiting average Lyapunov energy arising from random jerks. We also display hardware plots of the actual Phantom Omni Bundle angular position and velocity when the robot is driven by two torques, one a constant dc machine torque that balances the gravitational torque so that the angular positions of the links remain constant, two a jerk torque that comes from the hand operator. This jerk torque is generated by applying jerks at random time intervals. Alternately, if we are dealing with only a slave robot dynamics, the jerk torque comes from the environment in which the slave executes motion. The plots (Figs. 1, 2) of the angular position error q − q d show rapid fluctuations concentrated over small time intervals followed by gaps during which the angular positions are constant. The plots (Figs. 3, 4) of the corresponding angular velocity error show larger spikes in the oscillations followed by gaps. This is as expected since the rate of change in small discontinuities can be very large. There is a strong match between the hardware and software generated plots which confirms that our theory based on the Ito differential rule for Poisson processes is valid in practical situations involving jerk forces applied to dynamical systems. It should be noted that a general nonlinear dynamical system excited by jerk noise can be modeled using the following SDE [24] where N (t, dx) is a space-time Poisson random field with intensity λdtdF(x) and the integral above is taken over x ∈ E. Such a space-time Poisson random field can be generated by taking a sequence of i.i.d. random variables X n , n = 1, 2, . . . having common distribution F(x), a Poisson process N (t) independent of the X n s and defining for each Borel set A, where χ X k ∈A (ω) is the random variable that equals one if X k (ω) ∈ A and zero otherwise. The process X(t) constructed in this way gives jerks with an amplitude dependent on the current time and state and is nonzero only if at time t, the Poisson process N(t) makes a jump such that the amplitude of the that jump is g(t, X (t), X N (t)+1 ). This is the most general form of a SDE driven by a compound Poisson process and the solution is a Markov process.
It is an open problem as to how to simulate such a process using MATLAB. (72) suggests an example for computing approximately the mean, variance and higher moments of the output of an SDE driven by a Poisson field is as follows:

Remark 4 Equation
where EN (dt, dx) = λdtdF(x) and δ is a perturbation parameter introduced to show that the noise component is small. Let Then equating coefficients of δ 0 , δ 1 and δ 2 gives {X 0 (t)} is thus a non-random process,

Lévy disturbance rejection using nonlinear filtering theory
Consider the general Lévy-driven dynamical system as a simplified version of Eq. (74): where N is a Poisson noise with rate λ. We wish to eliminate the Poisson jerk N (t). For that, we must first measure ξ(t) or some output noise function of ξ(t) and then construct a real-time estimatorξ(t) of ξ(t) based on these noisy observations and then estimate the disturbance N (t) usinĝ or equivalently, To getξ(t), we consider a measurement output given by where v(t) is white Gaussian noise, i.e., v(t) = dB(t) dt where B(t) is standard vector-valued Brownian motion. The real-time filter can be derives as follows: For any function φ(ξ ) of ξ(t), define Then of the measurement up to time t. The functions F t (φ) and G t (φ) can be calculated as follows: Let where f (t) = (( f k (t))) d k=1 is a non-random vectorvalued function of the same dimension as z(t). Then C(t) is measurable with respect to z t and hence (using the MMSE property of π t (φ)). Thus or using Ito's formula, and since f (t) is arbitrary, we infer that These imply Here £ t is the generator of the Markov process ξ(t), i.e., These are the basic Kushner equations for dynamically estimating any function of the Markov process ξ(t) and as explained above can be applied to disturbance removal. This idea of constructing a nonlinear dynamical filter can be found in [25]. The quantum non-commutative case discussed in [25] is a generalization of the classical commutative situation which we have discussed here.

Experimentation and result
To validate our software simulation results, we have performed trajectory tracking experiments (without parametric uncertainties) on a Phantom Omni as shown in Fig. 11 by applying a desired machine torque in conjunction with a hand tremor torque and have plotted the angular position as a function of time. Figure 12 shows the schematic diagram of a two-link robot where Jerk acts and where non-random machine torque acts. The   [26]. Phantom Omni is a device which provides kinesthetic active force feedback on the x, y and z axes from the measurement of instantaneous position and velocity by use of joint optical encoders. IEEE.1394 Firewire port is used for the communication interface and allows realtime programming through class OpenHaptics based on C++. In this paper, real-time surface is implemented by MATLAB using Phansim toolbox [27]. Specifically, given a desired trajectory q d (t), we have computed the torque to be given by and have added to this torque a feedback component and a Lévy noise torque component w(t) given via the hand operator. The angular trajectory q(t) of the Omni Bundle satisfies the differential equation The actual parameters of the Phantom Omni are θ 0 = [m 1 , m 2 , l 1 , l 2 ] where m 1 = 0.035 kg, m 2 = 0.1 kg, l 1 = 0.135 m and l 2 = 0.135 m. The angular process q(t) traced out by the Omni Bundle has been measured and plotted along with the tracking error. Remarkable similarity with our software plots is seen which confirms the Poisson nature of the jerk noise. Figure 13 shows plots of the desired trajectory q d for both the links which is a constant and the actual trajectory q which is spiky. The difference process is e = q − q d , and its first component is precisely Fig. 14. This Fig. 14 shows a plot of the hardware robot trajectory tracking error for link one with time when no machine torque has been applied and only tremor torque via the hand operator has been applied in the form of jerks at random time intervals. The resulting trajectory is a continuous sequence of spikes at random intervals. Figure 15 is a plot of the error trajectory for link 2 which is similar to that of link one and a formula for this spike process can be derived by linearizing the robot SDE about the stable operating point as we saw earlier. In fact the error process satisfies a linear second-order differential equation with compound Poisson differential input and the solution for the error process is a stochastic integral with respect to the Poisson process which can be represented as a shot noise process: Again we observe a sequence of spikes clustered at random intervals but now the spikes assume even negative values owing to the compensation, i.e., subtraction of the mean. The clustering can be explained by the fact that in the process of discretization of the dynamical system differential equations for simulation, a finite time step size is used which means that within the time interval of the discretization, the spikes increase and decay. This is in contrast to the hardware simulation where the step size is nearly zero. Figure 2 displays the software simulated angular position trajectory for link two. This is similar to that of link 1. Figure 16 displays the angular trajectory of the first link when a background sinusoidal torque through machines is applied and in addition, a hand tremor torque is applied. The plot is a beautiful display of a sinusoidal process with random spikes. It confirms the linearization approximation that when the noise is small, we can treat it as a first-order perturbation and apply linearized perturbation theory to compute the effect of random tremor on the trajectory. If the noise amplitude is larger, then it can interfere with the background sinusoidal trajectory but we have not considered this case which would require higher-order perturbation theory. The error plot for this case is Fig. 17. It shows that the error consists principally of jerk noise which means that if the noise were absent then the non-random trajectory is accurately tracked. Of course, some small amplitude negative valued spikes are also observed in this plot and that can be considered as errors in the tracking of the non-random component of the trajectory.

Conclusion
In this paper, we have considered a two-link robot with its kinetic energy specified by a quadratic function of its link angular momenta having coefficients dependent on the angle between the two links. The potential energy is gravitational, and apart from this there is a contribution to the Hamiltonian coming from the torque processes applied to the links. We have first set up the Lagrangian for this system and have included a Lévy noise component in the torque in the classical equations of motion. This Lévy noise component is determined by the jerks appearing in the motors and in the hand operator. The torque has been generated from a desired angular trajectory by incorporating parametric uncertainties in its computation. The tracking error process satisfies a linear second-order differential equation with Lévy noise driving input. Further, the parametric error update equation also satisfies a linear SDE with the tracking error appearing linearly in it. We club the trajectory tracking error, the tracking error rate of change and the parametric error into a single eight-component vector, and this satisfies a linear stochastic differential equation with Lévy noise input. The Lévy noise has been generated as a compound Poisson process. Stability of this system has been analyzed using stochastic Lyapunov theory. Plots of the angular position and velocity error, the parametric estimation error as well as the Lyapunov energy and its mean value have been obtained using MATLAB as well as by actually driving a Phantom Omni Bundle with jerks applied by the human operator and good agreement with the theoretical results have been obtained. A crucial ingredient in the theoretical analysis is the Ito formula for a Poisson process, and we have shown how this leads to results that differ drastically from the case when noise is white Gaussian with the standard Ito formula for the Brownian motion process being applicable. Finally, we have illustrated how nonlinear filtering theory in the presence of measurement noise can be used for rejection of the Lévy noise component, i.e., we have addressed the disturbance rejection problem.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: General remarks on Lévy processdriven stochastic differential equations
Here, we explain how using time discretization, we can obtain the likelihood function for a stochastic delay differential equation driven by a Lévy process having a Gaussian component. We also indicate here how to compute the approximate probability density function for a Lévy process that is the sum of a Gaussian component and a stable process assuming that probability density function of both the processes are close to the delta-function. Using these formulas, it is a simple matter to formulate the likelihood function for parameter estimation in stochastic differential equation driven by Gaussian plus stable Lévy process [28,29]. Consider where T (t) is a time-varying delay and θ = θ 0 + δθ The parameter uncertainty δθ as well as the delay fluctuations δT (t) is to be estimated. Here, X (t) is a Lévy process of the form where Y i being iid random vectors. Then, where δ N (t) = N (t + δ) − N (t). Then, If O((δt) 2 ) is neglected, then this approximates to as δt → 0. Where 0 < K < ∞. The exact density of δ X (t) without making any approximations is obtained by Fourier inverting Eq. (105) w.r.t α and is given by This formula is exact if we assume that In general, the characteristic function of and hence the exact density of and f δ, We have approximately taken Equivalently the log likelihood function of (θ, r ) is Writing Ψ δ (x) = log f δ,X (x), this is where Δξ [n] = ξ [n + 1] − ξ [n]. If f δ,X has a white Gaussian component, then and L(θ, r |ξ [n], 0 ≤ n ≤ N ) where e −λ|α| s F −1 −−→ g s (x). Suppose where B(t) is Brownian motion and W (t) is Lévy process with stable law.
This is a convolution of two probability distributions, both of which are close to the δ-function. So If the Gaussian component is strong compared to the non-Gaussian component, then we have approximately An approximation to Eq. (100) can be expressed as and hence if f δ (x) is the density of δ X (t), we get the approximate likelihood function of {δT (t)} 0≤t≤T and δθ as 0≤t≤T maximizing this over {δT (t)|0 ≤ t ≤ T } and δθ, we can obtain their MLE's.

Appendix 2: Other methods for statistical inference for Lévy processes
Here, we explained how statistical inference, i.e., parameter estimation for a stochastic differential equation driven by a compound Poisson process, can be carried out based on measurements of the jump times, jump amplitudes and total number of jumps of the state process in a given time interval [0, T ]. We explain how the expression for the resulting likelihood function can be applied to a two-link robot with Lévy perturbation. We then generalize this expression for the likelihood function to Lévy processes expressible as the sum of a compound Poisson process and Brownian motion. The construction of the likelihood function is based on the fact that in between two successive jumps, the state process executes a continuous trajectory driven by white Gaussian noise for which the likelihood function is well known. Consider a stochastic differential equation of the form dξ(t) = f (t, ξ(t))dt + g(t, ξ(t))dV (t) where ξ(t) is an R n -valued random process and f : R + × R n ] → R n is a smooth map, g : R + × R n → R n×n is a smooth map and V (t) is an R n -valued Lévy process having a compound Poisson representation where N (t) is a Poisson process with rate λ and Y k , k = 1, 2, . . . is a sequence of independent identically distributed random vectors in R n independent of N (.). Let F Y denote the distribution of Y 1 and f Y its density. Then the moment generating function of V (t) is given by For our simulations, we have chosen Y k to be a N (0, R) random vector where R = E(Y 1 Y T 1 ) is an n×n positive definite matrix. Suppose the functions f and g depend on an unknown parameter vector θ to be estimated. This happens for example in robot where θ consists of the links and link masses. To construct the likelihood function for θ based on observations of ξ over the interval [0, T ], we take measurements of (a) the number of jumps in ξ over this interval, (b) the amplitude ξ just before and just after a jump, the jump times and the approximate jump duration dt. Since conditioned on the number of jumps of a Poisson process over an interval [0, T ] the jump times are uniformly distributed over [0, T ] [15], it follows that the likelihood function of θ given the measurements described is of the form p (t 1 , . . . , t n , N (T ) = n, ξ(t k −), ξ(t k +), k = 1, 2, . . . , n|θ ) The equivalent negative log likelihood function of θ is then given by −L(θ |ξ) = − n k=1 log |det (g (t k −, ξ(t k −), θ )) | −1 × f Y g (t k −, ξ(t k −), θ ) −1 × (Δξ(t k ) − f (t k −, ξ(t k −), θ ) dt))) The term involving f and dt can be ignored. If the contribution due to this term is to be taken into account, then we must have a Gaussian component in V (t), i.e., V (t) should be replaced by V (t) + σ B(t) where B(.) is Brownian motion. In this case, we also measure ξ between the successive jumps and then the negative log likelihood function has the form −L(θ|ξ) = − n k=1 log(|det (g(t k −, ξ(t k −), θ))| −1 × f Y (g(t k −, ξ(t k −), θ) −1 (Δξ(t k )) t, ξ(t), θ) T g(t, ξ(t)) −1 dξ(t) For our robot, the stochastic differential model can be taken as where with Y k ∈ R 2 , B(t) ∈ R 2 being standard twodimensional Brownian motion. In this case, ξ(t) = [q(t) T , q (t) T ] T and the matrix g(t, ξ(t), θ ) is singular. In this case, we replace the inverse of the g matrix by its pseudo-inverse.

Remark 5
We note that Eqs. (128) and (129) define the characteristic function of a compound Poisson process.
In [15], the Lévy-Khintchine formula log Ee ιαV (t) = tψ(α), (136a) where ν is a measure having certain kinds of special properties for a Lévy process (infinitely divisible) has been derived using limiting argument applied to the compound Poisson process. Thus, this formula can be derived from Eq. (129).

Appendix 3: Numerical simulations about Lévy processes and jerk noise
The procedure of simulating a Poisson process in MAT-LAB for a given rate λ has been carried out using the following program. Two methods are used, one is based on simulating independent Bernoulli random variables X k so that and we then add these X k 's, to generate the Poisson process. The second is based on simulating independent exponential random variables X 1 , X 2 , . . . with density λe −λx u(x) and noting that the Poisson process N (t) at time t is given by The MATLAB code for the second method (which has been used in our manuscript) is: