1. Introduction
The probability density is often used for a random variable to define its characteristic states and other properties. It is always considered that the probabilities always exists, at least in the weak sense. Probability computation is always centered around the mean values, the law of large numbers, and the central limit theorem. To calculate the very small probabilities or the probability of rare events, classical probability theory provides inconclusive outcomes. Therefore, unified theory of large deviations was introduced by Varadhan to calculate the probability of rare events. The theory of large deviations deals with rates at which probabilities of certain rare events decay as a natural parameter in the problem converges to a given values [Reference Varadhan9]. In robotics, trajectory tracking under noisy conditions is a complex problem. It is always difficult to predict the exit of actual trajectory following the desired trajectory under a small noise [Reference Gul, Zergeroglu, Tatlicioglu and Kilinc23]. The probability of this is very small, but increases as the noise input in increased. In this paper, large deviation theory is applied to the robot trajectory tracking problem to evaluate the rates at which the probability of exit of the trajectory tracking error from a threshold tends to zero in the almost zero noise limit. The method involves computing the LDP rate function of a mixture of Gaussian and Poisson process via the Legendre transform of the logarithmic moment generating function (MGF).
1.1. Related background
The prerequisite definitions and formulas are as follows: Consider a sum of random variables having the form
$S_n = \frac{1}{2} \sum_{i=0}^{n}X_i$
. For mutually independent and identically distributed (i.i.d.) random variables, the probability density function (pdf) of
$S_n$
in simple case is given by
$pS_n(s)$
. Thus the joint pdf of
$X_1,X_2,...,X_n$
factorizes as
$p(X_1,X_2,...,X_n) = \prod_{i=1}^{n} p(X_i)$
where
$p(X_i)$
is the pdf of
$X_i$
’s [Reference Touchette10].
-
1. In rough terms, the random variable
$S_n$
or its pdf
$p(S_n)$
satisfies the LDP if the following limit exists: (1)where I(s) is called rate function and
\begin{equation} \underset{n\rightarrow \infty}{\lim}-\frac{1}{n}\ln pS_n(s) = I(s) \end{equation}
$S_n$
denotes ensemble time average
$\frac{\tilde{S}_n}{n}$
.
-
2. The Gärter–Ellis (GE) theorem: It is an indirect method to establish LDPs of certain functions of
$S_n$
even when
$S_n$
is not a sum of i.i.d. random variables. It is based on the existence of following limiting logarithmic MGF: (2)where
\begin{equation} \lambda(k) = \underset{n\rightarrow \infty}{\lim}\frac{1}{n}\ln \mathbb{E}[e^{nkS_n}] \end{equation}
$\lambda(k)$
is known as scaled cumulant generating function,
$S_n$
is an arbitrary random variable and k is the real parameter. The GE theorem states that if
$\lambda(k)$
is differentiable then
$S_n$
satisfies the LDP and the rate function I(s) is given by Legendre–Fenchel transform of
$\lambda(k)$
: (3)where
\begin{equation} I(s)= \underset{k\epsilon\mathbb{R}}{\sup}\{ks-\lambda(k)\} \end{equation}
$\sup$
denotes the supremum.
-
3. Varadhan’s Theorem: This theorem is concerned with evaluation of functional expectation having the following form:
(4)where f is some function of
\begin{equation} W_n = \mathbb{E}\left[e^{nf(S_n)}\right]=\int_{\mathbb{R}} pS_n(s)e^{nf(s)}ds \end{equation}
$S_n$
satisfies an LDP with rate function I(s). We can write
$W_n\approx\int_{\mathbb{R}}e^{[f(s)-I(s)]}ds$
with sub-exponential corrections in n. By defining the following functional:
$\lambda[f]= \underset{n\rightarrow\infty}{\lim} \frac{1}{n}\ln W_n(f)$
using the LDP limit we obtain (5)
\begin{equation} \lambda(f)=\underset{n\epsilon \mathbb{R}}{\sup}\{f(s)-I(s) \} \end{equation}
-
4. The contraction principle: Let
$A_n$
be a families of random variables, which satisfies an LDP with rate function
$I_A(a)$
and
$B_n$
is the another families of random variable as a function of
$A_n$
that is
$B_n =f(A_n)$
. Now to find whether
$B_n$
satisfies an LDP and to find its rate function contraction principle is applied. The
$p_{B_n}(b)$
is thus given by
$\approx\exp(-n\underset{\{a:f(a)=b\}}{\inf}I_A(a))$
. This shows that
$p(B_n)$
satisfies an LDP with rate function
$I_B(b)= \underset{\{a:f(a)=b\}}{\inf}I_A(a)$
. This formula is known as contraction principle. -
5. Consider a stochastic differential equation (SDE) having the form:
(6)where F(t) is the random trajectory and
\begin{equation} \dot F(t)=f(F(t))+\sqrt{\epsilon}\xi(t) \end{equation}
$\xi(t)$
is white Gaussian noise. Now, it is interesting and important to find the pdf of the given random path
${F(t)}^T_{t=0} $
having duration T in the limit, when the noise power
$\epsilon$
becomes zero. The pdf can be defined using path integrals and it is denoted by
$p(\{F(t)\}_T^{t=0})$
. As
$\epsilon \rightarrow 0$
the random path appearing as the solutions of the SDE should converge in probability to the deterministic path
$F_d(t)$
solving the ordinary differential equation (ODE)
$\dot{F}_d(t)=f(F_d(t)),\ F_d(0)=0 $
, where
$\epsilon \rightarrow 0$
is the low noise limit or large deviation limit. The functional LDP characterization of the the probability that the random path
${F(t)}^T_{t=0} $
fluctuations deviates away from the deterministic path
$F_d{}$
(as
$\epsilon \rightarrow0$
) is given by
$p(\{F(t)\}_T^{t=0})\approx \exp\left({\frac{-I[F]}{\epsilon}}\right)$
, where
$I[F] = \int_{0}^{T}[\dot{F}(t)-f(F(t))^2]dt$
is the rate functional and [F] denotes the integer part of F. The minimum of the rate functional in the zero noise limit over a given set determine that trajectory of the system which contributes maximally to the exit probability. The contraction principle can be used to determine the pdf p(F,T) of the state F(T) reached after time T or in other words, the probability of reaching
$F(T) = F$
from
$F(0) = 0$
is provided by the path connecting these two endpoints having the largest probability. The LDP for the path eventually determines that path with a given subset of paths having the largest probability and approximates in the low noise limit the probability that the path falls in a set E by this largest probability (7)where the rate is denoted by
\begin{equation} p_\epsilon(E)\approx \exp\left(-\frac{1}{\epsilon}\underset{\xi\in E}{\inf I(\xi)}\right) = \underset{\xi \in E}{\sup}\exp\left(-\frac{1}{\epsilon}I(\xi)\right) \end{equation}
$\xi$
.
1.2. Related literature review
Trajectory tracking is a classical problem to test and investigate the proposed algorithm. In principle, the actual trajectory is made to follow the desired trajectory by applying constraints. These constraints may be controlled by appropriate control law like PD, PID, etc. The controller can be further divided into continuous and discrete. Although all controllers are strictly discrete as they are implemented on digital System on Chip (SoC) or microcontrollers, but if the time between samples are very small, then theoretically it is considered in continuous domain. Also, the control action can be broadly bifurcated into adaptive and nonadaptive control [Reference Alandoli and Lee22].
Adaptive sliding mode disturbance rejection control for robotic manipulators has been proposed in ref. [Reference Jing, Xu and Niu15]. Singularities are avoided by defining a terminal sliding mode surface. Uncertainty in the robot parameters are clubbed along with the disturbance torque.The trajectory tracking error is required to fall in a bounded interval, which changes with time.
Further, the dynamics of the transformation error is computed using, which the dynamics of a sliding mode surface is expressed as a function of the error. Finally, the dynamics of the sliding mode surface is used to obtain the control law for minimizing the tracking error and disturbance estimation error. The proof of convergence is based on Lyapunov function. This scheme fails in the presence of stochastic noise because the Lyapunov energy function is not bounded anymore with change in growth rate of the noise.
In ref. [Reference Lei and Chen16], an adaptive fault-tolerant control algorithm is designed for a space robot system with uncertain parameters. Using the Euler–Lagrange system, the authors first set up the basic second-order differential equation by taking model uncertainties into account as an additional disturbance effect on the dynamics. The trajectory tracking error is then defined as a linear combination of the error and its time derivative. The control law is obtained, which cause the error to converge to zero using Lyapunov energy iteration. If the controller coefficients suddenly change due to random fluctuations in disturbance above a given threshold, then the Lyapunov cannot guarantee its stability.
The differential equations for robots with unmeasurable angular velocity and multiple disturbances along with equation for trajectory is provided in ref. [Reference Zheng, Wang, Zhang and Wang17]. The control law is based on state feedback control and disturbance feed forward compensation. Finally, the trajectory error is minimized. However the fast convergence and asymptotic stability cannot be guaranteed.
In ref. [Reference Londhe, Mohan, Patre and Waghmare18], fuzzy controller has been used to control a vehicle manipulator system. The basic advantage of a fuzzy controller is that it replaces a nonlinear system by a fuzzy neural network (FNN) system. The FNN learns from iterative training, which is major drawback of the system. Also if rare events occur due to external disturbance then FNN fails.
Further, the LDP-related literature review is presented along with the dynamical systems. The empirical logarithmic MGF of a random walk has computed a function of the realization
$\omega\rightarrow\Omega$
and its convergence properties are used to deduce almost sure large deviation principles [Reference Barral and Loiseau1]. Local fluctuations of stochastic processes about the almost sure behavior can be studied using this method. We could also apply this technique to our robot dynamics
$q(t,\omega)$
by computing convergence properties of the empirical logarithmic MGF
$\Lambda_n^{\omega}(\lambda) = \frac{1}{n}\log(\frac{1}{k(n)}\int_{0}^{k(n)}exp(f(t)^T q(t,\omega)dt))$
to derive properties of local fluctuations of the trajectory around the almost sure trajectory, where
$k(n)\rightarrow \infty$
as
$n\rightarrow \infty$
.
Dynamical system with random perturbation coupled to transmutation process, which makes jumps has been presented in ref. [Reference Befekadu2]. Generalization of LDP results for diffusion exit from a domain has been applied to the diffusion–transmutation processes. Specifically, LDP results for the exit position and exit time have been obtained. We could also generalize our robot model by appending an additional SDE for the robot parameters causing the mode of the robot to make transmutations and then apply the LDP results of ref. [Reference Befekadu2]. By causing the robot parameters to satisfying another SDE driven by say jump processes, we are able to model our uncertain knowledge of the robot parameters.
A system of interacting diffusion (X(t)) is presented, when the dynamics contains a disorder random parameter k coming from the external magnetic filed [Reference Arous and Sortais3]. The diffusion satisfy a langevin SDE containing a potential term and linear interactions with neighbors characteristic of the Ising model. The LDP for the quenched system, that is the probability measure of the diffusions after averaging w.r.t the disorder parameter k is derived. We could also adopt this technique to our robot model by allowing the robot parameters to have a pdf then studying the quenched robot process using LDP methods, by quenching we mean the pdf of the robot trajectory after averaging w.r.t the pdf of the robot parameters.
Discrete time Gaussian process described by difference equation models like Autoregressive (AR) process, Autoregressive Fractionally Integrated Moving Average (ARFIMA) have been considered in ref. [Reference Massah, Nicol and Kantz4] and LDPs for their time averages has been derived sub-exponential for the LDP of long-range correlation and power law decay for intermittent maps have been reported. This time series theory can also be applied to our robot model by first linearizing the dynamics around the desired trajectories and then discretizing the second-order different equation to obtain an AR(2) vector valued time series model.
LDP for graphs has been described in ref. [Reference Varadhan5]. Specifically, a graph is defined by a random matrix f(x,y), where
$f(x,y)=1$
if vertex x connecting vertex y by an edge and 0 otherwise, p is the probability of connecting x and y. Then LDP results for the random matrix
$(f(x,y)), 1\leq y \leq N$
as the size N of the matrix
$\rightarrow 0$
has been surveyed in ref. [Reference Varadhan5]. If we have a system of N robots communicating with each other and p is the probability of the ith robot having a link with the jth robot, we require then to determine the LDP of the system of robot trajectories when
$N\rightarrow \infty$
. By a link, we mean that the ith robot gives feedback to the jth robot via a torque which forces the jth robot to follow the ith robot. In short, we require to formulate an LDP for a dynamical system moving on the vertices’s of a random graph.
A stationary stochastic process can be viewed as a dynamical system in sequence space having a shift-invariant measure. So all theorems for such dynamical systems like the ergodic theorem, LDP theory, etc., can also applied to stationary stochastic processes. In particular, we can use the results of ref. [Reference Young6] in our robot model after transients have died out and the robot trajectory in the absence of external nonrandom torque executes a stationary vector valued random process.
In ref. [Reference Chatterjee and Varadhan7], Varadhan has proved general theorems random walk in a random environment. Here, the random walk transition probability depends upon the chance outcome
$\omega$
in the environment with translations on the random walk lattice acting on this chance outcome in a special way. Large deviation rate functions have been computed for the quenched random walk, that is for the ensemble of probability measures on random walk path space indexed by the chance outcome of the environment. Further, the rate function for the random walk distribution averaged over the environment has also been computed. This can be applied to our own robot problem noting that
$\begin{bmatrix}q(t)&\dot{q}(t)\end{bmatrix}, t\geq 0$
is a Markov process in continuous time on
$\mathbb{R}^4$
, which when discretized in time and space gives a discrete time Markov chain
$\mathbb{Z}^4$
and the random environment
$\epsilon_0$
is taken as the random fluctuation and damping forces or the robot coming from the environment through which it moves. Translating the environment in one direction is equivalent to trajectory the space of the robot in the opposite direction, that is
$p(\tau_x,\omega,z) = p(\omega, z, x)$
where
$x,z = \begin{bmatrix}q\\[3pt] \dot{q}\end{bmatrix}$
and hence the LDP results of ref. [Reference Chatterjee and Varadhan7] are applicable to the robot problem.
In summary, the adaptive control methods heavily use Lyapunov method for stability purpose [Reference Izadbakhsh and Khorashadizadeh24–Reference Garcia, Castillo, Campos and Lozano26]. Although adaptive control with Lyapunov method may provide desired results, but in case of rare event it fails. The failure conditions are discussed in Section 5. Therefore, the failure outside the stability zone is worthy of investigation and LDP-based probability paves the path for investigation as discussed in the literature review. Although the probability to exit the stable zone would be small with the adaptive control.
1.3. Motivation
The LDP provides the probability of rare events. The rare events are always seen as the least unlikely of all the unlikely ways. In general, all the unwanted events like noise, disturbance, etc., are included during the system modeling, but still due to difference between mathematical models and actual measurements some events always occur, which cannot be mathematically modeled. Furthermore, it becomes necessary to calculate the probability of those unaccounted events or rare events. Therefore, in this paper, large deviations theory is applied to the ideal robotic model without any control method to calculate the probability of exit from the angular trajectory with mixture of Gaussian and Poisson noise. For example, to count the number of collisions from two different parallel random walks growth rate, free energy per monomer and mean value in equilibrium is needed. Further, the count can be represented in terms of probability using LDP. Similarly, the same LDP can be applied to find out the probability of exit from angular position data of robot in presence of small noise and further, we can develop various techniques to reduce the probability of exit.
1.4. Problem statement
The trajectory tracking problem is divided into two following parts: (1) Apply and derive LDP-based exit probability in case of rare event. (2) Design controller based on LDP. However, in this paper only the first part is considered.
The classical model [Reference Chen, Ballance, Gawthrop and O’Reilly12,Reference Agarwal and Parthasarathy13] for a robot is given by
\begin{equation}M(q)\ddot{q}+N(q,\dot{q})\dot{q} + G(q) = \tau(t) +d(t) +W(t)\end{equation}
where M(q) is the mass moment of inertia matrix,
$N(q,\dot{q})$
is the Coriolis and centrifugal force, G(q) is the gravity matrix,
$\tau(t)$
is the operator torque, d(t) is the disturbance, and W(t) is the noise. There is a general theorem of Levy–Khintchine [Reference Parthasarathy14], which states that any independent increment process (i.e., a process whose derivation is white non-Gaussian) can be represented as a continuous process. Further, any process in addition to the torque in Eq. (8) can be represented with any degree of accuracy as a superposition of Gaussian and Poisson process. The robot dynamical equation after neglecting the disturbance and in the presence of low-amplitude tremor noise comprising both a white Gaussian component and Poisson component is given by
\begin{equation}M(q(t))\ddot{q}(t) + N(q(t),\dot{q}(t)) =\tau(t) + (\sqrt{\epsilon}w_1(t)+ \epsilon w_2(t,\epsilon))\end{equation}
where
$ w_1(t) = \sigma \frac{dB(t)}{dt} $
and
$ w_2(t) = \sum_{k=1}^{p}c_k (\dot N_k(t)-\lambda_k(\epsilon)) $
with B(t) being standard d-dimensional Brownian motion and
$N_1(t), ...,N_2(t)$
independent Poisson processes with
$\epsilon$
dependent rates
$\lambda_{1}({\epsilon}),$
$\lambda_{2}({\epsilon}),...,\lambda_{p}({\epsilon})$
. The M(q(t)) is the mass moment of inertia matrix,
$ N(q(t),\dot{q}(t))$
is the Coriolis and centrifugal matrix added with gravity matrix that is
$N(q(t),\dot{q}(t))\dot{q}(t) + G(q)$
, and
$\tau(t)$
is the input torque. Further, in order to find the exit probability, few milestones given below are needed to achieved
-
1. Calculation of the logarithmic moment generating functional of the noise
$\sqrt{\epsilon}w_1 + \epsilon w_2(t, \epsilon),\ 0\leq t\leq T$
and then apply the contraction mapping principle to obtain the LDP rate functional of
$q(t),\ 0\leq t \leq T$
. -
2. Using the rate functional and LDP theory, calculation of the probability of exit from the stability zone;
$P\{\underset{0\leq t\leq T}{\max}||q(t)-q_d(t)||\geq\delta\}$
for small
$\epsilon$
where
$q_d(t)$
satisfies the noiseless
$\epsilon=0$
robot equation. -
3. Calculating approximate rate functional for the noise
$\sqrt{\epsilon}w_1 + \epsilon w_2(t, \epsilon),\ 0\leq t\leq,\ \epsilon \rightarrow 0$
in a quadratic from and also linearize the robot dynamics around the noiseless trajectory
$q_d$
and hence obtain approximate closed form formula from the probability of exit from the stability zone as
$\epsilon \rightarrow 0$
. -
4. Finally, computation of the error energy partition function over the time duration [0,T] using Varadhan’s variational theorem in large deviation theory.
1.5. Research highlights
In the proposed research work, the authors have investigated and developed a state-of-the-art technique using large deviation theory to accommodate the rare events. The noise crossing the desired variable boundary at any time t is considered as an rare event. The probability of occurrence of the rare event is derived for the robotic system. Further, the probability of occurrence of the rare event is used to estimate the probability of crossing the boundaries around the desired robotic angular trajectory. The following are the major highlights of the proposed research work.
-
1. The LDP is applied to the SDE that determine the dynamics of a robot with small-amplitude tremor torque comprising a white Gaussian and Poisson noise.
-
2. Using formula for the LDP rate function of Brownian motion and Poisson process and applying the contraction mapping method, the LDP rate function of the robot angular position process over a finite time duration [0,T] is determined.
-
3. Simplified approximate formulas for the rate function of the trajectories tracking error for the robot angular position process are obtained based on linearization of a dynamical system. The formulas for the logarithmic MGF of a mixture of Brownian and Poisson process and approximate calculation of the corresponding Legendre transform of the MGF are derived.
-
4. Varadhan’s integral variation formula for the low-temperature limit of the MGF of a functional family of random process is applied to calculate the low-temperature limit of the partition function of the average tracking error energy of the robot over a finite time duration [0,T].
-
5. The optimal trajectory error that maximizes the exponential integral (E-I) functional in Varadhan’s variation formula is shown to satisfy a fourth-order linear differential equation and hence the corresponding maximum values of E-I is obtained, which gives us the trajectory error partition function.
-
6. The superiority of the proposed LDP-based approach over the Lyapunov-based method is presented.
1.6. Organization of paper
Afterward, the paper is organized as follows: The application of LDP on d-link robot dynamics is presented in problem formulation Section 2. The calculation of rate function and calculation of LDP probabilities are derived in Sections 2.1 and 2.2. The error energy partition function is presented in Section 2.3. Stability analysis is presented in Section 3. The software–hardware validations and results are provided in Section 4. Further, the comparison between Lyapunov-based approach and LDP-based approach with respect to stability and computational complexity is presented in Sections 5.1 and 5.2 respectively. The extension of purposed approach to other robots is presented in Section 5.2. Furthermore, the conclusion is given in Section 6.
2. Problem Formulation
The dynamics of a d-link robot in the presence of white Gaussian noise and Poisson process is described by a system of coupled nonlinear SDE, where the noise amplitude is small. The large deviation theory is applied to calculate the probability of exact exit of the robot from a stability domain. Further, the approximate probability of trajectory exit from the stability domain within a finite interval of time [0,T] is calculated using the linearization process. These probabilities involve, computation of the LDP rate function of the robot state process. Furthermore, the approximate evaluation of the MGF of the average robot energy over a time interval [0,T] using Varadhan’s integral lemma combined with the variational calculus is presented. The block diagram for the problem formulation is provided in Fig. 1.
Block diagram showing problem formulation.
In the proposed work, we perform these computations in the general case where the noise has a white Gaussian component with Poisson component. The Poisson component is the correct model for jerk torque [Reference Rana, Gaur, Agarwal and Parthasarathy11]. Consider a d-link robot so that the angular position vector of its kind
$q(t) \in \mathbb{R}^d$
and angular velocity
$\dot{q}(t) \in \mathbb{R}^d$
satisfy following SDE:
\begin{align}\frac{dq(t)}{dt} = \dot{q}(t) \end{align}
\begin{align} M(q(t))\ddot{q}(t) + N(q(t),\dot{q}(t)) =\tau(t) + \left(\sqrt{\epsilon}w_1(t)+ \epsilon w_2(t,\epsilon)\right)\end{align}
or equivalently
\begin{equation}M(q(t))d\ddot{q}(t) + N(q(t),\dot{q}(t))dt =\tau(t)dt + \left(\sqrt{\epsilon}w_1(t)+ \epsilon w_2(t,\epsilon)\right)dt\end{equation}
where
$ w_1(t) = \sigma \frac{dB(t)}{dt} $
and
$ w_2(t) = \sum_{k=1}^{p}c_k (\dot N_k(t)-\lambda_k(\epsilon)),\ c_k\in \mathbb{R}^d$
with B(t) being standard d-dimensional Brownian motion and
$N_1(t), ...,N_2(t)$
independent Poisson processes with
$\epsilon$
dependent rates
$\lambda_{1}({\epsilon}), \lambda_{2}({\epsilon}),...,\lambda_{p}({\epsilon})$
. As
$\epsilon\rightarrow 0$
, the system becomes deterministic and we assume a stable solution
$q_0(t)$
to it, that is
\begin{equation}M(q_0(t))\ddot{{q}}_0(t) +N(q_0(t),\dot{q}_0(t))=\tau(t)\end{equation}
where
$\tau(t)$
is the sum of nonrandom machine torque and human operator torque, while
$\sqrt{\epsilon}w(t)$
is the motor noise with human hand tremor torque. The aim is to determine the any asymptotic probability of rare event
$\Big\{\underset{0\leq t\leq T}{\arg\max}||q(t)-q_0(t)||>\delta\Big\}$
as
$\epsilon \rightarrow 0$
using the LDP rate function of the
$w(t,\epsilon) =\sqrt{\epsilon}w_1(t)+ \epsilon w_2(t,\epsilon)$
. We shall use the notation
$w_2(t)$
interchangeably with
$w_2(t,\epsilon)$
. It should be noted that only with the
$\epsilon$
scaling defined above for the Gaussian and Poisson component
$\lambda_k(\epsilon)\epsilon \rightarrow \lambda_{k0}$
, we get meaningful LDP rate function as
$\epsilon \rightarrow 0$
. The limiting scaled logarithmic MGF
$\Lambda(f)$
of the noise
$\sqrt{\epsilon}w_1(t) + \epsilon w_2(t,\epsilon)$
is given by
\begin{align}\Lambda(f)&= \epsilon \log\mathbb{E}\Big(\exp\Big(\frac{1}{\epsilon}\int_{0}^{T}f(t)^T(\sqrt{\epsilon}w_1(t))+\epsilon w_2(t)dt\Big) \Big)\end{align}
\begin{align}&= \underset{\epsilon\rightarrow 0}{\lim}\epsilon \log\mathbb{E}\left[\exp\left(\frac{1}{\sqrt\epsilon}\int_{0}^{T}f(t)^T w_1(t)dt\right)+\exp\left(\int_{0}^{T}f(t)^T w_2(t)dt\right) \right]\end{align}
where f(t) is a parametric functional which indexes the MGF. It facilitates the calculation of the moments. Now, we note that
$\int_{0}^{T}f(t)^Tw_2(t)dt = \sum_{k}\int f(t)^T c_k(dN_k(t)-\lambda_k(\epsilon)dt)$
and
\begin{align} & \mathbb{E}\left[\exp\left(\frac{1}{\sqrt\epsilon}\int_{0}^{T}f(t)^TdB(t)+\sum_{k}\int_{0}^{T}(f(t)^T c_k)(dN_k((\lambda_k(\epsilon) dt))\right)\right]\nonumber\\[5pt]&\quad = \exp\left(\frac{1}{2}\int_{0}^{t}||f(t)||^2dt\right)\exp\left(\sum_{k}\lambda_k(\epsilon)\int_{0}^{T}\left(\exp(f(t)^T c_k)-1-f(t)^Tc_k\right)dt\right) \end{align}
So,
\begin{equation}\Lambda(f)= \underset{\epsilon\rightarrow 0}{\lim}\left\{\epsilon\left[\frac{1}{2\epsilon}\int_{0}^{T}||f(t)||^2+\sum_{k}\lambda_k(\epsilon)\int_{0}^{T}\left(\exp(c_k^Tf(t))-1-c_k^Tf(t)\right)dt\right]\right\}\end{equation}
assuming that
$\epsilon\lambda_k(\epsilon)\underset{\epsilon\rightarrow0}{\rightarrow}\lambda_{ko}$
we get
\begin{equation}\Lambda(f)=\left[\frac{1}{2}\int_{0}^{T}||f(t)||^2+\sum_{k}\lambda_{k0}\left(\int_{0}^{T}\exp(c_k^Tf(t))-1-c_k^Tf(t)\right)dt\right]\end{equation}
If a family of random processes satisfies an LDP, then a functional of this family will also satisfy an LDP with the rate function being determined by applying the contraction principle to the rate function of the original family [Reference Dembo and Zeitouni8] (see item 4 Section 1.1). Using this idea, we formulate an expression for the probability of the actual robot process to deviate from the desired noiseless process. The rate function for the noise is then Legendre transform of
$\Lambda(f)$
(see the Gärtner Ellis theorem in item 2 Section 1.1 Eq. (3) and [Reference Dembo and Zeitouni8]) and is calculated as follows:
\begin{align}I(\xi) &= \underset{f}{\sup}\left[\int_{0}^{T}f(t)^T\xi'(t)dt - \Lambda(f)\right]\nonumber\\[6pt]& = \underset{f}{\sup}\left[\int_{0}^{T}f(t)^T\xi'(t)dt-\frac{1}{2} \int_{0}^{T}||f(t)||^2 dt-\sum_{k=1}^{p}\lambda_{k0}\int_{0}^{T}\left(\exp\left(c_k^Tf(t)\right)-1-c_k^Tf(t)\right)dt\right] \end{align}
where
$\xi$
is the noise amplitude and
$I(\xi)$
is the rate function. Further, by the contraction principle, the exit probability is expressed as follows:
\begin{align}P \left\{\underset{0\leq t\leq T}{\max}||q(t)-q_0(t)||>\delta\right\}& \approx \exp\left(-\frac{1}{\epsilon}\inf \left\{I(\xi)\right\} | \xi(t)\right)\nonumber\\[3pt]& = M(q(t))\ddot{q}(t) + N(q(t),\dot{q}(t))-\tau(t),\ 0\leq t\leq T\,\ q\ \epsilon\ B_{q0}(\delta)^c \end{align}
where
$ B_{q0}(\delta) = \left\{\underset{0\leq t\leq T}{\max}||q(t)-q_d(t)||\leq\delta\right\}$
and c denotes the complement. Furthermore, the evaluation of the rate function in Eq. (19) is provided in Section 2.1.
2.1.
Evaluation of I(
$(\xi)$
)
Calculation of the rate function of a mixture of Gaussian and Poisson noise involves solving a highly nonlinear optimization problem. It is proved that the optimization hyperplane for the Legendre transform of the logarithmic MGF actually yields a global maxima. It is evaluated by the derivative involved in the Legendre transform and proving it to be negative definite everywhere. The rate function
$I(\xi)$
is calculated as follows:
\begin{align}I(\xi) &= \underset{f}{\sup}\left\{\int_{0}^{T}f(t)^T\xi'(t)dt-\Lambda(f)\right\}\end{align}
\begin{align}&= \underset{f}{\sup}\left[\int_{0}^{T}f(t)^T\xi'(t)dt -\frac{1}{2}\int_{0}^{T}f(t)^Tf(t)dt-\int_{0}^{T}\sum_{k=1}^{p}\lambda_{k0}\left(e^{c_k^Tf(t)}-1-c_k^Tf(t)\right)dt\right]\end{align}
\begin{align}&= \underset{f}{\sup}\left\{X(f)\right\}\end{align}
Setting the variational derivative w.r.t
$f(.)$
to zero that is
$\frac{\delta X(f)}{\dot\delta f(t)} = 0$
gives
\begin{equation}\xi'(t) -f(t) - \sum_{k=1}^{p}\lambda_{k0}c_k \left(e^{c_k^T f(t)}-1\right) =0\end{equation}
since
\begin{equation}\frac{\delta^2 X(f)}{\delta f(t)\delta f(s)} = \delta(t-s)\left\{-I_d - \sum_{k=1}^{p}\lambda_{k0}c_{k0}c_k^T e^{c_kf(t)} \right\}\end{equation}
is everywhere is negative definite, we are guaranteed that X will attain a global maximum, when f(t) satisfies (24). Further, high computational efforts are needed to calculate the exit probability from this method. Therefore, in Section 2.2 approximate calculation of exit probability is provided to reduce the computational effort.
2.2. Approximate calculation of LDP probabilities through linearization
Approximate evaluation of the rate functional of the trajectory tracking error process involves linearization of the dynamical system. The solution of the linearized system expresses the trajectory error process as a linear functional of the Gaussian and Poisson noise. The rate functional of the error process is directly calculated using the MGFs of Brownian motion and Poisson processes. Afterward, the rate function of the error process is approximately calculated by the logarithmic MGF of the Poisson component by a quadratic functional. The resulting approximate rate functional is then evaluated by quadratic optimization which yields a set of linear equations that are easily solved. Using the rate function, the deviation probability is evaluated using the KL spectral expansion of positive definite kernels. The approximate probability is calculated as follows: Consider Eq. (9)
\begin{equation}M(q(t))\ddot{q}(t) + N(q(t),\dot{q}(t)) =\tau(t) + \sqrt{\epsilon}w_1(t)+ \epsilon w_2(t,\epsilon)\end{equation}
let
$q{(t)} = q_0{(t)}+\delta q(t)$
, where
$q_0(t)$
is the noiseless trajectory. We need to evaluate
$P_r\{\underset{0\leq t\leq T}{\max}||\delta q(t)||>\delta\}$
, thus linearizing (26) gives
\begin{align}M(q_0(t))\delta\ddot{q}(t)& +\sum_{k=1}^{d}\frac{\partial M}{\partial q_k}\left(q_0(t)\delta q_k(t)\ddot{{q}}_o(t)\right)+\sum_{k}\frac{\partial N(q_0,\dot{q}_0)}{\partial q_k}\delta q_k\nonumber\\[5pt]& +\frac{\partial N}{\partial \dot{q}_k}(q_0,\dot{q}_k)\delta\dot{q}_k(t) = \sqrt{\epsilon}w_1+\epsilon w_2\end{align}
for
$q_0(t) = q_d(t)$
(i.e., the noiseless trajectory is the desired trajectory) the above equation can be expressed as
\begin{align}M(q_d(t))\delta\ddot{q}(t) & +\sum_{k=1}^{d}\frac{\partial M}{\partial q}(q_d(t)\delta(I_d\otimes\ddot{{q}}_d)\delta q(t) + \frac{\partial N}{\partial q}(q_d(t),\dot{q}_d(t))\delta q(t)\nonumber\\[5pt]& +\frac{\partial N}{\partial \dot{q}}(q_d(t),\dot{q}_d(t))\delta\dot{q}_d(t) = \sqrt{\epsilon}w_1(t)+\epsilon w_2(t,\epsilon)\end{align}
or
\begin{equation}\delta\ddot{q}(t)+F_1(t) \delta\dot{q}(t)+F_2(t) \delta q(t) = G(t) (\sqrt{\epsilon}w_1(t)+\epsilon w_2(t,\epsilon))\end{equation}
2.2.1. Evaluation of rate function with Gaussian Noise
If noise is purely Gaussian, that is
$w_2=0$
, then
$\delta q(t)$
is also a Gaussian process. Here,
$G(t)= M(q_d(t))^{-1}$
,
$F_1(t)=M(q_d(t))^{-1}\frac{\partial N}{\partial \dot{q}}(q_d(t),\dot{q}(t))$
, and
$F_2(t)=M(q_d(t))^{-1}\frac{\partial M}{\partial q}(q_d(t))(I_d\otimes\ddot{q}(t))$
so
\begin{equation}\frac{d}{dt}\begin{bmatrix}\delta q(t) \\[4pt] \delta \dot q(t)\end{bmatrix} = \begin{bmatrix}0 & \quad I_d\\[4pt]-F_2(t) & \quad -F_1(t)\end{bmatrix} \begin{bmatrix}\delta q(t) \\[4pt]\delta \dot{q} (t)\end{bmatrix} + \sqrt{\epsilon}\begin{bmatrix}0\\[4pt]G(t)\end{bmatrix}\sigma\frac{dB(t)}{dt}\end{equation}
now the solution to the above Eq. (30) is expressed as
\begin{equation}\begin{bmatrix}\delta q(t)\\[4pt]\delta \dot{q}(t)\end{bmatrix} = \sigma\sqrt{\epsilon}\int_{0}^{t}\Phi(t-s)\begin{bmatrix}0\\[4pt]G(s)\end{bmatrix}dB(s)\end{equation}
where
$\frac{\partial \Phi(t,s)}{\partial t} = \begin{bmatrix}0& \quad I_d\\[4pt] -F_2(t) & -F_1(t)\end{bmatrix}\Phi(t,s)$
for
$t> s$
and
$\Phi(s,s) = I_{2d}$
is the state transition matrix with initial condition when
$t=s$
. Further,
$\left\{\begin{bmatrix}\delta q(t)\\[4pt]\delta \dot{q}(t)\end{bmatrix} : t\geq 0\right\}$
is a Gaussian process with
$\delta q(t) = \sigma\sqrt{\epsilon}\int_{0}^{t}\Phi_{12}(t,s)G(s)dB(s)$
. Here we partition the
$2d\times 2d$
state transition matrix
$\Phi(t,s)$
into four
$d\times d$
blocks
$\Phi(t,s) = \begin{bmatrix}\Phi_{11}(t,s) & \quad\Phi_{12}(t,s)\\[4pt]\Phi_{21}(t,s) & \quad \Phi_{22}(t,s)\end{bmatrix}$
. Further, MGF function is expressed as follows:
\begin{align}\Lambda(f) &= \mathbb{E} \left\{\exp\left(\int_{0}^{T}f(t)^T\delta q(t)dt\right)\right\}\end{align}
\begin{align}&= \mathbb{E} \left\{\exp \sqrt{\epsilon}\int_{0\leq s\leq t\leq T}f(t)^T\Phi_{12}(t,s) G(s) dB(s)\right\}\end{align}
\begin{align}&= \exp\left(\frac{\sigma^2\epsilon}{2}\int_{0}^{T}\int\Big\|\int_{s}^{T}\Phi_{12}(t,s)^TG(s)f(t)dt\Big\|^2ds\right)\end{align}
so the rate function is
\begin{align}I(\xi) = \underset{f}{\sup}\Bigg\{\int_{0}^{T}f(t)^T\xi(t)dt-\frac{\sigma^2}{2}\underset{{0\leq s\leq t_1,t_2\leq T}}{\int} f(t_1)^T \Phi_{12}(t_1,s)G(s)G(s)^T f(t_2)^T \Phi_{12}(t_2,s)^Tdt_1dt_2ds\Bigg\}\end{align}
Further,
$\left\{\sigma^2\int^T_{\max(t_1,t_2)}\Phi_{12}(t_1,s)G(s)G(s)^T\Phi_{12}^T(t_2,s)ds\right\}$
in Eq. (35) denotes the autocorrelation of
$\xi(t)$
process, let it be denoted by
$R(t_1,t_2)$
. Afterward, the rate function is expressed as
\begin{align}I(\xi) &= \underset{f}{\sup}\left\{\int_{0}^{T}f(t)^T\xi(t)dt-\frac{1}{2}\int_{0}^{T}\int_{0}^{T}f(t_1)^T R(t_1,t_2) f(t_2)dt_1 dt_2\right\}\end{align}
\begin{align}&=\frac{1}{2}\int_{0}^{T}\int_{0}^{T}\xi(t_1)^TR^{-1}(t_1,t_2)\xi(t_2)dt_1 dt_2\end{align}
where
$\int_{0}^{T}R^{-1}(t_1,t_2)R(t_2,t_3)dt_2 = \delta(t_1-t_3)$
. Therefore, we have the approximate result that
$P_r\left\{\underset{0\leq t\leq T}{\max}||\delta q(t)||>\delta\right\} \approx \exp\left(\frac{-I_{\delta}}{\epsilon}\right)$
, where
\begin{equation}I_{\delta} = \underset{\underset{0\leq t\leq T}{\max||\xi(t)||>\delta}}{\inf} \frac{1}{2} \int_{0}^{t}\int_{0}^{t}\xi(t_1)^T R^{-1}(t_1,t_2)\xi(t_2)dt_1dt_2\end{equation}
Let
$R(t_1,t_2) = \left\{\sum_{k=1}^{\infty}\lambda_k\phi_k(t_1)\phi_k(t_2)^T ,\ 0\leq t_1,t_2\leq T \right\}$
be the KL spectral representation of
$R(t_1,t_2)$
. Therefore, the inverse of
$R^{-1}(t_1,t_2)$
is expressed as
$R^{-1}(t_1,t_2) = \sum_{k}\lambda_k^{-1}\phi_k(t_1)\phi_k(t_2)$
and the rate function is expressed as
\begin{equation}I(\xi) =-\frac{1}{2} \int_{0}^{T}\int_{0}^{T} \xi(t_1)^TR^{-1}(t_1,t_2)\xi(t_2)= \frac{1}{2}\sum_{k}\lambda_k^{-1}(<\xi,\phi_k>)^2\end{equation}
where
$<\xi,\phi_k> = \int_{0}^{T}\xi(t)^T\phi_k(t)dt$
and
$<\phi_k,\phi_j> = \delta_{k j}$
and
$ \sum_{k}|<\xi,\phi_k>|^2 = ||\xi||^2$
(attained when
$\xi = \delta\phi_{\max} = \delta \phi_{k0}$
). Thus, if
$\lambda_{\max} = \underset{k}{\max}\{\lambda_k\} = \lambda_{k0}$
say, then
$I_\delta = \frac{\delta^2}{2\lambda_{\max}}$
. Finally, the exit probability is expressed as
\begin{equation}P_r\left\{\underset{0\leq t\leq T}{\max} ||\delta q(t)||>\delta\right\} \approx \exp\left(\frac{-\delta^2}{2\epsilon\lambda_{\max}}\right)\end{equation}
Remark: If the unstable zone is defined in the
$L^2$
sense by
$\int_{0}^{T}||\delta(t)||^2 dt >\delta^2$
then its LDP probability evaluates to
$P_r\left\{\int_{0}^{T}||\delta q(t)||^2dt>\delta^2 \right\}\approx \exp(\frac{-I_{\delta}}{\epsilon})$
.
2.2.2. Evaluation of rate function with Poisson and Gaussian Noise
Further, if Poisson noise is also taken into account then, the linearized robot dynamic equation modifies to
\begin{equation}\delta\ddot{q}(t) + F_1(t)\delta\dot{q}(t) + F_2(t)\delta q(t) = G(t)(\sqrt{\epsilon}w_1(t)+\epsilon w_2(t,\epsilon))\end{equation}
Here, we assume only one Poisson component for simplicity, that is
$p=1$
and therefore Eq. (30) is expressed as
\begin{align}\delta \dot q(t) &= \int_{0}^{t}\Phi_{12}(t,s)G(s)(\sqrt{\epsilon} w_1(t) + \epsilon w_2(t))ds\end{align}
\begin{align}w_2(t) &= c\frac{dN(t)}{dt}\end{align}
where N(t) is Poisson with rate
$\lambda(\epsilon)=\frac{\lambda_0}{\epsilon}$
and c determines the strength of the jerk tremor torque in the robot links. To find the rate function, we need to find
$\mathbb{E} \Bigg\{\exp \left(\frac{1}{\epsilon}\int_{0}^{T}f(t)^T\delta q(t) dt\right)\Bigg\} $
and MGF as already denoted in Eq. (21), as follows:
\begin{multline}\mathbb{E} \Bigg\{\exp \left(\frac{1}{\epsilon}\int_{0}^{T}f(t)^T\delta q(t) dt\right)\Bigg\} = \exp \left(\frac{1}{2\epsilon} \int_{0}^{T}\int_{0}^{T}f(t_1)^T R(t_1,t_2) f(t_2) dt_1 dt_2\right)\\\exp\left(\lambda(\epsilon)\int \left(\exp\left(\int_{s}^{T}f(t)^T\Phi_{12}(t,s)G(s)cdt-1\right)ds\right)\right)\end{multline}
Further, the Gärtner–Ellis limiting logarithmic MGF is as follows:
\begin{align}\Lambda(f) = \frac{1}{2} \int_{0}^{T}\int_{0}^{T} f(t_1)^T R(t_1,t_2) f(t_2)dt_1 dt_2 + \lambda_0\int_{0}^{T}\left(\exp\left(\int_{s}^{T}f(t)^T\Phi_{12}(t,s)G(s)cdt\right)-1\right)ds\nonumber\\\end{align}
Finally, using Eqs. (44) and (45), the rate function can be expressed as
$I(\xi) = \underset{f}{\sup}({<}f,\xi'{>}-\Lambda(f))$
.
2.3. Error energy partition function
Further, to measure the distribution of error energy among the different trajectory error states, error energy partition function is computed in this section. The evaluation of error energy partition function provide a clear picture of error states with large energy difference and therefore, indication of trajectory exit. It is well known in statistical mechanics that if T is the temperature of a system and
$p_T(E)$
is the probability of the system being in the energy state E at temperature T, then
$\frac{p_T(E_1)}{p_T(E_2)} = \exp\left(\frac{E_2-E_1}{kT}\right)$
. Thus, it is interesting to compute the quantity
$\mathbb{E}\left[\exp\left(\frac{E(e)}{kT}\right)\right]$
in the low-temperature limit
$T\rightarrow 0$
. This statistical theory, using Varadhan’s integral lemma is used to find the small deviation in trajectory.
Varadhan’s integral lemma gives the evaluation of
$\underset{T\rightarrow0}{\lim}k T\log\mathbb{E}\left[\exp\left(\frac{E(e)}{kT}\right)\right]$
as
$\underset{q(.)}{\sup}[E(e)-I(q)]$
(see item 3 Section 1.1 (5)), where I(q) is the LDP rate function and E(e) is the difference of error energy between the two energy states, that is
$(E_2-E_1)$
.
Analogously, here
$e(t)=q_d(t) - q(t)$
is the deviation or trajectory error, where
$\left\{q(t) \text{ ,}0\leq t\leq \tau\right\}$
is the system trajectory and
$\left\{q_d(t)\text{,} 0\leq t\leq \tau\right\}$
is the desired trajectory. Further,
$E(e) = \int_{0}^{\tau} (C_1||q_d(t)-q(t)||^2 + C_2 ||\dot{q}_d(t)-\dot{q}(t)||^2)dt$
is the tracking error energy, where
$C_1, C_2$
are the position constants, which determine the relative contributions of the position and velocity error energies to the total error energy. For our linearized system
$e(t)=\delta q(t)$
and
$I(q)\rightarrow I(\delta q)$
. Now, the partition function is as follows:
\begin{align} \underset{T\rightarrow 0}{\lim}kT\log \mathbb{E}\left[\exp\left(\frac{E(e)}{kT}\right)\right] &= \underset{\delta q(.)}{\sup} \left[E\left\{\delta q\right\}-I(\delta q)\right]\\[5pt]&=\underset{\delta q}{\sup}\Bigg[\int_{0}^{\tau}(C_1||\delta q(t)||^2+C_2||\delta\dot{q}(t)||^2)dt\nonumber\end{align}
\begin{align} -\int_{0}^{\tau}\Lambda^{*}(G(t)^{-1}(\delta\ddot{q}(t)+F_1(t)\delta\dot{q}(t)+F_2(t)\delta q(t)))dt\Bigg]\end{align}
To fully evaluate the above Eq. (46),
$E({\delta q})$
and
$I(\delta q)$
needs to be evaluated first. Therefore, the rate function
$\int_{0}^{\tau}\Lambda^{*}(\dot\xi(t))dt$
is evaluated as follows (note:
$\dot\xi$
and
$\xi'$
are same and represents the derivative). The noise process
$\xi_\epsilon$
is expressed as
\begin{equation}\dot{\xi}_\epsilon(t) = \sqrt{\epsilon} w_1(t) + \epsilon w_2 (t,\epsilon) = \sqrt{\epsilon}\sigma \frac{dB(t)}{dt}+\epsilon c\left[\frac{dN(t)}{dt}-\lambda(\epsilon) \right]\end{equation}
with the
$\lambda(\epsilon)=\frac{\lambda_0}{\epsilon}$
.
Remark: The form of the noise as in Eq. (48) with
$\epsilon\rightarrow0$
has been assumed. This means that the noise is weak with the Poisson noise being weaker than the Gaussian noise by the factor
$\sqrt{\epsilon}$
. However, we assume that the Poisson rate
$\propto \frac{1}{\epsilon}$
, that is the rate grows as
$\epsilon \rightarrow 0$
. These assumptions are made so that the deviation probability can be well approximated using LDP rate function, which is easier to compute. The method works well when the noise amplitude is weak and the spike rate is large.
Further, the scaled logarithmic MGF is calculated as follows:
\begin{align}\Lambda_\epsilon &= \epsilon \log \left\{\frac{1}{\epsilon}\int_{0}^{\tau}f(t)^T\dot{\xi}_\epsilon(t)dt\right\}\exp\left\{\frac{\lambda_0}{\epsilon}\left(\int_{0}^{\tau}(\exp (f(t)^Tc)-c^Tf-1)dt\right)\right\}\\\nonumber\end{align}
\begin{align}&= \frac{\sigma^2}{2}\int_{0}^{\tau}||f(t)||^2 dt+\lambda_0\int_{0}^{\tau}\left(\exp(c^Tf(t))-c^Tf-1\right)dt\end{align}
with the corresponding rate function being
$\int_{0}^{t}\Lambda^*(\dot{\xi}(t))$
, where
\begin{equation}\Lambda^* (\dot{\xi}(t)) = \underset{f\in\mathbb{R}^d}{\sup} \left\{f^T\dot{\xi}(t)-\frac{\sigma^2}{2}||f||^2-\lambda_0\left(e^{c^Tf}-c^Tf-1\right)\right\}\end{equation}
Making the quadratic approximations
$(e^{c^{T}f}-1-c^Tf) \approx \frac{1}{2} (f^Tc c^Tf)$
, we get
$\Lambda^*(\xi')\approx \underset{f}{\sup}\{f^T\xi'- \frac{f^T}{2}(\sigma^2 I_d +\lambda_0 c c^T)f=\frac{1}{2}\xi'(\sigma^2 I_d +\lambda_0 c c^T) \}$
. The partition function is expressed as
\begin{multline}\underset{\epsilon\rightarrow 0}{\lim} \log \mathbb{E}\left\{\exp\left(\frac{1}{\epsilon}\int_{0}^{\tau}E(\delta q(t), \delta \dot{q}(t))dt\right)\right\}\\[6pt]=\underset{\delta q(t)\ 0\leq t\leq \tau}{\sup}\left\{\int_{0}^{\tau}E(\delta q(t),\delta\dot{q}(t))dt - \int_{0}^{\tau}\Lambda^*(G(t)^{-1}(\delta\ddot{q}(t))+F_1(t)\delta\dot{q}(t)+F_2(t)\delta q(t))dt\right\}\end{multline}
where
\begin{equation}E(\delta q(t),\delta\dot{q}(t)) = C_1 ||q_d(t)-q(t)||^2 + C_2 ||\dot{q}_d(t)-\dot{q}(t)||^2 = C_1||\delta q(t)||^2 + C_2||\delta q(t)||^2\end{equation}
and
$E(\delta q)= E(e) = \int_{0}^{T}E(\delta q(t),\delta\dot{q}(t))dt$
. Note that we are assuming that
$q_d(t)$
satisfies the noiseless robot differential equation
$ M(q_d(t))\ddot{q}_d + N(q_d(t),\dot{q}_d(t))-\tau(t) = 0$
. The optimum choice of the error trajectory
$\delta q(.)$
that maximizes the quantity within
$\{\}$
of Eq. (52) is given by the Euler–Lagrange equation
\begin{equation}-\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \delta\dot{q}} + \frac{d^2}{dt^2}\frac{\partial \mathcal{L}}{\partial \delta\ddot{q}} + \frac{\partial \mathcal{L}}{\partial \delta q}=0\end{equation}
where
\begin{align}&\mathcal{L}(t, \delta q(t), \delta\dot{q}(t), \delta\ddot{q}(t)) = E(\delta q(t),\delta \dot{q}(t)) - \Lambda^* (G(t)^{-1}(\delta\ddot{q}(t)+ F_1(t)\delta\dot{q}(t)+F_2(t)\delta q(t)))\\\nonumber\end{align}
\begin{align}&= c_1\delta q(t)^T \delta q(t) + c_2\delta \dot{q}(t)^T \delta \dot{q}(t) - \Lambda^* (G(t)^{-1}(\delta\ddot{q}(t)+ F_1(t)\delta\dot{q}(t)+F_2(t)\delta q(t)))\end{align}
Now,
\begin{equation}\mathcal{L}(t, \delta q(t), \delta\dot{q}(t), \delta\ddot{q}(t)) \approx c_1\delta q^T \delta q +c_2 \delta\dot{q}^T\delta\dot{q} -\frac{1}{2}(\eta - \lambda_0c)^T\left[\frac{\sigma^2}{2}*I_d + \lambda_0 c c^T\right]^{-1}(\eta-\lambda_0c)\end{equation}
where
$\eta = G(t)^{-1}[\delta\ddot{q} + F_1(t)\delta\dot{q}+F_2(t)\delta q]$
. The components of the Eq. (57) are as follows:
\begin{align}\frac{\partial \mathcal{L}}{\partial \delta q} &= 2c_1 \delta q- F_2^TG^{-T}\left(\frac{\sigma^2}{2}I_d+\lambda_0cc^T\right)^{-1}\eta\\[-6pt]\nonumber\end{align}
\begin{align}\frac{\partial \mathcal{L}}{\partial \delta \dot{q}} &= 2c_2 \delta \dot{q}- F_2^TG^{-T}\left(\frac{\sigma^2}{2}I_d+\lambda_0cc^T\right)^{-1}\eta\\[-6pt]\nonumber\end{align}
\begin{align}\frac{\partial \mathcal{L}}{\partial \delta \ddot{q}} &= -G^{-1}\left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\eta\\[-6pt]\nonumber\end{align}
\begin{align}\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \delta \dot{q}} &= 2c_2\delta \ddot{q}-(G^{-1}F_1)'^T\left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\eta-(G^{-1}F_1)^T\left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\dot{\eta}\\[-6pt]\nonumber\end{align}
\begin{align}\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \delta \ddot{q}} &= -(G^{-1})'^T\left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\eta-G^{-T}\left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\dot{\eta}\\[-6pt]\nonumber\end{align}
\begin{align}\frac{d^2}{dt^2}\frac{\partial \mathcal{L}}{\partial \delta \ddot{q}} & = -(G^{-1})''^T\left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\eta-G^{-T}\left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\ddot{\eta}- 2(G^{-1})'^T \left(\frac{\sigma^2}{2}I_d + \lambda_0cc^T\right)^{-1}\dot{\eta}\end{align}
where
$\dot{\eta}=(G^{-1})'[\delta\ddot{q}+F_1\delta\dot{q}+F_2\delta q]+G^{-1}[\delta\dddot{q}+F_1\delta\ddot{q}+F_2\delta \dot{q}+F_1'\delta\dot{q}+F_2'\delta q]$
3. Stability Analysis
The aim of this section is to use the standard perturbation theory for stochastic differential equations to derive an upper bound on the noise amplitude parameter
$\epsilon$
that has been used in our LDP-based calculation of the probability of deviation of the system trajectory from the desired trajectory by an amount greater than
$\delta$
. Let
$\delta q(t)=q(t)-q_0(t)$
, where
$q_0(t)$
is the desired nonrandom trajectory and q(t) is the true trajectory in the presence of noise. Further, the value of
$\epsilon$
is computed so that
\begin{equation}\underset{0\leq t\leq T}{\max}\mathbb{E}\left\{\parallel\delta q(t)\parallel^2\right\}<K\delta\end{equation}
where K is constant. Once this computation has been done, the same
$\epsilon$
can be used to compute the LDP probability
$\exp\left(-\epsilon^{-1}\underset{|\delta q|>\delta}{\min}I(\delta q)\right)$
. In other words, the noise amplitude bound is fixed first by the stability condition that the mean square value of the trajectory error is smaller than some constant multiple K of
$\delta$
. Afterward, the probability of deviation greater than
$\delta$
is computed using this noise amplitude bound. It is clear that if the multiple of
$\delta$
is made larger, then the bound on
$\epsilon$
will increase thereby causing the LDP probability to increase. More generally, we can choose an energy function
$V(\delta q)$
rather than
$|\delta q|^2$
and compute a bound on
$\epsilon$
so that
\begin{equation}\underset{0\leq t\leq T}{\max}\mathbb{E}\left\{V(\delta q(t))\right\}<K\delta\end{equation}
and then use this
$\epsilon$
in the proposed LDP calculation.
Linearizing the robot equation, we can write using first-order perturbation theory,
\begin{equation}\delta q''(t)+A(t)\delta q'(t)+B(t)\delta q(t)=\sqrt\epsilon C(t)w(t)\end{equation}
where
$A(t), B(t) \text{ and } C(t)$
are the coefficients and the noise is
$\sqrt\epsilon w(t)$
. Solving this stochastic differential equation using standard state variable methods gives us
\begin{equation}\delta q(t)=\sqrt\epsilon\int_0^t\Phi(t,s)w(s)ds\end{equation}
where
$\Phi(t,s)$
is the state transition matrix. Further,
$\mathbb E\left\{|\delta q(t)|^2\right\}$
is as follows:
\begin{equation}\mathbb E\left\{|\delta q(t)|^2\right\}=\epsilon\sigma^2\int_0^t|\Phi(t,s)|^2ds\end{equation}
where since the noise is white, we have written it as
\begin{equation}\mathbb E\left\{w(t)w(s)^T\right\}=\sigma^2\delta(t-s)\end{equation}
Note that this white noise can be Gaussian or Poissonian or a mixture of the two. We then get that the condition on
$\epsilon$
for the noise to be characterized as being weak is given by
\begin{equation}\epsilon\leq \frac{K\delta}{[\sigma^2\underset{0\leq t\leq T}{\max}\int_0^t|\Phi(t,s)|^2ds]}\end{equation}
More generally if
$V(.)$
is an energy function, then we can approximate
\begin{equation}V(\delta q)\approx V(0)+V'(0)\delta q+(1/2)\delta q^TV''(0)\delta q\end{equation}
and hence
\begin{equation}\mathbb E\{V(\delta q(t))\}\approx V(0)+\frac{1}{2}Tr\left(V''(0)\sigma^2\int_0^t\Phi(t,s)\Phi(t,s)^T\right)\end{equation}
Assuming
$V(0)=0$
, we then get the required bound on
$\epsilon$
as
\begin{equation}\epsilon\leq \frac{K\delta}{[\sigma^2\underset{0\leq t\leq T}{\max}Tr(V''(0)\int_0^t\Phi(t,s).\Phi(t,s)^Tds)]}\end{equation}
For Gaussian noise,
$V(.)$
is usually taken as a quadratic function while for Poissonian and other kinds of non-Gaussian noise, an appropriate choice for
$V(.)$
would be a polynomial function modulated by exponential functions. The choice of V would be dictated by the nature of the pdf of the noise process.
4. Implementation and Results
4.1. Software simulation and results
The presented problem formulation is implemented and validated on a two-link planner robot coded in the MATLAB software. The self-explanatory sequential steps for implementation are provided in Fig. 2 and Algorithm 1 with all variables defined in algorithm comments. The dynamics of the two-axis planner robot is taken as Eq. (9) [Reference Schilling27], where the mass moment if inertia matrix M(q) elements are given by
\begin{equation}\begin{bmatrix}((((1/3)m_1)+m_2)l_1^2+((1/3)m_2l_2^2)+m_2l_2l_1\cos(q2)) & \quad m_2(1/3l_2^2+1/2l_1l_2\cos(q2))\\[8pt]m_2(1/3l_2^2+1/2l_1l_2\cos(q2)) & \quad 1/3m_2l_2^2\\\end{bmatrix}\end{equation}
the Coriolis and centrifugal force matrix C with gravity matrix G that is
$(N(q,\dot{q})) = C+G$
is given by
\begin{equation}G = \begin{bmatrix}((1/2m_1+m_2)+l_1\cos(q1)+(1/2)m_2l_2\cos(q1+q2))g \\[8pt]1/2m_2l_2\cos(q1+q2)g\end{bmatrix}\end{equation}
and
\begin{equation}C = \begin{bmatrix}-m_2l_2l_1\sin(q2)\dot{q}_1\dot{q}_2 - 1/2m_2l_1l_2\sin(q2)\dot{q}_2^2\\[8pt]1/2m_2l_1l_2\sin(q2)\dot{q}_1^2\end{bmatrix}\end{equation}
The relative parameter and symbols are defined in Table I. The calculated probabilities are provided in Table II. There are three observations in Table II with three different noise percentages, that is 0.3%, 0.5%, and 0.8% w.r.t to the input trajectory signal. Figure 3 shows the noise free desired trajectory
$q_0(t) =q_d(t)$
. The actual trajectory created by using the inverse dynamics of robot along with inputs: the desired trajectory and Poisson noise is added to the torque. Afterward, from noisy torque the angular trajectory is generated and the Brownian noise is added to the generated trajectory, the final actual trajectory is shown in Fig. 4. The
$\delta q$
is shown in Fig. 5 along with its envelope. Figure 6 shows the desired trajectory along with the noisy trajectory with envelope, to give a pictorial representation of the boundary (or tube) of the exit probability. The width of envelopes are not constant because the envelope is made randomly proportional to the magnitude of the noise and the randomness is controlled by
$\delta$
(the boundary of the tube). The variable nature of envelope provides the rare event scenario and whenever the noisy trajectory crosses the tube/boundary it is counted as a rare event. Therefore, at some places it touches the trajectory and at some places it is above the noisy trajectory.
Parameters and variables.
Exit probability.
aFor reference, 0.2 probability means exit probability is 20%.
Block diagram showing proposed technique.
Trajectory without noise.
Trajectory with Brownian noise and Poisson process noise. Note that the Poisson process noisy trajectory is generated by adding Poisson Process to forward dynamics profile of the robot.
Plot showing
$\delta q$
with respective envelopes.
4.2. Hardware implementation and results
The hardware implementation is done on a commercially available Omni Bundle robot. It is a 6-DOF robot with three actuated and three non-actuated joints. Further, the Omni Bundle robot communicates with MATLAB software through USB to IEEE 802.3 LAN with the following softwares installed: QUARCTM, Open haptics, and GeomagicTM touch driver. The Omni Bundle robot is also connected to a DC power supply of 20V, 4.5A. The joints data are acquired through MATLAB for analysis. For implementation purpose, only two-link (2-DOF), that is link2 (q2) and link3 (q3) configuration, is used as shown in Fig. 7. The dynamics of robot with parameters is adequately presented and discussed in ref. [Reference Rana, Gaur, Agarwal and Parthasarathy11]. The inertia matrix used for Omni Bundle robot simulation is given as follows
$M =\begin{bmatrix}M11 & M12 & M13\\M21 & M22 & M23\\M31 & M32 & M33\\\end{bmatrix}$
, where
$M11 = K_1 + K_22 \cos(2\theta_2) + K_3 \cos(2\theta_3)+K_4 \cos(\theta_2) \sin (\theta_3)$
,
$M12 = K5 \sin\theta_2$
,
$M13 = 0$
,
$M21 = K5 \sin\theta_2$
,
$M22 = K6$
,
$M23 = -0.5K4 \sin \theta_2 - \theta_3$
,
$M31 = 0$
,
$M32 = -0.5K4 \sin \theta_2 - \theta_3$
,
$M33 = K7$
., where
$K_1= 0.00179819707554751$
,
$K_2 = 0.000864793119787878$
,
$K_3 = 0.000486674040957256$
,
$K_4 = 0.00276612067958414$
,
$K_5 = 0.000308649491069651$
,
$K_6 = 0.00252639617221043$
,
$K_7 = 0.000652944405770658$
,
$K_8 = 0.164158326503058$
,
$K_9 = 0.0940502380783103$
and
$K_{10} = 0.117294768011206$
[Reference Nygaard19]. The link lengths are
$L_1 = 0.132$
and
$L_2 = 0.132$
m. The link masses are
$M_1 = 0.035$
and
$M_2 = 0.1$
kg respectively. The rest of parameters of Omni Bundle are given in ref. [Reference Singla, Parthasarathy and Agarwal20].
Plot showing the desired trajectory and the actual noisy trajectory with respective envelopes
Hardware setup shown with PC connected with Omni Bundle q2 and q3.
Plot showing torque trajectory of two Omni bundle links q2 and q3.
The same Algorithm 1 (obviously with Omni Bundle dynamics) is used to implement on MATLAB Simulink
$^{\circledR}$
by using MATLAB functions. Afterward, the Simulink model is downloaded in Omni bundle to test the algorithm. The sensory data received form Omni Bundle is used to calculate the LDP exit probability of the trajectory. The torque added with Poisson noise is shown in Fig. 8 and the noise angular trajectory is shown in Fig. 9. The calculated
$\delta q$
for both the links is shown in Fig. 10. Afterward, the data are saved in the MATLAB workspace and the exit probability was found to be
$p(q1)=0.1216$
,
$p(q2)=0.1428$
. The Gaussian noise and Poisson noise are considered at 1% of the input trajectory amplitude. The exit probability increases exponentially as the noise percentage in trajectory is increased in the same pattern as demonstrated in Table II of software simulation.
The software and hardware simulation shows that it is important to choose
$\delta$
(width of tube) wisely. This is the only controllable parameter to classify the rare event, and changes made to
$\delta$
would directly result in increase or decrease in number of rare events. Further, the probabilities of Table II matches with the exponential increase of probability when noise is increased with constant
$\delta$
. This probability can be used to further develop control methods as the existing control methods does not include probabilistic approach.
Plot showing noisy angular trajectory after adding Gaussian noise.
Plot showing
$\delta q$
for both links of Omni Bundle.
5. Discussion
In this section, the advantages of LDP approach and extendibility of the proposed LDP method to other robots is discussed.
5.1. Lyapunov versus LDP approach
The following shows why LDP approach is preferred over Lyapunov energy method. Consider the linear SDE Eq. (29).
\begin{equation}\delta q''+F_1(t)\delta q'+F_2(t)\delta q=G(t)w(t)\end{equation}
In state variable notation, Eq. (29) can be written as
\begin{equation}\frac{d}{dt}\left(\begin{array}{c}\delta q_t\\[4pt]\delta\omega_t\end{array}\right)=\left(\begin{array}{c@{\quad}c}0&I\\[4pt]-F_2(t)&-F_1(t)\end{array}\right)\left(\begin{array}{c}\delta q(t)\\[4pt]\delta\omega(t)\end{array}\right)+\left(\begin{array}{c}0\\[5pt]G(t)\end{array}\right).dB(t)\end{equation}
The solution to Eq. (78) is given by
\begin{equation}\left(\begin{array}{c}\delta q_t\\[4pt]\delta\omega_t\end{array}\right)=\int_0^t\Phi(t,s)\left(\begin{array}{c}0\\[4pt]G(s)\end{array}\right).dB(s)\end{equation}
and hence the correlation matrix of the state vector is given by
\begin{align}R(t)&=\Bbb E\Bigg[\left(\begin{array}{c}\delta q_t\\[4pt]\delta\omega_t\end{array}\right).[\delta q_t^T,\delta\omega_t^T]\Bigg]\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\end{align}
\begin{align} =\int_0^t\Phi(t,s)\left(\begin{array}{c}0\\[5pt]G(s)\end{array}\right)\left[0,G(s)^T\right]\Phi(t,s)^Tds\nonumber\end{align}
\begin{align}=\int_0^t\left(\begin{array}{c@{\quad}c}\Phi_{12}(t,s)G(s)G(s)^T\Phi_{12}(t,s)^T&\Phi_{12}(t,s)G(s)G(s)^T\Phi_{22}(t,s)^T\\[5pt]\Phi_{22}(t,s)G(s)G(s)^T\Phi_{12}(t,s)^T&\Phi_{22}(t,s)G(s)G(s)^T\Phi_{22}(t,s)^T\end{array}\right)ds\\[-18pt]\nonumber\end{align}
Now let Q be a positive definite matrix. Define the stochastic Lyapunov energy at time t as
\begin{equation}V(t)=\frac{1}{2}[\delta q_t^T,\delta\omega_t^T]Q\left(\begin{array}{c}\delta q_t\\\delta\omega_t\end{array}\right)\end{equation}
Then, we get by applying the Ito differential rule:
\begin{align}\frac{d}{dt}\Bbb E(V(t))&=\Bbb E\left[\frac{dV(t)}{dt}\right]\\[-4pt]\nonumber\end{align}
\begin{align}&=\Bbb E\left([\delta q_t^T,\delta\omega_t^T]\left(\begin{array}{c@{\quad}c}0&-F_2^T\\[5pt]I&-F_1^T\end{array}\right)Q\left(\begin{array}{c}\delta q_t\\[4pt]\delta\omega_t\end{array}\right)\right)+(1/2)Tr\left(Q\left(\begin{array}{c}0\\[4pt]G(t)\end{array}\right)[0,G(t)^T]\right)\end{align}
Equation (85) evaluates to
\begin{align}\frac{d}{dt}\Bbb E(V(t)) = & -\int_0^tTr(F_2(t)Q_{21}\Phi_{12}(t,s)G(s)G(s)^T\Phi_{12}(t,s)^T)ds\nonumber\\[5pt] & -\int_0^tTr(F_1(t)Q_{22}\Phi_{22}(t,s)G(s)G(s)^T\Phi_{22}(t,s)^T)ds\nonumber\\[5pt] & -\int_0^tTr(F_2(t)Q_{22}\Phi_{22}(t,s)G(s)G(s)^T\Phi_{12}(t,s)^T)ds Tr(F_1(t)Q_{21}\Phi_{12}(t,s)G(s)G(s)^T\Phi_{21}(t,s)^T)ds\nonumber\\[5pt] & +(1/2)Tr(Q_{22}G(t)G(t)^T)\end{align}
If we assume that
$F_1,F_2$
are approximately constant matrices for large times and the eigenvalues of the matrix
\begin{equation}X=\left(\begin{array}{c@{\quad}c}0&I\\[5pt]-F_2(t)&-F_1(t)\end{array}\right)\end{equation}
all have negative real part, then using the eigen decomposition of this matrix X in the spectral representation form
$\sum_k\lambda_ku_kv_k^T$
, where
$u_k$
right eigen vector,
$v_k$
is the left eigen vector, both corresponds to the eigen value
$\lambda_k$
of X. Therefore, the state transition matrix is approximately given by
\begin{equation}\Phi(t,s)=\sum_kexp(\lambda_k(t-s))u_kv_k^T\end{equation}
and the above expression (86) for the rate of change of the average Lyapunov energy is approximately given by an expression of the form (for large times t)
\begin{equation}\frac{d}{dt}\Bbb E(V(t))=\sum_{k,m}c(k,m)(exp(\lambda_k+\lambda_m)t-1)+(1/2)Tr(Q_{22}G(t)G(t)^T)\end{equation}
It is easily seen that under the condition that
$F_1,F_2,G$
are approximately constant matrices, the first two terms on the rhs of the above expression behave with t as
$-(1-exp(\lambda_kt))/\lambda_k$
, which converges to a constant as
$t\rightarrow\infty$
provided that the eigenvalues of the matrix X all have negative real part. Thus, we have proved the result that the rate of increase of the Lyapunov energy cannot grow with time for large times, it can atmost be bounded in time. Further, if
$\epsilon> \mathbb{E}\{\delta q^2(t)\} $
then
\begin{equation}\frac{d}{dt}\mathbb{E}\{V(t)\} > 0\end{equation}
and the Lyapunov method will fail in this case. Therefore, the two major concerns in Lyapunov method are that: first, although the above analysis shows that the rate of increase of Lyapunov energy is bounded by a constant, however, the Lyapunov energy may grow linearly with time. Second, in sudden jerks in
$\delta q$
, that is
$\epsilon> \mathbb{E}\{\delta q^2(t)\} $
.
On the other hand if we use LDP, the rate function
\begin{equation}I(\delta q) =-\frac{1}{2} \int_{0}^{T}\int_{0}^{T} \delta q(t_1)^TR^{-1}(t_1,t_2)\delta q(t_2)\end{equation}
The probability of exit from the stability zone is given by
\begin{equation}Pr\left\{\underset{0\leq t\leq T}{\max}|\delta q|>\delta\right\} \approx \exp\left(-\frac{1}{2\epsilon}\left(\underset{||\delta q||>\delta}{\min}I(x)\right)\right)\end{equation}
where
$||\delta q|| = \underset{0\leq t\leq T}{\max}|\delta q(t)|$
. In this case the only priority is to minimize the probability, which can be easily achieved by choosing the maximum value out of the set of minimum values of rate function with the boundary conditions where the rare event occurs that is
$||\delta q|| >\delta$
. Therefore, we choose control coefficient K in a prescribed region R, so that this probability is minimum, that is
\begin{equation}\hat{K} = \underset{a \in R}{\arg\max}\left(\underset{||\delta q||>\delta}{\min}I(x)\right)\end{equation}
The Lyapunov method gives in general the result that the average Lyapunov energy can grow linearly with time, when the noise is present and hence, we can atmost design the controller, so that the rate of growth is small. The LDP method of controller design, on the other hand, guarantees that the probability of the error exceeding a given threshold over the entire duration of the trajectory is very small and hence the error energy with a large probability will stay bounded. Equivalently, with a large probability, the rate of error energy increase will be zero if we use LDP. Further, in Lyapunov energy method, for general boundaries, trajectory exit may happen with rare spikes in
$\delta q$
. However, the LDP method-based controller design mitigates any abrupt change in the trajectory of this sort as this theory at its core is developed considering the probabilities of rare events only. Therefore, it is suitable even if the boundaries has arbitrary shape, whereas the Lyapunov method is considered only for general boundary shapes.
5.2. Computational complexity
A comparison of the computational complexities required to assess the performance of the Lyapunov energy-based method and the LDP-based method. It is seen that in the Lyapunov-based method, the Lyapunov energy rate of change converges to nearly a positive constant dictated by noise terms in a time duration T, say. Therefore, the time taken by Eq. (89) to converge is expressed as follows:
\begin{equation}\frac{1}{2}\min_{k,m}|\operatorname{Re}(\lambda_k+\lambda_m)|=\frac{1}{2}\min_k|\operatorname{Re}(\lambda_k)|=T\end{equation}
If the number of computations required in computing the Lyuapunov energy from the state of the system per iteration is p and the time discretization step size is
$\Delta$
, then the total computational complexity required in computing all the values of the Lyapunov energy until its rate of change converges to a constant value is approximately pN where
$N=T/\Delta$
. On the other hand, we have noted that the LDP rate function of the perturbation in the angular position velocity process over the time interval [0,T] is computed as a quadratic form
$I(\xi)=(1/2)]\xi^TR^{-1}\xi$
where R is an
$N\times N$
correlation matrix with
$N=T/\Delta$
. The probability of exit from the stability zone over this time duration is according to LDP theory expressed as
\begin{equation}\exp\left(\left(\frac{-1}{\epsilon}\right)\underset{\parallel\xi\parallel>\delta}{\min}I(\xi)\right)\end{equation}
which evaluates to
\begin{equation}\exp\left(\frac{-\delta^2}{2\epsilon\lambda_m}\right)\end{equation}
where
$\lambda_m$
is the maximum eigenvalue of the positive definite matrix R. Note that here we are using the
$L^2$
-norm and not the
$L^{\infty}$
-norm for defining the stability zone. The computational complexity required in computing the eigenvalues of R is the same as that of finding the roots of the
$N{\rm th}$
degree polynomial
$\left\{\det(R-\lambda I)\right\}$
and this is quite large even when using fast numerical algorithms. However, using Rayleigh’s variational principle for Hermitian matrices, this complexity can be reduced. Further, the LDP method is preferred to overcome the limitations of Lyapunov method already discussed in Section 5.1, with less weightage for complexity.
5.3. Extension of proposed LDP method to other robots
The proposed method may be applied on any type of robot where the robot differential equation has the following generalized form:
\begin{equation}\dot{X} = F(X|\theta)+\sqrt{\epsilon}G(X|\theta)W(t)\end{equation}
where X is the actual trajectory,
$F(X|\theta)$
and
$G(X|\theta)$
denotes the dependence of X on
$\theta$
as a function, and
$\epsilon$
denotes the scaling factor of noise W(t). The objective is to determine the control parameters
$\theta$
so that if d(t) is the desired robot trajectory, then the error
$||X(t)-d(t)||,0\leq t\leq T$
is a minimum in some appropriate sense. We have interest to minimize the
\begin{align}\lim\limits_{T\rightarrow0}kT\ln \mathbb{E}\left\{\exp\left(\frac{1}{kT}\int_{0}^{T}||X(t)-d(t)||^2 dt\right)\right\}=\sup\limits_{X}\int_{0}^{T}||X(t)-d(t)||^2-I(X|\theta)\end{align}
where
$I(X|\theta)$
is the rate function of X(t) and
$\epsilon=kT$
. The supremum of (98) is evaluated using the Euler–Lagrange differential equation. In similar way the proposed theory is applied to any robot as in ref. [Reference Rana, Agarwal, Gaur and Parthasarathy21].
Statistically, Eq. (8) is a generic model and is ideally the same for any number links with any degree of freedom. The basic electrical elements of robots are nearly the same except for the kind of motion they provide. They all have angular position sensors, motors, rigid body (i.e., links), and a motion controller. The use of other robot would not affect the flexibility of the proposed method as in two-link 2-DOF robot, all matrix elements of Eq. (8) are present like Coriolis and centrifugal force matrix, which are absent in one-link robotic system. The hardware results are incorporated on a commercial 6-DOF robot(not planner). However, to test our algorithm we took only two links (2-DOF). Therefore, the results and methodology may be replicated and extended on any available robotic system, which satisfies the generalized Eq. (8).
6. Conclusion
The presented application of LDP theory for trajectory tracking problem has been derived and presented with conclusive results. The LDP probabilities exit from desired zone have been evaluated for low-amplitude noise given by a superposition of Gaussian and Poisson processes. Further, the novel error energy partition function is presented as fourth-order linear differential equation to give the exact time and condition of trajectory exit form stable domain. The condition for noise amplitude
$\epsilon$
to get a stable system is also derived in stability analysis. It is proved conclusively that the LDP, a probability-based approach is better than the energy-based Lyapunov approach. The successful simulation and implementation on hardware with small noise
$\epsilon\rightarrow 0$
provided desirable results. Further, the simulation and hardware results are in good agreement. The future scope of this work is to develop LDP-based probabilistic control methods.
Conflicts of Interest
The authors do not have any conflicts of interest.
