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

Nonlinear Dynamics, Mar 2017

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

http://link.springer.com/content/pdf/10.1007%2Fs11071-017-3471-8.pdf

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

Classical robots perturbed by Lévy processes: analysis and Lévy disturbance rejection methods Rohit Singla Harish Parthasarathy Vijyant Agarwal 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. Lévy noise; Stochastic differential equation (SDE); Lyapunov function; Phantom Omni; Time delay estimation; Likelihood function - 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. 2 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 |tx(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 log |tx(t)| , and we need to apply it to the twolitn→k∞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 dq(t ) = q (t )dt, dq (t ) = −M −1(q(t ))N (q(t ), q (t )) u∈E 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 Vt is Lévy process. This can be expressed in the form de(t ) = [9] studies a generalization of this linear SDE in which the differential process + M −1(q(t )) F (t, u, q, q )N (dt, du) = L1u(t, x ) + (L2u(t, x ))w(t, x ) 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 is replaced by the differential Gt d Wt where Wt is a Lévy process different from Vt . 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 L1 and L2 are partial differential operators in the variable x ∈ Rn and w(t, x ) is a white noise field, not necessarily Gaussian. Thus is a Lévy noise field, i.e., its increments L(tk+1 − tk . . .), k = 1, 2, . . . are independent fields for t1, t2, . . .. 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 M (q(t, x )) U (q(t, x ), q(t, y))d3x d3 y the result of applying the Euler–Lagrange equation to this action functional is − 2 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 K (q(t, x ), (t, x ), x , y)w(t, y)d3 y 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–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. 3 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 qd (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 ) = qd (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 . e(t ) = qd (t ) − q(t ), k=1 Yk = YN (t)+1 where {Yk } is an i.i.d. sequence in R2 and {N (t ) : t ≥ 0} is a Poisson process. The above model for the robot dynamics is derived from the Lagrangian L(q, q˙ , t ) = 21 q˙ T M (q)q˙ − V (q) + τ (t )T q with a damping torque N (q, q˙ ) and noise torque included in τ . The Euler–Lagrange equations [19] are given by This results in the given differential equation, where M (q) is the mass-moment of inertia matrix and 1 q T M (q)q˙ is the kinetic energy of the robot. V (q) 2 ˙ is the gravitational potential energy. K pand Kd 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) and q¨ = q¨d − e¨, or equivalently, e¨ + K pe + Kd e˙ = W (t, q(t ), q˙ (t ))δθ (t ) − M (q(t ), θ0)−1w(t ) W (t, q, q˙ ) = −M(q, θ0)−1 ∂ M ∂ N ∂θ (q, θ0) ⊗ q¨d (t ) + ∂θ (q, q˙ , θ0) An explicit computation for W is given in [1]. Further, we assume a parameter update equation where e˜ = v(t ) = Z N˜ (t)+1 with {Zk } iid in R4 and N˜ (.) a Poisson process independent of N (.). 3.1 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 V (t ) = 21 e˜T (t )Q1e˜(t ) + 2 δθ (t )T Q2δθ (t ) 1 Equation (16b) can be approximately expressed as de˜(t ) = − G(t )YN (t)+1d N (t ) 0 . M −1(qd (t )) This approximation involves neglecting terms of order (qd − 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 e˜(t ). To solve this, we would have to write G1(t, e˜(t )) = W0(t ) + k≥1 G2(t, e˜(t )) = G2,0(t ) + k≥1 where δ is a small perturbation parameter. We can then expand the solution m+r1+2r2+···+mrm=k,m≥1,r1,r2,...,rm≥0 × (I ⊗ e˜0(t)⊗r0 ⊗ e˜1(t)⊗r1 ⊗ . . . ⊗ e˜m (t)⊗rm )]dt m≥0,1+r1+2r2+···+mrm=k,r0,r1,...,rm≥0 × e˜0(t)⊗r0 ⊗ e˜1(t)⊗r1 ⊗ . . . ⊗ e˜m (t)⊗rm This gives us an iterative stochastic integral algorithm for calculating the processes e˜k (t ). The expressions will be in term of multiple stochastic integrals with respect to the processes k=1 Yk = The successive approximation ek (t ) will have mean varying as the kth moment of the Yk s and variance as the 2kth moment of the Yk s. Hence, our linearized approximation (21) will be good enough if the Yk s have small variances as compared to the initial error energy ||e˜(0)||2. Likewise Eq. (18) can be expressed as e˜(t ) = k=0 substituting this expression into the SDE and equating equal powers of δ, gives Let A = 1 dV (t ) = e˜TQ1de˜ + 2 de˜T Q1de˜ . Then, Ito’s formula gives = 21 e˜T + δθ T Q2(F e˜ + H δθ )dt + 21 T r Q2 Z N˜ (t)+1 Z NT˜ (t)+1d N˜ (t ) We assume that Q1W˜ + F T Q2 = 0 i.e., F = Q2−1W˜ T Q1. Then, dV (t) = 21 e˜T 1 + 2 T r GT Q1GYN (t)+1YNT(t)+1 dN (t) + 21 δθ T Q2 H + H T Q2 δθ dt + 21 T r Q2 Z N˜ (t)+1 Z NT˜ (t)+1 dN˜ (t) Assuming that EYk = 0, EZk = 0, Ed N (t ) = λdt and Ed N˜ (t ) = λ˜dt , we get with Re˜e˜(t ) = E(e˜(t )e˜(t )T ), Rθθ (t ) = E(δθ (t )δθ (t )T ) where RY = (Y1Y1T ), RZ = (Z1 Z1T ). We note more generally that if V (e˜, δθ ) is any function of e˜ and δθ , then G(t)⊗k YN⊗(kt)+1(−1)k dN (t) Z ⊗N˜ k(t)+1dN˜ (t) = E V (e˜, δθ ) ≤ a0||e˜||2 + b0||δθ ||2 ≤ K0V (e˜, δθ ) Suppose in addition, the partial derivatives of V (e˜, δθ ) satisfy the inequalities ≤ K1 a0||e˜||2 + b0||δθ ||2 1/2 (36a) ≤ K2 a0||e˜||2 + b0||δθ ||2 1/2 1 E 1 E ∂k V ∂e˜k where μk (Y ) = E(Y1⊗k ), μk (Z ) = E(Z1⊗k ). 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 RY , RZ . Remark 2 We can use (34) to arrive at an upper bound on EV (e˜(t ), δθ (t )) which is a step toward proving stability. Suppose that the Lyapunov function V (e˜, δθ ) satisfies the inequalities where it is assumed that ϕ (t ) ≤ a ϕ(t ) + bϕ(t ) Writing √ϕ(t ) = χ (t ) gives ≤ a (EV (e˜, δθ ))1/2 + bE V (e˜, δθ ) (39) e b2t −1 So we can put a upper bound on the rate of growth of V (e˜(t ), δθ (t )). Stability may be hard to prove but we can assert that log E[V (e˜(t ), δθ (t ))] grows at most linearly with t . 3.2 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 e˜(t ) (i.e., qd (t ), q˙d (t ), q(t ), q˙ (t ) or equivalently e(t ), e˙(t ), e¨(t )) and let fδt (w) denote the probability density of w(t )dt , i.e., of YN (t)+1d N (t ). Then from Eq. (16b), the approximate likelihood function of δθ (assuming that it is a constant uncertainty) is given by p(e(t ) : 0 ≤ t ≤ T |δθ ) = 0≤t≤T where δe(t ) = e(t + δ) − e(t ), δe˙(t ) = e˙(t + δ) − e˙(t ), δ being the time discretization step size and by ψ (t ) we actually mean ψ (kδ). Note E[eιαT YN(t)+1dN (t)] = λdt EeιαT Y1 + (1 − λdt ) where ΦY (α) is the characteristic function of Y1. Thus, we have an alternate description of the Likelihood function as a path integral: 0≤t≤T Alternately, if we measure e¨(tk ), e(tk ), e˙(tk ), k = 1, 2, . . . , N , then we can write the likelihood function for δθ in the form p (e(tk ), e˙(tk ), e¨(tk ), 1 ≤ k ≤ N |δθ ) N k=1 k=1 Feynman path integrals are used extensively in quantum physics to calculate probabilities of events by solving the Schrodinger equation. In classical probability, k=1 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 ot e : 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 k=1 We have dξ(t ) = YN (t)+1d N (t ) and so EeιαT dξ(t) = E eιαT Y1 λdt + (1 − λdt ) = 1 + λdt EeιαT Y1 − 1 From which we declare that ⎛ = exp ⎝ λ By infinite-dimensional Fourier inversion of the characteristic functional, we obtain the probability law as a path integral. 3.3 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. qd is known except for a time delay, i.e., the reference signal is qd (t − T ) where T = T0 + δT with T0 known but δT unknown and small. We write qd0(t ) in place of qd (t − T0). Then qd (t − T ) ≈ qd (t − T0) − q˙d (t − T0)δT Then, it follows from Eq. (17) that In the compound Poisson case of Eq. (53), i.e., when the noise is w(t ) = YN (t)+1 which is white but non-Gaussian, we have approximately where ΔN (t ) = N (t + Δt ) − N (t ) and E[eαT w(t)Δt ] = λΔt EeαT Y1 + (1 − λΔt ) so that the probability density function of w(t )Δt is given by fw(t)Δt (W ) = δ(W ) + λΔt ( fy1 (W ) − δ(W ) ... W (t ) = −M (q, θ0)−1 ∂∂Mθ (q, θ0) ⊗ (q¨d0 − q d0δT ) The probability density functional of {e(t ), 0 ≤ t ≤ T } is then approximately given by e¨(t ) + K pe(t ) + Kd e˙(t ) = (W0(t ) + δT W1(t ))δθ d N (t ) − YN (t)+1 dt If the noise were white Gaussian then the maximum likelihood estimator of δT , δθ would be (δTˆ , δθˆ) = arg min δT,δθ ||e¨(t ) + K pe(t ) + Kd e˙(t ) − W0(t )δθ − W1(t )δT δθ ||2dt We write δT δθ = δϕ and then setting the derivative of the above w.r.t δθ and δϕ to zero gives 0T W0(t)T W0(t)dt 0T W0(t)T W1(t)dt −1 0T W1(t)T W0(t)dt 0T W1(t)T W1(t)dt 0T W0(t)T (e¨(t) + K pe(t) + Kd e˙(t))dt 0T W1(t)T (e¨(t) + K pe(t) + Kd e˙(t))dt Δ f (e¨(t ) + K pe(t ) + Kd e˙(t ) t∈[0,T ] 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 k=1 where θ is the parametric uncertainty and τ is the time delay. Nk (.), 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 gk (t, X (t ), X (t − τ ), θ )d Nk (t ) k=1 and likewise for g˜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–23]. 4 Simulation and numerical results The differential equation that describes the tracking error evolution around a fixed point q(t ) = q0 = qd0∀t and hence q˙ (t ) = 0, q¨(t ) = 0∀t along with the parameter update equation e˜(t ) = has been carried out in our experiments. We assumed that Yk are constant vectors in R2 and Zk are constant vectors in R4. 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 Yk s are independent and identically distributed (i.i.d.) random vectors and also the Zk 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 , P = A1 = G1 = ⎝ −K p I2 −Kd I2 , A2 = G2 = 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 q0 arising from the Poisson jerk noise. Figure 9 shows Tracking Error Fig. 1 Tracking error fluctuations in first angular position due to Poisson noise around a fixed position Tracking Error Tracking Error Fig. 3 Tracking error fluctuations in first angular velocity due to Poisson noise around a fixed position and zero velocity Tracking Error Fig. 4 Tracking error fluctuations in second angular velocity due to Poisson noise around a fixed position and zero velocity Parameter Estimation Error Time (seconds) Fig. 5 Parametric error fluctuations in mass of link 1 (m1) estimate due to Poisson noise around a fixed parameter estimate a plot of the Lyapunov energy function as a function of time. This function is defined as 1 1 V (t ) = 2 e˜(t )T Q11e˜(t ) + 2 δθ (t )T Q22δθ (t ) Fig. 2 Tracking error fluctuations in second angular position due to Poisson noise around a fixed position 0.4 0.3 0.2 )g0.1 (K 0 ro-0.1 r rE-0.2 -0.3 -0.4 -0.5 0 ) (m 0 r ro-0.1 r E-0.2 Fig. 6 Parametric error fluctuations in mass of link 2 (m2) estimate due to Poisson noise around a fixed parameter estimate Fig. 9 Instantaneous Lyapunov tracking error cum parameter estimation error energy function Parameter Estimation Error Average Lyapunov Energy Function Parameter Estimation Error Lyapunov Energy Function yg30 r e nE20 y 5 g re 4 n E3 2 1 0 Fig. 7 Parametric error fluctuations in length of link 1 (l1) estimate due to Poisson noise around a fixed parameter estimate Parameter Estimation Error Time (seconds) Fig. 8 Parametric error fluctuations in length of link 2 (l2) estimate due to Poisson noise around a fixed parameter estimate 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 Fig. 10 Average Lyapunov tracking error cum parameter δ estimation error energy function the Lyapunov energy function. Further, the oscillations appear to be bounded which show that the system is in some sense stochastically stable. Figure 10 confirms this fact by showing the expected value of V (t ) against time. This plot is based on the exact formula E(V (t )) = T r (QE(ξ (t )ξ (t )T )) 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 − qd 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] d X (t ) = f (t, X (t ))dt + g(t, X (t ), x )N (dt, dx ) k=1 where N (t, dx ) is a space-time Poisson random field with intensity λdt d F (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 Xn, n = 1, 2, . . . having common distribution F (x ), a Poisson process N (t ) independent of the Xns and defining for each Borel set A, where χXk ∈A(ω) is the random variable that equals one if Xk (ω) ∈ 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. Remark 4 Equation (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: d X (t ) = f (t, X (t ))dt + δ g(t, X (t ), x )N (dt, dx ) where EN (dt, dx ) = λdt d F (x ) and δ is a perturbation parameter introduced to show that the noise component is small. Let X (t ) = X0(t ) + δ X1(t ) + δ2 X2(t ) + · · · . g (t, X0(t ), x )X1(t )N (dt, dx ) (76c) g(t, X0(t ), x )N (dt, dx ) X1(t ) = 0≤s≤t,x E 0≤s≤t 0≤s≤t,x E {X0(t )} is thus a non-random process, where v(t ) is white Gaussian noise, i.e., v(t ) = dBd(tt) where B(t ) is standard vector-valued Brownian motion. The real-time filter can be derives as follows: For any function φ (ξ ) of ξ(t ), define πt (φ) = E(φ (ξ(t ))|z(s), s ≤ t ) dπt (φ) = Ft (φ)dt + Gt (φ)T dz(t ) × g(s2, X0(s2), x2)ds1ds2dF(x1)dF(x2) where Ft (φ), Gt (φ) are measurable with respect to the σ − algeba 5 Lévy disturbance rejection using nonlinear filtering theory Consider the general Lévy-driven dynamical system as a simplified version of Eq. (74): dξ(t ) = f (t, ξ(t ))dt + g(t, ξ(t ))d N (t ) dC (t ) = C (t ) f (t )T dz(t ), C (0) = 1 of the measurement up to time t . The functions Ft (φ) and Gt (φ) can be calculated as follows: Let 0≤s≤t,x E 0≤s≤t,x E 0≤s≤t 0≤s≤t,x E Nˆ (t ) = or equivalently, d Nˆ (t ) = 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 ) using −1 × dξˆ (t ) − f t, ξˆ (t ) dt −1 × dξˆ (t ) − f t, ξˆ (t ) dt To get ξˆ (t ), we consider a measurement output given by where f (t ) = (( fk (t )))kd=1 is a non-random vectorvalued function of the same dimension as z(t ). Then C (t ) is measurable with respect to zt and hence E[(φ (ξ(t )) − πt (φ))C (t )] = 0∀t dE[(φ (ξ(t )) − πt (φ))C (t )] = 0 or using Ito’s formula, E[(dφ (ξ(t )) − dπt (φ))C (t )] + E[(φ (ξ(t )) − πt (φ))dC (t )] + E[(dφ (ξ(t )) − dπt (φ))dC (t )] = 0 and since f (t ) is arbitrary, we infer that E[(dφ (ξ(t )) − dπt (φ))C (t )] = 0, E[(φ (ξ(t )) − πt (φ)) f (t )T dz(t )C (t )] +E[(dφ (ξ(t )) − dπt (φ))( f (t )T dz(t ))C (t )] = 0 (91b) These imply E[dφ (ξ(t )) − dπt (φ)|zt ] = 0, E[(φ (ξ(t )) − πt (φ))dz(t )|zt ] + E[(dφ (ξ(t )) − dπt (φ))dz(t )|zt ] = 0 Equation (92a) gives E[(£t φ)(ξ(t))|zt ] − Ft (φ) − Gt (φ)T E[h(t, ξ(t))|zt ] = 0, (93a) while (92b) gives E[(φ (ξ(t )) − πt (φ))h(t, ξ(t ))|zt ] − σv2Gt (φ) = 0 (93b) πt (£t φ) = Ft (φ) + Gt (φ)T πt (ht ), πt (ht φ) − πt (φ)πt (ht ) = σv2Gt (φ) Here £t is the generator of the Markov process ξ(t ), i.e., φ (ξ )T f (t, ξ ) + λ(φ (ξ + g(t, ξ )) − φ (ξ )) (95) 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. 6 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 Fig. 11 Phantom Omni haptic device robot Fig. 12 Schematic of two-link robot displaying where jerk acts and where non-random torque acts, q1 angle between first link and ground plane, q2 angle between first link and second link experiment is run in real time on MATLAB Simulink with the Phansim toolbox. Phansim toolbox builds up an interface on the top of the OpenHaptics Toolkit. The user can input torque or force to the box and can read the joint angle, gimbal angle and gimbal position [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 qd (t ), we have computed the torque to be given by Fig. 13 Hardware robot plot for both links of classical Phantom angular trajectory with constant desired trajectory and jerk noise given by hand operator and have added to this torque a feedback component Phantom Trajectory Tracking τ f (t ) = M (qd (t )(K p(qd (t )−q(t )+Kd (qd (t )−q (t ))) (97) 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 M (q(t ))q (t ) + N (q(t ), q (t )) = M (qd (t ))(qd (t ) + K p(qd (t ) − q(t )) + Kd (qd (t ) − q (t ))) + N (qd (t ), qd (t )) + w(t ) (98) The actual parameters of the Phantom Omni are θ0 = [m1, m2, l1, l2] where m1 = 0.035 kg, m2 = 0.1 kg, l1 = 0.135 m and l2 = 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 qd for both the links which is a constant and the actual trajectory q which is spiky. The difference process is e = q − qd , 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: h(t , s)d N (s) = )d0.4 a r l(e0.3 g nA0.2 Fig. 14 Hardware robot plot for tracking error of link 1 with constant desired trajectory and jerk noise given by hand operator Time (seconds) Fig. 15 Hardware robot plot for tracking error of link 2 with constant desired trajectory and jerk noise given by hand operator where sk are the arrival times of the Poisson process N (t ) and h(t , s) = h(t − s) is the impulse response of the differential equation. If the pulse h(t ) is sharp, then sharp spikes appear while if this pulse is spread out and smooth, then smoothened spikes appear. Further, Fig. 1 shows a plot of the error trajectory for link one angular Trajectory Tracking ) 0.2 d ra0.1 ( lge 0 n A-0.1 Fig. 16 Hardware robot plot for link 1 of classical Phantom angular trajectory with sinusoidal desired trajectory and jerk noise given by hand operator position under the influence of compensated Poisson noise, i.e., Poisson noise with its mean subtracted out. 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. Fig. 17 Hardware robot plot for tracking error of link 1 with sinusoidal desired trajectory and jerk noise given by hand operator 7 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 dξ(t) = F(t, ξ(t), ξ(t − T )(t), θ )dt + G(t, ξ(t))d X (t) (100) 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 X (t ) = i=1 where Yi being iid random vectors. Then, d X (t ) = YN (t)+1d N (t ) δ X (t ) = X (t + δ) − X (t ) ≡ YN (t)+1δ N (t ) where δ N (t ) = N (t + δ) − N (t ). Then, n=0 If O((δt )2) is neglected, then this approximates to 1 − λδt + λδt ΦY (α) a result derived earlier. O((δt )2) means a term X (δt ) that satisfies 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 n≥0 This formula is exact if we assume that δ X (t ) = YN (t)+1 (N (t + δt ) − N (t )) In general, the characteristic function of i=N (t)+1 E exp {ια(X (t + δ) − X (t ))} n≥0 n≥0 fY∗n(x ) fY∗n(x ) = n≥0 e−λδ(λδ)n (2πnσY2 −21 exp −(x − nμY)2 n! 2nσY2 WehaveapproximatelytakenT(t) = T = rδ,ξ(nδ) = ξ[n] so that N−1 n=0 × (G(nδ,ξ[n])−1(ξ[n + 1] − ξ[n] − δF(nδ,ξ[n],ξ[n − r],θ))) where e−λ|α|s −F−−→1 gs(x). Suppose N−1 L(θ,r|ξ[n],0 ≤ n ≤ N) = log fδ,X(G(nδ,ξ[n])−1 Writing Ψδ(x) = log fδ,X(x), this is N−1 = Ψδ(G(nδ,ξ[n])−1(Δξ[n] where Δξ[n] = ξ[n + 1] − ξ[n]. If fδ,X has a white Gaussian component, then 1 N−1 (Δξ[n] − δF (nδ,ξ[n],ξ[n − r],θ))T × GRGT (nδ,ξ[n]) −1 N−1 n=0 − ηδ (Δξ[n] − δF (nδ,ξ[n],ξ[n − r],θ)) E eιαTYN(t)+1dN(t) ≈ λδΦY(α) + 1 − λδ Suppose ΦY(α) = −|α|s (Stable laws). e−λδ|α|s → δ−1/sgs δ−1/sx where B(t) is Brownian motion and W(t) is Lévy process with stable law. EeιαT X(δ) = exp −2δαT Rα exp(−λδ|α|s This is a convolution of two probability distributions, both of which are close to the δ-function. So F−1EeιαT X(δ) ≈ (2πδ)d/2|R|1/2 exp −21δ xT R−1x 1 If the Gaussian component is strong compared to the non-Gaussiancomponent, thenwehaveapproximately log fδ(x) ≈ −21δ xT R−1x + (2π)d/2δd/2−1s|R|1/2gs δ−s x 1 ≡ −21δ xT R−1x − ηδ(x) An approximation to Eq. (100) can be expressed as dξ(t) = F(t, ξ(t), ξ(t − T0), θ0)dt ∂ F − ∂ξ2 (t, ξ(t), ξ(t − T0), θ0)ξ (t − T0)δT (t)dt ∂ F + ∂θ (t, ξ(t), ξ(t − T0), θ0)δθ dt + G(t, ξ(t))d X (t) 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 |G(t, ξ(t ))|−1 fδ(G(t, ξ(t ))−1(δξ(t ) 0≤t≤T (t, ξ(t ), ξ(t − T0), θ0)ξ (t − T0)δT (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 Rn-valued random process and f : R+ × Rn] → Rn is a smooth map, g : R+ × Rn → Rn×n is a smooth map and V (t ) is an Rn-valued Lévy process having a compound Poisson representation V (t ) = k=1 where N (t ) is a Poisson process with rate λ and Yk , k = 1, 2, . . . is a sequence of independent identically distributed random vectors in Rn independent of N (.). Let FY denote the distribution of Y1 and fY its density. Then the moment generating function of V (t ) is given by exp(αT y) − 1 fY (y)dy For our simulations, we have chosen Yk to be a N (0, R) random vector where R = E(Y1Y1T ) 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 = exp(−λT ) ((λT )n/n! Πkn=1 fY g (tk −, ξ(tk −), θ )−1 × (Δξ(tk ) − f (tk −, ξ(tk −), θ ) dt)) The equivalent negative log likelihood function of θ is then given by log |det (g (tk −, ξ(tk −), θ )) |−1 k=1 × fY g (tk −, ξ(tk −), θ )−1 × (Δξ(tk ) − f (tk −, ξ(tk −), θ ) dt ))) (131) 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 log(|det (g(tk −, ξ(tk −), θ ))|−1 1 + 2 0 (− f (t, ξ(t), θ )T g(t, ξ(t))−1dξ(t) 1 + 2 f (t, ξ(t), θ )Tg(t, ξ(t), θ )−1 f (t, ξ(t), θ )dt) For our robot, the stochastic differential model can be taken as dq(t) = q (t)dt, + M(q(t))−1dV (t) V (t ) = k=1 with Yk ∈ R2, B(t ) ∈ R2 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 ψ (α), 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 MATLAB 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 Xk so that and we then add these Xk ’s, to generate the Poisson process. The second is based on simulating independent exponential random variables X1, X2, . . . with density λe−λx u(x ) and noting that the Poisson process N (t ) at time t is given by N (t ) = max n| Xk ≤ t k=1 The MATLAB code for the second method (which has been used in our manuscript) is: 1. Singla , R. , Agarwal , V. , Parthasarathy , H. : Statistical analysis of tracking and parametric estimation errors in a 2-link robot based on Lyapunov function . Nonlinear Dyn . ( 2015 ). doi:10.1007/s11071- 015 - 2151 -9 2. Rao , C.R.: Linear Statistical Inference and its Applications , 2nd edn. Wiley, New York ( 2008 ). doi:10.1002/ 9780470316436 3. Feller , W.: An Introduction to Probability Theory and its Applications . Vol. 2. Wiley, New York ( 1971 ). ISBN: 978- 0-471-25709-7 4. Xu , Y. , Wang , X.Y. , Zhang , H.Q. , Xu , W.: Stochastic stability for nonlinear systems driven by Lévy noise . Nonlinear Dyn . ( 2012 ). doi:10.1007/s11071- 011 - 0199 -8 5. Liu , D. , Wang , W. , Menaldi , J.L. : Almost sure asymptotic stabilization of differential equations with timer-varying delay by Lévy noise . Nonlinear Dyn . ( 2015 ). doi:10.1007/ s11071- 014 - 1653 -1 6. Applebaum , D. , Siakalli , M. : Stochastic stabilization of dynamical systems using Lévy noise . Stoch. Dyn. ( 2010 ). doi:10.1142/S0219493710003066 7. Applebaum , David: Lévy Processes and Stochastic Calculus , 2nd edn. Cambridge University Press, Cambridge ( 2009 ). ISBN-10: 0521738652 8. Albeverio , S. , Mandrekar , V. , Rudiger , B. : Existence of mild solutions for stochastic differential equations and semilinear equations with non-Gaussian Lévy noise . Stoch. Process. Appl . 119 , 835 - 863 ( 2009 ). doi:10.1016/j.spa. 2008 .03.006 9. Behme , A. , Lindner , A. , Maller , R.: Stationary solutions of the stochastic differential equation d Vt = Vt − dUt + d Lt with Lévy noise . Stoch. Process. Appl . 121 , 91 - 108 ( 2011 ). doi:10.1016/j.spa. 2010 .09.003 10. Mueller , C. : The heat equation with Lévy noise . Stoch. Process. Appl . 74 , 67 - 82 ( 1998 ) 11. Georgiou , P.G. , Kyriakakis , C. : Signal parameter estimation and localization via maximum likelihood using a sensor array in the presence of Lévy noise . Int. Conf. Signal Process . 1 , 382 - 385 ( 2002 ). doi:10.1109/ICOSP. 2002 .1181070 12. Singla , R. , Parthasarathy , H. , Agarwal , V. , Rana , R.: Feedback optimization problem for master-slave teleoperation tracking in the presence of random noise in dynamics and feedback . Nonlinear Dyn . ( 2016 ). doi:10.1007/ s11071- 016 - 2908 -9 13. Singla , R. , Agarwal , V. , Parthasarathy , H. : Nonlinear robot teleoperation with random fluctuations in the feedback controller . In: Proceedings of computer , communication and control IEEE, ( 2015 ). doi:10.1109/IC4.2015.7375514 14. Mariucci , E. : Asymptotic equivalence for pure jump Lévy processes with unknown Lévy density and Gaussian white noise . Stoch. Process. Appl. ( 2015 ). doi:10.1016/j.spa. 2015 . 09.009 15. Feller , W.: An Introduction to Probability Theory and Its Applications , vol. 2 , 2nd edn. Wiley, New York ( 1971 ). ISBN: 978-0-471-25709-7 16. Rosin´ski, J. : Lévy Processes, Simulation of Wiley StatsRef: Statistics Reference Online ( 2014 ) 17. Figueroa-López , J.E. : Jump-diffusion models driven by Lévy processes . Handbook of Computational Finance . Springer, Berlin ( 2012 ). doi:10.1007/978-3-642-17254-0-4 18. Hackmann , D. , Kuznetsov , A. : Approximating Lévy processes with completely monotone jumps . Ann. Appl. Probab. Inst. Math. Stat. ( 2016 ). doi:10.1214/ 14 - AAP1093 19. Goldstein , H. : Classical Mechanics, 3rd edn. Pearson Education , London ( 2011 ). ISBN-10: 8131758915 20. Buckwar , E. : Introduction to the numerical analysis of stochastic delay differential equations . J. Comput. Appl. Math. 125 , 297 - 307 ( 2000 ) 21. Mao , X. , Sabanis , S. : Numerical solutions of stochastic differential delay equations under local Lipschitz condition . J. Comput. Appl. Math. 151 , 215 - 227 ( 2003 ) 22. Prakasa Rao , B.L.S. : Parametric estimation for linear stochastic delay differential equations driven by fractional Brownian motion . Random Oper. Stoch. Equ . 16 ( 1 ), 27 - 38 ( 2008 ). doi:10.1515/ROSE. 2008 .003 23. Reib , M. , Reidle , M. , Van Gaans , O. : Delay differential equations driven by Lévy processes . Stoch. Process. Appl . 116 , 1409 - 1432 ( 2006 ) 24. Pugachev , V.S. , Sinitsyn, I.N. : Stochastic Differential Systems: Analysis and Filtering . Wiley, Chichester ( 1987 ). ISBN 0-471-91243-3 25. Gough, John, Kostler, Claus: Quantum filtering in coherent states . Commun. Stoch. Anal . 4 ( 4 ), 505 - 521 ( 2010 ) 26. Mohammadi , A. , Tavakoli , M. , Jazayeri , A. : PHANSIM: a simulink toolkit for the sensable Phantom haptic devices . In: Proceedings of the 23rd CANCAM , pp. 787 - 790 ( 2016 ). oai:CiteSeerX .psu: 10 .1.1.417.790 27. OpenHaptics Toolkit , The SenSable Technologies Inc., USA, http://www.sensable.com/productsopenhaptics toolkit.htm 28. Zanzotto , P.A. : On solutions of one-dimensional stochastic differential equations driven by stable Lévy motion . Stoch. Process. Appl . 68 , 209 - 228 ( 1997 ). doi:10.1016/ S0304-4149(97)00030- 6 29. Gupta , H.M. , Campanha , J.R. : The gradually truncated Lévy fight: stochastic process for complex systems . Phys. A 275 , 531 - 543 ( 2000 ). doi:10.1016/S0378-4371(99)00367- 2


This is a preview of a remote PDF: http://link.springer.com/content/pdf/10.1007%2Fs11071-017-3471-8.pdf

Classical robots perturbed by Lévy processes: analysis and Lévy disturbance rejection methods, Nonlinear Dynamics, 2017, DOI: 10.1007/s11071-017-3471-8